【数値解析】ルンゲクッタ法

分かっているもの

• $\displaystyle \frac{dy}{dx} = f(x, y)$ …dy/dx
• $y(0)$ …初期値
• $h$ …幅

求め方

1. $\displaystyle k_1 = f(x[n], y[n])$
2. $\displaystyle k_2 = f(x[n]+\frac{h}{2}, y[n]+\frac{h}{2}k_1)$
3. $\displaystyle k_3 = f(x[n] + \frac{h}{2}, y[n] + \frac{h}{2}k_2)$
4. $\displaystyle k_4 = f(x[n]+h, y[n]+hk_3)$
5. $\displaystyle y[n+1] = \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4)$

($k$は$n$が加算するごとに毎回計算する。)

ここで、$x[n] = hn$

サンプルプログラム(Ruby)

• $\displaystyle \frac{dy}{dx} = \cos{x}$
• $y(0) = 0$
• $h = 0.01$

を求めたい場合。

#!/usr/bin/ruby

def f t, y
    # 関数の微分
    Math::cos(t)
end

# 幅
h = 0.01

max = 1001
y = Array.new(max){nil}
t = Array.new(max){|i| i*h }
# 初期値
y[0] = 0.0

max.times do |n|
    k1 = f(t[n], y[n])
    k2 = f(t[n]+h/2, y[n]+h/2*k1)
    k3 = f(t[n]+h/2, y[n]+h/2*k2)
    k4 = f(t[n]+h, y[n]+h*k3)
    y[n+1] = y[n] + (k1 + 2*k2 + 2*k3 + k4)*h/6
end

max.times do |n|
    puts "#{t[n]} #{y[n]}"
end