Реализация некоторых численных методов
8.5.2.2 Метод Рунге-Кутта
Изложим идею метода на примере задачи Коши:




![]() |
( 8.10) |

Для удобства записи выражения (8.11) используем обозначение и замену переменной интегрирования
. Окончательно получим:
![]() |
( 8.12) |
В зависимости от способа вычисления интеграла в выражении (8.12) получают различные методы численного интегрирования обыкновенных дифференциальных уравнений.
Рассмотрим линейную комбинацию величин , которая будет являться аналогом квадратурной суммы и позволит вычислить приближенное значение приращения
:


Метод четвертого порядка для , являющийся аналогом широко известной в литературе четырехточечной квадратурной формулы "трех восьмых", имеет вид


Особо широко известно другое вычислительное правило типа Рунге-Кутта четвертого порядка точности:



Функция Maxima, реализующая метод Рунге-Кутта 4-го порядка, приведена в следующем примере (с печатью промежуточных результатов):
(%i1) rk4(rp,fun,y0,x0,xend,h):=block( [OK,n,h1,_x,_y,_k1,_k2,_k3,_k4,rez], _x:x0, _y:y0, rez:[_y], OK:-1, h1:h, n:length(_y), while OK<0 do ( if (_x+h1>=xend) then (h1:xend-_x, OK:1), _k1:makelist(float(h1*subst([fun[i]= float(_y[i]),x=float(_x)],rp[i])),i,1,n), _k2:makelist(float(h1*subst([fun[i]= float(_y[i]+_k1[i]/2),x=float(_x+h1/2)], rp[i])),i,1,n), _k3:makelist(float(h1*subst([fun[i]= float(_y[i]+_k2[i]/2),x=float(_x+h1/2)], rp[i])),i,1,n), _k4:makelist(float(h1*subst([fun[i]= float(_y[i]+_k3[i]),x=float(_x+h1)],rp[i])), i,1,n), _y1:makelist(float(_y[i]+ (_k1[i]+2*_k2[i]+2*_k3[i]+_k4[i])/6),i,1,n), rez:append(rez,[_y1]), print("x= ",_x, " y= ",_y), _x:_x+h1, _y:_y1 ), rez )$
Пример обращения к функции представлен следующей последовательностью команд (решалась та же система, что и при тестировании метода Эйлера):
(%i2) rk4([-2*y,-5*v,3*x],[y,v,z],[1,1,1],0,1,0.1); x= 0 y= [1,1,1] x= 0.1 y= [0.81873333333333,0.60677083333333,1.015] x= 0.2 y= [0.67032427111111,0.36817084418403,1.06] x= 0.3 y= [0.54881682490104,0.22339532993458,1.135] x= 0.4 y= [0.44933462844064,0.13554977050718,1.24] x= 0.5 y= [0.3678852381253,0.082247647208783,1.375] x= 0.6 y= [0.30119990729446,0.04990547343658,1.54] x= 0.7 y= [0.24660240409888,0.030281185705008,1.735] x= 0.8 y= [0.20190160831589,0.018373740284549,1.96] x= 0.9 y= [0.16530357678183,0.011148649703906,2.215] x= 1.0 y= [0.13533954843051,0.0067646754713805,2.5]