Family of implicit and explicit iterative methods
微分方程式のルンゲ・クッタ法の比較 （赤は正確な解） y ′ = sin ( t ) 2 ⋅ y {\displaystyle y'=\sin(t)^{2}\cdot y} 数値解析 では 、 ルンゲ・クッタ法 ( RUUNG -ə- KUUT -tah [1] 暗黙的および明示的な 反復法 のファミリーであり 同時非線形方程式 の近似解の 時間的離散化 に使用される オイラー 法 が含まれます 。 [2] カール・ルンゲ と ヴィルヘルム・クッタ によって開発されました 。
ルンゲ・クッタ法 古典的なルンゲ・クッタ法で使用される傾き ルンゲ・クッタ ファミリーの最も広く知られているメンバーは、一般に「RK4」、「古典的なルンゲ・クッタ法」、または単に「ルンゲ・クッタ法」と呼ばれています。
初期値問題を 次のように指定します 。
d y d t = f ( t , y ) , y ( t 0 ) = y 0 . {\displaystyle {\frac {dy}{dt}}=f(t,y),\quad y(t_{0})=y_{0}.} ここに 、時間 の未知の関数（スカラーまたはベクトル）があり 、これを近似したいと考えています。変化 率 はの関数であり、 自身 も の 関数であるとされています 。初期時刻における 対応する 値は です 。関数 と 初期条件 、 が与えられています。 y {\displaystyle y} t {\displaystyle t} d y d t {\displaystyle {\frac {dy}{dt}}} y {\displaystyle y} t {\displaystyle t} y {\displaystyle y} t 0 {\displaystyle t_{0}} y {\displaystyle y} y 0 {\displaystyle y_{0}} f {\displaystyle f} t 0 {\displaystyle t_{0}} y 0 {\displaystyle y_{0}}
ここで、ステップサイズ h > 0 を選択し、以下を定義します。
y n + 1 = y n + h 6 ( k 1 + 2 k 2 + 2 k 3 + k 4 ) , t n + 1 = t n + h {\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] を用いて
k 1 = f ( t n , y n ) , k 2 = f ( t n + h 2 , y n + h k 1 2 ) , k 3 = f ( t n + h 2 , y n + h k 2 2 ) , k 4 = f ( t n + h , y n + h k 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 によって指定された推定傾きの積です。 y n + 1 {\displaystyle y_{n+1}} y ( t n + 1 ) {\displaystyle y(t_{n+1})} y n + 1 {\displaystyle y_{n+1}} y n {\displaystyle y_{n}}
k 1 {\displaystyle k_{1}} は、 （ オイラー法 ） を使用した、区間の開始時の傾きです。 y {\displaystyle y} k 2 {\displaystyle k_{2}} は、および を使用した、区間の中点における傾きです 。 y {\displaystyle y} k 1 {\displaystyle k_{1}} k 3 {\displaystyle k_{3}} は再び中点の傾きですが、今度は とを使用します 。 y {\displaystyle y} k 2 {\displaystyle k_{2}} k 4 {\displaystyle k_{4}} は、および を使用した、区間の終わりの傾きです 。 y {\displaystyle y} k 3 {\displaystyle k_{3}} 4つの傾きを平均化する際に、中点の傾きに重みが置かれる。 が に依存しない場合 、つまり微分方程式が単積分と等価である場合、RK4は シンプソンの定理 である。 [5] f {\displaystyle f} y {\displaystyle y}
RK4 法は 4 次法であり、 局所的な切り捨て誤差は のオーダーである のに対し、 総累積誤差 は のオーダーで あることを意味します 。 O ( h 5 ) {\displaystyle O(h^{5})} O ( h 4 ) {\displaystyle O(h^{4})}
多くの実際のアプリケーションでは、関数 は から独立しており (いわゆる 自律システム 、または特に物理学では時間不変システム)、その増分はまったく計算されず、関数 に渡されず 、 の最終的な式のみが 使用されます。 f {\displaystyle f} t {\displaystyle t} f {\displaystyle f} t n + 1 {\displaystyle t_{n+1}}
陽的ルンゲ・クッタ法 陽的 ルンゲ・クッタ法は 、前述のRK4法の一般化です。次のように与えられます
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},} ここで [6]
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 ) , ⋮ k s = f ( t n + c s h , y n + ( a s 1 k 1 + a s 2 k 2 + ⋯ + a s , s − 1 k s − 1 ) h ) . {\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 ）、 bi （ i = 1, 2, ..., s ）、および c i （ i = 2, 3, ..., s ）を与える必要 が ある。行列[ a ij ]はルンゲ・クッタ行列 と呼ばれ 、 bi と c i は 重み 、 ノード と呼ばれる 。 [7]これらのデータは通常、 ブッチャー・タブロー （ ジョン・C・ブッチャー にちなんで） と呼ばれる記憶法で整理される。
0 {\displaystyle 0} c 2 {\displaystyle c_{2}} a 21 {\displaystyle a_{21}} c 3 {\displaystyle c_{3}} a 31 {\displaystyle a_{31}} a 32 {\displaystyle a_{32}} ⋮ {\displaystyle \vdots } ⋮ {\displaystyle \vdots } ⋱ {\displaystyle \ddots } c s {\displaystyle c_{s}} a s 1 {\displaystyle a_{s1}} a s 2 {\displaystyle a_{s2}} ⋯ {\displaystyle \cdots } a s , s − 1 {\displaystyle a_{s,s-1}} b 1 {\displaystyle b_{1}} b 2 {\displaystyle b_{2}} ⋯ {\displaystyle \cdots } b s − 1 {\displaystyle b_{s-1}} b s {\displaystyle b_{s}}
テイラー 級数 展開は、ルンゲ・クッタ法が矛盾しない条件が、
∑ i = 1 s b i = 1. {\displaystyle \sum _{i=1}^{s}b_{i}=1.} また、この手法に特定の位数 p を要求する場合、局所的な打ち切り誤差が O( h p +1 ) となるという付随的な要件もあります。これらは打ち切り誤差の定義自体から導き出せます。例えば、2段階法では、 b 1 + b 2 = 1、 b 2 c 2 = 1/2、 b 2 a 21 = 1/2 のとき、位数は 2 となります。 [8] なお、係数を決定するための一般的な条件は [8]です。
∑ j = 1 i − 1 a i j = c i for i = 2 , … , s . {\displaystyle \sum _{j=1}^{i-1}a_{ij}=c_{i}{\text{ for }}i=2,\ldots ,s.} しかし、この条件だけでは一貫性を保つには十分でも必要でもありません。 [9]
一般に、明示的 - 段ルンゲ・クッタ法が 次数 を持つ場合 、段数は を満たさなければならないことが証明でき 、 の場合 、 となる 。 [10]しかし、これらの境界が常に 厳密 であるかどうかはわかっていない 。場合によっては、境界を達成できないことが証明されている。例えば、Butcher は、 の場合、 段 を持つ明示的ルンゲ・クッタ法は存在しないことを証明した。 [11] Butcher はまた、 の場合 、段を持つ明示的ルンゲ・クッタ法は存在しないことを証明した 。 [12] しかし、一般に、明示的ルンゲ・クッタ法が 次数 を持つための 正確な最小段数が何であるかは未解決の問題である 。いくつかの既知の値は以下のとおりである。 [13] s {\displaystyle s} p {\displaystyle p} s ≥ p {\displaystyle s\geq p} p ≥ 5 {\displaystyle p\geq 5} s ≥ p + 1 {\displaystyle s\geq p+1} p > 6 {\displaystyle p>6} s = p + 1 {\displaystyle s=p+1} p > 7 {\displaystyle p>7} p + 2 {\displaystyle p+2} s {\displaystyle s} p {\displaystyle p}
p 1 2 3 4 5 6 7 8 min s 1 2 3 4 6 7 9 11 {\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]を参照のこと。 p = 1 , 2 , … , 6 {\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] この方法の主な利点は、ほとんどすべての誤差係数が一般的な方法よりも小さいことですが、時間ステップごとにわずかに多くのFLOP（浮動小数点演算）が必要です。そのブッチャー・タブローは
0 1/3 1/3 2/3 −1/3 1 1 1 −1 1 1/8 3/8 3/8 1/8
しかし、最も単純なルンゲ・クッタ法は（順方向） オイラー法 であり、式 で与えられます 。これは、1段階の唯一の一貫した陽的ルンゲ・クッタ法です。対応する表は y n + 1 = y n + h f ( t n , y n ) {\displaystyle y_{n+1}=y_{n}+hf(t_{n},y_{n})}
2段階の2次法 2段階の2次法の例としては、明示的中 点法 が挙げられます。
y n + 1 = y n + h f ( t n + 1 2 h , y n + 1 2 h f ( t n , y n ) ) . {\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] で与えられるそのような方法のファミリーが存在する。
y n + 1 = y n + h ( ( 1 − 1 2 α ) f ( t n , y n ) + 1 2 α f ( t n + α h , y n + α h f ( t n , y n ) ) ) . {\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 } ( 1 − 1 2 α ) {\displaystyle (1-{\tfrac {1}{2\alpha }})} 1 2 α {\displaystyle {\tfrac {1}{2\alpha }}}
この族では、 中点法 を与え 、 は Heun法 、 [5] を与え、 はRalston法を与える。 α = 1 2 {\displaystyle \alpha ={\tfrac {1}{2}}} α = 1 {\displaystyle \alpha =1} α = 2 3 {\displaystyle \alpha ={\tfrac {2}{3}}}
使用 例として、 ラルストン法 としても知られる、α = 2/3の2段階2次ルンゲ・クッタ法を考えてみましょう。これは次の表で与えられます
対応する方程式で
k 1 = f ( t n , y n ) , k 2 = f ( t n + 2 3 h , y n + 2 3 h k 1 ) , y n + 1 = y n + h ( 1 4 k 1 + 3 4 k 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}}} この方法は初期値問題を解くために使用される
d y d t = tan ( y ) + 1 , y 0 = 1 , t ∈ [ 1 , 1.1 ] {\displaystyle {\frac {dy}{dt}}=\tan(y)+1,\quad y_{0}=1,\ t\in [1,1.1]} ステップ サイズ h = 0.025 なので、この方法では 4 つのステップを実行する必要があります。
この方法は次のように進行します。
t 0 = 1 : {\displaystyle t_{0}=1\colon } y 0 = 1 {\displaystyle y_{0}=1} t 1 = 1.025 : {\displaystyle t_{1}=1.025\colon } y 0 = 1 {\displaystyle y_{0}=1} k 1 = 2.557407725 {\displaystyle k_{1}=2.557407725} k 2 = f ( t 0 + 2 3 h , y 0 + 2 3 h k 1 ) = 2.7138981400 {\displaystyle k_{2}=f(t_{0}+{\tfrac {2}{3}}h,\ y_{0}+{\tfrac {2}{3}}hk_{1})=2.7138981400} y 1 = y 0 + h ( 1 4 k 1 + 3 4 k 2 ) = 1.066869388 _ {\displaystyle y_{1}=y_{0}+h({\tfrac {1}{4}}k_{1}+{\tfrac {3}{4}}k_{2})={\underline {1.066869388}}} t 2 = 1.05 : {\displaystyle t_{2}=1.05\colon } y 1 = 1.066869388 {\displaystyle y_{1}=1.066869388} k 1 = 2.813524695 {\displaystyle k_{1}=2.813524695} k 2 = f ( t 1 + 2 3 h , y 1 + 2 3 h k 1 ) {\displaystyle k_{2}=f(t_{1}+{\tfrac {2}{3}}h,\ y_{1}+{\tfrac {2}{3}}hk_{1})} y 2 = y 1 + h ( 1 4 k 1 + 3 4 k 2 ) = 1.141332181 _ {\displaystyle y_{2}=y_{1}+h({\tfrac {1}{4}}k_{1}+{\tfrac {3}{4}}k_{2})={\underline {1.141332181}}} t 3 = 1.075 : {\displaystyle t_{3}=1.075\colon } y 2 = 1.141332181 {\displaystyle y_{2}=1.141332181} k 1 = 3.183536647 {\displaystyle k_{1}=3.183536647} k 2 = f ( t 2 + 2 3 h , y 2 + 2 3 h k 1 ) {\displaystyle k_{2}=f(t_{2}+{\tfrac {2}{3}}h,\ y_{2}+{\tfrac {2}{3}}hk_{1})} y 3 = y 2 + h ( 1 4 k 1 + 3 4 k 2 ) = 1.227417567 _ {\displaystyle y_{3}=y_{2}+h({\tfrac {1}{4}}k_{1}+{\tfrac {3}{4}}k_{2})={\underline {1.227417567}}} t 4 = 1.1 : {\displaystyle t_{4}=1.1\colon } y 3 = 1.227417567 {\displaystyle y_{3}=1.227417567} k 1 = 3.796866512 {\displaystyle k_{1}=3.796866512} k 2 = f ( t 3 + 2 3 h , y 3 + 2 3 h k 1 ) {\displaystyle k_{2}=f(t_{3}+{\tfrac {2}{3}}h,\ y_{3}+{\tfrac {2}{3}}hk_{1})} y 4 = y 3 + h ( 1 4 k 1 + 3 4 k 2 ) = 1.335079087 _ . {\displaystyle y_{4}=y_{3}+h({\tfrac {1}{4}}k_{1}+{\tfrac {3}{4}}k_{2})={\underline {1.335079087}}.}
数値解は下線部の値に対応します。
適応型ルンゲ・クッタ法 適応型ルンゲ・クッタ法は、単一のルンゲ・クッタステップの局所的な打ち切り誤差の推定値を生成するように設計されています。これは、次数 と次数 の2つの方法を用いることで行われます 。これらの方法は相互に絡み合っており、つまり共通の中間ステップを持っています。これにより、高次法を用いたステップと比較して、誤差の推定にかかる計算コストはほとんど、あるいは無視できるほどになります p {\displaystyle p} p − 1 {\displaystyle p-1}
積分中、ステップサイズは推定誤差がユーザー定義の閾値を下回るように調整されます。誤差が大きすぎる場合は、ステップサイズを小さくしてステップを繰り返します。誤差がはるかに小さい場合は、時間を節約するためにステップサイズを大きくします。これにより、（ほぼ）最適なステップサイズが得られ、計算時間を節約できます。さらに、ユーザーは適切なステップサイズを見つけるために時間を費やす必要がありません。
低次のステップは次のように与えられる。
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 {\displaystyle k_{i}}
e n + 1 = y n + 1 − y n + 1 ∗ = h ∑ i = 1 s ( b i − b i ∗ ) k i , {\displaystyle e_{n+1}=y_{n+1}-y_{n+1}^{*}=h\sum _{i=1}^{s}(b_{i}-b_{i}^{*})k_{i},} これは である 。この種の方法のブッチャー表は の値を与えるように拡張される ： O ( h p ) {\displaystyle O(h^{p})} b i ∗ {\displaystyle b_{i}^{*}}
0 c 2 {\displaystyle c_{2}} a 21 {\displaystyle a_{21}} c 3 {\displaystyle c_{3}} a 31 {\displaystyle a_{31}} a 32 {\displaystyle a_{32}} ⋮ {\displaystyle \vdots } ⋮ {\displaystyle \vdots } ⋱ {\displaystyle \ddots } c s {\displaystyle c_{s}} a s 1 {\displaystyle a_{s1}} a s 2 {\displaystyle a_{s2}} ⋯ {\displaystyle \cdots } a s , s − 1 {\displaystyle a_{s,s-1}} b 1 {\displaystyle b_{1}} b 2 {\displaystyle b_{2}} ⋯ {\displaystyle \cdots } b s − 1 {\displaystyle b_{s-1}} b s {\displaystyle b_{s}} b 1 ∗ {\displaystyle b_{1}^{*}} b 2 ∗ {\displaystyle b_{2}^{*}} ⋯ {\displaystyle \cdots } b s − 1 ∗ {\displaystyle b_{s-1}^{*}} b s ∗ {\displaystyle b_{s}^{*}}
ルンゲ ・クッタ・フェールベルグ法には、 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次の オイラー法 を組み合わせたものです。その拡張ブッチャー・タブローは次のとおりです
他の適応型ルンゲ・クッタ法には、 ボガッキ・シャンパイン法 （3次と2次）、 キャッシュ・カープ法 、 ドルマンド・プリンス法 （どちらも5次と4次）
があります
非合流型ルンゲ・クッタ法 ルンゲ・クッタ法は、 すべての が異なる場合、 非合流型 [21]であると言われます c i , i = 1 , 2 , … , s {\displaystyle c_{i},\,i=1,2,\ldots ,s}
ルンゲ・クッタ・ニストローム法 ルンゲ・クッタ・ニストロム法は、2階微分方程式に最適化されている特殊なルンゲ・クッタ法である。 [22] [23] 2階常微分方程式系に対する一般的な（つまり、明示的バージョンと暗黙的バージョンの両方を含む）ルンゲ・クッタ・ニストロム法
y ¨ i = f i ( y 1 , y 2 , … , y n ) {\displaystyle {\ddot {y}}_{i}=f_{i}(y_{1},y_{2},\ldots ,y_{n})}
順序 は 「ステージ」の数であり 、次の式で表されます。 s + 1 {\displaystyle s+1} s {\displaystyle s}
{ g i = y m + c i h y ˙ m + h 2 ∑ j = 1 s a i j f ( g j ) , i = 1 , 2 , … , s y m + 1 = y m + h y ˙ m + h 2 ∑ j = 1 s b ¯ j f ( g j ) y ˙ m + 1 = y ˙ m + h ∑ j = 1 s b j f ( g j ) {\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}}}
この方法の明示的なバージョンでは、の式における 和の 上限 を に置き換えることができる 。 [24] すなわち、和を に置き換えることができる 。これは、次の形式のブッチャー表を形成する。 s {\displaystyle s} ∑ j = 1 s a i j f ( g j ) {\displaystyle \sum _{j=1}^{s}a_{ij}f(g_{j})} g i {\displaystyle g_{i}} i − 1 {\displaystyle i-1} ∑ j = 1 i − 1 {\displaystyle \sum _{j=1}^{i-1}}
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 = c A b ¯ ⊤ b ⊤ {\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}}}
次のブッチャー テーブルには、2 つの 4 次明示的 RKN 法が示されています。
c i a i j 3 + 3 6 0 0 0 3 − 3 6 2 − 3 12 0 0 3 + 3 6 0 3 6 0 b i ¯ 5 − 3 3 24 3 + 3 12 1 + 3 24 b i 3 − 2 3 12 1 2 3 + 2 3 12 {\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}}} c i a i j 3 − 3 6 0 0 0 3 + 3 6 2 + 3 12 0 0 3 − 3 6 0 − 3 6 0 b i ¯ 5 + 3 3 24 3 − 3 12 1 − 3 24 b i 3 + 2 3 12 1 2 3 − 2 3 12 {\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つのスキームは、元の方程式が保存的な古典力学システムから導かれる場合、すなわち
f i ( x 1 , … , x n ) = ∂ V ∂ x i ( x 1 , … , x n ) {\displaystyle f_{i}(x_{1},\ldots ,x_{n})={\frac {\partial V}{\partial x_{i}}}(x_{1},\ldots ,x_{n})}
あるスカラー関数に対して 。 [25] V {\displaystyle V}
陰的ルンゲ・クッタ法 これまでに述べたルンゲ・クッタ法はすべて 陽的手法 です。陽的ルンゲ・クッタ法は、絶対安定領域が狭く、特に有界であるため、一般的に 硬い方程式 の解法には適していません 。 [26]この問題は、 偏微分方程式 の解法において特に重要です
陽的ルンゲ・クッタ法の不安定性は、陰的ルンゲ・クッタ法の開発の動機となった。陰的ルンゲ・クッタ法は以下の形をとる。
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 s a i j k j ) , i = 1 , … , s . {\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.} [27] 陽解法との違いは、陽解法では jの和が i -1 までしか上がらないことです。 [28] これはブッチャー・テーブルにも現れています。陽解法の 係数行列は下三角行列です。陰解法では、 jの和は s まで上がり 、係数行列は三角行列ではないため、ブッチャー・テーブルは次のようになります。 [18] a i j {\displaystyle a_{ij}}
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 ∗ = c A b T {\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}}} 行の説明については、上記の適応型ルンゲ・クッタ法を参照してください 。 b ∗ {\displaystyle b^{*}}
この違いの結果として、各ステップで代数方程式系を解かなければなりません。これにより計算コストが大幅に増加します。s段階の手法を用いて m個 の 要素を持つ微分方程式を解くと 、代数方程式系は ms個の要素を持つことになります。これは、暗黙的 線形多段階法 （常微分方程式のためのもう一つの大きな手法群） とは対照的です。暗黙的 s段階線形多段階法は、 m 個の要素のみを持つ代数方程式系を解く必要がある ため、ステップ数が増加しても方程式系のサイズは増加しません。 [29]
例 暗黙的ルンゲ・クッタ法の最も単純な例は 、 後退オイラー法です
y n + 1 = y n + h f ( t n + h , y n + 1 ) . {\displaystyle y_{n+1}=y_{n}+hf(t_{n}+h,\ y_{n+1}).\,} これに対するブッチャーのタブローは単純です：
1 1 1 {\displaystyle {\begin{array}{c|c}1&1\\\hline &1\\\end{array}}} このブッチャーのタブローは、次の式に対応しています。
k 1 = f ( t n + h , y n + h k 1 ) and y n + 1 = y n + h k 1 , {\displaystyle k_{1}=f(t_{n}+h,\ y_{n}+hk_{1})\quad {\text{and}}\quad y_{n+1}=y_{n}+hk_{1},} これを変形すると、上記の後退オイラー法の式が得られます。
暗黙的ルンゲ・クッタ法のもう一つの例は 台形則 です。そのブッチャー表は次のようになります。
0 0 0 1 1 2 1 2 1 2 1 2 1 0 {\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}}} 台形則は（その記事で議論されているように） 選点法の 一種である。すべての選点法は暗黙的ルンゲ・クッタ法であるが、すべての暗黙的ルンゲ・クッタ法が選点法であるわけではない。 [30]
ガウス ・ルジャンドル法は、 ガウス積分法 に基づく選点法の一種である。s 段階 のガウス・ルジャンドル法は、次数が 2 である （したがって、任意の高次数を持つ方法を構築できる）。 [31] 2段階（したがって次数が4）の方法は、ブッチャー・タブローを持つ。
1 2 − 1 6 3 1 4 1 4 − 1 6 3 1 2 + 1 6 3 1 4 + 1 6 3 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 {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}}} [29]
安定性 陰的ルンゲ・クッタ法が陽的ルンゲ・クッタ法よりも優れている点は、特に 硬い方程式 に適用した場合の安定性が高いことです。線形テスト方程式を考えてみましょう 。この方程式にルンゲ・クッタ法を適用すると、反復回数は となり 、 r は次のように与えられます y ′ = λ y {\displaystyle y'=\lambda y} y n + 1 = r ( h λ ) y n {\displaystyle y_{n+1}=r(h\lambda )\,y_{n}}
r ( z ) = 1 + z b T ( I − z A ) − 1 e = det ( I − z A + z e b T ) det ( I − z A ) , {\displaystyle r(z)=1+zb^{T}(I-zA)^{-1}e={\frac {\det(I-zA+zeb^{T})}{\det(I-zA)}},} [32] ここで eは 1のベクトルを表します。関数 rは 安定性関数 と呼ばれます 。 [33] 式から、 rは2つの s 次多項式の商であり 、その方法は s 段階です。陽解法は厳密に下三角行列 A を持ち、det( I − zA ) = 1であり、安定性関数は多項式であることを意味します。 [34]
線形検定方程式の数値解は、 z = h λ において | r ( z ) | < 1 のときゼロに減衰する。このような z の集合は 絶対安定領域 と呼ばれる 。特に、 Re( z ) < 0 となるすべての z が 絶対安定領域内にあるとき、この方法は 絶対安定 であると言われる。陽的ルンゲ・クッタ法の安定関数は多項式であるため、陽的ルンゲ・クッタ法はA安定にはなり得ない。 [34]
この方法が p 次数を持つ場合、安定関数は を満たす 。したがって、指数関数を最もよく近似する、与えられた次数の多項式の商を調べることは興味深い。これらは パデ近似として知られている。分子が m 次、分母が n次 のパデ近似が A安定であるための必要十分条件は、 m ≤ n ≤ m + 2である 。[35] r ( z ) = e z + O ( z p + 1 ) {\displaystyle r(z)={\textrm {e}}^{z}+O(z^{p+1})} z → 0 {\displaystyle z\to 0}
s 段階のガウス・ルジャンドル法は 2s の位数を持つため、その安定性関数は m = n = s のパデ近似となる 。したがって、この方法はA安定である。 [36]これは、A安定ルンゲ・クッタ法が任意の高位数を持つことができることを示している。対照的に、A安定 線形多段階法 の位数は 2を超えることはできない。 [37]
B安定性 微分方程式の解における A安定性の 概念は、線形自律方程式に関連しています 。Dahlquist (1963) は、単調性条件を満たす非線形システムに適用された場合の数値スキームの安定性の調査を提案しました。対応する概念は、 マルチステップ法（および関連するワンレッグ法）の場合は G安定性、ルンゲ・クッタ法の場合は B安定性 （Butcher, 1975）として定義されました。非線形システムに適用され 、が成り立つことが証明されるルンゲ・クッタ法は、この条件が 2つの数値解に対して 成り立つ場合、 B安定 と呼ばれます y ′ = λ y {\displaystyle y'=\lambda y} y ′ = f ( y ) {\displaystyle y'=f(y)} ⟨ f ( y ) − f ( z ) , y − z ⟩ ≤ 0 {\displaystyle \langle f(y)-f(z),\ y-z\rangle \leq 0} ‖ y n + 1 − z n + 1 ‖ ≤ ‖ y n − z n ‖ {\displaystyle \|y_{n+1}-z_{n+1}\|\leq \|y_{n}-z_{n}\|}
、、 をそれぞれ定義される 3つの 行列
とします 。ルンゲ・クッタ法は、行列 と がともに非負定値であるとき、代数的に安定であるといわれます [ 38 ] 。B 安定 性 の十分条件 [39] は、 と がともに 非負定値であることです。 B {\displaystyle B} M {\displaystyle M} Q {\displaystyle Q} s × s {\displaystyle s\times s} B = diag ( b 1 , b 2 , … , b s ) , M = B A + A T B − b b T , Q = B A − 1 + A − T B − A − T b b T A − 1 . {\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}}} B {\displaystyle B} M {\displaystyle M} B {\displaystyle B} Q {\displaystyle Q}
ルンゲ・クッタ4次法の導出 一般に、ルンゲ・クッタ法は 次のように記述できます。 s {\displaystyle s}
y t + h = y t + h ⋅ ∑ i = 1 s a i k i + O ( h s + 1 ) , {\displaystyle y_{t+h}=y_{t}+h\cdot \sum _{i=1}^{s}a_{i}k_{i}+{\mathcal {O}}(h^{s+1}),} ここで：
k i = ∑ j = 1 s β i j f ( k j , t n + α i h ) {\displaystyle k_{i}=\sum _{j=1}^{s}\beta _{ij}f(k_{j},\ t_{n}+\alpha _{i}h)} は、の階 微分 を評価することで得られる増分です y t {\displaystyle y_{t}} i {\displaystyle i}
我々は、上で説明したように、任意の区間の開始点、中点、終了点 で評価された一般的な公式を使用して、ルンゲ・クッタ4次法の 導出 [40] を展開します。したがって、次のように選択します。 s = 4 {\displaystyle s=4} ( t , t + h ) {\displaystyle (t,\ t+h)}
α i β i j α 1 = 0 β 21 = 1 2 α 2 = 1 2 β 32 = 1 2 α 3 = 1 2 β 43 = 1 α 4 = 1 {\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}}} そうで なければ、まず以下の量を定義します。 β i j = 0 {\displaystyle \beta _{ij}=0}
y t + h 1 = y t + h f ( y t , t ) y t + h 2 = y t + h f ( y t + h / 2 1 , t + h 2 ) y t + h 3 = y t + h f ( y t + h / 2 2 , t + h 2 ) {\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}}} ここで 、 次のように定義します。 y t + h / 2 1 = y t + y t + h 1 2 {\displaystyle y_{t+h/2}^{1}={\dfrac {y_{t}+y_{t+h}^{1}}{2}}} y t + h / 2 2 = y t + y t + h 2 2 . {\displaystyle y_{t+h/2}^{2}={\dfrac {y_{t}+y_{t+h}^{2}}{2}}.}
k 1 = f ( y t , t ) k 2 = f ( y t + h / 2 1 , t + h 2 ) = f ( y t + h 2 k 1 , t + h 2 ) k 3 = f ( y t + h / 2 2 , t + h 2 ) = f ( y t + h 2 k 2 , t + h 2 ) k 4 = f ( y t + h 3 , t + h ) = f ( y t + h k 3 , t + h ) {\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}}} また、前の関係については、 まで次の等式が成り立つことが示せます 。 ここで、 は の時間に関する 全微分です。 O ( h 2 ) {\displaystyle {\mathcal {O}}(h^{2})} k 2 = f ( y t + h / 2 1 , t + h 2 ) = f ( y t + h 2 k 1 , t + h 2 ) = f ( y t , t ) + h 2 d d t f ( y t , t ) k 3 = f ( y t + h / 2 2 , t + h 2 ) = f ( y t + h 2 f ( y t + h 2 k 1 , t + h 2 ) , t + h 2 ) = f ( y t , t ) + h 2 d d t [ f ( y t , t ) + h 2 d d t f ( y t , t ) ] k 4 = f ( y t + h 3 , t + h ) = f ( y t + h f ( y t + h 2 k 2 , t + h 2 ) , t + h ) = f ( y t + h f ( y t + h 2 f ( y t + h 2 f ( y t , t ) , t + h 2 ) , t + h 2 ) , t + h ) = f ( y t , t ) + h d d t [ f ( y t , t ) + h 2 d d t [ f ( y t , t ) + h 2 d d t f ( y t , t ) ] ] {\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}}} d d t f ( y t , t ) = ∂ ∂ y f ( y t , t ) y ˙ t + ∂ ∂ t f ( y t , t ) = f y ( y t , t ) y ˙ t + f t ( y t , t ) := y ¨ t {\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}} f {\displaystyle f}
今導出した式を使って一般的な式を表現すると、次のようになります。 y t + h = y t + h { a ⋅ f ( y t , t ) + b ⋅ [ f ( y t , t ) + h 2 d d t f ( y t , t ) ] + + c ⋅ [ f ( y t , t ) + h 2 d d t [ f ( y t , t ) + h 2 d d t f ( y t , t ) ] ] + + d ⋅ [ f ( y t , t ) + h d d t [ f ( y t , t ) + h 2 d d t [ f ( y t , t ) + h 2 d d t f ( y t , t ) ] ] ] } + O ( h 5 ) = y t + a ⋅ h f t + b ⋅ h f t + b ⋅ h 2 2 d f t d t + c ⋅ h f t + c ⋅ h 2 2 d f t d t + + c ⋅ h 3 4 d 2 f t d t 2 + d ⋅ h f t + d ⋅ h 2 d f t d t + d ⋅ h 3 2 d 2 f t d t 2 + d ⋅ h 4 4 d 3 f t d t 3 + O ( h 5 ) {\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}}}
これを の周り の テイラー級数 と比較すると次のようになります 。 y t + h {\displaystyle y_{t+h}} t {\displaystyle t} y t + h = y t + h y ˙ t + h 2 2 y ¨ t + h 3 6 y t ( 3 ) + h 4 24 y t ( 4 ) + O ( h 5 ) = = y t + h f ( y t , t ) + h 2 2 d d t f ( y t , t ) + h 3 6 d 2 d t 2 f ( y t , t ) + h 4 24 d 3 d t 3 f ( y t , 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}}}
係数に関する制約システムが得られます。
{ a + b + c + d = 1 1 2 b + 1 2 c + d = 1 2 1 4 c + 1 2 d = 1 6 1 4 d = 1 24 {\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}}} これを解くと 上記のようになります。 a = 1 6 , b = 1 3 , c = 1 3 , d = 1 6 {\displaystyle a={\frac {1}{6}},b={\frac {1}{3}},c={\frac {1}{3}},d={\frac {1}{6}}}
参照
注釈 ^ 「ルンゲ・クッタ法」 Dictionary.com 。 2021年 4月4日 閲覧 ^ DEVRIES, Paul L.; HASBUN, Javier E. 『計算物理学入門』第2版. Jones and Bartlett Publishers, 2011. p. 215. ^ プレス他。 2007、p. 908; Süli & Mayers 2003、p. 328 ^ ab 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 値を段階として用いている。 ^ ab Süli & Mayers 2003、p. 328 ^ プレス他 2007年、907ページ ^ イゼルレス 1996、38ページ ^ イゼルレス 1996、39ページ ^ 反例として、 と をランダムに選択した任意の明示的2段階ルンゲ・クッタ法を考えてみましょう 。 この 方法は矛盾がなく、（一般に）1次収束します。一方、 を とする1段階法は 矛盾しており、 が自明に成立しているにもかかわらず収束しません 。 b 1 = b 2 = 1 / 2 {\displaystyle b_{1}=b_{2}=1/2} c 1 {\displaystyle c_{1}} a 21 {\displaystyle a_{21}} b 1 = 1 / 2 {\displaystyle b_{1}=1/2} ∑ j = 1 i − 1 a i j = c i for i = 2 , … , s . {\displaystyle \sum _{j=1}^{i-1}a_{ij}=c_{i}{\text{ for }}i=2,\ldots ,s.} ^ ブッチャー 2008、187ページ ^ abc ブッチャー 1965年、408ページ ^ ブッチャー 1985 ^ ブッチャー 2008、187–196ページ ^ ブッチャー 1964 ^ カーティス 1970、268ページ ^ ハイラー、ノーセット、ワナー、1993、p. 179 ^ ブッチャー 1996、247ページ ^ ab Süli & Mayers 2003、p. 352 ^ Hairer、Nørsett & Wanner (1993、p. 138) は Kutta (1901) を参照。 ^ スーリ＆メイヤーズ 2003、327ページ ^ ランバート 1991, 278ページ ^ Dormand, JR; Prince, PJ (1978年10月). 「力学天文学における数値シミュレーションのための新しいルンゲ・クッタアルゴリズム」. 天体力学 . 18 (3): 223– 232. Bibcode :1978CeMec..18..223D. doi :10.1007/BF01230162. S2CID 120974351. ^ Fehlberg, E. (1974年10月). 一般2階微分方程式に対するステップサイズ制御による古典的7次、6次、5次ルンゲ・クッタ・ニストローム公式（報告書）（NASA TR R-432版）. マーシャル宇宙飛行センター、アラバマ州：アメリカ航空宇宙局. ^ ブッチャー 2008、94ページ ^ 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. ^ Süli & Mayers 2003、pp. 349–351 ^ イゼルレス、1996、p. 41; Süli & Mayers 2003、351–352 ページ ^ ブッチャー 2008、94ページ ^ ab Süli & Mayers 2003、p. 353 ^ イゼルレス 1996、43～44ページ ^ イゼルレス 1996、47ページ ^ ヘアラー&ワナー 1996年、40–41ページ ^ ヘアラー＆ワナー 1996年、40ページ ^ イゼルレス 1996、60ページ ^ イゼルレス 1996、62～63ページ ^ イゼルレス 1996、63ページ ^ この結果はDahlquist (1963)によるものです。 ^ ランバート 1991, 275ページ ^ ランバート 1991, 274ページ ^ 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 。 Cellier, 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年）、 微分方程式の数値解析入門 、 ケンブリッジ大学出版局 、 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版)、 ケンブリッジ大学出版局 、 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) , Journal of Mathematical Science & Mathematics Education , 7 (2): 1– 10 。 離散数学上級 ignou 参考書（コード：mcs033） ジョン・C・ブッチャー：「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 。 Curtis, AR (1970)、 「ステップごとに11回の関数評価を行う8次ルンゲ・クッタ過程」 、 Numerische Mathematik 、 16 (3): 268– 277、 doi :10.1007/BF02219778 。 クーパー, GJ; ヴァーナー, 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。
.svg/1280px-Comparison_of_the_Runge-Kutta_methods_for_the_differential_equation_(red_is_the_exact_solution).svg.png)
![{\displaystyle {\frac {dy}{dt}}=\tan(y)+1,\quad y_{0}=1,\ t\in [1,1.1]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/4a9f561fa8e63f2b86d6226df93dd25d29f76fc0)
![{\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}BA^{-T}bb^{T}A^{-1}.\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/fabfb494d0d855eb81e1ebf6bc78ebc21f82575d)
![{\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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/d7fab64b114fcd6b0ff7ccc1bb7d9754e3eace2d)
![{\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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/c5d443e3cc50f9e0636ccb9094ff5c4c304b455d)
![{\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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/d5fde6d4547411755a49dbe4fb9078743a9f1a08)
.svg){kind=link}
