ルンゲ＝クッタ法

ルンゲ＝クッタ法（Runge-Kutta method）とは、数値解析において常微分方程式の近似解を求める一連の方法である。この技法は1900年頃に数学者カール・ルンゲとマルティン・ヴィルヘルム・クッタ (英語: Martin Wilhelm 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）と推定された勾配の積を加えたものである。その勾配は次の4つの勾配の重み付け平均で求める。

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次の方法である。すなわち、厳密解と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_{s1}$	$a_{s2}$	$\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 + hf(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

誤差と収束性[編集]

関数 $f$ が $p$ 階連続偏微分可能であるならば、 $p$ 次RK法から得られた近似解と厳密解の局所誤差にはステップ幅をパラメータとした誤差限界が存在する。

しかし、この誤差限界は $p$ が高次になると意味を成さなくなる。そのため、ステップ幅の選定には様々な補外法が提案されている。

用法[編集]

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/3hf(t_0,u_1) = 1.042623462$
	$y_1=y_0+h(1/4 f(t_0, u_1) + 3/4 f(t_0+2/3h,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/3hf(t_1,u_1) = 1.113761467$
	$y_2=y_1+h(1/4 f(t_1, u_1) + 3/4 f(t_1+2/3h,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/3hf(t_2,u_1) = 1.194391125$
	$y_3=y_2+h(1/4 f(t_2, u_1) + 3/4 f(t_2+2/3h,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/3hf(t_3,u_1) = 1.290698676$
	$y_4=y_3+h(1/4 f(t_3, u_1) + 3/4 f(t_3+2/3h,u_2)=1.335079087$