ルンゲ・クッタ法

微分方程式のルンゲ・クッタ法の比較（赤は正確な解）${\displaystyle y'=\sin ^{2}(t)\cdot y}$

数値解析では、ルンゲ・クッタ法(英語: / ˈ 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 0}$
${\displaystyle c_{2}}$	${\displaystyle a_{21}}$
${\displaystyle c_{3}}$	${\displaystyle a_{31}}$	${\displaystyle a_{32}}$
${\displaystyle \vdots}$	${\displaystyle \vdots}$		${\displaystyle \ddots }$
${\displaystyle c_{s}}$	${\displaystyle a_{s1}}$	${\displaystyle a_{s2}}$	${\displaystyle \cdots }$	${\displaystyle a_{s,s-1}}$
	${\displaystyle b_{1}}$	${\displaystyle b_{2}}$	${\displaystyle \cdots }$	${\displaystyle b_{s-1}}$	${\displaystyle b_{s}}$

テイラー級数展開は、ルンゲ・クッタ法が矛盾しない条件が、

${\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 の明示的方法の例は、文献 [14] に記載されています。9^{段階の順序}7 の明示的方法^{[ 11 ]}と 11 段階の順序 8 の明示的方法^{[ 15 ]}も知られています。要約については 文献^{[ 16 ] [ 17 ]を参照してください。}${\displaystyle p=1,2,\ldots ,6}$

例

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

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

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

0
1/3	1/3
2/3	−1/3	1
1	1	−1	1
	1/8	3/8	3/8	1/8

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

0
	1

2段階の2次法

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

対応するタブローは

0
1/2	1/2
	0	1

^{中点法は、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 )}.}$

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

0
${\displaystyle \alpha }$	${\displaystyle \alpha }$
	${\displaystyle (1-{\tfrac {1}{2\alpha }})}$	${\displaystyle {\tfrac {1}{2\alpha }}}$

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

使用

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

	0
	2/3	2/3
		1/4	3/4

対応する方程式

${\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 つのステップを実行する必要があります。

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

${\displaystyle t_{0}=1\colon }$
	${\displaystyle y_{0}=1}$
${\displaystyle t_{1}=1.025\colon }$
	${\displaystyle y_{0}=1}$	${\displaystyle k_{1}=2.557407725}$	${\displaystyle k_{2}=f(t_{0}+{\tfrac {2}{3}}h,\ y_{0}+{\tfrac {2}{3}}hk_{1})=2.7138981400}$
	${\displaystyle y_{1}=y_{0}+h({\tfrac {1}{4}}k_{1}+{\tfrac {3}{4}}k_{2})={\underline {1.066869388}}}$
${\displaystyle t_{2}=1.05\colon }$
	${\displaystyle y_{1}=1.066869388}$	${\displaystyle k_{1}=2.813524695}$	${\displaystyle k_{2}=f(t_{1}+{\tfrac {2}{3}}h,\ y_{1}+{\tfrac {2}{3}}hk_{1})}$
	${\displaystyle y_{2}=y_{1}+h({\tfrac {1}{4}}k_{1}+{\tfrac {3}{4}}k_{2})={\underline {1.141332181}}}$
${\displaystyle t_{3}=1.075\colon }$
	${\displaystyle y_{2}=1.141332181}$	${\displaystyle k_{1}=3.183536647}$	${\displaystyle k_{2}=f(t_{2}+{\tfrac {2}{3}}h,\ y_{2}+{\tfrac {2}{3}}hk_{1})}$
	${\displaystyle y_{3}=y_{2}+h({\tfrac {1}{4}}k_{1}+{\tfrac {3}{4}}k_{2})={\underline {1.227417567}}}$
${\displaystyle t_{4}=1.1\colon }$
	${\displaystyle y_{3}=1.227417567}$	${\displaystyle k_{1}=3.796866512}$	${\displaystyle k_{2}=f(t_{3}+{\tfrac {2}{3}}h,\ y_{3}+{\tfrac {2}{3}}hk_{1})}$
	${\displaystyle y_{4}=y_{3}+h({\tfrac {1}{4}}k_{1}+{\tfrac {3}{4}}k_{2})={\underline {1.335079087}}.}$

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

暗黙的ルンゲ・クッタ法

明示的ルンゲ・クッタ法は、絶対安定領域が狭く、特に有界であるため、一般に硬い方程式の解法には適していません。 ^{[ 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 ]}これはブッチャー表にも現れている。陽解法の係数行列は下三角行列である。陰解法では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}\\\end{array}}={\begin{array}{c|c}\mathbf {c} &A\\\hline &\mathbf {b^{T}} \\\end{array}}}$

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

例

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

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

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

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

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

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

ガウス・ルジャンドル法は、ガウス積分法に基づく選点法の一種である。s段階のガウス・ルジャンドル法は、次数が2である（したがって、任意の高次数を持つ方法を構築できる）。^{[ 26 ]} 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}}}$ ^{[ 24 ]}

安定性

陽的ルンゲ・クッタ法に対する陰的ルンゲ・クッタ法の利点は、特に硬い方程式に適用した場合の安定性の高さである。線形テスト方程式 を考えてみよう。この方程式にルンゲ・クッタ法を適用すると、反復回数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)}},}$ ^{[ 27 ]}

ここで、eは1のベクトルを表す。関数rは安定性関数と呼ばれる。^{[ 28 ]}式から、rは2つのs次多項式の商であり、その方法はs段階からなる。陽解法ではAは厳密に下三角行列であるため、det( I − zA ) = 1となり、安定性関数は多項式となる。^{[ 29 ]}

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

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

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

適応ルンゲ・クッタ法

適応型法は、単一のルンゲ・クッタステップの局所打ち切り誤差の推定値を生成するように設計されています。これは、 次数と 次数 の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 つの方法があります。その拡張ブッチャー表は次のとおりです。

	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

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

	0
	1	1
		1/2	1/2
		1	0

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

非合流型ルンゲ・クッタ法

ルンゲ・クッタ法は、すべての が異なる場合、非合流性であると言われる^{[ 33 ]}。 ${\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つのルンゲ・クッタ行列と2組の重みを使用しますが、必要なノードは1組だけです。これにより、次の形式のブッチャー表が生成されます。 ${\displaystyle a_{ij},a'_{ij}}$ ${\displaystyle b_{i},b'_{i}}$ ${\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}}$

ルンゲ・クッタ・ニストロム法は、 と が厳密に下三角関数であるとき陽的であるとされ、この場合、の式における和は^{[ 36 ]}に置き換えられる。さらに、と の局所的打ち切り誤差が であるとき、ルンゲ・クッタ・ニストロム法は の順序であると言われる。 ${\displaystyle A,A'}$ ${\textstyle \sum _{j=1}^{s}}$ ${\displaystyle g_{i},g'_{i}}$ ${\textstyle \sum _{j=1}^{i-1}}$ ${\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})}$$

あるスカラー関数の場合。^{[ 37 ]}${\displaystyle V}$

B安定性

微分方程式の解における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つの行列 とする 。ルンゲ・クッタ法は、行列 と がともに非負定値であるとき、代数的に安定であると言われる[ 38 ] 。B^安定性の十分条件^{[ 39 ]}は、と がともに非負定値であることである。 ${\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}$

ルンゲ・クッタ4次法の導出

一般に、ルンゲ・クッタ法は次のように記述できます。 ${\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}}}$

参照

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

注記

1. 「ルンゲ・クッタ法」 . Dictionary.com . 2021年4月4日閲覧。
2. DEVRIES, Paul L.; HASBUN, Javier E. 『計算物理学入門』第2版. Jones and Bartlett Publishers, 2011. p. 215.
3. プレス他。 2007 年、p. 908; Süli & Mayers 2003、p. 328
4. Atkinson (1989 , p. 423)、Hairer, Nørsett & Wanner (1993 , p. 134)、Kaw & Kalu (2008 , §8.4)、Stoer & Bulirsch (2002 , p. 476)は、段階の定義において係数hを省略している。Ascher & Petzold (1998 , p. 81)、Butcher (2008 , p. 93)、Iserles (1996 , p. 38)は、y値を段階として用いている。
5. Süli & Mayers 2003、p. 328
6. プレス他 2007年、907ページ
7. イザールズ 1996、38ページ
8. イザールズ 1996、39ページ
9. 反例として、 と をランダムに選択した任意の明示的2段階ルンゲ・クッタ法を考えてみましょう。この方法は矛盾がなく、（一般に）1次収束します。一方、 を とする1段階法は矛盾しており、 が自明に成立しているにもかかわらず収束しません。${\displaystyle b_{1}=b_{2}=1/2}$ ${\displaystyle c_{1}}$ ${\displaystyle a_{21}}$ ${\displaystyle b_{1}=1/2}$ ${\displaystyle \sum _{j=1}^{i-1}a_{ij}=c_{i}{\text{ for }}i=2,\ldots ,s.}$
10. ブッチャー 2008、187ページ
11. ブッチャー 1965年、408ページ
12. ブッチャー 1985
13. ブッチャー 2008、187～196ページ
14. ブッチャー 1964
15. カーティス 1970、268ページ
16. ハイラー、ノーセット、ワナー、1993 年、p. 179
17. ブッチャー 1996、247ページ
18. Süli & Mayers 2003、p. 352
19. Hairer、Nørsett & Wanner (1993、p. 138) はKutta (1901)を参照。
20. Süli & Mayers 2003、p. 327
21. Süli & Mayers 2003、pp. 349–351
22. イゼルレス、1996 年、p. 41; Süli & Mayers 2003、pp. 351–352
23. ブッチャー 2008、94ページ
24. Süli & Mayers 2003、p. 353
25. イゼルレス 1996、43～44ページ
26. イザールズ 1996、47ページ
27. ヘアラー＆ワナー 1996、40–41 ページ
28. ヘアラー＆ワナー 1996、p. 40
29. イザールズ 1996、60ページ
30. イゼルレス 1996、62～63ページ
31. イゼルレス 1996、63ページ
32. この結果はDahlquist (1963)によるものです。
33. ランバート 1991、278ページ
34. Dormand, JR; Prince, PJ (1978年10月). 「力学天文学における数値シミュレーションのための新しいルンゲ・クッタアルゴリズム」.天体力学. 18 (3): 223– 232. Bibcode : 1978CeMec..18..223D . doi : 10.1007/BF01230162 . S2CID 120974351 . 
35. Fehlberg, E. (1974年10月). 一般2階微分方程式に対するステップサイズ制御による古典的7次、6次、5次ルンゲ・クッタ・ニストローム公式（報告書）（NASA TR R-432版）. マーシャル宇宙飛行センター、アラバマ州：アメリカ航空宇宙局.
36. ブッチャー 2008、94ページ
37. Qin, Meng-Zhao; Zhu, Wen-Jie (1991-01-01). 「2次常微分方程式に対する正準ルンゲ・クッタ・ニストローム（RKN）法」 . Computers & Mathematics with Applications . 22 (9): 85– 95. doi : 10.1016/0898-1221(91)90209-M . ISSN 0898-1221 . 
38. ランバート 1991、275ページ
39. ランバート 1991、274ページ
40. Lyu, Ling-Hsiao (2016年8月). 「付録C. 数値積分公式の導出」(PDF) .宇宙プラズマの数値シミュレーション(I)講義ノート. 国立中央大学宇宙科学研究所. 2022年4月17日閲覧。

参考文献

• 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。
• Iserles, Arieh (1996), 『微分方程式の数値解析入門』、ケンブリッジ大学出版局、Bibcode : 1996fcna.book.....I、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) 「Section 17.1 Runge-Kutta Method」、Numerical Recipes: The Art of Scientific Computing (第3版)、Cambridge University Press、ISBN 978-0-521-88068-8また、セクション17.2「ルンゲ・クッタ法の適応ステップサイズ制御」も参照してください。
• ステア、ジョセフ。 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)、数学科学・数学教育ジャーナル、7 (2): 1– 10。
• 離散数学上級 ignou 参考書（コード：mcs033）
• ジョン・C・ブッチャー：「Bシリーズ：数値解析法の代数的解析」、Springer（SSCM、第55巻）、ISBN 978-3030709556（2021年4月）。
• Butcher, 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。
• Curtis, AR (1970)、「ステップごとに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。
• ブッチャー, JC (1996)、「ルンゲ・クッタ法の歴史」、応用数値数学、20 (3): 247– 260、doi : 10.1016/0168-9274(95)00108-5。

外部リンク

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