ルンゲ＝クッタ法
ルンゲ＝クッタ法（）とは、数値解析において微分方程式の初期値問題に対して近似解を与える一連の方法である。この技法は1900年頃に数学者カール・ルンゲとによって発展された。
一連のルンゲ＝クッタ公式の中で最も広く知られているのが、古典的ルンゲ=クッタ法 (RK4、もしくは単に狭義の ルンゲ=クッタ法、) などと呼ばれる4次の公式である。
次の初期値問題を考える。
但し、 が近似的に求めたい未知関数であり、その における勾配は によって 及び の関数として与えられている。時刻 における初期値は で与えられている。
今、時刻 における値 = が既知のとき、十分に小さなステップ幅 に対して , を以下の式で与えると、 は の 4次精度の近似になっている。このステップを逐次的に繰り返すことによって、初期値 から任意の時刻 における近似値 が求められる。
ここで、
である。次の値 () は、現在の値 () に増分を加えたものであり、増分は勾配の推定値に間隔 を乗じたものになっている。勾配の推定値は、, ..., の4つの勾配の重み付け平均で求める。, ..., のそれぞれの勾配は、特定の に対する によって与えられ、以下のように解釈できる。
重み付き平均では、中央の勾配に対して大きな重みを用いる。シンプソン則を用いた平均と同等の形になる。
RK4は4次の方法である。厳密解とRK4のテイラー展開が4次の項まで一致し、1ステップの推定誤差は formula_5 のオーダーになる。目的の時刻の を求めるのに必要なステップ数は formula_6 になるので、全体の推定誤差は formula_7 になる。
前述の RK4 の一般化として、以下の形式を持つ 段のルンゲ゠クッタ法を構成することができる。整数 をそのルンゲ゠クッタ法の段数 (stage) という。
但し、
具体的なルンゲ゠クッタ公式は、解ができるだけ高い精度持つように適切に選ばれた係数 (ルンゲ゠クッタ行列)、 (重み)、 (節点) で指定される (formula_10)。特に formula_11 に対して formula_12 を満たす方法が広く用いられ、総称して 陽的ルンゲ゠クッタ法 (ERK、) と呼ぶ。そうでないものを 陰的ルンゲ゠クッタ法 (IRK、) と呼ぶ。
近似値 を から計算するときに発生する誤差の大きさが formula_13 のとき、そのルンゲ゠クッタ公式は 次精度を持つといい、 を次数 (または位数) と呼ぶ。 次のルンゲ゠クッタ公式は、誤差の大きさの条件に誤差の表式を代入し、係数の条件を求めることによって得られる。例えば、2段の陽的方法が2次精度を持つための係数に対する条件は、formula_14, formula_15, かつ formula_16 である。この条件を満たす範囲内で様々な方法を考え得る。
これらの係数を分かりやすく表記する方法として、以下のような形式のブッチャー配列 ("Butcher tableau") が知られている:
実際には、それぞれのルンゲ゠クッタ法について各要素に具体的な値を入れて用いる。陽的な方法ではルンゲ゠クッタ行列の上三角成分は常に 0 になるので表記を省略する。例えば RK4 はブッチャー配列を用いて以下のように表現される。
ルンゲ゠クッタ法が一貫しているためには、次の条件が満たされている必要がある。
ルンゲ＝クッタ法は、数値積分における求積法 (quadrature) と深く繋がる。時刻 での値から での値を求めるときの方程式は以下のように定める。
求積法は、与えられた区間での定積分の値を被積分関数の値の線型結合として近似する方法である。簡単のために、区間を とする。よって求積法の公式は
となる。ここで、 と は先に選ばれた定数であり、前述の重みと節点に対応する。上記式に対し、等式がすべての 次以下の多項式に成立するとき（すなわち、誤差が0のとき）、その求積法は 次精度であり、 を次数と呼ぶ。節点が 個の時、最大次数は 2 であり、その方法は 次ガウス・ルジャンドル公式と呼ばれる。
そして上記の方程式を積分形式に変形し、求積法を用いると次の公式となる。
formula_23 とおく。, あるいは を適切に（線型結合として）近似することでルンゲ＝クッタ法の公式となる。その上、係数をテイラー展開より正しく選択すると、方法の収束性も求積法の収束性より保証される。しかし、局所誤差のオーダーや上界は、方法によって大きく異なるので、方法別に計算しなければならない。
陽的ルンゲ゠クッタ法では formula_24 (formula_25) が満たされるので、ブッチャー配列は以下の形になる。
この時、各勾配 , ..., を以下のように逐次的に求めることができる。
陽的ルンゲ゠クッタ法が目的の精度 になる係数を持つためには、十分に大きな段数 が必要になる。少なくとも次数以上の段数が必要 (formula_27) であり、更に5次以上の場合には段数は次数よりも大きく取らなければならない (formula_28) ことが知られている。各次数 に対する最低段数 の一般式は未解決問題である。具体的な値が幾つか知られている:
最も単純なルンゲ゠クッタ法が 1段1次 の方法であり、一意に定まる。(前進) オイラー法と等価になる。
以下のブッチャー配列で表現される。
2段2次の方法は1パラメータの自由度 を持ち、以下のブッチャー配列で与えられる。
対応する公式は
である。formula_32 の場合が (または修正オイラー法, modified Euler method)
に対応し、以下のブッチャー配列で与えられる。
formula_34 の場合は （または改良オイラー法, improved Euler method)として知られ、ブッチャー配列は以下の通りである。
古典的ルンゲ゠クッタ法 (RK4) は既に述べた通り以下のブッチャー配列で与えられる:
修正版としてクッタの3/8公式 () が知られている。
他に使われている方法は、ルンゲ＝クッタ法のリストに参照。
前述の通り、2段の陽的方法が2次精度を持つための係数に対する条件は、formula_14, formula_15, かつ formula_16 である。例として、それらの条件の導出を見てみよう。
一般的に2段陽的ルンゲ＝クッタ法に対応する配列は以下の通りである(一貫性を用いて formula_38 がわかる)。
公式から formula_40 の に関するテイラー展開を考える
ここで、 は の に関する偏微分である。それを formula_42 に代入し、formula_43 を用いて整理すると
となる。同時に、方程式 formula_45 の に関する微分を取って を に置き換えると
となる。厳密解を とする。上記式を用いて の に関するテイラー展開を考える
2次精度を持つ条件は、局所誤差 formula_48 である。二つの展開式の係数を比較すると、その条件がformula_14, formula_15, かつ formula_16 であることがわかる。
ルンゲ＝クッタ法の局所誤差を精確に計算することが難しいので、実践では誤差を一定の範囲にコントロールするのが望ましいである。そのために開発されたのが 埋め込み型ルンゲ＝クッタ法（Embedded Runge-Kutta method）である。適応型ルンゲ＝クッタ法 (adaptive Runge–Kutta method) とも呼ばれる。線型多段法にもミルンデバイスと呼ばれる相似の方法が存在する。
埋め込み方法は陽的ルンゲ＝クッタ法二つを利用する（下の方法を上の方法に埋め込むようで埋め込み型）。ブッチャー配列の以下のように拡張する。
ここで、 は 次陽的方法に対応し、 は 次陽的方法に対応する。二つの方法は係数 と を共用する。正しく係数を選ぶと、二つの方法をともに収束させることができる。そのとき、埋め込み方法の時刻 での相対（局所）誤差 は次の公式で与えられる。
ルンゲ＝クッタ法の収束性から、相対誤差も0に収束することがわかる。埋め込みルンゲ＝クッタ法は、アルゴリズムを用いて一時刻ごと（自動的）に刻み幅 を調整し、誤差をコントロールすることができる（故に適応型）。よって絶対誤差を計算せずに刻み幅を正しく設定することができる。そのアルゴリズムは、大体以下の通りである。
誤差の許容値を とする。刻み幅 の近似値 と をルンゲ＝クッタ法で計算する。二つの値に対し、formula_53 が成立するとき、 が大きすぎて小さくする必要がある。よって新しい刻み幅 formula_54 と設定し再び近似値を計算する。代わりに formula_55 が成立するとき、 が小さすぎて大きくするほうが効率良いである。よって新しい刻み幅 formula_56 と設定し、次の時刻での近似値を計算する（この場合、誤差が許容値より小さいので今の近似値を二度と計算する必要はない）。二つの不等式ともが成立しないとき、そのままで次の時刻での近似値を計算するといいである。
ただし、使用したスカラー（1/2他）は方程や方法によって変更することも可能である。また、実践では を小さくしすぎないようにするため、下界 を先に設定することもある。formula_57 のときにアルゴリズムを停止し、代わりにもっと高い次数を持つ方法を使う。
もっとも単純な埋め込み型ルンゲ＝クッタ法は、ホイン法（2次）とオイラー法（1次）を組み合わせる方法である。以下の拡張ブッチャー配列で与えられる。
そして実践でよく使われているのが5次と4次の陽的方法を組み合わせる である。対応する拡張ブッチャー配列は以下の通りである。
他に使われている方法は、(次数3と2)、 と （両方次数5と4）である。
いままでに述べた方法はすべて陽公式である。陽的ルンゲ＝クッタ方法の絶対安定性領域（region of absolute stability）が小さくて有界であるため、硬い方程式の解を計算する場合、方法は不適切である。この問題は特に偏微分方程式の解を計算するときに重要である。
陽的方法の不安定さは陰的ルンゲ＝クッタ法の開発の動機となる。陰的ルンゲ＝クッタ法は上述の一般的な公式と同じく、以下の形で与えられる
但し、
であり、対応するルンゲ＝クッタ行列は厳密な下三角行列ではない。結果として、一時刻ごとに代数方程式系を解かなければならなくなる。応じて、計算コストもかなり上がる。 段法を使って 個の成分からなる微分方程式系（すなわち formula_60 のとき）に適用する場合に対応する代数方程式の数は となる。これは陰的線型多段法とも比較できる： 段陰的線型多段法を使って同じ方程式系に適用する場合に対応する代数方程式の数は であり、段数 に依存しない。
一番単純な陰的ルンゲ＝クッタ法は、後退オイラー法
である。対応するブッチャー配列も単純に以下の通りである。
他の例として、台形公式と呼ばれる方法は以下の公式で与えられる。
対応する配列は以下の通りである。
その対応は、formula_65 を覚えてて、公式を以下のように書き換えると明らかになる。
台形公式は、選点法である。選点法はすべて陰的ルンゲ＝クッタ法であるけど、陰的ルンゲ＝クッタ法がすべて選点法であるわけではない。
選点法の中では、ガウス求積に基づいたが一番次数の高い方法である。 段ガウス・ルジャンドル法の次数は となる（よって、任意に高い次数を持つ方法を構造できるようになる）。例として、2段ガウス・ルジャンドル法は次の配列で与えられる。
数値分析における安定性は、それぞれ異なる定義が複数存在する。ルンゲ＝クッタ法の安定性を反映できる概念は、主に以下の二つである。
線型テスト方程式 formula_68 を考える。一つのルンゲ＝クッタ法を使ってこの方程式に適用すると formula_69 となる。ここで、
は安定性関数と呼ばれる 上の有理関数である（ はすべての成分が 1 のベクトルである）。 段法の場合、行列式の展開によって は二つの 次多項式の商となる。陽的方法の場合、対応するルンゲ＝クッタ行列が狭義下三角行列であるため、formula_71 がわかる。したがって、 は 次多項式となる。
ルンゲ＝クッタ法による上記の方程式の数値解がゼロに減衰する条件が formula_72 (formula_73) である。上記の条件を満たす の集合は方法に対する 絶対安定性領域 （region of absolute stability） である。 特に、絶対安定性領域が左複素数平面を含むとき、その方法は A-安定 (A-stable) という。陽的ルンゲ＝クッタ法は、安定性関数が多項式であるため、A-安定な方法になれない。
もっと一般的に、方法の次数が のときに、安定性関数に対し formula_74 が成立する。ゆえに、指数関数 の、与えられた次数の多項式の商からなる有理関数の中での最善近似が重要だと考えられる。 そのような有理関数はと呼ばれる。分子の次数が で分母の次数が のパデ近似式がA-安定性の条件 formula_72, formula_76 を満足するための必要十分条件は、formula_77 である。
以上のことから、陰的ルンゲ＝クッタ法は陽的方法より優れる安定性を持つこともわかる。
A-安定性という概念は線型自励方程式 formula_68 を考える結果である。 は、数値的方法を使って、とある単調性条件を満たす非線型方程式系に適用するときの安定性を主張した。対応する安定性は、線型多段法の場合に G-安定性 （G-stability）で、ルンゲ＝クッタ法の場合に B-安定性 (B-stability) と呼ぶ。 一般的なテスト方程 formula_80 と単調性条件 formula_81 を満足する を考える。もしすべての formula_82 に対し、 不等式 formula_83 が成立するとき、そのルンゲ＝クッタ法は B-安定 という。 ここで、 と はそれぞれの初期値に対する数値解である。 formula_84 に方法を適用することで、B-安定性はA-安定性より強い条件であることがわかる。