ルンゲ＝クッタ法

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

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

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

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

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

ここで、

である。次の値 () は、現在の値 () に増分を加えたものであり、増分は勾配の推定値に間隔 を乗じたものになっている。勾配の推定値は、, ..., の4つの勾配の重み付け平均で求める。, ..., のそれぞれの勾配は、特定の に対する によって与えられ、以下のように解釈できる。 は区間の最初 における勾配である (オイラー法で用いる勾配に一致する)。 は区間の中央 における勾配の近似値である (中点法で用いる勾配)。計算に用いる中央の の値は、初期位置の勾配 を用いてオイラー法で推定する。 は区間の中央における勾配のもう一つの近似値である。中央の を の値から推定して用いる。 は区間の最後 における勾配の近似値であり、 の値から推定された最後の点の の値を用いる。 重み付き平均では、中央の勾配に対して大きな重みを用いる。シンプソン則を用いた平均と同等の形になる.。

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

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

但し; .、

(文献によって、等価であるが上と異なる定義の仕方をしているものがあることに注意する), , and leave out the factor h in the definition of the stages. , and use the y values as stages.。 具体的なルンゲ゠クッタ公式は、解ができるだけ高い精度を持つように適切に選ばれた係数 (ルンゲ゠クッタ行列)、 (重み)、 (節点) で指定される ().。特に に対して を満たす方法が広く用いられ、総称して 陽的ルンゲ゠クッタ法 (ERK、) と呼ぶ。そうでないものを 陰的ルンゲ゠クッタ法 (IRK、) と呼ぶ。

近似値 を から計算するときに発生する誤差の大きさが のとき、そのルンゲ゠クッタ公式は 次精度を持つといい、 を次数 (または位数) と呼ぶ。 次のルンゲ゠クッタ公式は、誤差の大きさの条件に誤差の表式を代入し、係数の条件を求めることによって得られる。例えば、2段の陽的方法が2次精度を持つための係数に対する条件は、, , かつ である.。この条件を満たす範囲内で様々な方法を考え得る。

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

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

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

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

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

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

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

とおく。 あるいは を適切に（線型結合により）近似することでルンゲ＝クッタ法の公式が得られる。さらに、テイラー展開より係数を正しく選択すれば、方法の収束性は求積法の収束性から保証される。しかし、局所誤差のオーダーや上界は、方法によって大きく異なるので、方法別に計算しなければならない。

陽的ルンゲ゠クッタ法 陽的ルンゲ゠クッタ法では () が満たされるので、ブッチャー配列は以下の形になる。 0 この時、各勾配 , ..., を以下のように逐次的に求めることができる。

(文献によって、等価であるが上と異なる定義の仕方をしているものがあることに注意する)。 陽的ルンゲ゠クッタ法が目的の次数 を実現するためには十分な段数 が必要になる。少なくとも次数以上の段数が必要 () であり、更に5次以上の場合には段数を次数より大きく取らねばならない () ことが知られている。各次数 に対してそれを実現するための最小の段数 を与える一般式は未解決問題である。具体的な値が幾つか知られている:

各段数 に対して達成が可能な最大の次数 は以下の表になる。

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

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

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

である。 の場合が 中点法 (または修正オイラー法, modified Euler method)

に対応し、以下のブッチャー配列で与えられる。 01/2 1/2 0 1 の場合は （または改良オイラー法, improved Euler method)として知られ、ブッチャー配列は以下の通りである。 0 1 11/21/2

4段4次の公式 古典的ルンゲ゠クッタ法 (RK4) は既に述べた通り以下のブッチャー配列で与えられる: 01/21/21/201/210011/61/31/31/6 修正版としてクッタの3/8公式 () が知られている refer to .。 01/31/32/3 -1/3111 -111/83/83/81/8

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

計算例：2段2次陽的方法の条件の導出 前述の通り、2段の陽的方法が2次精度を持つための係数に対する条件は、, , かつ である。例として、それらの条件の導出を見てみよう。

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

公式から の に関するテイラー展開を考える

ここで、 は の に関する偏微分である。それを に代入し、 を用いて整理すると

となる。同時に、方程式 の に関する微分を取って を に置き換えると

となる。厳密解を とする。上記式を用いて の に関するテイラー展開を考える

2次精度を持つ条件は、局所誤差 である。二つの展開式の係数を比較すると、その条件が, , かつ であることがわかる。

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

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

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

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

誤差の許容値を とする。刻み幅 の近似値 と をルンゲ＝クッタ法で計算する。二つの値に対し、 が成立するとき、 が大きすぎて小さくする必要がある。よって新しい刻み幅 と設定し再び近似値を計算する。代わりに が成立するとき、 が小さすぎて大きくするほうが効率が良い。よって新しい刻み幅 と設定し、次の時刻での近似値を計算する（この場合、誤差が許容値より小さいので今の近似値を二度と計算する必要はない）。二つの不等式がどちらも成立しないとき、そのままの刻み幅で次の時刻での近似値を計算してよい。

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

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

0 1 1 1/2 1/2 1 0

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

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）が小さくて有界であるため、硬い方程式の解を計算する場合、方法は不適切である。この問題は特に偏微分方程式の解を計算するときに重要である。

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

但し、

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

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

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

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

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

その対応は、 を覚えてて、公式を以下のように書き換えると明らかになる。

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

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

安定性 数値分析における安定性は、それぞれ異なる定義が複数存在する。ルンゲ＝クッタ法の安定性を反映できる概念は、主に以下の二つである。 A-安定性 線型テスト方程式 を考える。一つのルンゲ＝クッタ法を使ってこの方程式に適用すると となる。ここで、

は安定性関数と呼ばれる 上の有理関数である（ はすべての成分が 1 のベクトルである）。 段法の場合、行列式の展開によって は二つの 次多項式の商となる。陽的方法の場合、対応するルンゲ＝クッタ行列が狭義下三角行列であるため、 がわかる。したがって、 は 次多項式となる.。

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

もっと一般的に、方法の次数が のときに、安定性関数に対し が成立する。ゆえに、指数関数 の、与えられた次数の多項式の商からなる有理関数の中での最善近似が重要だと考えられる。そのような有理関数はパデ近似と呼ばれる。分子の次数が で分母の次数が のパデ近似式がA-安定性の条件 , はそのような有理関数を A-acceptable と呼ぶ。 を満足するための必要十分条件は、 である。

段ガウス・ルジャンドル法の次数は前述通り 2 である。よって安定性関数の分子と分母は同じく 次多項式となる。すなわち、 である。上記の条件が満たされるので、ガウス・ルジャンドル法はA-安定であることがわかる。故に任意に高い次数を持つ、A-安定なルンゲ＝クッタ法が存在する。比べて、A-安定性を持つ線型多段法の次数は2以下であるこの結果（時々にsecond Dahlquist barrier）は、 によるものである。。

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

B-安定性 A-安定性という概念は線型自励方程式 を考える結果である。 は、数値的方法を使って、とある単調性条件を満たす非線型方程式系に適用するときの安定性を主張した。対応する安定性は、線型多段法の場合に G-安定性 （G-stability）で、ルンゲ＝クッタ法の場合に B-安定性 (B-stability) と呼ぶ。 一般的なテスト方程式 と単調性条件 を満足する を考える。もしすべての に対し、 不等式 が成立するとき、そのルンゲ＝クッタ法は B-安定 という。 ここで、 と はそれぞれの初期値に対する数値解である。 に方法を適用することで、B-安定性はA-安定性より強い条件であることがわかる。

脚注

参考文献 . . . . . . . . . . . Also, Section 17.2. Adaptive Stepsize Control for Runge-Kutta. . .

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

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

関連項目

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