ルンゲ・クッタ法

数値解析では、ルンゲ・クッタ法(英語: / ˈ r ʊ ŋ ə ˈ k ʊ t ɑː / ^ⓘ RUUNG-ə-KUUT-tah^[1]、暗黙的および明示的な反復法のファミリーであり同時非線形方程式の近似解の時間的離散化に使用されるオイラー法。^[2]これらの方法は、1900年頃にドイツの数学者カール・ルンゲとヴィルヘルム・クッタ。

ルンゲ・クッタ法ファミリーの最も広く知られているメンバーは、一般に「RK4」、「古典的なルンゲ・クッタ法」、または単に「ルンゲ・クッタ法」と呼ばれています。

初期値問題を次のように指定します。

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

ここでは、近似したい時間 の未知の関数 (スカラーまたはベクトル) があります。が変化する速度 は、および自身の関数であるとされています。 初期時間では、対応する値は です。 関数と初期条件、 が与えられています。 ${\displaystyle y}$${\displaystyle t}$${\displaystyle {\frac {dy}{dt}}}$${\displaystyle y}$${\displaystyle t}$${\displaystyle y}$${\displaystyle t_{0}}$${\displaystyle y}$${\displaystyle y_{0}}$${\displaystyle f}$ ${\displaystyle t_{0}}$${\displaystyle y_{0}}$

ここで、ステップサイズh > 0 を選択し、次のように定義します。

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

n = 0, 1, 2, 3, ...の場合、 ^[3]

${\displaystyle {\begin{aligned}k_{1}&=\ f(t_{n},y_{n}),\\k_{2}&=\ f\!\left(t_{n}+{\frac {h}{2}},y_{n}+h{\frac {k_{1}}{2}}\right),\\k_{3}&=\ f\!\left(t_{n}+{\frac {h}{2}},y_{n}+h{\frac {k_{2}}{2}}\right),\\k_{4}&=\ f\!\left(t_{n}+h,y_{n}+hk_{3}\right).\end{aligned}}}$

（注：上記の式は、異なるテキストでは異なるが同等の定義を持っています。^[4]）

ここではの RK4 近似値であり、次の値 ( ) は現在の値 ( ) と4 つの増分の加重平均によって決定されます。ここで、各増分は、区間のサイズhと、微分方程式の右側の 関数fによって指定される推定傾きの積です。${\displaystyle y_{n+1}}$${\displaystyle y(t_{n+1})}$${\displaystyle y_{n+1}}$${\displaystyle y_{n}}$

• ${\displaystyle k_{1}}$は、 （オイラー法）を使用した区間の開始時の傾きです。${\displaystyle y}$
• ${\displaystyle k_{2}}$は、およびを使用した区間の中点における傾きです。${\displaystyle y}$${\displaystyle k_{1}}$
• ${\displaystyle k_{3}}$は再び中間点の傾きですが、今度はと を使用します。${\displaystyle y}$${\displaystyle k_{2}}$
• ${\displaystyle k_{4}}$は、およびを使用した、区間の終了時の傾きです。${\displaystyle y}$${\displaystyle k_{3}}$

4つの傾きを平均化する際に、中点の傾きに大きな重みが与えられます。 がと独立であり、微分方程式が単純な積分と等価である場合、RK4はシンプソンの定理です。^[5]${\displaystyle f}$${\displaystyle y}$

RK4 法は 4 次法であり、局所的な切り捨て誤差はのオーダーで あるのに対し、 総累積誤差は のオーダーであることを意味します。 ${\displaystyle O(h^{5})}$${\displaystyle O(h^{4})}$

多くの実際のアプリケーションでは、関数は から独立しており(いわゆる自律システム、または特に物理学における時間不変システム)、その増分はまったく計算されず、関数 に渡されず、 の最終的な式のみが使用されます。 ${\displaystyle f}$${\displaystyle t}$${\displaystyle f}$${\displaystyle t_{n+1}}$

明示的ルンゲ・クッタ法は、上述のRK4法の一般化であり、次のように表される。

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

ここで^[6]

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

（注：上記の式は、テキストによっては異なるが同等の定義を持つ場合があります。^[4]）

特定の方法を指定するには、整数s (ステージ数)、係数a _ij (1 ≤ j < i ≤ sの場合)、b _i ( i = 1, 2, ..., s の場合)、c _i ( i = 2, 3, ..., sの場合) を指定する必要があります。行列 [ a _ij ] はルンゲ・クッタ行列と呼ばれ、b _iとc _{i は}重み、ノードと呼ばれます。^{[7]これらのデータは通常、}ブッチャー・タブロー(ジョン・C・ブッチャーにちなんで)と呼ばれる記憶術で配置されます。

テイラー級数展開によれば、ルンゲ・クッタ法が矛盾しないのは、

${\displaystyle \sum _{i=1}^{s}b_{i}=1.}$

また、この方法に特定の次数p を要求する場合、局所的打ち切り誤差が O( h ^{p +1} )になるという付随要件もあります。これらは、打ち切り誤差の定義自体から導き出すことができます。たとえば、2 段階法では、b ₁ + b ₂ = 1、b ₂ c ₂ = 1/2、b ₂ a ₂₁ = 1/2 の場合、次数が 2 になります。^[8]係数を決定するための一般的な条件は、^[8]です。

${\displaystyle \sum _{j=1}^{i-1}a_{ij}=c_{i}{\text{ for }}i=2,\ldots ,s.}$

しかし、この条件だけでは一貫性を保つには十分でも必要でもありません。 ^[9]

一般に、明示的段ルンゲ・クッタ法が 次数 を持つ場合、段数は を満たす必要があり、の場合 であることが証明できます。^{[10]ただし、これらの境界がすべての場合に}明確であるかどうか はわかっていません。 場合によっては、境界を達成できないことが証明されています。 たとえば、Butcher は、 に対して、段数を持つ明示的方法は存在しないことを証明しました。 ^[11] Butcher はまた、 に対して、段数 を持つ明示的ルンゲ・クッタ法は存在しないことを証明しました。^[12]ただし、一般に、明示的ルンゲ・クッタ法が 次数 を持つための正確な最小段数が何であるかは未解決の問題です。 既知の値には次のものがあります。^[13]${\displaystyle s}$${\displaystyle p}$${\displaystyle s\geq p}$${\displaystyle p\geq 5}$${\displaystyle s\geq p+1}$${\displaystyle p>6}$${\displaystyle s=p+1}$${\displaystyle p>7}$${\displaystyle p+2}$${\displaystyle s}$${\displaystyle p}$

${\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}}}$

上記の証明可能な境界は、これらの順序について既に知られている方法よりも少ない段階を必要とする順序の方法を見つけることができないことを意味します。Butcher の研究は、7 次と 8 次の方法にはそれぞれ最小で 9 段階と 11 段階があることも証明しています。^{[11] [12]} 7 段階の順序 6 の明示的方法の例は、Ref. ^[14]にあります。9段階の順序 7 の明示的方法^[11]と 11 段階の順序 8 の明示的方法^[15]も知られています。要約については、 Ref. ^{[16] [17]を参照してください。}${\displaystyle p=1,2,\ldots ,6}$

RK4法はこの枠組みに当てはまる。その表は^[18]

ルンゲ・クッタ法のわずかなバリエーションも1901年にクッタによって考案され、3/8ルールと呼ばれています。^[19]この方法の主な利点は、ほとんどすべての誤差係数が一般的な方法よりも小さいことですが、時間ステップごとにわずかに多くのFLOP（浮動小数点演算）が必要です。そのブッチャー・タブローは

しかし、最も単純なルンゲ・クッタ法は（順方向）オイラー法であり、式 で与えられる。これは、1段階の唯一の一貫した明示的ルンゲ・クッタ法である。対応する表は、 ${\displaystyle y_{n+1}=y_{n}+hf(t_{n},y_{n})}$

2段階の2次法の例は、明示的な中点法によって提供されます。

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

対応するタブローは

中点法は、2段階の2次ルンゲ・クッタ法の唯一の方法ではない。このような方法は、αによってパラメータ化され、式^{[20]で与えられるファミリーが存在する。}

${\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 )}.}$

そのブッチャータブローは

この族では、中点法を与え、はHeun法、^[5]であり、はRalston法である。 ${\displaystyle \alpha ={\tfrac {1}{2}}}$${\displaystyle \alpha =1}$${\displaystyle \alpha ={\tfrac {2}{3}}}$

例として、α = 2/3の2段階2次ルンゲ・クッタ法（ラルストン法とも呼ばれる）を考えてみましょう。これは次の表で与えられます。

対応する方程式

${\displaystyle {\begin{aligned}k_{1}&=f(t_{n},\ y_{n}),\\k_{2}&=f(t_{n}+{\tfrac {2}{3}}h,\ y_{n}+{\tfrac {2}{3}}hk_{1}),\\y_{n+1}&=y_{n}+h\left({\tfrac {1}{4}}k_{1}+{\tfrac {3}{4}}k_{2}\right).\end{aligned}}}$

この方法は初期値問題を解くために使用される

${\displaystyle {\frac {dy}{dt}}=\tan(y)+1,\quad y_{0}=1,\ t\in [1,1.1]}$

ステップサイズh = 0.025 なので、この方法では 4 つのステップを実行する必要があります。

方法は次のように進行します。

数値解は下線部の値に対応します。

適応型法は、単一のルンゲ・クッタステップのローカル打ち切り誤差の推定値を生成するように設計されています。これは、次数と次数 の2 つの方法を使用して行われます。これらの方法は相互に絡み合っており、つまり、共通の中間ステップがあります。これにより、高次法のステップと比較して、誤差の推定にかかる計算コストは​​ほとんどまたは無視できるほどになります。 ${\displaystyle p}$${\displaystyle p-1}$

積分中、ステップ サイズは、推定誤差がユーザー定義のしきい値を下回るように調整されます。誤差が大きすぎる場合は、ステップ サイズを小さくしてステップを繰り返します。誤差がはるかに小さい場合は、時間を節約するためにステップ サイズを大きくします。これにより、(ほぼ) 最適なステップ サイズが得られ、計算時間が節約されます。さらに、ユーザーは適切なステップ サイズを見つけるために時間を費やす必要がありません。

低次のステップは次のように与えられる。

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

ここで、高階法の場合と同じです。すると、誤差は ${\displaystyle k_{i}}$

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

これは です。この種の方法の Butcher タブローは の値を与えるように拡張されます。 ${\displaystyle O(h^{p})}$${\displaystyle b_{i}^{*}}$

ルンゲ・クッタ・フェルベルグ法には、 5 次と 4 次の 2 つの方法があります。その拡張されたブッチャー表は次のとおりです。

しかし、最も単純な適応型ルンゲ・クッタ法では、次数 2 のホイン法と次数 1 の オイラー法を組み合わせる。その拡張ブッチャー表は次のようになる。

その他の適応型ルンゲ・クッタ法には、ボガッキ・シャンピン法(次数 3 および 2)、キャッシュ・カープ法、ドルマンド・プリンス法(いずれも次数 5 および 4) があります。

ルンゲ・クッタ法は、すべてのが異なっている場合、非合流性であると言われる ^[21]。 ${\displaystyle c_{i},\,i=1,2,\ldots ,s}$

ルンゲ・クッタ・ニストローム法は、2次微分方程式に最適化された特殊なルンゲ・クッタ法である。^{[22] [23]} 2次常微分方程式系に対する一般的なルンゲ・クッタ・ニストローム法

$${\displaystyle {\ddot {y}}_{i}=f_{i}(y_{1},y_{2},\ldots ,y_{n})}$$

注文はフォームで ${\displaystyle s}$

$${\displaystyle {\begin{cases}g_{i}=y_{m}+c_{i}h{\dot {y}}_{m}+h^{2}\sum _{j=1}^{s}a_{ij}f(g_{j}),&i=1,2,\ldots ,s\\y_{m+1}=y_{m}+h{\dot {y}}_{m}+h^{2}\sum _{j=1}^{s}{\bar {b}}_{j}f(g_{j})\\{\dot {y}}_{m+1}={\dot {y}}_{m}+h\sum _{j=1}^{s}b_{j}f(g_{j})\end{cases}}}$$

これは、ブッチャーテーブルの形状を形成します

$${\displaystyle {\begin{array}{c|cccc}c_{1}&a_{11}&a_{12}&\dots &a_{1s}\\c_{2}&a_{21}&a_{22}&\dots &a_{2s}\\\vdots &\vdots &\vdots &\ddots &\vdots \\c_{s}&a_{s1}&a_{s2}&\dots &a_{ss}\\\hline &{\bar {b}}_{1}&{\bar {b}}_{2}&\dots &{\bar {b}}_{s}\\&b_{1}&b_{2}&\dots &b_{s}\end{array}}={\begin{array}{c|c}\mathbf {c} &\mathbf {A} \\\hline &\mathbf {\bar {b}} ^{\top }\\&\mathbf {b} ^{\top }\end{array}}}$$

次の Butcher テーブルには、2 つの 4 次明示的 RKN 法が示されています。

$${\displaystyle {\begin{array}{c|ccc}c_{i}&&a_{ij}&\\{\frac {3+{\sqrt {3}}}{6}}&0&0&0\\{\frac {3-{\sqrt {3}}}{6}}&{\frac {2-{\sqrt {3}}}{12}}&0&0\\{\frac {3+{\sqrt {3}}}{6}}&0&{\frac {\sqrt {3}}{6}}&0\\\hline {\overline {b_{i}}}&{\frac {5-3{\sqrt {3}}}{24}}&{\frac {3+{\sqrt {3}}}{12}}&{\frac {1+{\sqrt {3}}}{24}}\\\hline b_{i}&{\frac {3-2{\sqrt {3}}}{12}}&{\frac {1}{2}}&{\frac {3+2{\sqrt {3}}}{12}}\end{array}}}$$ $${\displaystyle {\begin{array}{c|ccc}c_{i}&&a_{ij}&\\{\frac {3-{\sqrt {3}}}{6}}&0&0&0\\{\frac {3+{\sqrt {3}}}{6}}&{\frac {2+{\sqrt {3}}}{12}}&0&0\\{\frac {3-{\sqrt {3}}}{6}}&0&-{\frac {\sqrt {3}}{6}}&0\\\hline {\overline {b_{i}}}&{\frac {5+3{\sqrt {3}}}{24}}&{\frac {3-{\sqrt {3}}}{12}}&{\frac {1-{\sqrt {3}}}{24}}\\\hline b_{i}&{\frac {3+2{\sqrt {3}}}{12}}&{\frac {1}{2}}&{\frac {3-2{\sqrt {3}}}{12}}\end{array}}}$$

これら2つのスキームは、元の方程式が保存的な古典力学システムから導かれる場合、すなわち、

$${\displaystyle f_{i}(x_{1},\ldots ,x_{n})={\frac {\partial V}{\partial x_{i}}}(x_{1},\ldots ,x_{n})}$$

あるスカラー関数に対して^[24]${\displaystyle V}$

これまで述べたルンゲ・クッタ法はすべて陽的解法である。陽的ルンゲ・クッタ法は絶対安定領域が小さく、特に有界であるため、一般に硬い方程式の解法には適していない。 ^{[25]この問題は}偏微分方程式 の解法において特に重要である。

陽的ルンゲ・クッタ法の不安定性は、陰的ルンゲ・クッタ法の開発の動機となった。陰的ルンゲ・クッタ法は、次の形式を持つ。

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

どこ

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

明示的方法との違いは、明示的方法ではjの合計がi − 1までしか上がらないことです。これはブッチャー表にも現れています。明示的方法の係数行列は下三角です。暗黙的方法では、 jの合計はsまでで、係数行列は三角ではないため、ブッチャー表は次のようになります^[18]。${\displaystyle a_{ij}}$

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

行の説明については、上記の適応型ルンゲ・クッタ法を参照してください。 ${\displaystyle b^{*}}$

この違いの結果、各ステップで代数方程式系を解かなければなりません。これにより、計算コストが大幅に増加します。s 段階の手法を使用してm個の要素を持つ微分方程式を解くと、代数方程式系にはms 個の要素が含まれます。これは、暗黙的線形多段階法(常微分方程式のもう 1 つの大きな手法ファミリー)とは対照的です。暗黙的sステップ線形多段階法では、 m個の要素のみを持つ代数方程式系を解く必要があるため、ステップ数が増えてもシステムのサイズは増加しません。^[27]

暗黙的ルンゲ・クッタ法の最も単純な例は、後退オイラー法である。

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

これに対するブッチャーのタブローは単純です:

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

このブッチャーのタブローは、次の式に対応しています。

${\displaystyle k_{1}=f(t_{n}+h,\ y_{n}+hk_{1})\quad {\text{and}}\quad y_{n+1}=y_{n}+hk_{1},}$

これを変形すると、上記の後退オイラー法の式が得られます。

暗黙的ルンゲ・クッタ法のもう1つの例は台形則です。そのブッチャー・タブローは次のようになります。

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

台形則は選点法である（その記事で議論されているように）。選点法はすべて暗黙的ルンゲ・クッタ法であるが、すべての暗黙的ルンゲ・クッタ法が選点法であるわけではない。^[28]

ガウス・ルジャンドル法は、ガウス積分法に基づく選点法のファミリーを形成します。 s段階のガウス・ルジャンドル法の次数は2です（したがって、任意の高次数を持つ方法を構築できます）。^[29] 2段階（したがって次数は 4）の方法にはブッチャー・タブローがあります。

${\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}}\\&{\frac {1}{2}}+{\frac {1}{2}}{\sqrt {3}}&{\frac {1}{2}}-{\frac {1}{2}}{\sqrt {3}}\end{array}}}$ ^[27]

明示的ルンゲ・クッタ法に対する暗黙的ルンゲ・クッタ法の利点は、特に硬い方程式に適用した場合に安定性が高いことです。線形テスト方程式を考えてみましょう。この方程式にルンゲ・クッタ法を適用すると、反復rが次のよう に表されます。${\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)}},}$ ^[30]

ここで、eは1のベクトルを表す。関数rは安定性関数と呼ばれる。^[31]この式から、この方法がs段階を持つ場合、 rは2つのs次多項式の商であることがわかる。明示的方法では、厳密に下三角行列Aを持ち、これはdet( I − zA )=1であり、安定性関数が多項式であることを意味する。 ^[32]

線形検定方程式の数値解は、z = h λ で | r ( z ) | < 1の場合にゼロに減少する。このようなzの集合は絶対安定領域と呼ばれる。特に、 Re( z ) < 0となるすべてのz が絶対安定領域内にある場合、この方法は絶対安定であると言われる。陽的ルンゲ・クッタ法の安定性関数は多項式であるため、陽的ルンゲ・クッタ法が A 安定になることはない。^[32]

方法の次数がpの場合、安定関数は を満たす。したがって、指数関数を最もよく近似する、与えられた次数の多項式の商を調べることは興味深い。これらはパデ近似として知られている。分子がm次、分母がn 次であるパデ近似は、 m ≤ n ≤ m + 2の場合にのみ A 安定である。 ^[33]${\displaystyle r(z)={\textrm {e}}^{z}+O(z^{p+1})}$${\displaystyle z\to 0}$

s段階のガウス・ルジャンドル法は次数が2 sなので、その安定性関数はm = n = sのパデ近似となる。したがって、この方法はA安定である。^{[34]これは、A安定ルンゲ・クッタ法が任意の高次数を持つことができることを示している。対照的に、A安定}線形多段階法の次数は2を超えることはできない。^[35]

微分方程式の解の A 安定性の概念は、線形自律方程式 と関連しています。Dahlquist ( 1963) は、単調性条件を満たす非線形システムに適用した場合の数値スキームの安定性の調査を提案しました。対応する概念は、マルチステップ法 (および関連するワンレッグ法) の場合はG 安定性、ルンゲ・クッタ法の場合はB 安定性(Butcher、1975) として定義されました。を証明する非線形システム に適用されたルンゲ・クッタ法は、この条件が2 つの数値解を 意味する場合、 B 安定と呼ばれます。${\displaystyle y'=\lambda y}$${\displaystyle y'=f(y)}$${\displaystyle \langle f(y)-f(z),\ y-z\rangle \leq 0}$${\displaystyle \|y_{n+1}-z_{n+1}\|\leq \|y_{n}-z_{n}\|}$

、およびを、によって定義される3つの行列 とします。 ルンゲ・クッタ法は、行列およびが両方とも非負定値である場合に、代数的に安定しているといわれます^{[36] 。B}安定性の十分条件^[37]は次のとおりです。およびは非負定値です。 ${\displaystyle B}$${\displaystyle M}$${\displaystyle Q}$${\displaystyle s\times s}$$${\displaystyle {\begin{aligned}B&=\operatorname {diag} (b_{1},b_{2},\ldots ,b_{s}),\\[4pt]M&=BA+A^{T}B-bb^{T},\\[4pt]Q&=BA^{-1}+A^{-T}B-A^{-T}bb^{T}A^{-1}.\end{aligned}}}$$${\displaystyle B}$${\displaystyle M}$${\displaystyle B}$${\displaystyle Q}$

一般に、ルンゲ・クッタ法は次のように記述できます。 ${\displaystyle s}$

${\displaystyle y_{t+h}=y_{t}+h\cdot \sum _{i=1}^{s}a_{i}k_{i}+{\mathcal {O}}(h^{s+1}),}$

どこ：

${\displaystyle k_{i}=\sum _{j=1}^{s}\beta _{ij}f(k_{j},\ t_{n}+\alpha _{i}h)}$

は、次の微分値を評価して得られる増分です。 ${\displaystyle y_{t}}$${\displaystyle i}$

上で説明したように、任意の区間の開始点、中点、終了点で評価された一般式を使用して、ルンゲ・クッタ4次法の導出^[38]を展開します。したがって、次のように選択します。 ${\displaystyle s=4}$${\displaystyle (t,\ t+h)}$

${\displaystyle {\begin{aligned}&\alpha _{i}&&\beta _{ij}\\\alpha _{1}&=0&\beta _{21}&={\frac {1}{2}}\\\alpha _{2}&={\frac {1}{2}}&\beta _{32}&={\frac {1}{2}}\\\alpha _{3}&={\frac {1}{2}}&\beta _{43}&=1\\\alpha _{4}&=1&&\\\end{aligned}}}$

そうでなければ、次の量を定義することから始めます。 ${\displaystyle \beta _{ij}=0}$

${\displaystyle {\begin{aligned}y_{t+h}^{1}&=y_{t}+hf\left(y_{t},\ t\right)\\y_{t+h}^{2}&=y_{t}+hf\left(y_{t+h/2}^{1},\ t+{\frac {h}{2}}\right)\\y_{t+h}^{3}&=y_{t}+hf\left(y_{t+h/2}^{2},\ t+{\frac {h}{2}}\right)\end{aligned}}}$

ここで、 次のように定義します。 ${\displaystyle y_{t+h/2}^{1}={\dfrac {y_{t}+y_{t+h}^{1}}{2}}}$${\displaystyle y_{t+h/2}^{2}={\dfrac {y_{t}+y_{t+h}^{2}}{2}}.}$

${\displaystyle {\begin{aligned}k_{1}&=f(y_{t},\ t)\\k_{2}&=f\left(y_{t+h/2}^{1},\ t+{\frac {h}{2}}\right)=f\left(y_{t}+{\frac {h}{2}}k_{1},\ t+{\frac {h}{2}}\right)\\k_{3}&=f\left(y_{t+h/2}^{2},\ t+{\frac {h}{2}}\right)=f\left(y_{t}+{\frac {h}{2}}k_{2},\ t+{\frac {h}{2}}\right)\\k_{4}&=f\left(y_{t+h}^{3},\ t+h\right)=f\left(y_{t}+hk_{3},\ t+h\right)\end{aligned}}}$

また、前の関係式に対して、次の等式が まで成り立つことが示せます。 ここで、 はの時間に関する 全微分です。${\displaystyle {\mathcal {O}}(h^{2})}$$${\displaystyle {\begin{aligned}k_{2}&=f\left(y_{t+h/2}^{1},\ t+{\frac {h}{2}}\right)=f\left(y_{t}+{\frac {h}{2}}k_{1},\ t+{\frac {h}{2}}\right)\\&=f\left(y_{t},\ t\right)+{\frac {h}{2}}{\frac {d}{dt}}f\left(y_{t},\ t\right)\\k_{3}&=f\left(y_{t+h/2}^{2},\ t+{\frac {h}{2}}\right)=f\left(y_{t}+{\frac {h}{2}}f\left(y_{t}+{\frac {h}{2}}k_{1},\ t+{\frac {h}{2}}\right),\ t+{\frac {h}{2}}\right)\\&=f\left(y_{t},\ t\right)+{\frac {h}{2}}{\frac {d}{dt}}\left[f\left(y_{t},\ t\right)+{\frac {h}{2}}{\frac {d}{dt}}f\left(y_{t},\ t\right)\right]\\k_{4}&=f\left(y_{t+h}^{3},\ t+h\right)=f\left(y_{t}+hf\left(y_{t}+{\frac {h}{2}}k_{2},\ t+{\frac {h}{2}}\right),\ t+h\right)\\&=f\left(y_{t}+hf\left(y_{t}+{\frac {h}{2}}f\left(y_{t}+{\frac {h}{2}}f\left(y_{t},\ t\right),\ t+{\frac {h}{2}}\right),\ t+{\frac {h}{2}}\right),\ t+h\right)\\&=f\left(y_{t},\ t\right)+h{\frac {d}{dt}}\left[f\left(y_{t},\ t\right)+{\frac {h}{2}}{\frac {d}{dt}}\left[f\left(y_{t},\ t\right)+{\frac {h}{2}}{\frac {d}{dt}}f\left(y_{t},\ t\right)\right]\right]\end{aligned}}}$$$${\displaystyle {\frac {d}{dt}}f(y_{t},\ t)={\frac {\partial }{\partial y}}f(y_{t},\ t){\dot {y}}_{t}+{\frac {\partial }{\partial t}}f(y_{t},\ t)=f_{y}(y_{t},\ t){\dot {y}}_{t}+f_{t}(y_{t},\ t):={\ddot {y}}_{t}}$$${\displaystyle f}$

今導出した式を使って一般式を表現すると次のようになります。$${\displaystyle {\begin{aligned}y_{t+h}={}&y_{t}+h\left\lbrace a\cdot f(y_{t},\ t)+b\cdot \left[f(y_{t},\ t)+{\frac {h}{2}}{\frac {d}{dt}}f(y_{t},\ t)\right]\right.+\\&{}+c\cdot \left[f(y_{t},\ t)+{\frac {h}{2}}{\frac {d}{dt}}\left[f\left(y_{t},\ t\right)+{\frac {h}{2}}{\frac {d}{dt}}f(y_{t},\ t)\right]\right]+\\&{}+d\cdot \left[f(y_{t},\ t)+h{\frac {d}{dt}}\left[f(y_{t},\ t)+{\frac {h}{2}}{\frac {d}{dt}}\left[f(y_{t},\ t)+\left.{\frac {h}{2}}{\frac {d}{dt}}f(y_{t},\ t)\right]\right]\right]\right\rbrace +{\mathcal {O}}(h^{5})\\={}&y_{t}+a\cdot hf_{t}+b\cdot hf_{t}+b\cdot {\frac {h^{2}}{2}}{\frac {df_{t}}{dt}}+c\cdot hf_{t}+c\cdot {\frac {h^{2}}{2}}{\frac {df_{t}}{dt}}+\\&{}+c\cdot {\frac {h^{3}}{4}}{\frac {d^{2}f_{t}}{dt^{2}}}+d\cdot hf_{t}+d\cdot h^{2}{\frac {df_{t}}{dt}}+d\cdot {\frac {h^{3}}{2}}{\frac {d^{2}f_{t}}{dt^{2}}}+d\cdot {\frac {h^{4}}{4}}{\frac {d^{3}f_{t}}{dt^{3}}}+{\mathcal {O}}(h^{5})\end{aligned}}}$$

これをの周りのテイラー級数と比較すると次のようになります。${\displaystyle y_{t+h}}$${\displaystyle t}$$${\displaystyle {\begin{aligned}y_{t+h}&=y_{t}+h{\dot {y}}_{t}+{\frac {h^{2}}{2}}{\ddot {y}}_{t}+{\frac {h^{3}}{6}}y_{t}^{(3)}+{\frac {h^{4}}{24}}y_{t}^{(4)}+{\mathcal {O}}(h^{5})=\\&=y_{t}+hf(y_{t},\ t)+{\frac {h^{2}}{2}}{\frac {d}{dt}}f(y_{t},\ t)+{\frac {h^{3}}{6}}{\frac {d^{2}}{dt^{2}}}f(y_{t},\ t)+{\frac {h^{4}}{24}}{\frac {d^{3}}{dt^{3}}}f(y_{t},\ t)\end{aligned}}}$$

係数に対する制約のシステムが得られます。

${\displaystyle {\begin{cases}&a+b+c+d=1\\[6pt]&{\frac {1}{2}}b+{\frac {1}{2}}c+d={\frac {1}{2}}\\[6pt]&{\frac {1}{4}}c+{\frac {1}{2}}d={\frac {1}{6}}\\[6pt]&{\frac {1}{4}}d={\frac {1}{24}}\end{cases}}}$

これを解くと上記のようになります。 ${\displaystyle a={\frac {1}{6}},b={\frac {1}{3}},c={\frac {1}{3}},d={\frac {1}{6}}}$

• オイラー法
• ルンゲ・クッタ法の一覧
• 常微分方程式の数値解析法
• ルンゲ・クッタ法 (SDE)
• 一般線形法
• リー群積分器

• Runge、Carl David Tolmé (1895)、「Über die numerische Auflösung von Differentialgleichungen」、Mathematische Annalen、46 (2)、Springer : 167–178、doi :10.1007/BF01446807、S2CID  119924854。
• Kutta、Wilhelm (1901)、「Beitrag zur näherungsweisen Integration totaler Differentialgleichungen」、数学と物理学に関する研究、46 : 435–453。
• Ascher, Uri M.; Petzold, Linda R. (1998)、「常微分方程式と微分代数方程式のコンピュータ手法」、フィラデルフィア：産業応用数学協会、ISBN 978-0-89871-412-8。
• アトキンソン、ケンドール A. (1989)、数値解析入門（第 2 版）、ニューヨーク：ジョン ワイリー アンド サンズ、ISBN 978-0-471-50023-0。
• ブッチャー、ジョン C. (1963 年 5 月)、「ルンゲ・クッタ積分過程の研究のための係数」、オーストラリア数学会誌、3 (2): 185–201、doi : 10.1017/S1446788700027932。
• ブッチャー、ジョン C. (1964 年 5 月)、「高次のルンゲ・クッタ過程について」、オーストラリア数学会誌、4 (2): 179–194、doi : 10.1017/S1446788700023387
• ブッチャー、ジョン C. (1975)、「暗黙的ルンゲ・クッタ法の安定性特性」、BIT、15 (4): 358–361、doi :10.1007/bf01931672、S2CID  120854166。
• ブッチャー、ジョン C. (2000)、「20 世紀における常微分方程式の数値解析法」、J. Comput. Appl. Math.、125 (1–2): 1–29、Bibcode :2000JCoAM.125....1B、doi : 10.1016/S0377-0427(00)00455-6。
• ブッチャー、ジョン C. (2008)、常微分方程式の数値解析法、ニューヨーク：ジョン ワイリー アンド サンズ、ISBN 978-0-470-72335-7。
• セリエ、F. Kofman, E. (2006)、連続システム シミュレーション、Springer Verlag、ISBN 0-387-26102-8。
• Dahlquist, Germund (1963)、「線形多段階法の特殊な安定性問題」、BIT、3 :27–43、doi :  10.1007 /BF01963532、hdl : 10338.dmlcz/103497、ISSN0006-3835  、S2CID120241743。
• フォーサイス、ジョージ E.; マルコム、マイケル A.; モーラー、クリーブ B. (1977)、数学的計算のためのコンピュータ手法、プレンティス・ホール（第6章を参照）。
• ハイラー、エルンスト。ノーセット、シベール・ポール。 Wanner、Gerhard (1993)、常微分方程式の解法 I: Nonstiff 問題、ベルリン、ニューヨーク: Springer-Verlag、ISBN 978-3-540-56670-0。
• ヘアラー、エルンスト、ワナー、ゲルハルト（1996）、常微分方程式の解法 II：スティフおよび微分代数問題（第 2 版）、ベルリン、ニューヨーク：シュプリンガー・フェアラーク、ISBN 978-3-540-60452-5。
• イザールズ、アリエ（1996）、微分方程式の数値解析入門、ケンブリッジ大学出版局、ISBN 978-0-521-55655-2。
• ランバート、JD (1991)、常微分システムの数値解析法。初期値問題、John Wiley & Sons、ISBN 0-471-92990-5
• Kaw, Autar; Kalu, Egwu (2008)、数値計算法とその応用 (第 1 版)、autarkaw.com。
• Press, William H.; Teukolsky, Saul A .; Vetterling, William T.; Flannery, Brian P. (2007)、「セクション 17.1 ルンゲ・クッタ法」、Numerical Recipes: The Art of Scientific Computing (第 3 版)、ケンブリッジ大学出版局、ISBN 978-0-521-88068-8また、セクション 17.2「ルンゲ・クッタ法の適応ステップサイズ制御」も参照してください。
• ステア、ジョセフ。 Bulirsch、Roland (2002)、数値解析入門(第 3 版)、ベルリン、ニューヨーク: Springer-Verlag、ISBN 978-0-387-95452-3。
• スーリ、エンドレ、メイヤーズ、デイヴィッド（2003）、数値解析入門、ケンブリッジ大学出版局、ISBN 0-521-00794-1。
• Tan, Delin; Chen, Zheng (2012)、「第 4 次ルンゲ・クッタ法の一般公式について」(PDF)、数学科学と数学教育ジャーナル、7 (2): 1–10。
• 高度な離散数学 ignou 参考書 (コード - mcs033)
• John C. Butcher：「Bシリーズ：数値手法の代数的解析」、Springer（SSCM、第55巻）、ISBN 978-3030709556（2021年4月）。 
• ブッチャー、JC (1985)、「10 段階 8 次明示的ルンゲ・クッタ法の非存在」、BIT 数値数学、25 (3): 521–540、doi :10.1007/BF01935372。
• ブッチャー、JC (1965)、「ルンゲ・クッタ法の達成可能な順序について」、計算数学、19 (91): 408–417、doi :10.1090/S0025-5718-1965-0179943-X。
• カーティス、AR（1970）、「ステップごとに11回の関数評価を伴う8次ルンゲ・クッタ過程」、Numerische Mathematik、16（3）：268–277、doi：10.1007 / BF02219778。
• クーパー、GJ; ヴァーナー、JH (1972)、「高次の明示的ルンゲ・クッタ法」、SIAM 数値解析ジャーナル、9 (3): 389–405、doi :10.1137/0709037。
• ブッチャー、JC（1996）、「ルンゲ・クッタ法の歴史」、応用数値数学、20（3）：247–260、doi：10.1016/0168-9274（95）00108-5。

• 「ルンゲ・クッタ法」、数学百科事典、EMS Press、2001 [1994]
• ルンゲ・クッタ 4 次法
• Matlab での Tracker コンポーネント ライブラリの実装 — に 32 個の埋め込み Runge Kutta アルゴリズムRungeKStep、に 24 個の埋め込み Runge-Kutta Nyström アルゴリズムRungeKNystroemSStep、に 4 個の汎用 Runge-Kutta Nyström アルゴリズムを実装しますRungeKNystroemGStep。