ルンゲ＝クッタ法

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

古典的ルンゲ＝クッタ法

一連のルンゲ＝クッタ公式の中で最も広く知られているのが、古典的ルンゲ=クッタ法 (RK4、もしくは単に狭義の ルンゲ=クッタ法、テンプレート:Lang-en-short) などと呼ばれる4次の公式である。

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

$y^{\prime }=f(t,y),y(t_0)=y_0.$

但し、テンプレート:Math が近似的に求めたい未知関数であり、その テンプレート:Mvar における勾配は テンプレート:Math によって テンプレート:Mvar 及び テンプレート:Math の関数として与えられている。時刻 テンプレート:Mvar における初期値は テンプレート:Mvar で与えられている。

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

$\begin{array}{cc}{y}_{n+1}& =y_n+\frac{h}{6}(k_1+2k_2+2k_3+k_4),\\ t_{n+1}& =t_n+h.\end{array}$

ここで、

$\begin{array}{cc}{k}_1& =f(t_n,y_n),\\ k_2& =f(t_n+\frac{h}{2},y_n+\frac{h}{2}{k}_1),\\ k_3& =f(t_n+\frac{h}{2},y_n+\frac{h}{2}{k}_2),\\ k_4& =f(t_n+h,y_n+h{k}_3)\end{array}$

であるテンプレート:Sfnテンプレート:Sfn。次の値 (テンプレート:Math) は、現在の値 (テンプレート:Mvar) に増分を加えたものであり、増分は勾配の推定値に間隔 テンプレート:Mvar を乗じたものになっている。勾配の推定値は、テンプレート:Mvar, ..., テンプレート:Mvar の4つの勾配の重み付け平均で求める。テンプレート:Mvar, ..., テンプレート:Mvar のそれぞれの勾配は、特定の テンプレート:Math に対する テンプレート:Mvar によって与えられ、以下のように解釈できる。

• テンプレート:Mvar は区間の最初 テンプレート:Mvar における勾配である (オイラー法で用いる勾配に一致する)。
• テンプレート:Mvar は区間の中央 テンプレート:Math における勾配の近似値である (中点法で用いる勾配)。計算に用いる中央の テンプレート:Mvar の値は、初期位置の勾配 テンプレート:Mvar を用いてオイラー法で推定する。
• テンプレート:Mvar は区間の中央における勾配のもう一つの近似値である。中央の テンプレート:Mvar を テンプレート:Mvar の値から推定して用いる。
• テンプレート:Mvar は区間の最後 テンプレート:Math における勾配の近似値であり、テンプレート:Mvar の値から推定された最後の点の テンプレート:Mvar の値を用いる。

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

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

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

ルンゲ゠クッタ法

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

$y_{n+1}=y_n+h\sum _{i=1}^s{b}_i{k}_i,$

但し、

$k_i=f(t_n+h{c}_i,y_n+h\sum _{j=1}^s{a}_{ij}{k}_j),i=1,\dots ,s.$

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

具体的なルンゲ゠クッタ公式は、解ができるだけ高い精度を持つように適切に選ばれた係数 テンプレート:Mvar (ルンゲ゠クッタ行列)、テンプレート:Mvar (重み)、テンプレート:Mvar (節点) で指定される ($i,j\le s$)。特に $i\le j$ に対して $a_{ij}=0$ を満たす方法が広く用いられ、総称して 陽的ルンゲ゠クッタ法 (ERK、テンプレート:Lang-en-short) と呼ぶ。そうでないものを 陰的ルンゲ゠クッタ法 (IRK、テンプレート:Lang-en-short) と呼ぶ。

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

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

$\begin{array}{cc}\begin{array}{ccccc}{c}_1& a_{11}& a_{12}& \cdots & a_{1s}\\ c_2& a_{21}& a_{22}& \cdots & a_{2s}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ c_s& a_{s1}& a_{s2}& \cdots & a_{ss}\\ & b_1& b_2& \cdots & b_s\end{array}& =\begin{array}{cc}𝐜& A\\ & {𝐛}^{\mathrm{T}}\end{array}\end{array}$

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

$\begin{array}{c}\begin{array}{cc}0& \\ 1/2& 1/2\\ 1/2& 0& 1/2\\ 1& 0& 0& 1\\ & 1/6& 1/3& 1/3& 1/6\end{array}\end{array}$

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

$\sum _{j=1}^s{a}_{ij}=c_i,i=1,\dots ,s.$

収束性

ルンゲ＝クッタ法は、数値積分による求積法 (quadrature) と深く関係にある。時刻 テンプレート:Mvar での値から テンプレート:Math での値を求めるときの方程式が以下であるとする。

$y^{\prime }=f(t,y),y(t_n)=y_n.$

求積法は、与えられた区間での定積分の値を被積分関数の値の線型結合により近似する方法である。簡単のために、区間を テンプレート:Math とする。すると求積法の公式は

$\underset{0}{\overset{1}{\int }}g(t)dt\approx \sum _{i=1}^s{b}_{i}g(c_i)$

となる。ここで、テンプレート:Mvar と テンプレート:Mvar はあらかじめ選ばれた定数であり、前述の重みと節点に対応する。上記の式に対し、すべての テンプレート:Math 次以下の多項式 g に対して等式が成立するとき（すなわち、誤差が0のとき）には、その求積法は テンプレート:Mvar 次精度であり、テンプレート:Mvar を次数と呼ぶテンプレート:Sfn。節点が テンプレート:Mvar 個の時、最大次数は 2テンプレート:Mvar であり、その方法は テンプレート:Mvar 次ガウス・ルジャンドル公式と呼ばれる。

そこで上記の常微分方程式を積分形式に変形して、求積法を被積分関数 $g(t)=f(t,y(t))$ に対して適用すると以下の近似式が得られる。

$y(t_{n+1})\approx y_n+\underset{y_n}{\overset{y_{n+1}}{\int }}f(t,y(t))dt=y_n+h\underset{0}{\overset{1}{\int }}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))$

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

陽的ルンゲ゠クッタ法

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

	0
	$c_2$	$a_{21}$
	$c_3$	$a_{31}$	$a_{32}$
	$⋮$	$⋮$		$\ddots$
	$c_s$	$a_{s1}$	$a_{s2}$	$\cdots$	$a_{s,s-1}$
		$b_1$	$b_2$	$\cdots$	$b_{s-1}$	$b_s$

この時、各勾配 テンプレート:Math, ..., テンプレート:Mvar を以下のように逐次的に求めることができるテンプレート:Sfn。

$\begin{array}{cc}{k}_1& =f(t_n,y_n),\\ k_2& =f(t_n+c_{2}h,y_n+a_{21}h{k}_1),\\ k_3& =f(t_n+c_{3}h,y_n+a_{31}h{k}_1+a_{32}h{k}_2),\\ & ⋮\\ k_s& =f(t_n+c_{s}h,y_n+a_{s1}h{k}_1+a_{s2}h{k}_2+\cdots +a_{s,s-1}h{k}_{s-1}).\end{array}$

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

陽的ルンゲ゠クッタ法が目的の次数 テンプレート:Mvar を実現するためには十分な段数 テンプレート:Mvar が必要になる。少なくとも次数以上の段数が必要 ($s\ge p$) であり、更に5次以上の場合には段数を次数より大きく取らねばならない ($s>p$) ことが知られているテンプレート:Sfn。各次数 テンプレート:Mvar に対してそれを実現するための最小の段数 テンプレート:Mvar を与える一般式は未解決問題である。具体的な値が幾つか知られているテンプレート:Sfn:

$\begin{array}{ccccccccc}p& 1& 2& 3& 4& 5& 6& 7& 8\\ \min s& 1& 2& 3& 4& 6& 7& 9& 11\end{array}$

各段数 テンプレート:Mvar に対して達成が可能な最大の次数 テンプレート:Mvar は以下の表になる。

$\begin{array}{cccccccccccc}s& 1& 2& 3& 4& 5& 6& 7& 8& 9& 10& 11\\ \max p& 1& 2& 3& 4& 4& 5& 6& 6& 7& 7& 8\end{array}$

オイラー法、ホイン法、古典的ルンゲ＝クッタ法 (RK4) の相対誤差の比較。初期値 $y(0)=0$ のもとでの常微分方程式 $\frac{dy}{dt}=\cos (y)$ の数値解の $t=1$ での値をステップサイズの関数として対数プロットした。各手法がそれぞれ1次精度、2次精度、4次精度であることに対応して、両対数グラフ上で直線の傾きがそれぞれ 1, 2, 4 でステップサイズが小さくなるにつれて誤差が減少している。

1段1次の公式

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

$y_{n+1}=y_n+hf(t_n,y_n).$

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

	0
		1

2段2次の公式

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

	0
	$\alpha$	$\alpha$
		$(1-\frac{1}{2\alpha })$	$\frac{1}{2\alpha }$

対応する公式はテンプレート:Sfn

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

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

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

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

	0
	1/2	1/2
		0	1

$\alpha =1$ の場合は テンプレート:仮リンク（または改良オイラー法, improved Euler method)として知られ、ブッチャー配列は以下の通りである。

	0
	1	1
		1/2	1/2

4段4次の公式

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

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

修正版としてクッタの3/8公式 (テンプレート:Lang-en-short) が知られている。

	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次精度を持つための係数に対する条件は、$b_1+b_2=1$, $b_2{c}_2=1/2$, かつ $b_2{a}_{21}=1/2$ である。例として、それらの条件の導出を見てみよう。

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

$\begin{array}{ccc}0& & \\ c_2& a_{21}& \\ & b_1& b_2\end{array}$

公式から $k_2=f(t_n+c_{2}h,y_n+h{a}_{21}f(t_n,y_n))$ の テンプレート:Math に関するテイラー展開を考える

$k_2=f(t_n,y_n)+c_{2}h{f}_t(t_n,y_n)+h{a}_{21}f(t_n,y_n)f_y(t_n,y_n)+O(h^2)$

ここで、テンプレート:Mvar は テンプレート:Mvar の テンプレート:Mvar に関する偏微分である。それを $y_{n+1}=y_n+h(b_1{k}_1+b_2{k}_2)$ に代入し、$k_1=f(t_n,y_n)$ を用いて整理すると

$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)$

となる。同時に、方程式 $y^{\prime }(t)=f(t,y)$ の テンプレート:Mvar に関する微分を取って テンプレート:Mvar を テンプレート:Math に置き換えると

$y^{″}(t)=f_t(t,y)+f(t,y)f_y(t,y)$

となる。厳密解を テンプレート:Math とする。上記式を用いて テンプレート:Math の テンプレート:Mvar に関するテイラー展開を考える

$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次精度を持つ条件は、局所誤差 $e_{n+1,h}=y(t_{n+1})-y_{n+1}=O(h^3)$ である。二つの展開式の係数を比較すると、その条件が$b_1+b_2=1$, $b_2{c}_2=1/2$, かつ $b_2{a}_{21}=1/2$ であることがわかる。

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

ルンゲ＝クッタ法の局所誤差を正確に求めることは難しいので、計算を行なう場合には誤差が一定の範囲に収まるように制御できることが望ましい。そのために開発されたのが 埋め込み型ルンゲ＝クッタ法（Embedded Runge-Kutta method）である。適応型ルンゲ＝クッタ法 (adaptive Runge–Kutta method) とも呼ばれる。線型多段法にもミルンデバイスと呼ばれる相似の方法が存在するテンプレート:Sfn。

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

	0
	$c_2$	$a_{21}$
	$c_3$	$a_{31}$	$a_{32}$
	$⋮$	$⋮$		$\ddots$
	$c_s$	$a_{s1}$	$a_{s2}$	$\cdots$	$a_{s,s-1}$
		$b_1$	$b_2$	$\cdots$	$b_{s-1}$	$b_s$
		$b_1^{*}$	$b_2^{*}$	$\cdots$	$b_{s-1}^{*}$	$b_s^{*}$

ここで、テンプレート:Mvar は テンプレート:Mvar 次陽的方法に対応し、テンプレート:Math は テンプレート:Math 次陽的方法に対応する。二つの方法は係数 テンプレート:Mvar と テンプレート:Mvar を共用する。係数を正しく選べば、二つの方法をともに収束させることができる。そのとき、埋め込み方法の時刻 テンプレート:Mvar での相対（局所）誤差 テンプレート:Mvar は次の公式で与えられる。

$e_n=y_n-y_n^{*}=h\sum _{i=1}^s(b_i-b_i^{*})k_i=O(h^p)$

ルンゲ＝クッタ法の収束性から、相対誤差も0に収束することがわかる。埋め込み型ルンゲ＝クッタ法は、アルゴリズムを用いて一時刻ごと（自動的）に刻み幅 テンプレート:Mvar を調整することで誤差を制御できる（故に適応型と呼ばれる）テンプレート:Sfn。よって絶対誤差を求めなくても刻み幅を適切に制御することができる。そのアルゴリズムは、大体以下の通りである。

誤差の許容値を テンプレート:Mvar とする。刻み幅 テンプレート:Mvar の近似値 テンプレート:Mvar と テンプレート:Math をルンゲ＝クッタ法で計算する。二つの値に対し、$\Vert e_{n+1}\Vert \ge h\delta$ が成立するとき、テンプレート:Mvar が大きすぎて小さくする必要がある。よって新しい刻み幅 $h^{\prime }=\frac{1}{2}h$ と設定し再び近似値を計算する。代わりに $\Vert e_{n+1}\Vert \le \frac{1}{10}h\delta$ が成立するとき、テンプレート:Mvar が小さすぎて大きくするほうが効率が良い。よって新しい刻み幅 $h^{\prime }=2h$ と設定し、次の時刻での近似値を計算する（この場合、誤差が許容値より小さいので今の近似値を二度と計算する必要はない）。二つの不等式がどちらも成立しないとき、そのままの刻み幅で次の時刻での近似値を計算してよい。

ただし、使用したスカラー（1/2他）の値は方程式や方法によって変更することもできる。また、実際の計算では テンプレート:Mvar が小さくなり過ぎることを避けるために、下界 テンプレート:Mvar を先に設定することもある。$h<h_{\mathrm{m}\mathrm{i}\mathrm{n}}$ のときにはアルゴリズムを停止して、別のより次数の高い方法を使う。

例

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

	0
	1	1
		1/2	1/2
		1	0

実際の計算でよく用いられているのは、5次と4次の陽的方法を組み合わせる ルンゲ＝クッタ＝フェールベルグ法 である。対応する拡張ブッチャー配列は以下の通りであるテンプレート:Sfn。

	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

使われている他の方法は、テンプレート:仮リンク(次数3と2)、テンプレート:仮リンク と ドルマン＝プリンス法（どちらの方法も次数の組合せは5と4）である。

陰的ルンゲ＝クッタ法

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

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

$y_{n+1}=y_n+h\sum _{i=1}^s{b}_i{k}_i,$

但し、

$k_i=f(t_n+h{c}_i,y_n+h\sum _{j=1}^s{a}_{ij}{k}_j),i=1,\dots ,s.$

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

例

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

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

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

$\begin{array}{cc}1& 1\\ & 1\end{array}$

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

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

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

$\begin{array}{ccc}0& 0& 0\\ 1& \frac{1}{2}& \frac{1}{2}\\ & \frac{1}{2}& \frac{1}{2}\end{array}$

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

$y_{n+1}=y_n+\frac{1}{2}f(t_n,y_n)+\frac{1}{2}f(t_n+h,y_n+\frac{1}{2}{k}_1+\frac{1}{2}{k}_2)$

台形公式は、選点法である。選点法はすべて陰的ルンゲ＝クッタ法であるけど、陰的ルンゲ＝クッタ法がすべて選点法であるわけではないテンプレート:Sfn。

選点法の中では、ガウス求積に基づいたテンプレート:仮リンクが一番次数の高い方法である。テンプレート:Mvar 段ガウス・ルジャンドル法の次数は テンプレート:Math となる（よって、任意に高い次数を持つ方法を構造できるようになる）。例として、2段ガウス・ルジャンドル法は次の配列で与えられる。

$\begin{array}{ccc}\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}\\ & \frac{1}{2}& \frac{1}{2}\end{array}$

安定性

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

A-安定性

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

$r(z)=1+z{b}^T(I-zA{)}^{-1}e=\frac{\det (I-zA+ze{b}^T)}{\det (I-zA)}$ テンプレート:Sfn

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

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

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

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

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

B-安定性

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

参考文献

• テンプレート:Citation.
• テンプレート:Citation.
• テンプレート:Citation.
• テンプレート:Citation.
• テンプレート:Citation.
• テンプレート:Citation.
• テンプレート:Citation.
• テンプレート:Citation.
• テンプレート:Citation.
• テンプレート:Citation.
• テンプレート:Citation. Also, Section 17.2. Adaptive Stepsize Control for Runge-Kutta.
• テンプレート:Citation.
• テンプレート:Citation.

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

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

関連項目

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