ルンゲ＝クッタ法

出典: フリー百科事典『ウィキペディア（Wikipedia）』

ルンゲ＝クッタ法（Runge-Kutta method）とは、数値解析において常微分方程式の近似解を求める一連の方法である。この技法は1900年頃に数学者C. RungeとM. W. Kuttaによって発展された。

[編集] 古典的な4次ルンゲ＝クッタ法

一般に用いられているルンゲ＝クッタ法は4次のルンゲ＝クッタ法（RK4）と呼ばれるものである。

初期値問題を次のように設定する。

$y' = f(t, y), \quad y(t_0) = y_0$

次に、RK4ではこの問題に対して次式を与える。

$y_{n+1} = y_n + {h \over 6} (k_1 + 2k_2 + 2k_3 + k_4)$

ここで、

$k_1 = f \left( t_n, y_n \right)$

$k_2 = f \left( t_n + {h \over 2}, y_n + {h \over 2} k_1 \right)$

$k_3 = f \left( t_n + {h \over 2}, y_n + {h \over 2} k_2 \right)$

$k_4 = f \left( t_n + h, y_n + hk_3 \right)$

である。

このようにして、次のステップの値（y_n+1）が、現在の値（y_n）とその間隔（h）の積から見積もられた勾配から求められる。その勾配は重み付け平均である。

k₁は初期値における勾配である。
k₂は区間の中央における勾配であり、勾配k₁を用いてt_n + h/2におけるyの値をオイラー法により決定したものである。

k₃は区間の中央における勾配を再計算したものであり、k₂の値から決められたyの値を用いる。
k₄は区間の最後における勾配であり、k₃の値から決められたyの値を用いる。

これら4つの平均を取るには、中央の勾配に対して大きな重み付けを与える。

$\mbox{slope} = \frac{k_1 + 2k_2 + 2k_3 + k_4}{6}.$

RK4は4次の方法であり、全体の推定誤差がh⁴のオーダーになることを意味する。

[編集] イクスプリシットルンゲ＝クッタ法

イクスプリシット（またはエクスプリシット）ルンゲ＝クッタ法は前述したRK4の一般化である。次式で与えられる。

$y_{n+1} = y_n + h\sum_{i=1}^s b_i k_i,$

ここで、

$k_1 = f(t_n, y_n), \,$

$k_2 = f(t_n+c_2h, y_n+a_{21}hk_1), \,$

$k_3 = f(t_n+c_3h, y_n+a_{31}hk_1+a_{32}hk_2), \,$

$\vdots$

$k_s = f(t_n+c_sh, y_n+a_{s1}hk_1+a_{s2}hk_2+\cdots+a_{s,s-1}hk_{s-1}).$

（注：これらの式は差分を含んでいるが、異なる書式による等価の定義である。）

個々の問題としては、整数s（段数）と、係数a_ij (for 1 ≤ j < i ≤ s), b_i (for i = 1, 2, ..., s) および c_i (for i = 2, 3, ..., s)を与える必要がある。それらのデータはRunge-Kutta tableauとして、記憶を助ける配置として知られている。:

0
$c 2$	$a 21$
$c 3$	$a 31$	$a 32$
$\vdots$	$\vdots$		$\ddots$
$c s$	$a s 1$	$a s 2$	$\cdots$	$a s, s - 1$
	$b 1$	$b 2$	$\cdots$	$b s - 1$	$b s$

ルンゲ＝クッタ法は次のように一貫したものである、

$\sum_{j=1}^{i-1} a_{ij} = c_i\ \mathrm{for}\ i=2, \ldots, s.$

私達があるオーダーpを持つ方法を必要とする場合、切断誤差が O(h^p+1)となる付随する条件がある。切断誤差自体の定義からこれらを得ることができる。例えば、2次の方法はb₁ + b₂ = 1, b₂c₂ = 1/2, and b₂a₂₁ = 1/2であれば2次のオーダーを持つ。

[編集] 例

RK4をこの枠組みに当てはめると、その tableauは:

0
1/2	1/2
1/2	0	1/2
1	0	0	1
	1/6	1/3	1/3	1/6

ところで、最も単純なルンゲ＝クッタ法は（前方）オイラー法であり、公式 $y n + 1 = y n + h f (t n, y n)$ で与えられる。これは1段のイクスプリシット・ルンゲ＝クッタ法から構成される。対応するtableauは:

	0
		1

2次の方法による例は中間値法により与えられる。

$y_{n+1} = y_n + hf\left(t_n+\frac{h}{2},y_n+\frac{h}{2}f(t_n, y_n)\right).$

対応するtableauは:

0
1/2	1/2
	0	1

[編集] 用法

2段のイクスプリシット・ルンゲ＝クッタ法の例を示す。

0
2/3	2/3
	1/4	3/4

次の初期値問題を解くものとする。

$y' = (\tan{y})+1,\quad y(1)=1,\ t\in [1, 1.1]$

ステップの値をh=0.025.とする。

tableauは、方法の定義のもとで同等な対応をする方程式をもたらす。:

$u_1 = y_n \,$

$u_2 = y_n + 2/3hf(t_n, u_1) \,$

$y_{n+1} = y_n + h(1/4f(t_n,u_1)+3/4f(t_n+2/3h,u_2))\,$

$t 0 = 1$
	$y 0 = 1$
$t 1 = 1.025$
	$u 1 = y 0 = 1$	$f (t 0, u 1) = 2.557407725$	$u 2 = y 0 + 2 / 3 h f (t 0, u 1) = 1.042623462$
	$y 1 = y 0 + h (1 / 4 f (t 0, u 1) + 3 / 4 f (t 0 + 2 / 3 h, u 2) = 1.066869388$
$t 2 = 1.05$
	$u 1 = y 1 = 1.066869388$	$f (t 1, u 1) = 2.813524695$	$u 2 = y 1 + 2 / 3 h f (t 1, u 1) = 1.113761467$
	$y 2 = y 1 + h (1 / 4 f (t 1, u 1) + 3 / 4 f (t 1 + 2 / 3 h, u 2) = 1.141332181$
$t 3 = 1.075$
	$u 1 = y 2 = 1.141332181$	$f (t 2, u 1) = 3.183536647$	$u 2 = y 2 + 2 / 3 h f (t 2, u 1) = 1.194391125$
	$y 3 = y 2 + h (1 / 4 f (t 2, u 1) + 3 / 4 f (t 2 + 2 / 3 h, u 2) = 1.227417567$
$t 4 = 1.1$
	$u 1 = y 3 = 1.227417567$	$f (t 3, u 3) = 3.796866512$	$u 2 = y 3 + 2 / 3 h f (t 3, u 1) = 1.290698676$
	$y 4 = y 3 + h (1 / 4 f (t 3, u 1) + 3 / 4 f (t 3 + 2 / 3 h, u 2) = 1.335079087$