метод Руеге-Кутты для автономной системы ОДУ

[tomsksmile]

Народ, поможите
как будет выглядеть метод Рунге-Кутты 4 порядка для решения автономной системы ОДУ? Или объясните, плиз, как перейти от неавтономной системы к автономной и как при этом модифицируется алгоритм вычислений (в котором явно абсциссы не присутствуют).

[Demo_S]

у-у:)
а чем отличается автономная система от неавтономной?

сам метод я  делал, только вот не знаю, для какой системы...

[shadeofgray]

Цитата tomsksmile, 03.07.03, 07:10:44

как будет выглядеть метод Рунге-Кутты 4 порядка для решения автономной системы ОДУ?

Так же, как и для неавтономной, только правая часть (функция F) теперь не зависит от х. Код приведен ниже (для автономной и неавтономной). Фактически, нет необходимости менять что-то в алгоритме вычисления, достаточно ввести функцию F, принимающую кроме аргумента y ещё и х, но зависящую только от y. Т.е. х будет фиктивным параметром.

Цитата tomsksmile, 03.07.03, 07:10:44

Или объясните, плиз, как перейти от неавтономной системы к автономной и как при этом модифицируется алгоритм вычислений (в котором явно абсциссы не присутствуют).

Есть один способ. Пусть система имеет n неизвестных y₁..y_n. Введем новую переменную y_n+1, такую, что y_n+1'=1. Заменим в правой части системы х на y_n+1. Теперь у системы увеличилась размерность на 1, но она явно не зависит от времени - хотя фактически, ничего не изменилось, кроме способа записи.
В общем случае, другого способа, не повышающего размерность, не существует.


Вот пример на псевдокоде. Перевести вручную на что хочешь не составит труда. Процедура решает систему из n уравнений с n неизвестными вида y_i'=F_i(x,y), где производная берется по x. Шаг интегрирования (делается один шаг) равен h. Начальное состояние находится в y, там же находится состояние на следующем шаге.

В псевдокоде используется вызов процедуры SetBounds, который устанавливает размер массивов так, чтобы в них были элементы с номерами от 1 до n.

Это код для решения неавтономной системы

<br>function F(I:Integer; X:Real; Y:array of Real):Real;<br>begin<br>...<br>end;<br><br>procedure SystemRungeKuttStep(<br>    const   x   :   Real;<br>    const   h   :   Real;<br>    const   n   :   Integer;<br>    var     y   :   array of Real);<br>var<br>   I        :   Integer;<br>   yt       :   array of Real;<br>   k1       :   array of Real;<br>   k2       :   array of Real;<br>   k3       :   array of Real;<br>   k4       :   array of Real;<br>begin<br>    SetBounds(yt,   [1,N]);<br>    SetBounds(k1,   [1,N]);<br>    SetBounds(k2,   [1,N]);<br>    SetBounds(k3,   [1,N]);<br>    SetBounds(k4,   [1,N]);<br><br>    for I:=1 to N do<br>        k1[i] := h*F(i,x,y);<br><br>    for i := 1 to n do<br>        yt[i] := y[i]+0.5*k1[i];<br>    for I:=1 to N do<br>        k2[i] := h*F(i,x+h*0.5,yt);<br><br>    for i := 1 to n do<br>        yt[i] := y[i]+ 0.5*k2[i];<br>    for I:=1 to N do<br>        k3[i] := h*F(i,x+h*0.5,yt);<br><br>    for i := 1 to n do<br>        yt[i] := y[i] + k3[i];<br>    for I:=1 to N do<br>        k4[i] := h*F(i,x+h,yt);<br><br>    for i := 1 to n do<br>        y[i] := y[i] + (k1[i]+2.0*k2[i]+2.0*k3[i]+k4[i])/6;<br>end;<br>

Это код для решения автономной системы

<br>function F(I:Integer;Y:array of Real):Real;<br>begin<br>...<br>end;<br><br>procedure SystemRungeKuttStep(<br>    const   h   :   Real;<br>    const   n   :   Integer;<br>    var     y   :   array of Real);<br>var<br>   I        :   Integer;<br>   yt       :   array of Real;<br>   k1       :   array of Real;<br>   k2       :   array of Real;<br>   k3       :   array of Real;<br>   k4       :   array of Real;<br>begin<br>    SetBounds(yt,   [1,N]);<br>    SetBounds(k1,   [1,N]);<br>    SetBounds(k2,   [1,N]);<br>    SetBounds(k3,   [1,N]);<br>    SetBounds(k4,   [1,N]);<br><br>    for I:=1 to N do<br>        k1[i] := h*F(i,y);<br><br>    for i := 1 to n do<br>        yt[i] := y[i]+0.5*k1[i];<br>    for I:=1 to N do<br>        k2[i] := h*F(i,yt);<br><br>    for i := 1 to n do<br>        yt[i] := y[i]+ 0.5*k2[i];<br>    for I:=1 to N do<br>        k3[i] := h*F(i,yt);<br><br>    for i := 1 to n do<br>        yt[i] := y[i] + k3[i];<br>    for I:=1 to N do<br>        k4[i] := h*F(i,yt);<br><br>    for i := 1 to n do<br>        y[i] := y[i] + (k1[i]+2.0*k2[i]+2.0*k3[i]+k4[i])/6;<br>end;<br>

[shadeofgray]

Цитата Demo_S, 03.07.03, 10:07:02

а чем отличается автономная система от неавтономной?

в автономной системе правая часть не содержит времени (в явной форме).