x_i = х₀ + ih  и y_i = y(x_i),   где   i = 0, 1, 2, ...,

      x_i - узлы сетки,

      y_i - значение интегральной функции в узлах.

      Иллюстрации к решению приведены на рисунке 2.

      Проведем  прямую АВ через точку (x_i, y_i) под углом α. При этом tg α = f(x_i, y_i)

      В соответствий с геометрическим смыслом  задачи, прямая АВ является касательной  к интегральной функции. Произведем замену точки интегральной функции точкой, лежащей на касательной АВ.

      Тогда y_i+1 = y_i + Δy

      Из  прямоугольного треугольника ABC

      Приравняем  правые части tg α = f(x_i, y_i) и . Получим

      Отсюда  Δу = h ∙ f(x_i, y_i).

      Подставим в это выражение формулу y_i+1 = y_i + Δy, а затем преобразуем его. В результате получаем формулу расчета очередной точки интегральной функции:

       . 
 
 
 
 
 
 
 
 
 
 

        
 
 
 

      Рисунок 2. Метод Эйлера.

 

      Из  формулы  видно, что для расчета каждой следующей точки интегральной функции необходимо знать значение только одной предыдущей точки. Таким образом, зная начальные условия, можно построить интегральную кривую на заданном промежутке.

      Блок-схема  процедуры решения дифференциального  уравнения методом Эйлера приведена на рисунке 3.

      

F(x, у) - заданная функция – должна

быть описана  отдельно.

Входные параметры:

Х0, XK—начальное и конечное

значения независимой переменной;

Y0 – значение y₀ из начального условия

y(x₀) = y₀;

N - количество отрезков разбиения;

Выходные параметры:

У - массив значений искомого решения

в узлах сетки. 
 

Рисунок 3. Блок-схема процедуры решения дифференциального уравнения методом Эйлера. 

      Метод Эйлера - один из простейших методов  численного решения обыкновенных дифференциальных уравнений. Но существенным его недостатком является большая погрешность вычислений. На рисунке 2 погрешность вычислений для i-гo шага обозначена ε. С каждым шагом погрешность вычислений увеличивается.

      2. 5. Метод Эйлера модифицированный

      Для уменьшения погрешности вычислений часто используется модифицированный метод Эйлера. Этот метод имеет так же следующие названия: метод Эйлера-Коши или метод Рунге-Кутта второго порядка точности.

      Пусть дано дифференциальное уравнение первого  порядка 

      

      с начальным условием:

      

      Выберем шаг h и введём обозначения:

      x_i = x₀ + ih  и y_i = y(x_i),   где   i = 0, 1, 2, ...,

      x_i  - узлы сетки,

      y_i - значение интегральной функции в узлах.

      При использовании модифицированного  метода Эйлера шаг h делится на два отрезка.

      Иллюстрации к решению приведены на рисунке 4. 

        
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

      Рисунок 4. Метод Эйлера модифицированный 

      Проведем  решение в несколько этапов:

1. Обозначим точки: А(х_i, y_i,), C(x_i + h/2, y_i + h/2 ∙ f(x_i, y_i)) и B(x_i+1, y_i+1);
2. Через точку А проведем прямую под углом α, где tg α = f(x_i, y_i);
3. На этой прямой найдем точку С(х_i + h/2, y_i + h/2 ∙ f(x_i, y_i));
4. Через точку С проведем прямую под углом α1, где tg α1 = f(x_i + h/2,y_i + h/2 ∙ f(x_i, y_i));
5. Через точку А проведем прямую, параллельную последней прямой;
6. Найдем  точку B(x_i+1, y_i+1).   Будем считать   B(x_i+1, y_i+1)  решением дифференциального уравнения при х = x_i+1;
7. После проведения вычислений, аналогичных вычислениям, описанным в методе Эйлера, получим формулу для определения значения у_i+1:

      y_i+1 = y_i + h ∙ f(x_i + h/2, y_i + h/2 ∙ f(x_i, y_i)).

      Модифицированный  метод Эйлера дает меньшую погрешность. На рисунке 4 это хорошо видно. Так величина εl характеризует погрешность метода Эйлера, а ε - погрешность метода Эйлера модифицированного.

      Блок-схема  процедуры решения дифференциального  уравнения методом Эйлера модифицированным приведена на рисунке 5. 

F(x, у) - заданная функция - должна

быть описана отдельно.

Входные параметры: 
Х0, XК - начальное и конечное

значения независимой

переменной;

Y0 – значение y₀ из начального условия

y(x₀)=y₀;

N - количество отрезков разбиения;

Выходные  параметры:

Y - массив значений искомого решения

в узлах сетки. 
 

      Рисунок 5. Блок-схема процедуры решения дифференциального уравнения методом Эйлера модифицированным. 

      2.6 Метод Рунге-Кутта 4 порядка

      Пусть дано дифференциальное уравнение первого  порядка    с начальным условием y(x₀)=y_0.   Выберем шаг h и введем обозначения:

      x_i = x₀ + ih  и y_i = y(x_i),   где   i = 0, 1, 2, ... .

     Аналогично  описанному выше методу производится решение

дифференциального уравнения. Отличие состоит в  делении шага на 4 части.

          Согласно методу Рунге-Кутта четвертого  порядка, последовательные значения  y_i   искомой функции y определяются по формуле:

      y_i+1 = y_i +∆y_i                 где i = 0, 1, 2 ...

              ∆y=(k1+2*k2+2*k3+k4)/6                          

                    

      a числа k1, k2 ,k3, k4 на каждом шаге вычисляются по формулам:

      k1  = h*f(x_i, y_i )

      k2  = f (x_i +h/2, y_i +k1 /2)*h

      k3  = F(x_i +h/2, y_i +k2 /2)*h

      k4  = F(x_i +h, y_i +k3 )*h

      Это явный четырехэтапный метод 4 порядка  точности. 
 
 

      Блок-схема  процедуры решения дифференциального  уравнения методом Рунге-Кутта приведена на рисунке 6.

F(x, у) - заданная функция - должна

быть описана отдельно.

Входные параметры: 
Х0, XК - начальное и конечное

значения независимой

переменной;

Y0 – значение y₀ из начального условия

y(x₀)=y₀;

N - количество отрезков разбиения;

Выходные  параметры:

Y - массив значений искомого решения

в узлах сетки.

        

      

      

      

        

        

      

        
 
 
 
 

      2. 7. Решение поставленной задачи методами Эйлера и Рунге-Кутта

      2. 7. 1. Метод Эйлера

      1. Строим оси координат;

      2. Отмечаем A(0; 0) – первую точку интегральной кривой;

      3. Ищем угол наклона касательной к графику в точке A:

      

      

      

        

      4. Строим касательную l₀ в точке А под углом α₀;

      5. Находим х₁ по формуле: x_i = х₀ + ih, где h – шаг интегрирования

      x₁ = 0 + 1 · 0,1 = 0,1;

6. Проводим прямую x = x₁ = 0,1 до пересечения с прямой l₀, отмечаем точку B(x₁; y₁);

      7. Ищем y точки B:

      Из  прямоугольного треугольника ABC ,

      Δy = y₁ – y₀,

      Δx = x₁ – x₀ = h,

      f(x₀; y₀) = (y₁ – y₀)/h =>

      y₁ = y₀ + h · (f(x₀; y₀)) = 0 + 0,1 · f(0;0) = 0 + 0,1 · 2 = 0,2

      Следовательно, точка B имеет координаты (0,1; 0,2). 
 

      2.7.2. Метод Рунге-Кутта 4 порядка

      1. Строим оси координат;

      2. Отмечаем А(0; 0) – первую точку интегральной кривой;

      3. Ищем угол наклона касательной  к графику в точке A:

        
 

        

      

      4. Строим касательную l₀ в точке А под углом α₀;

      5. Находим х₁ по формуле: x_i = х₀ + ih

      x₁ = 0 + 1 · 0,1 = 0,1;

5. Находим по формулам:

     k1=0,1·f(0,0)=0,1·2=0,2

     k2=0,1· f(0,1+0,1/2;0+0,2/2)= 0,2

     k3=0,1· f(0,1+0,1/2;0+0,2/2)= 0,2

       k4=0,1· f(0,1+0,1;0+0,2)= 0,2

       ∆y₁=(0,2+2·0,2+2·0,2+0,2)/6=0,2

       ∆y₂=0+0,2=0,2

     Следовательно, следующая точка графика решения  имеет координаты (0,1; 0,2)  
 

      . Алгоритм решения  задачи

      3. 1. Алгоритмы подпрограмм

      3. 1. 1. Подпрограмма  метода Эйлера 
 

        
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

      3.1.2 Подпрограмма метода Рунге-Кутта 4 порядка