Численные решения обыкновенных дифференциальных уравнений

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

      Если  учесть четыре члена ряда, то получим  формулу Рунге-Кутта:

      y_i=y_i-1+1/6 ^. (K₁^(i-1)+2 ^. K₂^(i-1)+2 ^. K₃^(i-1)+K₄^(i-1)),

где K_j^(i-1) определяется из формул :    

      K₁^(i-1)=h ^. f(x_i-1;y_i-1);

     K₂^(i-1)=h ^. f(x_i-1+h/2;y_i-1+K₁^(i-1)/2);

     K₃^(i-1)=h ^. f(x_i-1+h/2;y_i-1+K₂^(i-1)/2) ;

      K₄^(i-1)=h ^. f(x_i-1+h;y_i-1+K₃^(i-1)).

      Для решения обыкновенного дифференциального  уравнения n-го порядка необходимо сначала свести это уравнение к системе уравнений 1-го порядка. Начальными условиями такой системы будут : 

      y(x₀)=y⁰;

     y₁(x₀)=y₁⁰;

    y₂(x₀)=y₂⁰;

     . . . . . . . .

     y_n-1(x₀)=y⁰_n-1. 

      Решение находится с помощью простого обобщения метода Рунге-Кутта для уравнения 1-го порядка. Но теперь на каждом шаге находится значение не одной y функции , а n функций y;y₁;y₂;...;y_n-1.

      Пусть дано дифференциальное уравнение второго  порядка y”=f(x, y, y’) с начальными условиями y(x₀)=y₀ ; y;(x₀)=z₀.

      Представим  сначала его в виде системы уравнений

  y’=z

            z’=f(x, y,z).

с начальными условиями  y(x₀)=y₀ ; z(x₀)=z₀

      Так как получилось два уравнения, то теперь будут и две правые части:

f₁(x, y, z) ;f₂(x, y, z). На каждом шаге теперь будем пользоваться формулами :

      y_i=y_i-1+1/6 ^. (K1₁^(i-1)+2 ^. K1₂^(i-1)+2 ^. K1₃^(i-1)+K1₄^(i-1));

      z_i=z_i-1+1/6 ^. (K2₁^(i-1)+2 ^. K2₂^(i-1)+2 ^. K2₃^(i-1)+K2₄^(i-1)),

где коэффициенты K1_j соответствуют функции y, а коэффициенты K2_j – функции z. Они вычисляются по следующим формулам :

      K1₁^(i-1)=h ^. f₁(x_i-1;y_i-1;z_i-1);

      K1₂^(i-1)=h ^. f₁(x_i-1+h/2;y_i-1+1/2 ^. K1₁^(i-1);z_i-1+1/2 ^. K2₁^(i-1));

      K2₁^(i-1)=h ^. f₂(x_i-1;y_i-1;z_i-1);

      K2₂^(i-1)=h ^. f₂(x_i-1+h/2;y_i-1+1/2 ^. K1₁^(i-1);z_i-1+1/2 ^. K2₁^(i-1));

и так далее  для K1₃; K1₄; K2₃; K2₄. 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

2 Решение обыкновенного дифференциального уравнения первого порядка.

      2.1 Аналитическое решение. 

cosy-y’√x=0,         [3.8;4.8],       y₀(3.8)=10.2

Данное уравнение  представляет собой уравнение с разделяющимися переменными.

cosy-y’√x=0,

cosy-(dy∕dx)√x=0,/×(dx∕cosy)

dx∕√x-dy∕cosy=0,

dx/√x=dy/cosy,

∫dx/√x=∫dy/cosy.

Вычислим отдельно левую часть:

∫dx/√x=2√x

Вычислим  левую часть:

∫dy/cosy=ln│tg(y∕2+π/4)│

В итоге у  нас получилось:

ln│tg(y∕2+π/4)│= 2√x+ ln│c│

tg(y∕2+π/4)=ce^2√x

y∕2+π/4=arctg (ce^2√x)+πn,n€z.

y= -π/2 + arctg (ce^2√x)+2 πn, n€z.

-π/2 + arctg (ce^2√3,8)+2 πn=10,2,

arctg (ce^2√3,8)=1/2(10,2+π/2-2π),при n=2,

ce^2√3,8= tg(-0,76),

c=tg(-0,76)/e^2√3,8.

Далее на отрезке[3.8;4.8] вычислим приближенное значение y(x). Результаты, сделанные в ручную, поместим в таблицу 1. 

2.2 Решение методом  Эйлера.

      С помощью метода Эйлера найти приближенно  решение уравнения  на cosy-y’√x=0, отрезке [3,8;4,8] с шагом h=0,1 при начальных условиях y₀(3.8)=10.2

      Начиная с начальной точки (x₀;y₀) и применяя расчетную формулу метода  y_i = y_i-1 +h ^. f(x_i-1;y_i-1) , будем последовательно получать значения приближенного решения :

            y₁=y₀+h ^. f(x₀;y₀)

            y₂=y₁+h ^. f(x₁;y₁) и так далее.

Покажем на примере как применять данные формулы:

f(x₀;y₀)=f(3.8;10.2)=cos10.2∕√3.8=0.5048

y₁= y₀+h ^. f(x₀;y₀)

y₁=10.2+0.1×0.5048=10.251 

f(x₁;y₁)=f(3.9;10.251)=cos10.251∕√3.9=0.4983

y₂=y₁+h ^. f(x₁;y₁)

y₂=10.251+0.1×0.4983=10.299 и так далее.

Расчетная таблица 

 

      Программа метода на языке Qbasic :

CLS

h = .1: y= 10.2: x = 3.8

PRINT #1, USING "x=#.#   y=#.#####"; x; y

FOR x =3.8 TO 4.8 STEP .1

F=COS(y)/SQR(x)

Y= Y + H * F

PRINT #1, USING "x=#.#    y=#.#####"; x; y

NEXT x

END 

Результаты счета программы:

x=3.8    y=10.16336

x=3.9   y=10.25048

x=4.0    y=10.29865

x=4.1   y=10.34765

x=4.2    y=10.39658

x=4.3    y=10.44457

x=4.4    y=10.49199

x=4.5    y=10.53872

x=4.6    y=10.58542

x=4.7    y=10.63091

x=4.8    y=10.67583

   

Как видно результаты ручного счета почти полностью совпали с ответами, полученными в программе.

                                                                   
 

                    
 
 

Блок-схема  метода Эйлера 

начало

x0=3.8:h= .1:y0=10.2:x=x0:y=y0:х=х0

печать x, у

                F=f(x,y)

 

                                                             y=y+F*h

 

 печать x,y

x=x+h

                                                                          x<x2                  да                                                                                                                                                                                                                 

                                                                   

                                                                                  нет

                                                                        конец                        

                                                                                                               

Здесь  f(x,y) = cosy/√x 

2.3 Решение методом  Рунге-Кутта.

      Решить  дифференциальное уравнение cosy-y’√x=0, на отрезке [3.8;4.8], с начальными условиями y₀(3.8)=10.2.

      Начиная от точки (x₀;y₀), взятой из начального условия, подставляя координаты x₀ и y₀ вместо x_i-1, y_i-1 соответственно в формулы для получения коэффициентов K_j^(i-1) , получаем их значения. Затем, подставляя найденные коэффициенты в формулу Рунге-Кутта для нахождения y_i , найдем новое значение y₁ , соответствующее значению аргумента x₁=x₀+h и так далее.

                  Расчетная таблица.