3.3. Метод рунге-кутта 2-го порядка

Пусть имеем дифференциальное уравнение y'=f(x,y) с начальными условиями у(x₀)=y₀.

Ищем решение на отрезке [x₀,x_n].

Пусть имеем точку (x_m,y_m), принадлежащему искомому решению. Для того, чтобы найти следующую точку проведем касательную к кривой в точке (x_m,y_m) до пересечения с прямой x=x_m+0.5, где x_m+0.5=x_m+h/2. Тогда получим координату (по формуле Эйлера)

y_m+0.5=y_m+h/2·f(x_m,y_m).

Теперь найдём тангенс угла наклона касательной в т.В (x_m+0.5,y_m+0.5), (прямая L). Через точку А проведём прямую L'||L. Ординату точки пересечения прямых L' и x=x_m+1 возьмём в качестве y_m+1.

Таким образом

y_m+0.5= y_m+h/2·f(x_m,y_m), x_m+0.5=x_m+h/2

y_m+1= y_m+h·f(x_m+0.5,y_m+0.5), x_m+1=x_m+h

Для системы дифференциальных уравнений

y_i'= f_i(x,y₁,y₂,…y_k),

y_i(x₀)= y_i⁰, i= 1,2,…k

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

y_i^m+0.5= y_i^m+h/2·f_i(x_m,y₁^m,y₂^m,…y_k^m), x_m+0.5= x_m+h/2;

y_i^m+1= y_i^m+h·f_i(x_m+0.5,y₁^m+0.5,y₂^m+0.5,…y_k^m+0.5);

x_m+1= x_m+h.

ПРОГРАММА РУНГЕ-КУТТА 2-ГО ПОРЯДКА

Program Runge_Kutte_Method;

Uses Crt;

Type Vec=Array[1..6] of Real;

Var y, y0, k, f:Vec;

i, KolKom, KolRec, KolForPrint, count:Byte;

x, xk, h, s:Real;

Fout:text;

NameOut:string[12];

Procedure f_dir(y:vec; var f:vec);

var r:vec;

begin

r[1]:=k[1]*y[1];

r[2]:=k[2]*y[2];

r[3]:=k[3]*y[2];

r[4]:=k[4]*y[3];

r[5]:=k[5]*y[2];

f[1]:=-r[1]+r[5] ;

f[2]:=r[1]-r[2]-r[3]+r[4]-r[5];

f[3]:=r[2]-r[4];

f[4]:=r[3]

end;

BEGIN

ClrScr;

x:=0; xk:=15; h:=0.1;count:=0;s:=0;

write('введите количество компонентов ');readln(KolKom);

For i:=1 to KolKom do

begin

write(i,’введите концентрации компонентов ');readln(y[i])

end;

write('введите количество химических реакций ');readln(KolRec);

For i:=1 to KolRec do

begin

write(i,'введте константы химических реакции');readln(k[i])

end;

write('введите число точек для вывода значения на экран ');

readln(KolForPrint);

write('введите имя файлов для результатов');

readln(NameOut);

ClrScr;

assign(fout,NameOut);

rewrite(fout);

write(' время ');write(fout,' время ');

For i:=1 to KolKom do begin

write(' y[',i,']');

write(fout,' y[',i,']'); end;

writeln(' сумма '); writeln(fout,' сумма ');

write(x:5:1);write(fout,x:5:1);

for i:=1 to KolKom do

begin s:=s+y[i];write(y[i]:8:5);write(fout,y[i]:8:5) end;

writeln(s:8:2);writeln(fout,s:8:2);

repeat

f_dir(y,f);

For i:=1 to KolKom do

begin

y0[i]:=y[i];

y[i]:=y0[i]+0.5*h*f[i];

end;

f_dir(y,f);

For i:=1 to KolKom do y[i]:=y0[i]+h*f[i];

x:=x+h;

count:=count+1;

if count=round(xk/h/KolForPrint) then

begin

write(x:5:1);write(fout,x:5:1);s:=0;

for i:=1 to KolKom do

begin

s:=s+y[i];write(y[i]:8:5);write(fout,y[i]:8:5);

end;

writeln(s:8:2);writeln(fout,s:8:2);count:=0

end;

until x>xk;

close(fout);

readkey

end.