ルンゲ=クッタ法 は、以下の形の常微分方程式 の初期値問題 の解を数値で近似計算する方法 である。
y ′ = f ( t , y ) , y ( t 0 ) = y 0 {\displaystyle y'=f(t,y),\;y(t_{0})=y_{0}} 一般的に、ルンゲ=クッタ法は以下の形で与えられる。
y n + 1 = y n + h ∑ i = 1 s b i k i {\displaystyle y_{n+1}=y_{n}+h\sum _{i=1}^{s}b_{i}k_{i}\,} ただし、
k i = f ( t n + c i h , y n + h ∑ j = 1 i − 1 a i j k j ) , {\displaystyle k_{i}=f\left(t_{n}+c_{i}h,y_{n}+h\sum _{j=1}^{i-1}a_{ij}k_{j}\right),} 以下のリストで記述するすべての計算方法は、それに対応するブッチャー配列で与えられる。ある一つの方法に対する係数をブッチャー配列で以下の形で表す。
c 1 a 11 a 12 … a 1 s c 2 a 21 a 22 … a 2 s ⋮ ⋮ ⋮ ⋱ ⋮ c s a s 1 a s 2 … a s s b 1 b 2 … b s {\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}}} また、陽的ルンゲ=クッタ法に対応するルンゲ=クッタ行列は狭義の下三角行列であるので上三角成分の表記は省略される。
オイラー法は1次の方法である。 安定性と精度が低いため、オイラー法は入門の例でしか使われない(実際は大規模な計算では今でも使われることがある)。
0 1 {\displaystyle {\begin{array}{c|c}0&\\\hline &1\\\end{array}}} (陽的)中点法 は2段2次の方法である(陰的中点法も参照)。
0 1 / 2 1 / 2 0 1 {\displaystyle {\begin{array}{c|cc}0&&\\1/2&1/2&\\\hline &0&1\\\end{array}}} ホイン法 (英語版 ) も2段2次の方法である (陽的台形公式としても知られている)。
0 1 1 1 / 2 1 / 2 {\displaystyle {\begin{array}{c|cc}0&&\\1&1&\\\hline &1/2&1/2\\\end{array}}} Ralston法 は2段2次の方法のうちで局所誤差の上界が最小のものである[1] 。
0 2 / 3 2 / 3 1 / 4 3 / 4 {\displaystyle {\begin{array}{c|cc}0&&\\2/3&2/3&\\\hline &1/4&3/4\\\end{array}}} 0 x x 1 − 1 2 x 1 2 x {\displaystyle {\begin{array}{c|ccc}0&&\\x&x&\\\hline &1-{\frac {1}{2x}}&{\frac {1}{2x}}\\\end{array}}} 0 1 / 2 1 / 2 1 − 1 2 1 / 6 2 / 3 1 / 6 {\displaystyle {\begin{array}{c|ccc}0&&&\\1/2&1/2&&\\1&-1&2&\\\hline &1/6&2/3&1/6\\\end{array}}} [2] 0 1 / 2 1 / 2 1 / 2 0 1 / 2 1 0 0 1 1 / 6 1 / 3 1 / 3 1 / 6 {\displaystyle {\begin{array}{c|cccc}0&&&&\\1/2&1/2&&&\\1/2&0&1/2&&\\1&0&0&1&\\\hline &1/6&1/3&1/3&1/6\\\end{array}}} この方法は、上記の古典的方法と同じ論文で提出されたが[3] 、古典的方法に比べるとあまり用いられていない。
0 1 / 3 1 / 3 2 / 3 − 1 / 3 1 1 1 − 1 1 1 / 8 3 / 8 3 / 8 1 / 8 {\displaystyle {\begin{array}{c|cccc}0&&&&\\1/3&1/3&&&\\2/3&-1/3&1&&\\1&1&-1&1&\\\hline &1/8&3/8&3/8&1/8\\\end{array}}} 埋め込み型の方法はルンゲ=クッタ法の局所誤差を推定するために開発された方法である。それらの方法は誤差を制御するために刻み幅を調整する。
埋め込み型方法に対応するブッチャー配列は以下のように与えられる。
c 1 a 11 a 12 … a 1 s c 2 a 21 a 22 … a 2 s ⋮ ⋮ ⋮ ⋱ ⋮ c s a s 1 a s 2 … a s s b 1 b 2 … b s b 1 ∗ b 2 ∗ … b s ∗ {\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}}} ここで、上側の段の係数 bi は p 次陽的方法に対応するものであり、下側の段の係数 b * i は p -1 次陽的方法に対応するものである。
この方法は、2次のホイン法と1次のオイラー法を組み合わせる方法であり、もっとも単純な埋め込み型方法である。
0 1 1 1 / 2 1 / 2 1 0 {\displaystyle {\begin{array}{c|cc}0&\\1&1\\\hline &1/2&1/2\\&1&0\end{array}}} フェールベルグ法は3段で、次数2と1の方法を用いる[4] 。
0 1 / 2 1 / 2 1 1 / 256 255 / 256 1 / 512 255 / 256 1 / 512 1 / 256 255 / 256 0 {\displaystyle {\begin{array}{c|ccc}0&\\1/2&1/2\\1&1/256&255/256\\\hline &1/512&255/256&1/512\\&1/256&255/256&0\end{array}}} Bogacki–Shampine法 (英語版 ) は4段で、次数3と2の方法を用いる[5] 。MATLABのコマンド ode23 はこの方法の実装である[6] 。
0 1 / 2 1 / 2 3 / 4 0 3 / 4 1 2 / 9 1 / 3 4 / 9 2 / 9 1 / 3 4 / 9 0 7 / 24 1 / 4 1 / 3 1 / 8 {\displaystyle {\begin{array}{c|ccc}0&\\1/2&1/2\\3/4&0&3/4\\1&2/9&1/3&4/9\\\hline &2/9&1/3&4/9&0\\&7/24&1/4&1/3&1/8\end{array}}} ルンゲ=クッタ=フェールベルグ法 は5段で、次数5と4の方法を用いる[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 16 / 135 0 6656 / 12825 28561 / 56430 − 9 / 50 2 / 55 25 / 216 0 1408 / 2565 2197 / 4104 − 1 / 5 0 {\displaystyle {\begin{array}{c|cccccc}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\\\hline &16/135&0&6656/12825&28561/56430&-9/50&2/55\\&25/216&0&1408/2565&2197/4104&-1/5&0\end{array}}} Cash-Karp法 (英語版 ) はフェールベルグの最初のアイディアを変型した方法である[7] 。フェールベルグの方法と同じく5段で、次数5と4の方法を用いる。
0 1 / 5 1 / 5 3 / 10 3 / 40 9 / 40 3 / 5 3 / 10 − 9 / 10 6 / 5 1 − 11 / 54 5 / 2 − 70 / 27 35 / 27 37 / 378 0 250 / 621 125 / 594 0 512 / 1771 2825 / 27648 0 18575 / 48384 13525 / 55296 277 / 14336 1 / 4 {\displaystyle {\begin{array}{c|cccccc}0&\\1/5&1/5\\3/10&3/40&9/40\\3/5&3/10&-9/10&6/5\\1&-11/54&5/2&-70/27&35/27\\\hline &37/378&0&250/621&125/594&0&512/1771\\&2825/27648&0&18575/48384&13525/55296&277/14336&1/4\end{array}}} ドルマン=プリンス法 は6段で、次数5と4の方法を用いる[8] 。MATLABのコマンド ode45 はこの方法を実装したものである[6] 。
0 1 / 5 1 / 5 3 / 10 3 / 40 9 / 40 4 / 5 44 / 45 − 56 / 15 32 / 9 8 / 9 19372 / 6561 − 25360 / 2187 64448 / 6561 − 212 / 729 1 9017 / 3168 − 355 / 33 46732 / 5247 49 / 176 − 5103 / 18656 35 / 384 0 500 / 1113 125 / 192 − 2187 / 6784 11 / 84 0 5179 / 57600 0 7571 / 16695 393 / 640 − 92097 / 339200 187 / 2100 1 / 40 {\displaystyle {\begin{array}{c|ccccccc}0&\\1/5&1/5\\3/10&3/40&9/40\\4/5&44/45&-56/15&32/9\\8/9&19372/6561&-25360/2187&64448/6561&-212/729\\1&9017/3168&-355/33&46732/5247&49/176&-5103/18656\\\hline &35/384&0&500/1113&125/192&-2187/6784&11/84&0\\&5179/57600&0&7571/16695&393/640&-92097/339200&187/2100&1/40\end{array}}} 上述の埋め込み型方法と同じく、ブッチャー配列に二つの方法が含まれている場合、表で重ねて書かれた下側の段の係数の方法が誤差をコントロールするためのものとなる。
後退オイラー法 は1次の方法である。 この方法は、偏微分方程式である線型拡散方程式の時間方向の離散化に用いた場合には無条件に安定で非振動的な方法である。
1 1 1 {\displaystyle {\begin{array}{c|c}1&1\\\hline &1\\\end{array}}} 陰的中点法は2次方法である。選点法 であり、以下のガウス・ルジャンドル法 (英語版 ) の最も簡単な場合である。
1 / 2 1 / 2 1 {\displaystyle {\begin{array}{c|c}1/2&1/2\\\hline &1\end{array}}} これらの方法はガウス求積 法に基づいた方法であり、高い次数を持つ(s 段ガウス・ルジャンドル法の次数は 2s である)。
4次の方法は以下のブッチャー配列で与えられる[9] 。
1 2 − 3 6 1 4 1 4 − 3 6 1 2 + 3 6 1 4 + 3 6 1 4 1 2 1 2 1 2 + 1 2 3 1 2 − 1 2 3 {\displaystyle {\begin{array}{c|cc}{\frac {1}{2}}-{\frac {\sqrt {3}}{6}}&{\frac {1}{4}}&{\frac {1}{4}}-{\frac {\sqrt {3}}{6}}\\{\frac {1}{2}}+{\frac {\sqrt {3}}{6}}&{\frac {1}{4}}+{\frac {\sqrt {3}}{6}}&{\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}}} さらに6次の方法に対応する配列は以下で与えられる[9] 。
1 2 − 15 10 5 36 2 9 − 15 15 5 36 − 15 30 1 2 5 36 + 15 24 2 9 5 36 − 15 24 1 2 + 15 10 5 36 + 15 30 2 9 + 15 15 5 36 5 18 4 9 5 18 − 5 6 8 3 − 5 6 {\displaystyle {\begin{array}{c|ccc}{\frac {1}{2}}-{\frac {\sqrt {15}}{10}}&{\frac {5}{36}}&{\frac {2}{9}}-{\frac {\sqrt {15}}{15}}&{\frac {5}{36}}-{\frac {\sqrt {15}}{30}}\\{\frac {1}{2}}&{\frac {5}{36}}+{\frac {\sqrt {15}}{24}}&{\frac {2}{9}}&{\frac {5}{36}}-{\frac {\sqrt {15}}{24}}\\{\frac {1}{2}}+{\frac {\sqrt {15}}{10}}&{\frac {5}{36}}+{\frac {\sqrt {15}}{30}}&{\frac {2}{9}}+{\frac {\sqrt {15}}{15}}&{\frac {5}{36}}\\\hline &{\frac {5}{18}}&{\frac {4}{9}}&{\frac {5}{18}}\\&-{\frac {5}{6}}&{\frac {8}{3}}&-{\frac {5}{6}}\end{array}}} Lobatto法[10] は主に IIIA、 IIIB と IIIC(古典的な文献によって、記号 I と II は二種類のRadau法にのみ使われる)と呼ばれる三種類の方法を指している。方法の名称は Rehuel Lobatto にちなむ。それらの方法はすべて陰的であり、次数 2s -2 を持ち、係数に対し条件 c 1 = 0 と cs = 1 を満たす。
Lobatto IIIA法 はコロケーション法である。
2次の方法は陰的台形公式として知られ[11] 、以下の配列で与えられる。
0 0 0 1 1 / 2 1 / 2 1 / 2 1 / 2 {\displaystyle {\begin{array}{c|cc}0&0&0\\1&1/2&1/2\\\hline &1/2&1/2\\\end{array}}} さらに4次の方法は以下の配列で与えられる[12] 。
0 0 0 0 1 / 2 5 / 24 1 / 3 − 1 / 24 1 1 / 6 2 / 3 1 / 6 1 / 6 2 / 3 1 / 6 {\displaystyle {\begin{array}{c|ccc}0&0&0&0\\1/2&5/24&1/3&-1/24\\1&1/6&2/3&1/6\\\hline &1/6&2/3&1/6\\\end{array}}} これらの方法はどれもA-安定であるが、L-安定やB-安定ではない。
Lobatto IIIB法 はコロケーション法ではないけど、非連続的コロケーション法として見ることができる。
2次の方法は以下の配列で与えられる[13] 。
0 1 / 2 0 1 1 / 2 0 1 / 2 1 / 2 {\displaystyle {\begin{array}{c|cc}0&1/2&0\\1&1/2&0\\\hline &1/2&1/2\\\end{array}}} さらに4次の方法は以下の配列であたえられる[12] 。
0 1 / 6 − 1 / 6 0 1 / 2 1 / 6 1 / 3 0 1 1 / 6 5 / 6 0 1 / 6 2 / 3 1 / 6 {\displaystyle {\begin{array}{c|ccc}0&1/6&-1/6&0\\1/2&1/6&1/3&0\\1&1/6&5/6&0\\\hline &1/6&2/3&1/6\\\end{array}}} これらの方法はどれもA-安定であるが、L-安定やB-安定ではない。
Lobatto IIIC法 も非連続的コロケーション法である。
2次の方法は以下の配列で与えられる[11] 。
0 1 / 2 − 1 / 2 1 1 / 2 1 / 2 1 / 2 1 / 2 1 0 {\displaystyle {\begin{array}{c|cc}0&1/2&-1/2\\1&1/2&1/2\\\hline &1/2&1/2\\&1&0\\\end{array}}} さらに4次の方法は以下の配列で与えられる[12] 。
0 1 / 6 − 1 / 3 1 / 6 1 / 2 1 / 6 5 / 12 − 1 / 12 1 1 / 6 2 / 3 1 / 6 1 / 6 2 / 3 1 / 6 {\displaystyle {\begin{array}{c|ccc}0&1/6&-1/3&1/6\\1/2&1/6&5/12&-1/12\\1&1/6&2/3&1/6\\\hline &1/6&2/3&1/6\\\end{array}}} これらの方法は、すべてL-安定であり、さらに代数的安定(よってB-安定)でもある。そのため、硬い方程式 に対する適切な方法である。
Lobatto IIIC*法[13] は文献によって、Lobatto III法[14] 、ブッチャーのLobatto法やLobatto IIIC法としても知られている。
2次の方法は上述の陽的台形公式にあたり(よって陰的方法ではない)、以下の配列で与えられる[11] 。
0 0 0 1 1 0 1 / 2 1 / 2 {\displaystyle {\begin{array}{c|cc}0&0&0\\1&1&0\\\hline &1/2&1/2\\\end{array}}} さらに4次の方法は以下の配列で与えられる[12] 。
0 0 0 0 1 / 2 1 / 4 1 / 4 0 1 0 1 0 1 / 6 2 / 3 1 / 6 {\displaystyle {\begin{array}{c|ccc}0&0&0&0\\1/2&1/4&1/4&0\\1&0&1&0\\\hline &1/6&2/3&1/6\\\end{array}}} これらの方法は、どれもA-安定でも、L-安定でも、B-安定でもない。
上述のLobatto法の係数はすべて一意に定められるため、方法の線型結合も考えられる[13] 。一般的に、3つの実数パラメータ ( α A , α B , α C ) {\displaystyle (\alpha _{A},\alpha _{B},\alpha _{C})} からなるLobatto係数を以下のようにする。
a i , j ( α A , α B , α C ) = α A a i , j A + α B a i , j B + α C a i , j C + α C ∗ a i , j C ∗ {\displaystyle a_{i,j}(\alpha _{A},\alpha _{B},\alpha _{C})=\alpha _{A}a_{i,j}^{A}+\alpha _{B}a_{i,j}^{B}+\alpha _{C}a_{i,j}^{C}+\alpha _{C*}a_{i,j}^{C*}} 但し、
α C ∗ = 1 − α A − α B − α C {\displaystyle \alpha _{C*}=1-\alpha _{A}-\alpha _{B}-\alpha _{C}} ここで、a A i,j はLobatto IIIA法に対するルンゲ=クッタ行列である。
αA = 2 、αB = 2 、αC = -1 のとき、対応する方法はLobatto IIID法であり、Lobatto IIINW法とも呼ばれる。
2次の方法は以下の配列であたえられる。
0 1 / 2 1 / 2 1 − 1 / 2 1 / 2 1 / 2 1 / 2 {\displaystyle {\begin{array}{c|cc}0&1/2&1/2\\1&-1/2&1/2\\\hline &1/2&1/2\\\end{array}}} さらに4次の方法は以下の配列であたえられる。
0 1 / 6 0 − 1 / 6 1 / 2 1 / 12 5 / 12 0 1 1 / 2 1 / 3 1 / 6 1 / 6 2 / 3 1 / 6 {\displaystyle {\begin{array}{c|ccc}0&1/6&0&-1/6\\1/2&1/12&5/12&0\\1&1/2&1/3&1/6\\\hline &1/6&2/3&1/6\\\end{array}}} これらの方法もすべてL-安定であり、さらに代数的安定(よってB-安定)である。
Radau法[10] は次数 2s -1 を持ち、A-安定である。しかし、方法に対する計算コストが高い以上、方法の次数が落ちるという恐れもある。
Radau IA法の係数 ci は方程式
P s ( 2 x − 1 ) + P s − 1 ( 2 x − 1 ) = 0 {\displaystyle P_{s}(2x-1)+P_{s-1}(2x-1)=0} の解である。ここで、Ps は s 次ルジャンドル多項式 である。
3次の方法は以下の配列で与えられる[15] 。
0 1 / 4 − 1 / 4 2 / 3 1 / 4 5 / 12 1 / 4 3 / 4 {\displaystyle {\begin{array}{c|cc}0&1/4&-1/4\\2/3&1/4&5/12\\\hline &1/4&3/4\\\end{array}}} さらに5次の方法は以下の配列で与えられる[15] 。
0 1 9 − 1 − 6 18 − 1 + 6 18 3 5 − 6 10 1 9 11 45 + 7 6 360 11 45 − 43 6 360 3 5 + 6 10 1 9 11 45 + 43 6 360 11 45 − 7 6 360 1 9 4 9 + 6 36 4 9 − 6 36 {\displaystyle {\begin{array}{c|ccc}0&{\frac {1}{9}}&{\frac {-1-{\sqrt {6}}}{18}}&{\frac {-1+{\sqrt {6}}}{18}}\\{\frac {3}{5}}-{\frac {\sqrt {6}}{10}}&{\frac {1}{9}}&{\frac {11}{45}}+{\frac {7{\sqrt {6}}}{360}}&{\frac {11}{45}}-{\frac {43{\sqrt {6}}}{360}}\\{\frac {3}{5}}+{\frac {\sqrt {6}}{10}}&{\frac {1}{9}}&{\frac {11}{45}}+{\frac {43{\sqrt {6}}}{360}}&{\frac {11}{45}}-{\frac {7{\sqrt {6}}}{360}}\\\hline &{\frac {1}{9}}&{\frac {4}{9}}+{\frac {\sqrt {6}}{36}}&{\frac {4}{9}}-{\frac {\sqrt {6}}{36}}\\\end{array}}} Radau IIA法の係数 ci は方程式
P s ( 2 x − 1 ) − P s − 1 ( 2 x − 1 ) = 0 {\displaystyle P_{s}(2x-1)-P_{s-1}(2x-1)=0} の解である。
3次の方法は以下の配列で与えられる[15] 。
1 / 3 5 / 12 − 1 / 12 1 3 / 4 1 / 4 3 / 4 1 / 4 {\displaystyle {\begin{array}{c|cc}1/3&5/12&-1/12\\1&3/4&1/4\\\hline &3/4&1/4\\\end{array}}} さらに5次の方法は以下の配列で与えられる[11] 。
2 5 − 6 10 11 45 − 7 6 360 37 225 − 169 6 1800 − 2 225 + 6 75 2 5 + 6 10 37 225 + 169 6 1800 11 45 + 7 6 360 − 2 225 − 6 75 1 4 9 − 6 36 4 9 + 6 36 1 9 4 9 − 6 36 4 9 + 6 36 1 9 {\displaystyle {\begin{array}{c|ccc}{\frac {2}{5}}-{\frac {\sqrt {6}}{10}}&{\frac {11}{45}}-{\frac {7{\sqrt {6}}}{360}}&{\frac {37}{225}}-{\frac {169{\sqrt {6}}}{1800}}&-{\frac {2}{225}}+{\frac {\sqrt {6}}{75}}\\{\frac {2}{5}}+{\frac {\sqrt {6}}{10}}&{\frac {37}{225}}+{\frac {169{\sqrt {6}}}{1800}}&{\frac {11}{45}}+{\frac {7{\sqrt {6}}}{360}}&-{\frac {2}{225}}-{\frac {\sqrt {6}}{75}}\\1&{\frac {4}{9}}-{\frac {\sqrt {6}}{36}}&{\frac {4}{9}}+{\frac {\sqrt {6}}{36}}&{\frac {1}{9}}\\\hline &{\frac {4}{9}}-{\frac {\sqrt {6}}{36}}&{\frac {4}{9}}+{\frac {\sqrt {6}}{36}}&{\frac {1}{9}}\\\end{array}}} Bogacki, Przemyslaw; Shampine, Lawrence F. (1989), “A 3(2) pair of Runge–Kutta formulas”, Applied Mathematics Letters 2 (4): 321–325, doi :10.1016/0893-9659(89)90079-7 , ISSN 0893-9659 . Butcher, John C. (2008), Numerical Methods for Ordinary Differential Equations , New York: John Wiley & Sons , ISBN 978-0-470-72335-7 . Cash, J. R.; Karp, Alan H. (1990). “A variable order Runge-Kutta method for initial value problems with rapidly varying right-hand sides”. ACM Transactions on Mathematical Software (ACM) 16 (3): 201-222. doi :10.1145/79505.79507 . . Dormand, J. R.; Prince, P. J. (1980), “A family of embedded Runge-Kutta formulae”, Journal of Computational and Applied Mathematics 6 (1): 19–26, doi :10.1016/0771-050X(80)90013-3 . Fehlberg, Erwin (1970). “Klassische Runge-Kutta-Formeln vierter und niedrigerer Ordnung mit Schrittweiten-Kontrolle und ihre Anwendung auf Wärmeleitungsprobleme”. Computing (Springer-Verlag) 6 (1): 61-67. doi :10.1007/BF02241732 . . Hairer, Ernst; Nørsett, Syvert Paul; Wanner, Gerhard (1993), Solving ordinary differential equations I: Nonstiff problems , Berlin, New York: Springer-Verlag , ISBN 978-3-540-56670-0 . Iserles, Arieh (2008), A First Course in the Numerical Analysis of Differential Equations (Second Edition) , Cambridge University Press, ISBN 978-0-521-73490-5 . Jay, Laurent O.. “Lobatto Methods ”. 2016年12月31日 閲覧。 Kutta, Martin Wilhelm (1901), “Beitrag zur näherungsweisen Integration totaler Differentialgleichungen”, Zeitschrift für Mathematik und Physik 46 : 435–453 . Moler, Cleve (2014年). “Ordinary Differential Equation Solvers ODE23 and ODE45 ”. 2016年12月31日 閲覧。 Ralston, Anthony (1962). “Runge-Kutta Methods with Minimum Error Bounds”. Mathematics of Computation 16 (80): 431-437. doi :10.2307/2003133 . .