Реализация некоторых численных методов
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]
![y(x + h) = y(x) + \int\limits_x^{x + h} {f\,[t,y(t)]dt},](/sites/default/files/tex_cache/7c95a7b85a6bd92e2861bb535ea3ddca.png)
![\Delta y=h\int\limits_{0}^{1}{f[x}+\alpha
\,h,y(x+\alpha \,h)]d\alpha \,](/sites/default/files/tex_cache/38101b3b01a678117748ec9c54359abd.png)