ルンゲクッタ法（メモ）

matlab

常微分方程式\reverse\opaque\frac{dy}{dx}=f(x, y)，初期値\reverse\opaque x_0, y_0が与えられたとき，ルンゲクッタ法で計算することができる．
\reverse\opaque\frac{dy}{dx} = f(x_{n-1}, y_{n-1})
\reverse\opaque k_1=\Delta x f(x_{n-1}, y_{n-1})
\reverse\opaque k_2=\Delta x f(x_{n-1}+\frac{\Delta x}{2}, y_{n-1}+\frac{k_1}{2})
\reverse\opaque k_3=\Delta x f(x_{n-1}+\frac{\Delta x}{2}, y_{n-1}+\frac{k_2}{2})
\reverse\opaque k_4=\Delta x f(x_{n-1}+\Delta x, y_{n-1}+k_3)
\reverse\opaque y_n = y_{n-1} +\left( k_1+2k_2+2k_3+k_4\right)/6

\reverse\opaque\Delta xはxの刻み幅で計算したい値までを繰り返して計算する．n階の微分方程式もn個の１階微分方程式に分けることでルンゲクッタ法で計算できる．
単振動を考えるとき，\reverse\opaque \frac{d^2x}{dt^2}=-xという式が与えられたとする．この式は
\reverse\opaque \frac{dv}{dt}=-x=f_1(t_{n-1},v_{n-1},x_{n-1})
\reverse\opaque \frac{dx}{dt}=v=f_2(t_{n-1},v_{n-1},x_{n-1})
という2式で表せ，初期値が与えら得たとき，
\reverse\opaque v_n = v_{n-1} +\left( k_1+2k_2+2k_3+k_4\right)/6
\reverse\opaque x_n = x_{n-1} +\left( l1+2l_2+2l_3+l_4\right)/6
という形で求める．したがって
\reverse\opaque k_1=\Delta t f_1(t_{n-1},v_{n-1}, x_{n-1})
\reverse\opaque l_1=\Delta t f_2(t_{n-1},v_{n-1}, x_{n-1})
\reverse\opaque k_2=\Delta t f_1(t_{n-1}+\frac{\Delta t}{2},v_{n-1}+\frac{k_1}{2}, x_{n-1}+\frac{l_1}{2})
\reverse\opaque l_2=\Delta t f_2(t_{n-1}+\frac{\Delta t}{2},v_{n-1}+\frac{k_1}{2}, x_{n-1}+\frac{l_1}{2})
\reverse\opaque k_3=\Delta t f_1(t_{n-1}+\frac{\Delta t}{2},v_{n-1}+\frac{k2}{2}, x_{n-1}+\frac{l_2}{2})
\reverse\opaque l_3=\Delta t f_2(t_{n-1}+\frac{\Delta t}{2},v_{n-1}+\frac{k2}{2}, x_{n-1}+\frac{l_2}{2})
\reverse\opaque k_4=\Delta t f_1(t_{n-1}+\Delta t,v_{n-1}+k_3, x_{n-1}+l_3)
\reverse\opaque l_4=\Delta t f_2(t_{n-1}+\Delta t,v_{n-1}+k_3 t, x_{n-1}+l_3)
から計算することができる．
次のようなカオスを表す3つの微分方程式の場合でも同じ．
\reverse\opaque \frac{dx}{dt}=y=f_1(t_{n-1},x_{n-1},y_{n-1},z_{n-1})
\reverse\opaque \frac{dy}{dt}=-x+yz=f_2(t_{n-1},x_{n-1},y_{n-1},z_{n-1})
\reverse\opaque \frac{dz}{dt}=1-y^2=f_3(t_{n-1},x_{n-1},y_{n-1},z_{n-1})
\reverse\opaque \left( x_0,y_0,z_0 \right) = \left( 1, -0.01,0.2\right)

matlabの場合，ソルバーがあるのでそれを利用する．上のカオスの例を解くには，まず，1階の微分方程式を入力したmファイルを作る．

function  dy = chaos(t,y);
dy=zeros(3,1);
dy(1) = y(2);
dy(2) = -y(1)+ y(2)*y(3);
dy(3) = 1- y(2)^2;

次に，ソルバーを利用したコードを書く．

options = odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4 1e-5]);
[t,y] = ode45('chaos',[0 100], [-1 -0.01 0.2], options);

結果は↓．



参考