ルンゲ・クッタ法

ルンゲ・クッタ法を用いて、微分方程式の初期値問題を解く。

微分方程式と初期条件が次式で与えられたとする。




このとき、4次のルンゲ・クッタ法により次式が与えられる。


ここで h は刻み幅、










である。初期条件x₀，y₀からx₁，y₁が求まり、x₁，y₁からx₂，y₂が求まる・・・というふうに計算を繰り返して数値解を得る。

例えば、微分方程式 f(x,y) = y-3x と、初期条件 y(0)=1 が与えられたとすると、

#coding: utf-8

def f(x, y):
    return (y-3*x)

#初期条件
x=0.
y=1.

#増加幅
h=0.1

while x <= 1:
    print x,y

    k1 = f(x, y)
    k2 = f(x+h/2, y+h/2*k1)
    k3 = f(x+h/2, y+h/2*k2)
    k4 = f(x+h, y+h*k3)
    k  = (k1+2*k2+2*k3+k4)/6
    y = y+k*h
    x+=h

計算すると、

x   y
0.0 1.0
0.1 1.08965833333
0.2 1.1571948583
0.3 1.20028300587
0.4 1.21635151984
0.5 1.20255872281
0.6 1.15576407582
0.7 1.07249674681
0.8 0.948920873415
0.9 0.78079717244
1.0 0.56344051173

と求まる。

ちなみに解析解は y=3x-2e^x+3 である。グラフにプロットするとこんな感じ。良い精度で近似できていることがわかる。