[微分方程式] ルンゲクッタ法

4次のルンゲクッタ法

初期条件 ($ x_0, y_0 $) 付き微分方程式 : $ \frac{dy}{dx} = f(x,y) $ の精度よい解法

原理

まずは差分化 :

$$ \frac{y_{i+1}-y_i}{\Delta x} = f(x_i, y_i) $$

ここで、オイラー法では $ y_{i+1} = y_i + f \Delta x $ だった。

ルンゲクッタ法では次のようにする。

$$ y_{i+1} = y_i + \frac{1}{6} (k_1 + 2k_2 + 2k_3 + k_4) $$

ここで

$$ k_1 = \Delta x f(x_i, y_i) $$ $$ k_2 = \Delta x f(x_i + \frac{\Delta x}{2}, y_i+\frac{k_1}{2}) $$ $$ k_3 = \Delta x f(x_i + \frac{\Delta x}{2}, y_i+\frac{k_2}{2}) $$ $$ k_4 = \Delta x f(x_i+\Delta x, y_i+k_3) $$

詳しい原理は Wiki様 を参考に。

例題

微分方程式 : $ \frac{dy}{dx} = 2x - 2xy $ を初期条件 $x_0=0$, $y_0=3$ で解け。

計算範囲は $x=0-3$ として、$\Delta x$ が0.1, 0.01 の場合で試して理論解との誤差がどうなるか調べよ。

理論解は $y_{true} =1+2e^{-x^2}$ である。

結果

PROGRAM RungeKutta
IMPLICIT NONE
    REAL(8) :: k1, k2, k3, k4, x, y, dx=1d-1    ! 1d-2
    INTEGER :: i, Nmax=30     ! 300

    x = 0.d0	! initial value
    y = 3.d0

    DO i = 1, Nmax
    	WRITE(10,*) x, y, 1.d0+2.d0*exp(-x**2)
    	x = x+dx
    	k1 = dx * f(x        , y        )	! ルンゲクッタの成分
    	k2 = dx * f(x+dx/2.d0, y+k1/2.d0)
    	k3 = dx * f(x+dx/2.d0, y+k2/2.d0)
    	k4 = dx * f(x+dx     , y+k3     )
    	y = y + (k1+2.d0*k2+2.d0*k3+k4)/6.d0
    END DO

CONTAINS
FUNCTION f(x,y)		! 微分方程式の右辺
REAL(8) :: x, y, f
	f = 2.d0*x - 2.d0*x*y
END
END PROGRAM RungeKutta

係数は足し合わせたあとに dx をかけたほうがよかったかも。

$\Delta x = 0.1$ (上)と$\Delta x = 0.01$ (下)