ルンゲ・クッタ法

数値解析では、ルンゲ・クッタ法(英語: / ˈ 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}$それ自体。最初の時点で$t0$対応する${\displaystyle y}$値は$\displaystyle y_0}}$. 機能${\displaystyle f}$そして初期条件$t0$、 $\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}+k_{1}{\frac {h}{2}}\right),\\k_{3}&=\ f\!\left(t_{n}+{\frac {h}{2}},y_{n}+k_{2}{\frac {h}{2}}\right),\\k_{4}&=\ f\!\left(t_{n}+h,y_{n}+hk_{3}\right).\end{aligned}}}$

（注：上記の数式は、文献によって定義は異なるが、同等の意味を持つ。^{[ 4 ]}）

ここ$\displaystyle y_{n+1}}$RK4近似は${\displaystyle y(t_{n+1})}$、次の値（$\displaystyle y_{n+1}}$）は現在の価値によって決定されます（$\displaystyle y_n}}$) に加えて、4 つの増分の加重平均を加えます。ここで、各増分は、区間の大きさhと、微分方程式の右辺の関数fによって指定される推定勾配の積です。

• ${\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つの傾斜を平均化する場合、中間点の傾斜に重みが大きくなります。${\displaystyle f}$独立している${\displaystyle y}$微分方程式が単純な積分と等価である場合、RK4 はシンプソンの公式である。^{[ 5 ]}

RK4法は4次法であり、局所的な打ち切り誤差はオーダーである。${\displaystyle O(h^{5})}$累積誤差の合計は、 ${\displaystyle O(h^{4})}$。

多くの実用的なアプリケーションでは、${\displaystyle f}$独立している${\displaystyle t}$（いわゆる自律系、または時間不変系、特に物理学において）であり、その増分は全く計算されず、関数に渡されない。${\displaystyle f}$最終的な式は$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 ]}これらのデータは通常、ブッチャー表( John C. Butcherにちなんで)として知られる記憶術で整理されます。

テイラー級数展開により、ルンゲ・クッタ法が一致性を持つのは、

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

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

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

しかし、この条件だけでは、一貫性の十分条件でも必要条件でもない。 ^{[ 9 ]}

一般的に、明示的な${\displaystyle s}$-段階ルンゲ・クッタ法は次数を持つ${\displaystyle p}$すると、ステージ数は次の条件を満たす必要があることが証明できる。${\displaystyle s\geq p}$そしてもし${\displaystyle p\geq 5}$、 それから${\displaystyle s\geq p+1}$^{[ 10 ]}しかし、これらの境界がすべての場合に厳密 であるかどうかは不明である。場合によっては、境界が達成できないことが証明されている。例えば、Butcher は、の場合に次のことを証明した。${\displaystyle p>6}$明示的なメソッドはありません${\displaystyle s=p+1}$段階。^{[ 11 ]}ブッチャーはまた、${\displaystyle p>7}$明示的なルンゲ・クッタ法はありません${\displaystyle p+2}$段階。^{[ 12 ]}しかし一般的には、正確な最小段階数がいくつであるかは未解決の問題である。${\displaystyle s}$明示的なルンゲ・クッタ法の次数が${\displaystyle p}$既知の値には以下のようなものがあります: ^{[ 13 ]}

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

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

例として、α = 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つのステップを実行する必要があります。

この方法は以下の手順で進められます。

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

明示的ルンゲ・クッタ法は、絶対安定領域が小さいため、一般的に硬い方程式の解法には適していません。特に、その領域は有界です。^{[ 21 ]}この問題は、偏微分方程式 の解法において特に重要です。

明示的ルンゲ・クッタ法の不安定性は、暗黙的方法の開発を促す。暗黙的ルンゲ・クッタ法は次の形式をとる。

${\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.}$^{[ 22 ]}

明示的な方法との違いは、明示的な方法では、jに関する和はi − 1までしか取られないという点である。 ^{[ 23 ]}これはブッチャー表にも表れる。係数行列${\displaystyle a_{ij}}$明示的方法の場合、係数行列は下三角行列になります。暗黙的方法では、jに関する和はsまでになり、係数行列は厳密には三角行列ではないため、 ^{[ 18 ]}の形式のブッチャー表が得られます。

${\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}\\\end{array}}={\begin{array}{c|c}\mathbf {c} &A\\\hline &\mathbf {b^{T}} \\\end{array}}}$

この違いの結果、各ステップで代数方程式系を解く必要が生じます。これにより計算コストが大幅に増加します。m個の成分を持つ微分方程式を解くためにs段階の方法を使用する場合、代数方程式系はms個の成分を持ちます。これは、暗黙的線形多段階法（常微分方程式のもう1つの大きな方法群）とは対照的です。暗黙的s段階線形多段階法では、 m個の成分のみを持つ代数方程式系を解く必要があるため、ステップ数が増加してもシステムのサイズは増加しません。^{[ 24 ]}

適応法は、単一のルンゲ・クッタステップの局所打ち切り誤差の推定値を生成するように設計されています。これは、次数を持つ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},}$

それは${\displaystyle O(h^{p})}$この種の方法のブッチャー表は、次の値を与えるように拡張されます。${\displaystyle b_{i}^{*}}$:

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

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

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

その他の適応型ルンゲ・クッタ法としては、ボガッキ・シャンパイン法（次数3および2）、キャッシュ・カープ法、ドーマンド・プリンス法（いずれも次数5および4）などがある。

^ルンゲ・クッタ法は^、すべて^の${\displaystyle c_{i},\,i=1,2,\ldots ,s}$それらは異なる。

ルンゲ・クッタ・ニューストロム（RKN）法は、ルンゲ・クッタ法と同じ原理に基づく手法群ですが、2階の初期値問題^{[ 34 ] [ 35 ]}を対象としており、したがって、次の形式の問題に対応します 。

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

導関数は2つあり、近似値も2つあるため、ルンゲ・クッタ・ニューストロム法では2つのルンゲ・クッタ行列を使用します。${\displaystyle a_{ij},a'_{ij}}$、そして2組の重り${\displaystyle b_{i},b'_{i}}$しかし、それでも必要なノードのセットは1つだけです${\displaystyle c_{i}}$これにより、次の形式のブッチャーテーブルが生成されます 。

$${\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 &a'_{11}&a'_{12}&\dots &a'_{1s}\\&a'_{21}&a'_{22}&\dots &a'_{2s}\\&\vdots &\vdots &\ddots &\vdots \\&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} &\mathbf {A} \\\hline &\mathbf {A'} \\\hline &\mathbf {b} ^{\top }\\&\mathbf {b'} ^{\top }\end{array}}}$$

近似は、${\displaystyle t_{n}}$、 と${\displaystyle y_{n}}$近似値${\displaystyle y(t_{n})}$そして${\displaystyle y'_{n}}$近似値${\displaystyle {\frac {dy}{dt}}(t_{n})}$近似値${\displaystyle y_{n+1},y'_{n+1}}$で${\displaystyle t_{n+1}=t_{n}+h}$以下の連立方程式の解は次のとおりです 。

$${\displaystyle {\begin{cases}g_{i}=y_{n}+c_{i}hy'_{n}+h^{2}\sum _{j=1}^{s}a_{ij}f(t_{n}+c_{j}h,g'_{j},g_{j}),&i=1,2,\ldots ,s\\g'_{i}=y'_{n}+h\sum _{j=1}^{s}a'_{ij}f(t_{n}+c_{j}h,g'_{j},g_{j}),&i=1,2,\ldots ,s\\\\y_{n+1}=y_{n}+hy'_{n}+h^{2}\sum _{j=1}^{s}b_{j}f(t_{n}+c_{j}h,g'_{j},g_{j})\\y'_{n+1}=y'_{n}+h\sum _{j=1}^{s}b'_{j}f(t_{n}+c_{j}h,g'_{j},g_{j})\end{cases}}}$$

どこ${\displaystyle g_{i},g'_{i}}$は中間近似値である${\displaystyle y}$そして${\displaystyle {\frac {dy}{dt}}}$値を扱うことは厳密に同等です${\displaystyle k_{j}=f(t_{n}+c_{j}h,g'_{j},g_{j})}$どこで${\displaystyle g_{j},g'_{j}}$彼らの公式に置き換えられた代わりに、${\displaystyle g_{i},g'_{i}}$これは、以前にルンゲ・クッタ法で行ったことと同様ですが、この方法の方がシステムを記述しやすくなっています。

ルンゲ・クッタ・ニューストロム法は、以下の条件を満たす場合に明示的であると言われる。${\displaystyle A,A'}$厳密には下三角行列であり、この場合、和は${\textstyle \sum _{j=1}^{s}}$表現において${\displaystyle g_{i},g'_{i}}$、は以下に置き換えられる場合があります${\textstyle \sum _{j=1}^{i-1}}$^{[ 36 ]}さらに、ルンゲ・クッタ・ニューストロム法は次数であると言われています。${\displaystyle p}$両方の局所的切り捨て誤差が${\displaystyle y_{n+1},y'_{n+1}}$は${\displaystyle O(h^{p+1})}$。

関数が${\displaystyle f}$検討対象の初期値問題の${\displaystyle {\frac {dy}{dt}}}$中間値を近似する必要はありません${\displaystyle g'_{i}}$近似値を計算するには、重み${\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} &\mathbf {A} \\\hline &\mathbf {b} ^{\top }\\&\mathbf {b'} ^{\top }\end{array}}}$$

この特殊なケースは、一般的にルンゲ・クッタ・ニューストロム法が達成できる以上の高次の次数を実現できるため、特に興味深い。例えば、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 b_{i}&{\frac {3-2{\sqrt {3}}}{12}}&{\frac {1}{2}}&{\frac {3+2{\sqrt {3}}}{12}}\\\hline b'_{i}&{\frac {5-3{\sqrt {3}}}{24}}&{\frac {3+{\sqrt {3}}}{12}}&{\frac {1+{\sqrt {3}}}{24}}\\\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 b_{i}&{\frac {3+2{\sqrt {3}}}{12}}&{\frac {1}{2}}&{\frac {3-2{\sqrt {3}}}{12}}\\\hline b'_{i}&{\frac {5+3{\sqrt {3}}}{24}}&{\frac {3-{\sqrt {3}}}{12}}&{\frac {1-{\sqrt {3}}}{24}}\\\end{array}}}$$

これらの2つのスキームは、元の方程式が保存的な古典力学システムから導出されている場合、つまり

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

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

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

させて${\displaystyle B}$、${\displaystyle M}$そして${\displaystyle Q}$3つ${\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}}}$$ ルンゲ・クッタ法は、行列が代数的に安定していると言われている^{[ 38 ]。}${\displaystyle B}$そして${\displaystyle M}$両方とも非負定値である。B安定性^{[ 39 ]}の十分条件は次のとおりである。${\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次法の導出^{[ 40 ]を展開する。}${\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}}}$上記のとおりです。

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

• 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", Zeitschrift für Mathematik und Physik, 46: 435–453.
• Ascher, Uri M.; Petzold, Linda R. (1998), Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations, Philadelphia: Society for Industrial and Applied Mathematics, ISBN 978-0-89871-412-8.
• 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. (1963年5月)、「ルンゲ・クッタ積分過程の研究のための係数」、Journal of the Australian Mathematical Society、3 (2): 185–201、doi : 10.1017/S1446788700027932。
• Butcher, John C. (1964年5月)、「高次のルンゲ・クッタ過程について」、Journal of the Australian Mathematical Society、4 (2): 179–194、doi : 10.1017/S1446788700023387
• Butcher, John C. (1975)、「陰的ルンゲ・クッタ法の安定性特性」、BIT、15 (4): 358–361、doi : 10.1007/bf01931672、S2CID 120854166 。
• Butcher, John 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、ISSN 0006-3835、S2CID 120241743  。
• フォーサイス、ジョージ・E.、マルコム、マイケル・A.、モーラー、クリーブ・B. (1977)、『数理計算のためのコンピュータ手法』、プレンティス・ホール（第6章を参照）。
• ハイラー、エルンスト。ノーセット、シベール・ポール。 Wanner、Gerhard (1993)、常微分方程式の解法 I: Nonstiff 問題、ベルリン、ニューヨーク: 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 (1996), 『微分方程式の数値解析入門』、ケンブリッジ大学出版局、Bibcode : 1996fcna.book.....I、ISBN 978-0-521-55655-2。
• ランバート、JD (1991)、常微分方程式系の数値解法：初期値問題、ジョン・ワイリー・アンド・サンズ、ISBN 0-471-92990-5
• Kaw, Autar; Kalu, Egwu (2008), Numerical Methods with Applications (1st  ed.), 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 版)、Cambridge University Press、ISBN 978-0-521-88068-8また、セクション 17.2. Runge-Kutta の適応ステップサイズ制御。
• ステア、ジョセフ。 Bulirsch、Roland (2002)、数値解析入門(第 3 版)、ベルリン、ニューヨーク: Springer-Verlag、ISBN 978-0-387-95452-3。
• Süli, Endre; Mayers, David (2003),数値解析入門, Cambridge University Press , ISBN 0-521-00794-1。
• Tan, Delin; Chen, Zheng (2012)、「4次ルンゲ・クッタ法の一般公式について」（PDF）、Journal of Mathematical Science & Mathematics Education、7（2）：1–10。
• IGNOU高度離散数学参考書（コード：MCS033）
• ジョン・C・ブッチャー著「Bシリーズ ：数値解析の代数的分析」、シュプリンガー（SSCM、第55巻）、ISBN 978-3030709556（2021年4月）
• Butcher, JC (1985)、「10段階8次陽的ルンゲ・クッタ法の非存在」、BIT Numerical Mathematics、25 (3): 521–540、doi : 10.1007/BF01935372。
• Butcher, JC (1965)、「ルンゲ・クッタ法の達成可能な次数について」、Mathematics of Computation、19 (91): 408–417、doi : 10.1090/S0025-5718-1965-0179943-X。
• Curtis, AR (1970)、「1ステップあたり11回の関数評価を行う8次ルンゲ・クッタ法」、Numerische Mathematik、16 (3): 268–277、doi : 10.1007/BF02219778。
• Cooper, GJ; Verner, JH (1972)、「高次の明示的ルンゲ・クッタ法のいくつか」、SIAM Journal on Numerical Analysis、9 (3): 389–405、Bibcode : 1972SJNA....9..389C、doi : 10.1137/0709037。
• Butcher, JC (1996)、「ルンゲ・クッタ法の歴史」、応用数値数学、20 (3): 247–260、doi : 10.1016/0168-9274(95)00108-5。