ルンゲ＝クッタ法

数値解析においてルンゲ＝クッタ法（英: Runge–Kutta method）とは、初期値問題に対して近似解を与える常微分方程式の数値解法に対する総称である。この技法は1900年頃に数学者カール・ルンゲとマルティン・クッタによって発展を見た。

古典的ルンゲ＝クッタ法

一連のルンゲ＝クッタ公式の中で最も広く知られているのが、古典的ルンゲ=クッタ法 (RK4、もしくは単に狭義の ルンゲ=クッタ法、英: the (classical) Runge–Kutta method) などと呼ばれる4次の公式である。

次の初期値問題を考える。

${\displaystyle y'=f(t,y),\quad y(t_{0})=y_{0}.}$

但し、y(t) が近似的に求めたい未知関数であり、その t における勾配は f(t, y) によって t 及び y(t) の関数として与えられている。時刻 t₀ における初期値は y₀ で与えられている。

今、時刻 t_n における値 y_n = y(t_n) が既知のとき、十分に小さなステップ幅 h に対して y_n+1, t_n+1 を以下の式で与えると、y_n+1 は y(t_n+1) の 4次精度の近似になっている。このステップを逐次的に繰り返すことによって、初期値 y₀ から任意の時刻 t_n における近似値 y_n が求められる。

${\displaystyle {\begin{aligned}y_{n+1}&=y_{n}+{h \over 6}(k_{1}+2k_{2}+2k_{3}+k_{4}),\\t_{n+1}&=t_{n}+h.\end{aligned}}}$

ここで、

${\displaystyle {\begin{aligned}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)\end{aligned}}}$

である^[1][2]。次の値 (y_n+1) は、現在の値 (y_n) に増分を加えたものであり、増分は勾配の推定値に間隔 h を乗じたものになっている。勾配の推定値は、k₁, ..., k₄ の4つの勾配の重み付け平均で求める。k₁, ..., k₄ のそれぞれの勾配は、特定の (t, y) に対する f によって与えられ、以下のように解釈できる。

• k₁ は区間の最初 t_n における勾配である (オイラー法で用いる勾配に一致する)。
• k₂ は区間の中央 t_n + h/2 における勾配の近似値である (中点法で用いる勾配)。計算に用いる中央の y の値は、初期位置の勾配 k₁ を用いてオイラー法で推定する。
• k₃ は区間の中央における勾配のもう一つの近似値である。中央の y を k₂ の値から推定して用いる。
• k₄ は区間の最後 t_n + h における勾配の近似値であり、k₃ の値から推定された最後の点の y の値を用いる。

重み付き平均では、中央の勾配に対して大きな重みを用いる。シンプソン則を用いた平均と同等の形になる^[3]。

${\displaystyle {\mbox{slope}}={\frac {k_{1}+2k_{2}+2k_{3}+k_{4}}{6}}.}$

RK4は4次の方法である。厳密解とRK4のテイラー展開が4次の項まで一致し、1ステップの推定誤差は ${\displaystyle O(h^{5})}$ のオーダーになる。目的の時刻の y を求めるのに必要なステップ数は ${\displaystyle O(h^{-1})}$ になるので、全体の推定誤差は ${\displaystyle O(h^{4})}$ になる。

ルンゲ゠クッタ法

前述の RK4 の一般化として、以下の形式を持つ s 段のルンゲ゠クッタ法を構成することができる。整数 s をそのルンゲ゠クッタ法の段数 (stage) という。

${\displaystyle y_{n+1}=y_{n}+h\sum _{i=1}^{s}b_{i}k_{i},}$

但し^[4]、

${\displaystyle k_{i}=f\left(t_{n}+hc_{i},y_{n}+h\sum _{j=1}^{s}a_{ij}k_{j}\right),\quad i=1,\ldots ,s.}$

(文献によって、等価であるが上と異なる定義の仕方をしているものがあることに注意する)^[5]。

具体的なルンゲ゠クッタ公式は、解ができるだけ高い精度を持つように適切に選ばれた係数 a_ij (ルンゲ゠クッタ行列)、b_i (重み)、c_j (節点) で指定される (${\displaystyle i,j\leq s}$)^[6]。特に ${\displaystyle i\leq j}$ に対して ${\displaystyle a_{ij}=0}$ を満たす方法が広く用いられ、総称して 陽的ルンゲ゠クッタ法 (ERK、英: explicit Runge–Kutta methods) と呼ぶ。そうでないものを 陰的ルンゲ゠クッタ法 (IRK、英: implicit Runge–Kutta methods) と呼ぶ。

近似値 y_n を y_n+1 から計算するときに発生する誤差の大きさが ${\displaystyle O(h^{p+1})}$ のとき、そのルンゲ゠クッタ公式は p 次精度を持つといい、p を次数 (または位数) と呼ぶ。p 次のルンゲ゠クッタ公式は、誤差の大きさの条件に誤差の表式を代入し、係数の条件を求めることによって得られる。例えば、2段の陽的方法が2次精度を持つための係数に対する条件は、${\displaystyle b_{1}+b_{2}=1}$, ${\displaystyle b_{2}c_{2}=1/2}$, かつ ${\displaystyle b_{2}a_{21}=1/2}$ である^[7]。この条件を満たす範囲内で様々な方法を考え得る。

これらの係数を分かりやすく表記する方法として、以下のような形式のブッチャー配列 (Butcher tableau) が知られている^[8]:

${\displaystyle {\begin{aligned}{\begin{array}{c|ccccc}c_{1}&a_{11}&a_{12}&\cdots &a_{1s}\\c_{2}&a_{21}&a_{22}&\cdots &a_{2s}\\\vdots &\vdots &\vdots &\ddots &\vdots \\c_{s}&a_{s1}&a_{s2}&\cdots &a_{ss}\\\hline &b_{1}&b_{2}&\cdots &b_{s}\end{array}}&={\begin{array}{c|c}\mathbf {c} &A\\\hline &\mathbf {b} ^{\mathrm {T} }\end{array}}\end{aligned}}}$

実際には、それぞれのルンゲ゠クッタ法について各要素に具体的な値を入れて用いる。陽的な方法ではルンゲ゠クッタ行列の上三角成分は常に 0 になるので表記を省略する。例えば RK4 はブッチャー配列を用いて以下のように表現される。

${\displaystyle {\begin{aligned}{\begin{array}{c|ccccc}0&\\1/2&1/2\\1/2&0&1/2\\1&0&0&1\\\hline &1/6&1/3&1/3&1/6\end{array}}\end{aligned}}}$

ルンゲ゠クッタ法が一貫しているためには、次の条件が満たされている必要がある。

${\displaystyle \sum _{j=1}^{s}a_{ij}=c_{i},\;i=1,\ldots ,s.}$

収束性

ルンゲ＝クッタ法は、数値積分における求積法 (quadrature) と深く繋がる。時刻 t_n での値から t_n+1 = t_n + h での値を求めるときの方程式は以下のように定める。

${\displaystyle y'=f(t,y),\;y(t_{n})=y_{n}.}$

求積法は、与えられた区間での定積分の値を被積分関数の値の線型結合として近似する方法である。簡単のために、区間を [0, 1] とする。よって求積法の公式は

${\displaystyle \int _{0}^{1}f(t)dt\approx \sum _{i=1}^{s}b_{i}f(c_{i})}$

となる。ここで、b_i と c_i は先に選ばれた定数であり、前述の重みと節点に対応する。上記式に対し、等式がすべての p − 1 次以下の多項式に成立するとき（すなわち、誤差が0のとき）、その求積法は p 次精度であり、p を次数と呼ぶ^[9]。節点が s 個の時、最大次数は 2s であり、その方法は s 次ガウス・ルジャンドル公式と呼ばれる。

そして上記の方程式を積分形式に変形し、求積法を用いると次の公式となる^[6]。

${\displaystyle y(t_{n+1})=y_{n}+\int _{y_{n}}^{y_{n+1}}f(t,y(t))dt=y_{n}+h\int _{0}^{1}y(t_{n}+h\tau ,y(t_{n}+h\tau ))d\tau =y_{n}+h\sum _{i=1}^{s}b_{i}f(t_{n}+c_{i}h,y(t_{n}+c_{i}h))}$

${\displaystyle k_{i}=f(t_{n}+c_{i}h,y(t_{n}+c_{i}h))}$ とおく。k_i あるいは y(t_n + c_i h) を適切に（線型結合として）近似することでルンゲ＝クッタ法の公式となる。その上、係数をテイラー展開より正しく選択すると、方法の収束性も求積法の収束性より保証される。しかし、局所誤差のオーダーや上界は、方法によって大きく異なるので、方法別に計算しなければならない。

陽的ルンゲ゠クッタ法

陽的ルンゲ゠クッタ法では ${\displaystyle a_{ij}=0}$ (${\displaystyle i\leq j}$) が満たされるので、ブッチャー配列は以下の形になる。

	0
	${\displaystyle c_{2}}$	${\displaystyle a_{21}}$
	${\displaystyle c_{3}}$	${\displaystyle a_{31}}$	${\displaystyle a_{32}}$
	${\displaystyle \vdots }$	${\displaystyle \vdots }$		${\displaystyle \ddots }$
	${\displaystyle c_{s}}$	${\displaystyle a_{s1}}$	${\displaystyle a_{s2}}$	${\displaystyle \cdots }$	${\displaystyle a_{s,s-1}}$
		${\displaystyle b_{1}}$	${\displaystyle b_{2}}$	${\displaystyle \cdots }$	${\displaystyle b_{s-1}}$	${\displaystyle b_{s}}$

この時、各勾配 k₁, ..., k_s を以下のように逐次的に求めることができる^[10]。

${\displaystyle {\begin{aligned}k_{1}&=f(t_{n},y_{n}),\\k_{2}&=f(t_{n}+c_{2}h,y_{n}+a_{21}hk_{1}),\\k_{3}&=f(t_{n}+c_{3}h,y_{n}+a_{31}hk_{1}+a_{32}hk_{2}),\\&\vdots \\k_{s}&=f(t_{n}+c_{s}h,y_{n}+a_{s1}hk_{1}+a_{s2}hk_{2}+\cdots +a_{s,s-1}hk_{s-1}).\end{aligned}}}$

(文献によって、等価であるが上と異なる定義の仕方をしているものがあることに注意する)^[5]。

陽的ルンゲ゠クッタ法が目的の精度 p になる係数を持つためには、十分に大きな段数 s が必要になる。少なくとも次数以上の段数が必要 (${\displaystyle s\geq p}$) であり、更に5次以上の場合には段数は次数よりも大きく取らなければならない (${\displaystyle s>p}$) ことが知られている^[11]。各次数 p に対する最低段数 s の一般式は未解決問題である。具体的な値が幾つか知られている^[12]:

${\displaystyle {\begin{array}{c|cccccccc}p&1&2&3&4&5&6&7&8\\\hline \min s&1&2&3&4&6&7&9&11\end{array}}}$

1段1次の公式

最も単純なルンゲ゠クッタ法が 1段1次 の方法であり、一意に定まる。(前進) オイラー法と等価になる。

${\displaystyle y_{n+1}=y_{n}+hf(t_{n},y_{n}).}$

以下のブッチャー配列で表現される。

	0
		1

2段2次の公式

2段2次の方法は1パラメータの自由度 α を持ち、以下のブッチャー配列で与えられる。

	0
	${\displaystyle \alpha }$	${\displaystyle \alpha }$
		${\displaystyle (1-{\tfrac {1}{2\alpha }})}$	${\displaystyle {\tfrac {1}{2\alpha }}}$

対応する公式は^[14]

${\displaystyle y_{n+1}=y_{n}+h{\bigl (}(1-{\tfrac {1}{2\alpha }})f(t_{n},y_{n})+{\tfrac {1}{2\alpha }}f(t_{n}+\alpha h,y_{n}+\alpha hf(t_{n},y_{n})){\bigr )}.}$

である。${\displaystyle \alpha =1/2}$ の場合が 中点法 (または修正オイラー法, modified Euler method)

${\displaystyle y_{n+1}=y_{n}+hf\left(t_{n}+{\frac {h}{2}},y_{n}+{\frac {h}{2}}f(t_{n},y_{n})\right)}$

に対応し、以下のブッチャー配列で与えられる。

	0
	1/2	1/2
		0	1

${\displaystyle \alpha =1}$ の場合は ホイン法（英語版）（または改良オイラー法, improved Euler method)として知られ、ブッチャー配列は以下の通りである^[3]。

	0
	1	1
		1/2	1/2

4段4次の公式

古典的ルンゲ゠クッタ法 (RK4) は既に述べた通り以下のブッチャー配列で与えられる^[8]:

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

修正版としてクッタの3/8公式 (英: Kutta's 3/8-rule) が知られている^[15]。

	0
	1/3	1/3
	2/3	-1/3	1
	1	1	-1	1
		1/8	3/8	3/8	1/8

他に使われている方法は、ルンゲ＝クッタ法のリストに参照。

計算例：2段2次陽的方法の条件の導出

前述の通り、2段の陽的方法が2次精度を持つための係数に対する条件は、${\displaystyle b_{1}+b_{2}=1}$, ${\displaystyle b_{2}c_{2}=1/2}$, かつ ${\displaystyle b_{2}a_{21}=1/2}$ である^[7]。例として、それらの条件の導出を見てみよう。

一般的に2段陽的ルンゲ＝クッタ法に対応する配列は以下の通りである(一貫性を用いて ${\displaystyle c_{1}=0}$ がわかる)。

${\displaystyle {\begin{array}{c|cc}0&&\\c_{2}&a_{21}&\\\hline &b_{1}&b_{2}\\\end{array}}}$

公式から ${\displaystyle k_{2}=f(t_{n}+c_{2}h,y_{n}+ha_{21}f(t_{n},y_{n}))}$ の f(t_n,y_n) に関するテイラー展開を考える

${\displaystyle k_{2}=f(t_{n},y_{n})+c_{2}hf_{t}(t_{n},y_{n})+ha_{21}f(t_{n},y_{n})f_{y}(t_{n},y_{n})+O(h^{2})}$

ここで、f_t は f の t に関する偏微分である。それを ${\displaystyle y_{n+1}=y_{n}+h(b_{1}k_{1}+b_{2}k_{2})}$ に代入し、${\displaystyle k_{1}=f(t_{n},y_{n})}$ を用いて整理すると

${\displaystyle y_{n+1}=y_{n}+h(b_{1}+b_{2})f(t_{n},y_{n})+h^{2}b_{2}(c_{2}f_{t}(t_{n},y_{n})+a_{21}f(t_{n},y_{n})f_{y}(t_{n},y_{n}))+O(h^{3})}$

となる。同時に、方程式 ${\displaystyle y'(t)=f(t,y)}$ の t に関する微分を取って y' を f(t,y) に置き換えると

${\displaystyle y''(t)=f_{t}(t,y)+f(t,y)f_{y}(t,y)}$

となる。厳密解を y(t) とする。上記式を用いて y(t_n+1) の t_n に関するテイラー展開を考える

${\displaystyle y(t_{n+1})=y(t_{n})+hf(t_{n},y_{n})+{\frac {1}{2}}h^{2}(f_{t}(t_{n},y_{n})+f(t_{n},y_{n})f_{y}(t_{n},y_{n}))+O(h^{3})}$

2次精度を持つ条件は、局所誤差 ${\displaystyle e_{n+1,h}=y(t_{n+1})-y_{n+1}=O(h^{3})}$ である。二つの展開式の係数を比較すると、その条件が${\displaystyle b_{1}+b_{2}=1}$, ${\displaystyle b_{2}c_{2}=1/2}$, かつ ${\displaystyle b_{2}a_{21}=1/2}$ であることがわかる。

埋め込み型ルンゲ＝クッタ法

ルンゲ＝クッタ法の局所誤差を精確に計算することが難しいので、実践では誤差を一定の範囲にコントロールするのが望ましい。そのために開発されたのが 埋め込み型ルンゲ＝クッタ法（Embedded Runge-Kutta method）である。適応型ルンゲ＝クッタ法 (adaptive Runge–Kutta method) とも呼ばれる。線型多段法にもミルンデバイスと呼ばれる相似の方法が存在する^[16]。

埋め込み方法は陽的ルンゲ＝クッタ法二つを利用する（下の方法を上の方法に埋め込むように見えるので埋め込み型と呼ばれる）。ブッチャー配列の以下のように拡張する。

	0
	${\displaystyle c_{2}}$	${\displaystyle a_{21}}$
	${\displaystyle c_{3}}$	${\displaystyle a_{31}}$	${\displaystyle a_{32}}$
	${\displaystyle \vdots }$	${\displaystyle \vdots }$		${\displaystyle \ddots }$
	${\displaystyle c_{s}}$	${\displaystyle a_{s1}}$	${\displaystyle a_{s2}}$	${\displaystyle \cdots }$	${\displaystyle a_{s,s-1}}$
		${\displaystyle b_{1}}$	${\displaystyle b_{2}}$	${\displaystyle \cdots }$	${\displaystyle b_{s-1}}$	${\displaystyle b_{s}}$
		${\displaystyle b_{1}^{*}}$	${\displaystyle b_{2}^{*}}$	${\displaystyle \cdots }$	${\displaystyle b_{s-1}^{*}}$	${\displaystyle b_{s}^{*}}$

ここで、b_i は p 次陽的方法に対応し、b *
i  は p-1 次陽的方法に対応する。二つの方法は係数 a_ij と c_i を共用する。正しく係数を選ぶと、二つの方法をともに収束させることができる。そのとき、埋め込み方法の時刻 t_n での相対（局所）誤差 e_n は次の公式で与えられる。

${\displaystyle e_{n}=y_{n}-y_{n}^{*}=h\sum _{i=1}^{s}(b_{i}-b_{i}^{*})k_{i}=O(h^{p})}$

ルンゲ＝クッタ法の収束性から、相対誤差も0に収束することがわかる。埋め込みルンゲ＝クッタ法は、アルゴリズムを用いて一時刻ごと（自動的）に刻み幅 h を調整し、誤差をコントロールすることができる（故に適応型と呼ばれる）^[17]。よって絶対誤差を計算せずに刻み幅を正しく設定することができる。そのアルゴリズムは、大体以下の通りである。

誤差の許容値を δ とする。刻み幅 h の近似値 y_n+1 と y *
n+1  をルンゲ＝クッタ法で計算する。二つの値に対し、${\displaystyle \|e_{n+1}\|\geq h\delta }$ が成立するとき、h が大きすぎて小さくする必要がある。よって新しい刻み幅 ${\displaystyle h'={\frac {1}{2}}h}$ と設定し再び近似値を計算する。代わりに ${\displaystyle \|e_{n+1}\|\leq {\frac {1}{10}}h\delta }$ が成立するとき、h が小さすぎて大きくするほうが効率が良い。よって新しい刻み幅 ${\displaystyle h'=2h}$ と設定し、次の時刻での近似値を計算する（この場合、誤差が許容値より小さいので今の近似値を二度と計算する必要はない）。二つの不等式ともが成立しないとき、そのままの刻み幅で次の時刻での近似値を計算してよい。

ただし、使用したスカラー（1/2他）は方程や方法によって変更することも可能である。また、実践では h を小さくしすぎないようにするため、下界 h_min を先に設定することもある。${\displaystyle h<h_{\mathrm {min} }}$ のときにアルゴリズムを停止し、代わりにもっと高い次数を持つ方法を使う。

例

もっとも単純な埋め込み型ルンゲ＝クッタ法は、ホイン法（2次）とオイラー法（1次）を組み合わせる方法である。以下の拡張ブッチャー配列で与えられる。

	0
	1	1
		1/2	1/2
		1	0

そして実践でよく使われているのが5次と4次の陽的方法を組み合わせる ルンゲ＝クッタ＝フェールベルグ法 である。対応する拡張ブッチャー配列は以下の通りである^[18]。

	0
	1/4	1/4
	3/8	3/32	9/32
	12/13	1932/2197	−7200/2197	7296/2197
	1	439/216	−8	3680/513	-845/4104
	1/2	−8/27	2	−3544/2565	1859/4104	−11/40
		16/135	0	6656/12825	28561/56430	−9/50	2/55
		25/216	0	1408/2565	2197/4104	−1/5	0

他に使われている方法は、Bogacki–Shampine法（英語版）(次数3と2)、Cash–Karp法（英語版） と ドルマン＝プリンス法（両方次数5と4）である。

陰的ルンゲ＝クッタ法

いままでに述べた方法はすべて陽公式である。陽的ルンゲ＝クッタ方法の絶対安定性領域（region of absolute stability）が小さくて有界であるため^[19]、硬い方程式の解を計算する場合、方法は不適切である。この問題は特に偏微分方程式の解を計算するときに重要である。

陽的方法の不安定さは陰的ルンゲ＝クッタ法の開発の動機となる。陰的ルンゲ＝クッタ法は上述の一般的な公式と同じく、以下の形で与えられる

${\displaystyle y_{n+1}=y_{n}+h\sum _{i=1}^{s}b_{i}k_{i},}$

但し、

${\displaystyle k_{i}=f\left(t_{n}+hc_{i},y_{n}+h\sum _{j=1}^{s}a_{ij}k_{j}\right),\quad i=1,\ldots ,s.}$

であり^[4]、対応するルンゲ＝クッタ行列は厳密な下三角行列ではない。結果として、一時刻ごとに代数方程式系を解かなければならなくなる。応じて、計算コストもかなり上がる。s 段法を使って m 個の成分からなる微分方程式系（すなわち ${\displaystyle {\vec {y}}=(y_{1},\ldots ,y_{m})}$ のとき）に適用する場合に対応する代数方程式の数は ms となる。これは陰的線型多段法とも比較できる：s 段陰的線型多段法を使って同じ方程式系に適用する場合に対応する代数方程式の数は m であり、段数 s に依存しない^[20]。

例

一番単純な陰的ルンゲ＝クッタ法は、後退オイラー法

${\displaystyle y_{n+1}=y_{n}+hf(t_{n+1},y_{n+1})}$

である。対応するブッチャー配列も単純に以下の通りである。

${\displaystyle {\begin{array}{c|c}1&1\\\hline &1\\\end{array}}}$

他の例として、台形公式と呼ばれる方法は以下の公式で与えられる。

${\displaystyle y_{n+1}=y_{n}+{\frac {1}{2}}h(f(t_{n},y_{n})+f(t_{n+1},y_{n+1}))}$

対応する配列は以下の通りである。

${\displaystyle {\begin{array}{c|cc}0&0&0\\1&{\frac {1}{2}}&{\frac {1}{2}}\\\hline &{\frac {1}{2}}&{\frac {1}{2}}\\\end{array}}}$

その対応は、${\displaystyle y_{n+1}=y_{n}+{\frac {1}{2}}k_{1}+{\frac {1}{2}}k_{2}}$ を覚えてて、公式を以下のように書き換えると明らかになる。

${\displaystyle y_{n+1}=y_{n}+{\frac {1}{2}}f(t_{n},y_{n})+{\frac {1}{2}}f\left(t_{n}+h,y_{n}+{\frac {1}{2}}k_{1}+{\frac {1}{2}}k_{2}\right)}$

台形公式は、選点法である。選点法はすべて陰的ルンゲ＝クッタ法であるけど、陰的ルンゲ＝クッタ法がすべて選点法であるわけではない^[21]。

選点法の中では、ガウス求積に基づいたガウス・ルジャンドル法（英語版）が一番次数の高い方法である。s 段ガウス・ルジャンドル法の次数は 2s となる（よって、任意に高い次数を持つ方法を構造できるようになる）^[22]。例として、2段ガウス・ルジャンドル法は次の配列で与えられる。

${\displaystyle {\begin{array}{c|cc}{\frac {1}{2}}-{\frac {1}{6}}{\sqrt {3}}&{\frac {1}{4}}&{\frac {1}{4}}-{\frac {1}{6}}{\sqrt {3}}\\{\frac {1}{2}}+{\frac {1}{6}}{\sqrt {3}}&{\frac {1}{4}}+{\frac {1}{6}}{\sqrt {3}}&{\frac {1}{4}}\\\hline &{\frac {1}{2}}&{\frac {1}{2}}\\\end{array}}}$ ^[22]

安定性

数値分析における安定性は、それぞれ異なる定義が複数存在する。ルンゲ＝クッタ法の安定性を反映できる概念は、主に以下の二つである。

A-安定性

線型テスト方程式 ${\displaystyle y'=\lambda y}$ を考える。一つのルンゲ＝クッタ法を使ってこの方程式に適用すると ${\displaystyle y_{n+1}=r(h\lambda )y_{n}}$ となる。ここで、

${\displaystyle r(z)=1+zb^{T}(I-zA)^{-1}e={\frac {\det(I-zA+zeb^{T})}{\det(I-zA)}}}$ ^[23]

は安定性関数と呼ばれる C 上の有理関数である（e はすべての成分が 1 のベクトル^{[要曖昧さ回避]}である）^[24]。s 段法の場合、行列式の展開によって r(z) は二つの s 次多項式の商となる。陽的方法の場合、対応するルンゲ＝クッタ行列が狭義下三角行列であるため、${\displaystyle \det(I-zA)=1}$ がわかる。したがって、r(z) は s 次多項式となる^[25]。

ルンゲ＝クッタ法による上記の方程式の数値解がゼロに減衰する条件が ${\displaystyle |r(z)|<1}$ (${\displaystyle z=h\lambda }$) である。上記の条件を満たす z の集合は方法に対する 絶対安定性領域 （region of absolute stability） である。 特に、絶対安定性領域が左複素数平面を含むとき、その方法は A-安定 (A-stable) という。陽的ルンゲ＝クッタ法は、安定性関数が多項式であるため、A-安定な方法になれない^[25]。

もっと一般的に、方法の次数が p のときに、安定性関数に対し ${\displaystyle r(z)=e^{z}+O(z^{p+1})}$ が成立する。ゆえに、指数関数 e^z の、与えられた次数の多項式の商からなる有理関数の中での最善近似が重要だと考えられる。そのような有理関数はパデ近似と呼ばれる。分子の次数が m で分母の次数が n のパデ近似式がA-安定性の条件 ${\displaystyle |r(z)|<1}$, ${\displaystyle \mathrm {Re} (z)\leq 0}$^[26] を満足するための必要十分条件は、${\displaystyle m\leq n\leq m+2}$ である^[27]。

s 段ガウス・ルジャンドル法の次数は前述通り 2s である。よって安定性関数の分子と分母は同じく s 次多項式となる。すなわち、 ${\displaystyle m=n=s}$ である。上記の条件が満たされるので、ガウス・ルジャンドル法はA-安定であることがわかる^[28]。故に任意に高い次数を持つ、A-安定なルンゲ＝クッタ法が存在する。比べて、A-安定性を持つ線型多段法の次数は2以下である^[29]。

以上のことから、陰的ルンゲ＝クッタ法は陽的な方法よりも優れた安定性を持つこともわかる。

B-安定性

A-安定性という概念は線型自励方程式 ${\displaystyle y'=\lambda y}$ を考える結果である。 ダールキスト（英語版）は、数値的方法を使って、とある単調性条件を満たす非線型方程式系に適用するときの安定性を主張した。対応する安定性は、線型多段法の場合に G-安定性 （G-stability）で、ルンゲ＝クッタ法の場合に B-安定性 (B-stability)^[30] と呼ぶ。 一般的なテスト方程式 ${\displaystyle y'=f(x,y)}$ と単調性条件 ${\displaystyle \langle f(x,y)-f(x,z),y-z\rangle \leq 0}$ を満足する f を考える。もしすべての ${\displaystyle h\geq 0}$ に対し、 不等式 ${\displaystyle \|y_{n+1}-z_{n+1}\|\leq \|y_{n}-z_{n}\|}$ が成立するとき、そのルンゲ＝クッタ法は B-安定 という。 ここで、y_n と z_n はそれぞれの初期値に対する数値解である。 ${\displaystyle f(x,y)=\lambda y}$ に方法を適用することで、B-安定性はA-安定性より強い条件であることがわかる^[31]。

脚注

参考文献

• Atkinson, Kendall A. (1989), An Introduction to Numerical Analysis (2nd ed.), New York: John Wiley & Sons, ISBN 978-0-471-50023-0 .
• Butcher, John C. (May 1963), “Coefficients for the study of Runge-Kutta integration processes”, Journal of the Australian Mathematical Society 3 (2): 185–201, doi:10.1017/S1446788700027932, http://journals.cambridge.org/action/displayAbstract?fromPage=online&aid=4922056 .
• Butcher, John C. (1975), “A stability property of implicit Runge-Kutta methods”, BIT 15: 358–361, doi:10.1007/bf01931672 .
• Butcher, John C. (2008), Numerical Methods for Ordinary Differential Equations, New York: John Wiley & Sons, ISBN 978-0-470-72335-7 .
• Dahlquist, Germund (1963), “A special stability problem for linear multistep methods”, BIT 3: 27–43, doi:10.1007/BF01963532, ISSN 0006-3835 .
• Hairer, Ernst; Nørsett, Syvert Paul; Wanner, Gerhard (1993), Solving ordinary differential equations I: Nonstiff problems, Berlin, New York: Springer-Verlag, ISBN 978-3-540-56670-0 .
• Hairer, Ernst; Wanner, Gerhard (1996), Solving ordinary differential equations II: Stiff and differential-algebraic problems (2nd ed.), Berlin, New York: Springer-Verlag, ISBN 978-3-540-60452-5 .
• Iserles, Arieh (2008), A First Course in the Numerical Analysis of Differential Equations (Second Edition), Cambridge University Press, ISBN 978-0-521-73490-5 .
• Kaw, Autar; Kalu, Egwu (2008), Numerical Methods with Applications (1st ed.), autarkaw.com, http://numericalmethods.eng.usf.edu/topics/textbook_index.html .
• Kutta, Martin Wilhelm (1901), “Beitrag zur näherungsweisen Integration totaler Differentialgleichungen”, Zeitschrift für Mathematik und Physik 46: 435–453 .
• Press, William H.; Flannery, Brian P.; Teukolsky, Saul A.; Vetterling, William T. (2007), “Section 17.1 Runge-Kutta Method”, Numerical Recipes: The Art of Scientific Computing (3rd ed.), Cambridge University Press, ISBN 978-0-521-88068-8, http://apps.nrbook.com/empanel/index.html#pg=907 . Also, Section 17.2. Adaptive Stepsize Control for Runge-Kutta.
• Stoer, Josef; Bulirsch, Roland (2002), Introduction to Numerical Analysis (3rd ed.), Berlin, New York: Springer-Verlag, ISBN 978-0-387-95452-3 .
• Süli, Endre; Mayers, David (2003), An Introduction to Numerical Analysis, Cambridge University Press, ISBN 0-521-00794-1 .

• John C. Butcher: "B-Series : Algebraic Analysis of Numerical Methods", Springer（SSCM, volume 55), ISBN 978-3030709556 (April, 2021).

• Runge-Kutta 法についてのノート (早川尚男) # 公式を導出する例を示している。

関連項目

• ルンゲ＝クッタ法のリスト
• 常微分方程式の数値解法
• 線型多段法