3.4. Метод рунге-кутта 4-го порядка
Этот метод один из самых распространённых методов интегрирования дифференциальных уравнений.
Для одиночного дифференциального уравнения y'= f(x,y), y(x0)=y0 расчётные формулы имеют следующий вид:
ym+1= ym+h/6(k1+2k2+2k3+k4),
где k1= f(xm,ym);
k2= f(xm+,ym+);
k3= f(xm+,ym+);
k4= f(xm+h,ym+hk3);
xm+1= xm+h.
Для системы дифференциальных уравнений
yi'= fi(x,y1,y2,…yk), yi(x0)=yi0, i=1,2,..k,
расчетные формулы запишутся следующим образом:
yim+1= yim+(k1i+2k2i+2k3i+k4i), i=1,2,…k,
где k1i= fi(xm,y1m,y2m,…ykm);
k2i= fi(xm+h/2,y1m+k11,y2m+k12,…ykm+k1k);
k3i= fi(xm+,y1m+k21,y2m+k22,…ykm+k2k);
k4i= fi(xm+h,y1m+hk31,y2m+hk32,…ykm+hk3k).
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.
МОДЕЛИРОВАНИЕ КИНЕТИКИ ХИМИЧЕСКОЙ РЕАКЦИИ
Пример.
Задано:
Схема механизма химической реакции
К1 К2
А В С
К5К3К4
Д
Константы скоростей отдельных стадий реакции
К1 = 0,76 с-1
К2 = 0,9 с-1
К3 = 0,5 с-1
К4 = 0,45 с-1
К5 = 0,6 с-1
Начальные концентрации компонентов
СА(0) = 0,7 моль/л
СВ(0) = 0
СС(0) = 0
СД(0) = 0
Продолжительность реакции 15 секунд
Метод численного решения – Эйлера.
Выполнение работы:
Система дифференциальных уравнений, представляющая кинетическую модель данной химической реакции:
Расчетные формулы метода Эйлера:
CAm+1= CAm+h(-K1CAm+K5CBm);
CBm+1= CBm+h(-K2CBm+K1CAm+K4CCm-K5CBm-K3CBm);
CCm+1= CCm+h(K2CBm-K4CCm);
CDm+1= CDm+h·K3CBm;
tm+1= tm+h.
Результаты численного решения системы дифференциальных уравнений на калькуляторе. h = 0.2; 5 шагов по времени.
Пример вычислений:
СA1 = 0.7+0.2(-0.76*0.7+0.6*0)=0.5936;
CB1 = 0+0.2(-0.2*0+0.76*0.7+0.45*0-0.6*0-0.5*0)=0.01064
CC1 = 0+0.2(0.9*0-0.45*0)=0;
CD1=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.