4.3.2. Метод Рунге-Кутта

Дальнейшее улучшение точности решения ОДУ первого порядка возможно за счет увеличения точности приближенного вычисления интеграла в выражении.

Воспользовавшись формулой Симпсона, можно получить еще более точную формулу для решения ОДУ первого порядка – широко используемого в вычислительной практике метода Рунге-Кутта.

В формуле Симпсона для приближенного вычисления определенного интеграла используются значения подинтегрального выражения в трех точках. В интеграле их всего две, поэтому введем дополнительную точку в середине отрезка [x_i+1 x_i].

,

тогда можно переписать так:

Полученное выражение является неявным, так как в правой части содержатся еще не определенные значения функции y_i+h/2 и y_i+1. Чтобы воспользоваться этой формулой, надо использовать некоторое приближение для вычисления этих значений

При использовании различных методов приближенного вычисления этих величин, получаются выражения для методов Рунге-Кутта различного порядка точности.

Алгоритм Рунге-Кутта третьего порядка (погрешность порядка h³):

где

Алгоритм Рунге-Кутта четвертого порядка (погрешность порядка h⁴):

где

Пример 4.3.2 Рассчитать кинетику гомогенной химической реакции, используя метод Рунге-Кутта.

VAR ca, cb, cc, cd, ca0, cb0, cc0, cd0, k1, k2, k3, k4, t0, tk, h:real;

{Ca, Cb, Cc, Cd - Концентрации веществ до реакции }

{k1, k2, k3, k4 - Скорости для каждой из реакций}

{t0, tk - Начальное и конечное время реагирования}

{h - шаг}

x, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, q13, q14, q15, q16:real;

{Qn - Коэффициенты Рунге-Кутты}

f1, f2:text;

BEGIN

assign(f1,'runge.in');{Файл с данными}

assign(f2,'runge.out'); {Файл с результатами }

reset(f1); read(f1,ca0,cb0,cc0,cd0,k1,k2,k3,k4,t0,tk,h);

close(f1);

cc:=cc0;ca:=ca0;cb:=cb0;cd:=cd0; x:=t0;

rewrite(f2);

Writeln(f2,' Метод Рунге - Кутта');

writeln(f2,' Время Концентрация ');

writeln(f2,'A B С D ');

REPEAT

q1:=h*((-1)*(k1+k4)*ca0+(k3*cd0));

q2:=h*((-1)*(k1+k4)*(ca0+q1/2)+(k3*(cd0+q1/2)));

q3:=h*((-1)*(k1+k4)*(ca0+q2/2)+(k3*(cd0+q2/2)));

q4:=h*((-1)*(k1+k4)*(ca0+q3)+(k3*(cd0+q3)));

q5:=h*((k1*ca0)-(k2*cb0));

q6:=h*((k1*(ca0+q5/2))-k2*(cb0+q5/2)));

q7:=h*((k1*(ca0+q6/2))-(k2*(cb0+q6/2)));

q8:=h*((k1*(ca0+q7))-(k2*(cb0+q7)));

q9:=h*(k2*cb0);

q10:=h*(k2*(cb0+q9/2));

q11:=h*(k2*(cb0+q10/2));

q12:=h*(k2*(cb0+q10));

q13:=h*((k4*ca0)-(k3*cb0));

q14:=h*((k4*(ca0+q13/2))-(k3*(cb0+q13/2)));

q15:=h*((k4*(ca0+q14/2))-(k3*(cb0+q14/2)));

q16:=h*((k4*(ca0+q15))-(k3*(cb0+q15)));

ca:=ca0+((q1+q2+q3+q4)/6);

cb:=cb0+((q5+q6+q7+q8)/6);

cc:=cc0+((q9+q10+q11+q12)/6);

cd:=cd0+((q13+q14+q15+q16)/6);

writeln(f2,' ',x:8:4,' ',ca:8:4,' ',cb:8:4,' ',cc:8:4,' ',cd:8:4,' ');

ca0:=ca; cb0:=cb; cc0:=cc; cd0:=cd;

x:=x+h;

UNTIL x>=tk;

writeln(f2,' ',x:8:4,' ',ca:8:4,' ',cb:8:4,' ',cc:8:4,' ',cd:8:4,' ');

close(f2);

END.