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.

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

Этот метод один из самых распространённых методов интегрирования дифференциальных уравнений.

Для одиночного дифференциального уравнения y'= f(x,y), y(x₀)=y₀ расчётные формулы имеют следующий вид:

y_m+1= y_m+h/6(k₁+2k₂+2k₃+k₄),

где k₁= f(x_m,y_m);

k₂= f(x_m+,y_m+);

k₃= f(x_m+,y_m+);

k₄= f(x_m+h,y_m+hk₃);

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+1= y_i^m+(k_1i+2k_2i+2k_3i+k_4i), i=1,2,…k,

где k_1i= f_i(x_m,y₁^m,y₂^m,…y_k^m);

k_2i= f_i(x_m+h/2,y₁^m+k₁₁,y₂^m+k₁₂,…y_k^m+k_1k);

k_3i= f_i(x_m+,y₁^m+k₂₁,y₂^m+k₂₂,…y_k^m+k_2k);

k_4i= f_i(x_m+h,y₁^m+hk₃₁,y₂^m+hk₃₂,…y_k^m+hk_3k).

PROGRAM RYNGE_KYTT;

USES CRT;

VAR

X,X0,Y0,XMAX:REAL;

YT,DY,K1,K2,K3,K4:REA;

H:REAL;

I,N:LONGINT;

FUNCTION F(X,Y:REAL):REAL;

BEGIN

F:=(-2*Y/X+EXP(-X*X)/X;

END;

BEGIN

WRITE(‘H=’); READLN(H);

WRITE(‘X0=’); READLN(X0);

WRITE(‘XMAX=’); READLN(XMAX);

WRITE(‘Y0=’); READLN(Y0);

N:=ROUND((XMAX-X0)/H);

YT:=-EXP(-X0*X0)/(2*X0*X0)+15/(X0*X0);

WRITELN(‘X0=’,X0:8:4,’ Y0=YT= ’,YT:8:4);

FOR I:=1 TO N DO

BEGIN

K1:=F(X0,Y0);

K2:=F(X0+H/2,Y0+H*K1/2);

K3:=F(X0+H/2,Y0+H*K2/2);

K4:=F(X0+H,Y0+H*K3);

Y0:=Y0+H/6*(K1+2*K2+2*K3+K4);

X0:=X0+H;

YT=-EXP(-X0*X0)/(2*X0*X0)+15/(X0*X0);

DY:=ABS(YT-Y0);

WRITELN(‘X=’,X0:8:4,’YT=’,YT:8:4,’YR=’,YO:8:4,’DY=’,DY:10:8);

END;

READKEY;

END.

Program Runge_Kutte_4_Method;

Uses Crt;

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

Var y, y0, k, ka, 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]);

y0[i]:=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

ka[i]:=h*f[i];

y[i]:=y0[i] + ka[i]/2;

end;

f_dir(y,f);

For i:=1 to KolKom do

begin

ka[i]:=ka[i] + 2*h*f[i];

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

end;

f_dir(y,f);

For i:=1 to KolKom do

begin

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

end;

f_dir(y,f);

For i:=1 to KolKom do

begin

y[i]:=y0[i] + (ka[i] + h*f[i])/6;

y0[i]:=y[i]

end;

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.

4. МОДЕЛИРОВАНИЕ КИНЕТИКИ ХИМИЧЕСКОЙ РЕАКЦИИ

Пример.

Задано:

1. Схема механизма химической реакции

К₁ К₂

А В С

К₅К₃К₄

Д

2. Константы скоростей отдельных стадий реакции

К₁ = 0,76 с^-1

К₂ = 0,9 с^-1

К₃ = 0,5 с^-1

К₄ = 0,45 с^-1

К₅ = 0,6 с^-1

3. Начальные концентрации компонентов

С_А(0) = 0,7 моль/л

С_В(0) = 0

С_С(0) = 0

С_Д(0) = 0

4. Продолжительность реакции 15 секунд

5. Метод численного решения – Эйлера.

Выполнение работы:

1. Система дифференциальных уравнений, представляющая кинетическую модель данной химической реакции:

2. Расчетные формулы метода Эйлера:

C_A^m+1= C_A^m+h(-K₁C_A^m+K₅C_B^m);

C_B^m+1= C_B^m+h(-K₂C_B^m+K₁C_A^m+K₄C_C^m-K₅C_B^m-K₃C_B^m);

C_C^m+1= C_C^m+h(K₂C_B^m-K₄C_C^m);

C_D^m+1= C_D^m+h·K₃C_B^m;

t^m+1= t^m+h.

3. Результаты численного решения системы дифференциальных уравнений на калькуляторе. h = 0.2; 5 шагов по времени.

Пример вычислений:

С_A1 = 0.7+0.2(-0.76*0.7+0.6*0)=0.5936;

C_B1 = 0+0.2(-0.2*0+0.76*0.7+0.45*0-0.6*0-0.5*0)=0.01064

C_C1 = 0+0.2(0.9*0-0.45*0)=0;

C_D1=0+0.2(0.5*0)=0;

t=0+0.2=0.2

Аналогично вычисляем С_А2, С_В2, С_С2, С_Д2,…С_А5, С_В5, С_С5, С_Д5.

Таблица результатов расчёта

№	t	CA	CB	CC	CD
0	0	0,7	0	0	0
1	0,2	0,5936	0,1064	0	0
2	0,4	0,5161	0,1241	0,0192	0,0106
3	0,6	0,4562	0,1726	0,0452	0,026
4	0,8	0,4076	0,177	0,0722	0,0433
5	1	0,3669	0,1747	0,0976	0,061

Графики кинетических кривых

4. ПАСКАЛЬ-программа решения системы дифференциальных уравнений

PROGRAM EILER;

USES CRT;

VAR

C, Z: ARRAY[1..4] OF REAL;

I: INTEGER;

T, TMAX: REAL;

BEGIN

CLRSCR;

WRITE (‘TMAX=’); READLN (TMAX);

T:=0;

WRITE (‘ВВЕДИТЕ С1, С2, С3, С4’);

READLN(C[1]);

READLN(C[2]);

READLN(C[3]);

READLN(C[4]);

WRITELN (‘T C1 C2 C3 C4’);

WHILE T<=TMAX DO

BEGIN

Z[1]:=C[1]+0.2*(-0.76*C[1]+0.6C[2]);

Z[2]:=C[2]+0.2*(-0.9*C[2]+0.76*C[1]+0.45*C[3]-0.6*C[2]-0.5*C[2]);

Z[3]:=C[3]+0.2*(0.9*C[2]-0.45*C[3]);

Z[4]:=C[4]+0.2*0.5*C[2];

T:=T+0.2;

WRITE (T:5:1);

FOR I:=1 TO 4 DO

WRITE (Z[1]:10:5);

WRITELN;

FOR I:=1 TO 4 DO

C[I]:=Z[I];

END;

READKEY;

END.