ルンゲ＝クッタ法のソースを表示

[[数値解析]]において'''ルンゲ＝クッタ法'''（{{lang-en-short|Runge&ndash;Kutta method}}）とは、[[初期値問題]]に対して近似解を与える[[常微分方程式の数値解法]]に対する総称である。この技法は1900年頃に数学者[[カール・ルンゲ]]と[[マルティン・クッタ]]によって発展を見た。

== 古典的ルンゲ＝クッタ法 ==
一連のルンゲ＝クッタ公式の中で最も広く知られているのが、'''古典的ルンゲ=クッタ法''' ('''RK4'''、もしくは単に狭義の '''ルンゲ=クッタ法'''、{{lang-en-short|the (classical) Runge&ndash;Kutta method}}) などと呼ばれる4次の公式である。

次の[[初期値問題]]を考える。
: <math>  y' = f(t, y), \quad y(t_0) = y_0. </math>
但し、{{math|''y''(''t'')}} が近似的に求めたい未知関数であり、その {{mvar|t}} における勾配は {{math|''f''(''t'', ''y'')}} によって {{mvar|t}} 及び {{math|''y''(''t'')}} の関数として与えられている。時刻 {{mvar|t{{sub|0}}}} における初期値は {{mvar|y{{sub|0}}}} で与えられている。

今、時刻 {{mvar|t{{sub|n}}}} における値 {{mvar|y{{sub|n}}}} = {{math|''y''(''t{{sub|n}}'')}} が既知のとき、十分に小さなステップ幅 {{mvar|h}} に対して {{mvar|y{{sub|n+1}}}}, {{mvar|t{{sub|n+1}}}} を以下の式で与えると、{{mvar|y{{sub|n+1}}}} は {{math|''y''(''t{{sub|n+1}}'')}} の 4次精度の近似になっている。このステップを逐次的に繰り返すことによって、初期値 {{mvar|y{{sub|0}}}} から任意の時刻 {{mvar|t{{sub|n}}}} における近似値 {{mvar|y{{sub|n}}}} が求められる。
: <math>\begin{align}
  y_{n+1} &= y_n + {h \over 6} (k_1 + 2k_2 + 2k_3 + k_4), \\
  t_{n+1} &= t_n + h.
\end{align}</math>
ここで、
: <math>\begin{align} k_1 &= f \left( t_n, y_n \right), \\
k_2 &= f \left( t_n + {h \over 2}, y_n + {h \over 2} k_1 \right), \\
k_3 &= f \left( t_n + {h \over 2}, y_n + {h \over 2} k_2 \right), \\
k_4 &= f \left( t_n + h, y_n + hk_3 \right) \end{align}</math>
である{{sfn|Press|Teukolsky|Vetterling|Flannery|2007|p=908}}{{sfn|Süli|Mayers|2003|p=328}}。次の値 ({{math|''y''{{sub|''n''+1}}}}) は、現在の値 ({{mvar|y{{sub|n}}}}) に増分を加えたものであり、増分は勾配の推定値に間隔 {{mvar|h}} を乗じたものになっている。勾配の推定値は、{{mvar|k{{sub|1}}}}, ..., {{mvar|k{{sub|4}}}} の4つの勾配の重み付け平均で求める。{{mvar|k{{sub|1}}}}, ..., {{mvar|k{{sub|4}}}} のそれぞれの勾配は、特定の {{math|(''t'', ''y'')}} に対する {{mvar|f}} によって与えられ、以下のように解釈できる。
* {{mvar|k{{sub|1}}}} は区間の最初 {{mvar|t{{sub|n}}}} における勾配である ([[オイラー法]]で用いる勾配に一致する)。
* {{mvar|k{{sub|2}}}} は区間の中央 {{math|''t{{sub|n}}'' + ''h''/2}} における勾配の近似値である ([[中点法]]で用いる勾配)。計算に用いる中央の {{mvar|y}} の値は、初期位置の勾配 {{mvar|k{{sub|1}}}} を用いて[[オイラー法]]で推定する。
* {{mvar|k{{sub|3}}}} は区間の中央における勾配のもう一つの近似値である。中央の {{mvar|y}} を {{mvar|k{{sub|2}}}} の値から推定して用いる。
* {{mvar|k{{sub|4}}}} は区間の最後 {{math|''t{{sub|n}}'' + ''h''}} における勾配の近似値であり、{{mvar|k{{sub|3}}}} の値から推定された最後の点の {{mvar|y}} の値を用いる。
重み付き平均では、中央の勾配に対して大きな重みを用いる。[[シンプソン則]]を用いた平均と同等の形になる<ref name="Süli 2003 328">{{harvnb|Süli|Mayers|2003|p=328}}.</ref>。
: <math>\mbox{slope} = \frac{k_1 + 2k_2 + 2k_3 + k_4}{6}.</math>
RK4は4次の方法である。厳密解とRK4の[[テイラー展開]]が4次の項まで一致し、1ステップの推定誤差は <math>O(h^5)</math> の[[ランダウの記号|オーダー]]になる。目的の時刻の {{mvar|y}} を求めるのに必要なステップ数は <math>O(h^{-1})</math> になるので、全体の推定誤差は <math>O(h^4)</math> になる。

== ルンゲ゠クッタ法 ==
前述の RK4 の一般化として、以下の形式を持つ {{mvar|s}} 段のルンゲ゠クッタ法を構成することができる。整数 {{mvar|s}} をそのルンゲ゠クッタ法の'''段数''' (stage) という。
: <math> y_{n+1} = y_n + h\sum_{i=1}^s b_i k_i, </math>
但し<ref name=":0">{{harvnb|Iserles|2008|p=41}}; {{harvnb|Süli|Mayers|2003|pp=351–352}}.</ref>、
: <math>k_i = f\left(t_n + hc_i, y_n + h\sum_{j=1}^s a_{ij} k_j\right), \quad i = 1,\ldots,s.</math>
: (文献によって、等価であるが上と異なる定義の仕方をしているものがあることに注意する)<ref name=notation>{{harvtxt|Atkinson|1989|p=423}}, {{harvtxt|Hairer|Nørsett|Wanner|1993|p=134}}, {{harvtxt|Kaw|Kalu|2008|loc=§8.4}} and {{harvtxt|Stoer|Bulirsch|2002|p=476}} leave out the factor ''h'' in the definition of the stages. {{harvtxt|Ascher|Petzold|1998|p=81}}, {{harvtxt|Butcher|2008|p=93}} and {{harvtxt|Iserles|2008|p=38}} use the ''y'' values as stages.</ref>。
具体的なルンゲ゠クッタ公式は、解ができるだけ高い精度を持つように適切に選ばれた係数 {{mvar|a{{sub|ij}}}} ('''ルンゲ゠クッタ行列''')、{{mvar|b{{sub|i}}}} ('''重み''')、{{mvar|c{{sub|j}}}} ('''節点''') で指定される (<math> i,j \le s </math>)<ref name="Iserles200838">{{harvnb|Iserles|2008|p=38}}.</ref>。特に <math> i \le j</math> に対して <math> a_{ij} = 0</math> を満たす方法が広く用いられ、総称して '''陽的ルンゲ゠クッタ法''' (ERK、{{lang-en-short|explicit Runge&ndash;Kutta methods}}) と呼ぶ。そうでないものを '''陰的ルンゲ゠クッタ法''' (IRK、{{lang-en-short|implicit Runge&ndash;Kutta methods}}) と呼ぶ。

近似値 {{mvar|y{{sub|n}}}} を {{mvar|y{{sub|n+1}}}} から計算するときに発生する誤差の大きさが <math>O(h^{p+1})</math> のとき、そのルンゲ゠クッタ公式は {{mvar|p}} 次精度を持つといい、{{mvar|p}} を'''次数''' (または'''位数''') と呼ぶ。{{mvar|p}} 次のルンゲ゠クッタ公式は、誤差の大きさの条件に誤差の表式を代入し、係数の条件を求めることによって得られる。例えば、2段の陽的方法が2次精度を持つための係数に対する条件は、<math> b_1 + b_2 =1 </math>, <math> b_2c_2 = 1/2 </math>, かつ <math> b_2a_{21} = 1/2</math> である<ref name="Iserles200839">{{harvnb|Iserles|2008|p=39}}.</ref>。この条件を満たす範囲内で様々な方法を考え得る。

これらの係数を分かりやすく表記する方法として、以下のような形式の'''ブッチャー配列''' (''Butcher tableau'') が知られている<ref name="Süli2003352">{{harvnb|Süli|Mayers|2003|p=352}}.</ref>:
: <math>\begin{align}
  \begin{array}{c|ccccc}
  c_1 & a_{11} & a_{12} & \cdots & a_{1s} \\
  c_2 & a_{21} & a_{22} & \cdots & a_{2s} \\
  \vdots & \vdots & \vdots & \ddots & \vdots \\
  c_s & a_{s1} & a_{s2} & \cdots & a_{ss} \\ \hline
      & b_1 & b_2 & \cdots & b_s
  \end{array} &= \begin{array}{c|c} \mathbf{c} & A \\ \hline & \mathbf{b}^\mathrm{T} \end{array}
\end{align}</math>
実際には、それぞれのルンゲ゠クッタ法について各要素に具体的な値を入れて用いる。陽的な方法ではルンゲ゠クッタ行列の上三角成分は常に 0 になるので表記を省略する。例えば RK4 はブッチャー配列を用いて以下のように表現される。
: <math>\begin{align}
  \begin{array}{c|ccccc}
  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}
\end{align}</math>
ルンゲ゠クッタ法が'''一貫'''しているためには、次の条件が満たされている必要がある。
: <math>\sum_{j=1}^{s} a_{ij} = c_i, \; i=1, \ldots, s.</math>

=== 収束性 ===
ルンゲ＝クッタ法は、数値積分による[[求積法]] (quadrature) と深く関係にある。時刻 {{mvar|t{{sub|n}}}} での値から {{math|1=''t{{sub|n+1}}'' = ''t{{sub|n}}'' + ''h''|''t{{sub|n+1}}'' = }} での値を求めるときの方程式が以下であるとする。
: <math>y' = f(t,y), \; y(t_n) = y_n.</math>
求積法は、与えられた区間での[[定積分]]の値を被積分関数の値の[[線型結合]]により近似する方法である。簡単のために、区間を {{math|[0, 1]}} とする。すると求積法の公式は
: <math>\int_0^1 g(t)dt \approx \sum_{i=1}^{s} b_i g(c_i)</math>
となる。ここで、{{mvar|b{{sub|i}}}} と {{mvar|c{{sub|i}}}} はあらかじめ選ばれた定数であり、前述の重みと節点に対応する。上記の式に対し、すべての {{math|''p'' − 1}} 次以下の多項式 g に対して等式が成立するとき（すなわち、誤差が0のとき）には、その求積法は {{mvar|p}} 次精度であり、{{mvar|p}} を次数と呼ぶ{{sfn|Iserles|2008|p=34}}。節点が {{mvar|s}} 個の時、最大次数は 2{{mvar|s}} であり、その方法は {{mvar|s}} 次[[ガウス求積|ガウス・ルジャンドル公式]]と呼ばれる。

そこで上記の常微分方程式を積分形式に変形して、求積法を被積分関数 <math>g(t)=f(t,y(t))</math> に対して適用すると以下の近似式が得られる<ref name = "Iserles200838"/>。
: <math>y(t_{n+1}) \approx y_n + \int_{y_n}^{y_{n+1}} f(t,y(t)) dt = y_n + h\int_0^1 y(t_n+h\tau, y(t_n+h\tau))d\tau = y_n + h\sum_{i=1}^s b_if(t_n + c_i h, y(t_n + c_i h)) </math>
<math> k_i = f(t_n + c_i h, y(t_n + c_i h)) </math> とおく。{{mvar|k{{sub|i}}}} あるいは {{math|''y''(''t{{sub|n}}'' + ''c{{sub|i}}'' ''h'')}} を適切に（線型結合により）近似することでルンゲ＝クッタ法の公式が得られる。さらに、テイラー展開より係数を正しく選択すれば、方法の収束性は求積法の収束性から保証される。しかし、局所誤差のオーダーや上界は、方法によって大きく異なるので、方法別に計算しなければならない。

== 陽的ルンゲ゠クッタ法 ==
陽的ルンゲ゠クッタ法では <math>a_{ij}=0</math> (<math> i \le j </math>) が満たされるので、ブッチャー配列は以下の形になる。
{| cellspacing="0px" cellpadding="3px"
| width="20px" |
| style="border-right:1px solid;" |0
|-
|
| style="border-right:1px solid;" |<math> c_2 </math>
|<math> a_{21} </math>
|-
|
| style="border-right:1px solid;" |<math> c_3 </math>
|<math> a_{31} </math>
|<math> a_{32} </math>
|-
|
| style="border-right:1px solid;" |<math> \vdots </math>
|<math> \vdots </math>
|
|<math> \ddots </math>
|-
|
| style="border-right:1px solid; border-bottom:1px solid;" |<math> c_s </math>
| style="border-bottom:1px solid;" |<math> a_{s1} </math>
| style="border-bottom:1px solid;" |<math> a_{s2} </math>
| style="border-bottom:1px solid;" |<math> \cdots </math>
| style="border-bottom:1px solid;" |<math> a_{s,s-1} </math>
| style="border-bottom:1px solid;" |
|-
|
| style="border-right:1px solid;" |
|<math> b_1 </math>
|<math> b_2 </math>
|<math> \cdots </math>
|<math> b_{s-1} </math>
|<math> b_s </math>
|}
この時、各勾配 {{math|''k''{{sub|1}}}}, ..., {{mvar|k{{sub|s}}}} を以下のように逐次的に求めることができる{{sfn|Press|Teukolsky|Vetterling|Flannery|2007|p=907}}。
: <math>\begin{align} k_1 &= f(t_n, y_n), \\
k_2 &= f(t_n+c_2h, y_n+a_{21}hk_1), \\
k_3 &= f(t_n+c_3h, y_n+a_{31}hk_1+a_{32}hk_2), \\
&\vdots\\
k_s &= f(t_n+c_sh, y_n+a_{s1}hk_1+a_{s2}hk_2+\cdots+a_{s,s-1}hk_{s-1}). \end{align}</math>
: (文献によって、等価であるが上と異なる定義の仕方をしているものがあることに注意する)<ref name=notation />。
陽的ルンゲ゠クッタ法が目的の次数 {{mvar|p}} を実現するためには十分な段数 {{mvar|s}} が必要になる。少なくとも次数以上の段数が必要 (<math> s \ge p </math>) であり、更に5次以上の場合には段数を次数より大きく取らねばならない (<math> s > p </math>) ことが知られている{{sfn|Butcher|2008|p=187}}。各次数 {{mvar|p}} に対してそれを実現するための最小の段数 {{mvar|s}} を与える一般式は未解決問題である。具体的な値が幾つか知られている{{sfn|Butcher|2008|pp=187–196}}:
: <math>
\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} </math>

各段数 {{mvar|s}} に対して達成が可能な最大の次数 {{mvar|ｐ}} は以下の表になる。
: <math>
\begin{array}{c|cccccccccc}
s & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 & 10 & 11 \\
\hline
\max p & 1 & 2 & 3 & 4 & 4 & 5 & 6 & 6 & 7 & 7 & 8

\end{array} </math>


[[File:Error of ODE Solvers (Japanese label).svg|thumb|320px|オイラー法、ホイン法、古典的ルンゲ＝クッタ法 (RK4) の相対誤差の比較。初期値 <math>y ( 0 ) = 0</math> のもとでの常微分方程式 <math>\frac{ d y }{ d t } = \cos ( y )</math> の数値解の <math>t = 1</math> での値をステップサイズの関数として対数プロットした。各手法がそれぞれ1次精度、2次精度、4次精度であることに対応して、両対数グラフ上で直線の傾きがそれぞれ 1, 2, 4 でステップサイズが小さくなるにつれて誤差が減少している<ref>{{Cite book |和書 |author=齊藤宣一 |title=数値解析 （共立講座 数学探検 17） |publisher=共立出版 |date=2017 |isbn=9784320992740 |page=107}}</ref>。]]

=== 1段1次の公式 ===
最も単純なルンゲ゠クッタ法は 1段1次 の方法で一意に定まり、'''(前進) オイラー法'''に一致する。
: <math> y_{n+1} = y_n + hf(t_n,y_n).</math>
以下のブッチャー配列で表現される。
{| cellspacing="3px" cellpadding="3px"
| width="20px" |
| style="border-right:1px solid; border-bottom:1px solid;" width="10px" | 0
| style="border-bottom:1px solid;" width="10px" | 
|-
|
| style="border-right:1px solid;" |
| 1
|}

=== 2段2次の公式 ===
2段2次の方法は1パラメータの自由度 {{mvar|α}} を持ち、以下のブッチャー配列で与えられる。
{| cellspacing="3px" cellpadding="3px"
| width="20px" |
| style="border-right:1px solid;" |0
|-
|
| style="border-right:1px solid; border-bottom:1px solid;" |<math>\alpha</math>
| style="border-bottom:1px solid;text-align:center;" |<math>\alpha</math>
| style="border-bottom:1px solid;text-align:center;" |
|-
|
| style="border-right:1px solid;" |
|<math>(1-\tfrac1{2\alpha})</math>
|<math>\tfrac1{2\alpha}</math>
|}
対応する公式は{{sfn|Süli|Mayers|2003|p=327}}
: <math> y_{n+1} = y_n + h\bigl( (1-\tfrac1{2\alpha}) f(t_n, y_n) + \tfrac1{2\alpha} f(t_n + \alpha h, y_n + \alpha h f(t_n, y_n)) \bigr). </math>
である。<math> \alpha = 1/2 </math> の場合が [[中点法]] (または'''修正オイラー法''', modified Euler method)
: <math> y_{n+1} = y_n + hf\left(t_n+\frac{h}{2},y_n+\frac{h}{2}f(t_n, y_n)\right)</math>
に対応し、以下のブッチャー配列で与えられる。
{| cellspacing="0px" cellpadding="3px"
| width="20px" |
| style="border-right:1px solid;" |0
|-
|
| style="border-right:1px solid; border-bottom:1px solid;" |1/2 
| style="border-bottom:1px solid;" |1/2
| style="border-bottom:1px solid;" |
|-
|
| style="border-right:1px solid;" | 
|0 
|1
|}
<math> \alpha = 1 </math> の場合は {{仮リンク|ホイン法|en|Heun's method}}（または'''改良オイラー法''', improved Euler method)として知られ、ブッチャー配列は以下の通りである<ref name="Süli 2003 328"/>。
{| cellspacing="0px" cellpadding="3px"
| width="20px" |
| style="border-right:1px solid;" | 0
|-
|
| style="border-right:1px solid; border-bottom:1px solid;" | 1
| style="border-bottom:1px solid;" | 1
| style="border-bottom:1px solid;" |
|-
|
| style="border-right:1px solid;" |
|1/2
|1/2
|}

=== 4段4次の公式 ===
古典的ルンゲ゠クッタ法 (RK4) は既に述べた通り以下のブッチャー配列で与えられる<ref name="Süli2003352"/>:
{| cellspacing="0px" cellpadding="3px"
| width="20px" |
| style="border-right:1px solid;" |0
|-
|
| style="border-right:1px solid;" |1/2
|1/2
|-
|
| style="border-right:1px solid;" |1/2
|0
|1/2
|-
|
| style="border-right:1px solid; border-bottom:1px solid;" |1
| style="border-bottom:1px solid;" |0
| style="border-bottom:1px solid;" |0
| style="border-bottom:1px solid;" |1
| style="border-bottom:1px solid;" |
|-
|
| style="border-right:1px solid;" |
|1/6
|1/3
|1/3
|1/6
|}
修正版として'''クッタの3/8公式''' ({{lang-en-short|Kutta's 3/8-rule}}) が知られている<ref>{{harvtxt|Hairer|Nørsett|Wanner|1993|p=138}} refer to {{harvtxt|Kutta|1901}}.</ref>。
{| cellspacing="0px" cellpadding="3px"
| width="20px" |
| style="border-right:1px solid;" |0
|-
|
| style="border-right:1px solid;" |1/3
|1/3
|-
|
| style="border-right:1px solid;" |2/3
| -1/3
|1
|-
|
| style="border-right:1px solid; border-bottom:1px solid;" |1
| style="border-bottom:1px solid;" |1
| style="border-bottom:1px solid;" | -1
| style="border-bottom:1px solid;" |1
| style="border-bottom:1px solid;" |
|-
|
| style="border-right:1px solid;" |
|1/8
|3/8
|3/8
|1/8
|}

他に使われている方法は、[[ルンゲ＝クッタ法のリスト]]に参照。

== 計算例：2段2次陽的方法の条件の導出 ==
前述の通り、2段の陽的方法が2次精度を持つための係数に対する条件は、<math> b_1 + b_2 =1 </math>, <math> b_2c_2 = 1/2 </math>, かつ <math> b_2a_{21} = 1/2</math> である<ref name="Iserles200839" />。例として、それらの条件の導出を見てみよう。

一般的に2段陽的ルンゲ＝クッタ法に対応する配列は以下の通りである(一貫性を用いて <math>c_1 = 0</math> がわかる)。 

: <math>
\begin{array}{c|cc}
0 & &\\
c_2 & a_{21}&\\
\hline
  &  b_1&b_2\\
\end{array}
</math>

公式から <math>k_2 = f(t_n + c_2h, y_n + ha_{21} f(t_n,y_n))</math> の {{math|''f''(''t{{sub|n}}'',''y{{sub|n}}'')}} に関するテイラー展開を考える
: <math> k_2 = f(t_n,y_n) + c_2h f_t(t_n,y_n) + ha_{21}f(t_n,y_n)f_y(t_n,y_n) + O(h^2) </math>
ここで、{{mvar|f{{sub|t}}}} は {{mvar|f}}  の {{mvar|t}} に関する偏微分である。それを <math> y_{n+1} = y_n + h(b_1k_1 + b_2k_2) </math> に代入し、<math> k_1 = f(t_n,y_n) </math> を用いて整理すると
: <math> y_{n+1} = y_n + h(b_1 + b_2)f(t_n,y_n) + h^2b_2(c_2f_t(t_n,y_n) + a_{21}f(t_n,y_n)f_y(t_n,y_n)) + O(h^3) </math>
となる。同時に、方程式 <math> y'(t) = f(t,y) </math> の {{mvar|t}} に関する微分を取って {{mvar|y'}} を {{math|''f''(''t'',''y'')}} に置き換えると
: <math> y''(t) = f_t(t,y) + f(t,y)f_y(t,y) </math>
となる。厳密解を {{math|''y''(''t'')}} とする。上記式を用いて {{math|''y''(''t{{sub|n+1}}'')}} の {{mvar|t{{sub|n}}}} に関するテイラー展開を考える
: <math> y(t_{n+1}) = y(t_n) + hf(t_n,y_n) + \frac{1}{2}h^2(f_t(t_n,y_n) + f(t_n,y_n)f_y(t_n,y_n)) + O(h^3) </math>
2次精度を持つ条件は、局所誤差 <math> e_{n+1,h} = y(t_{n+1}) - y_{n+1} = O(h^3) </math> である。二つの展開式の係数を比較すると、その条件が<math> b_1 + b_2 =1 </math>, <math> b_2c_2 = 1/2 </math>, かつ <math> b_2a_{21} = 1/2</math> であることがわかる。

== 埋め込み型ルンゲ＝クッタ法 ==
ルンゲ＝クッタ法の局所誤差を正確に求めることは難しいので、計算を行なう場合には誤差が一定の範囲に収まるように制御できることが望ましい。そのために開発されたのが '''埋め込み型ルンゲ＝クッタ法'''（Embedded Runge-Kutta method）である。'''適応型ルンゲ＝クッタ法''' (adaptive Runge–Kutta method) とも呼ばれる。線型多段法にもミルンデバイスと呼ばれる相似の方法が存在する{{sfn|Iserles|2008|p=107}}。

埋め込み型の方法では二つの陽的ルンゲ＝クッタ法を用いる（下の方法を上の方法に埋め込むように見えるので埋め込み型と呼ばれる）。ブッチャー配列を以下のように拡張する。

{| cellpadding=3px cellspacing=0px
|width="20px"| || style="border-right:1px solid;" | 0
|-
||| style="border-right:1px solid;" | <math> c_2 </math> || <math> a_{21} </math>
|-
||| style="border-right:1px solid;" | <math> c_3 </math> || <math> a_{31} </math> || <math> a_{32} </math>
|-
||| style="border-right:1px solid;" | <math> \vdots </math> || <math> \vdots </math> || || <math> \ddots </math>
|-
||| style="border-right:1px solid; border-bottom:1px solid;" | <math> c_s </math>
| style="border-bottom:1px solid;" | <math> a_{s1} </math>
| style="border-bottom:1px solid;" | <math> a_{s2} </math>
| style="border-bottom:1px solid;" | <math> \cdots </math>
| style="border-bottom:1px solid;" | <math> a_{s,s-1} </math> || style="border-bottom:1px solid;" |
|-
||| style="border-right:1px solid;" | || <math> b_1 </math> || <math> b_2 </math> || <math> \cdots </math> || <math> b_{s-1} </math> || <math> b_s </math>
|-
||| style="border-right:1px solid;" | || <math> b^*_1 </math> || <math> b^*_2 </math> || <math> \cdots </math> || <math> b^*_{s-1} </math> || <math> b^*_s </math>
|}
ここで、{{mvar|b{{sub|i}}}} は {{mvar|p}} 次陽的方法に対応し、{{math|{{SubSup|''b''|''i''|*}}}} は {{math|''p''-1}} 次陽的方法に対応する。二つの方法は係数 {{mvar|a{{sub|ij}}}} と {{mvar|c{{sub|i}}}} を共用する。係数を正しく選べば、二つの方法をともに収束させることができる。そのとき、埋め込み方法の時刻 {{mvar|t{{sub|n}}}} での相対（局所）誤差 {{Mvar|e{{sub|n}}}} は次の公式で与えられる。
: <math> e_{n} = y_n - y^*_n = h \sum_{i=1}^s (b_i - b^*_i)k_i = O(h^p) </math>
ルンゲ＝クッタ法の収束性から、相対誤差も0に収束することがわかる。埋め込み型ルンゲ＝クッタ法は、アルゴリズムを用いて一時刻ごと（自動的）に刻み幅 {{mvar|h}} を調整することで誤差を制御できる（故に適応型と呼ばれる）{{sfn|Iserles|2008|p=109}}。よって絶対誤差を求めなくても刻み幅を適切に制御することができる。そのアルゴリズムは、大体以下の通りである。

誤差の許容値を {{mvar|δ}} とする。刻み幅 {{mvar|h}} の近似値 {{mvar|y{{sub|n+1}}}} と {{math|{{SubSup|''y''|''n''+1|*|}}}} をルンゲ＝クッタ法で計算する。二つの値に対し、<math> \|e_{n+1}\| \ge h\delta </math> が成立するとき、{{mvar|h}} が大きすぎて小さくする必要がある。よって新しい刻み幅 <math>h' = \frac{1}{2}h</math> と設定し再び近似値を計算する。代わりに <math> \|e_{n+1}\| \le \frac{1}{10}h\delta </math> が成立するとき、{{mvar|h}} が小さすぎて大きくするほうが効率が良い。よって新しい刻み幅 <math>h' = 2h</math> と設定し、次の時刻での近似値を計算する（この場合、誤差が許容値より小さいので今の近似値を二度と計算する必要はない）。二つの不等式がどちらも成立しないとき、そのままの刻み幅で次の時刻での近似値を計算してよい。

ただし、使用したスカラー（1/2他）の値は方程式や方法によって変更することもできる。また、実際の計算では {{mvar|h}} が小さくなり過ぎることを避けるために、下界 {{mvar|h{{sub|min}}}} を先に設定することもある。<math> h < h_{\mathrm{min}} </math> のときにはアルゴリズムを停止して、別のより次数の高い方法を使う。 

=== 例 ===
もっとも単純な埋め込み型ルンゲ＝クッタ法は、ホイン法（2次）とオイラー法（1次）を組み合わせる方法である。以下の拡張ブッチャー配列で与えられる。

{| cellpadding=3px cellspacing=0px
|width="20px"| || style="border-right:1px solid;" | 0
|-
||| style="border-right:1px solid; border-bottom:1px solid;" | 1 || style="border-bottom:1px solid;" | 1 || style="border-bottom:1px solid;" |
|-
||| style="border-right:1px solid;" | || 1/2 || 1/2
|-
||| style="border-right:1px solid;" | || 1 || 0
|}

実際の計算でよく用いられているのは、5次と4次の陽的方法を組み合わせる [[ルンゲ＝クッタ＝フェールベルグ法]] である。対応する拡張ブッチャー配列は以下の通りである{{sfn|Hairer|Nørsett|Wanner|1993|p=177}}。

{| cellpadding=3px cellspacing=0px
|width="20px"| || style="border-right:1px solid;" | 0
|-
||| style="border-right:1px solid;" | 1/4 || 1/4
|-
||| style="border-right:1px solid;" | 3/8 || 3/32 || 9/32
|-
||| style="border-right:1px solid;" | 12/13 || 1932/2197 || −7200/2197 || 7296/2197
|-
||| style="border-right:1px solid;" | 1 || 439/216 || −8 || 3680/513 || -845/4104
|-
||| style="border-right:1px solid; border-bottom:1px solid;" | 1/2 || style="border-bottom:1px solid;" | −8/27 || style="border-bottom:1px solid;" | 2 || style="border-bottom:1px solid;" | −3544/2565 || style="border-bottom:1px solid;" | 1859/4104 || style="border-bottom:1px solid;" | −11/40 || style="border-bottom:1px solid;" |
|-
||| style="border-right:1px solid;" | || 16/135 || 0 || 6656/12825 || 28561/56430 || −9/50 || 2/55
|-
||| style="border-right:1px solid;" | || 25/216 || 0 || 1408/2565 || 2197/4104 || −1/5 || 0
|}
使われている他の方法は、{{仮リンク|Bogacki–Shampine法|en|Bogacki–Shampine method}}(次数3と2)、{{仮リンク|Cash–Karp法|en|Cash–Karp method}} と [[ドルマン＝プリンス法]]（どちらの方法も次数の組合せは5と4）である。

== 陰的ルンゲ＝クッタ法 ==
いままでに述べた方法はすべて陽公式である。陽的ルンゲ＝クッタ方法の絶対安定性領域（region of absolute stability）が小さくて有界であるため{{sfn|Süli|Mayers|2003|pp=349–351}}、[[硬い方程式]]の解を計算する場合、方法は不適切である。この問題は特に[[偏微分方程式]]の解を計算するときに重要である。

陽的方法の不安定さは陰的ルンゲ＝クッタ法の開発の動機となる。陰的ルンゲ＝クッタ法は上述の一般的な公式と同じく、以下の形で与えられる
: <math> y_{n+1} = y_n + h\sum_{i=1}^s b_i k_i, </math>
但し、
: <math>k_i = f\left(t_n + hc_i, y_n + h\sum_{j=1}^s a_{ij} k_j\right), \quad i = 1,\ldots,s.</math>
であり<ref name=":0" />、対応するルンゲ＝クッタ行列は厳密な[[下三角行列]]ではない。結果として、一時刻ごとに[[代数方程式|代数方程式系]]を解かなければならなくなる。応じて、計算コストもかなり上がる。{{Mvar|s}} 段法を使って {{Mvar|m}} 個の成分からなる微分方程式系（すなわち <math>\vec{y} = (y_1, \ldots, y_m)</math> のとき）に適用する場合に対応する代数方程式の数は {{Mvar|ms}} となる。これは陰的[[線型多段法]]とも比較できる：{{Mvar|s}} 段陰的線型多段法を使って同じ方程式系に適用する場合に対応する代数方程式の数は {{Mvar|m}} であり、段数 {{Mvar|s}} に依存しない<ref name="Süli 2003 353">{{harvnb|Süli|Mayers|2003|p=353}}.</ref>。

=== 例 ===
一番単純な陰的ルンゲ＝クッタ法は、[[後退オイラー法]]
: <math>y_{n+1} = y_n + hf(t_{n+1},y_{n+1})</math>
である。対応するブッチャー配列も単純に以下の通りである。
: <math>
\begin{array}{c|c}
1 & 1 \\
\hline
  & 1 \\
\end{array}
</math>
他の例として、[[台形公式]]と呼ばれる方法は以下の公式で与えられる。
: <math> y_{n+1} = y_n + \frac{1}{2}h(f(t_n,y_n) + f(t_{n+1},y_{n+1})) </math>
対応する配列は以下の通りである。
: <math>
\begin{array}{c|cc}
0 & 0& 0\\
1 & \frac{1}{2}& \frac{1}{2}\\
\hline
  &  \frac{1}{2}&\frac{1}{2}\\
\end{array}
</math> 
その対応は、<math> y_{n+1} = y_n + \frac{1}{2}k_1 + \frac{1}{2}k_2 </math> を覚えてて、公式を以下のように書き換えると明らかになる。
: <math> y_{n+1} = y_n + \frac{1}{2}f(t_n,y_n) + \frac{1}{2}f\left(t_n + h, y_n + \frac{1}{2}k_1 + \frac{1}{2}k_2\right) </math>
台形公式は、[[選点法]]である。選点法はすべて陰的ルンゲ＝クッタ法であるけど、陰的ルンゲ＝クッタ法がすべて選点法であるわけではない{{sfn|Iserles|2008|pp=43–44}}。

選点法の中では、[[ガウス求積]]に基づいた{{仮リンク|ガウス・ルジャンドル法|en|Gauss–Legendre method}}が一番次数の高い方法である。{{Mvar|s}} 段ガウス・ルジャンドル法の次数は {{math|2''s''}} となる（よって、任意に高い次数を持つ方法を構造できるようになる）<ref name = "Iserles200847">{{harvnb|Iserles|2008|p=47}}.</ref>。例として、2段ガウス・ルジャンドル法は次の配列で与えられる。
: <math>
\begin{array}{c|cc}
\frac12 - \frac16 \sqrt3 & \frac14                  & \frac14 - \frac16 \sqrt3 \\
\frac12 + \frac16 \sqrt3 & \frac14 + \frac16 \sqrt3 & \frac14 \\
\hline
                         & \frac12                  & \frac12 \\
\end{array}
</math> <ref name = "Iserles200847"/>

== 安定性 == 
数値分析における[[数値的安定性|安定性]]は、それぞれ異なる定義が複数存在する。ルンゲ＝クッタ法の安定性を反映できる概念は、主に以下の二つである。
=== A-安定性 ===
線型テスト方程式 <math> y' = \lambda y </math> を考える。一つのルンゲ＝クッタ法を使ってこの方程式に適用すると <math> y_{n+1} = r(h\lambda)y_n </math> となる。ここで、
: <math> r(z) = 1 + zb^T(I-zA)^{-1}e = \frac{\det (I-zA+zeb^T)}{\det(I-zA)} </math> {{sfn|Hairer|Wanner|1996|pp=40–41}}
は安定性関数と呼ばれる {{math|'''C'''}} 上の[[有理関数]]である（{{Mvar|e}} はすべての成分が 1 の[[ベクトル]]{{要曖昧さ回避|date=2021年7月}}である）{{sfn|Hairer|Wanner|1996|p=40}}。{{Mvar|s}} 段法の場合、[[行列式]]の展開によって {{Math|''r''(''z'')}} は二つの {{Mvar|s}} 次[[多項式]]の商となる。陽的方法の場合、対応するルンゲ＝クッタ行列が狭義下三角行列であるため、<math>\det(I-zA) = 1</math> がわかる。したがって、{{Math|''r''(''z'')}} は {{Mvar|s}} 次多項式となる<ref name="Iserles200860">{{harvnb|Iserles|2008|p=60}}.</ref>。

ルンゲ＝クッタ法による上記の方程式の数値解がゼロに減衰する条件が <math>|r(z)| < 1</math> (<math>z = h\lambda</math>) である。上記の条件を満たす {{Mvar|z}} の[[集合]]は方法に対する '''絶対安定性領域''' （region of absolute stability） である。 特に、絶対安定性領域が左[[複素数平面]]を含むとき、その方法は '''A-安定''' (A-stable) という。陽的ルンゲ＝クッタ法は、安定性関数が多項式であるため、A-安定な方法になれない<ref name="Iserles200860"/>。

もっと一般的に、方法の次数が {{Mvar|p}} のときに、安定性関数に対し <math>r(z) = e^z + O(z^{p+1})</math> が成立する。ゆえに、[[指数関数]] {{Math|''e''{{sup|''z''}}}} の、与えられた次数の多項式の商からなる有理関数の中での最善近似が重要だと考えられる。そのような有理関数は[[パデ近似]]と呼ばれる。分子の次数が {{Mvar|m}} で分母の次数が {{Mvar|n}} のパデ近似式がA-安定性の条件 <math>|r(z)| < 1</math>, <math> \mathrm{Re}(z) \le 0</math><ref>{{harvtxt|Iserles|2008}} はそのような有理関数を A-acceptable と呼ぶ。</ref> を満足するための[[必要十分条件]]は、<math>m \le n \le m+2</math> である{{sfn|Iserles|2008|pp=62–63}}。

{{Mvar|s}} 段ガウス・ルジャンドル法の次数は前述通り 2{{Mvar|s}} である。よって安定性関数の分子と分母は同じく {{Mvar|s}} 次多項式となる。すなわち、 <math>m = n = s</math> である。上記の条件が満たされるので、ガウス・ルジャンドル法はA-安定であることがわかる{{sfn|Iserles|2008|p=63}}。故に任意に高い次数を持つ、A-安定なルンゲ＝クッタ法が存在する。比べて、A-安定性を持つ線型多段法の次数は2以下である<ref>この結果（時々にsecond Dahlquist barrier）は、 {{harvtxt|Dahlquist|1963}} によるものである。</ref>。

以上のことから、陰的ルンゲ＝クッタ法は陽的な方法よりも優れた安定性を持つこともわかる。

=== B-安定性 ===
A-安定性という概念は線型[[自励系|自励方程式]] <math> y' = \lambda y </math> を考える結果である。 {{仮リンク|ダールキスト|en|Germund Dahlquist}}は、数値的方法を使って、とある単調性条件を満たす非線型方程式系に適用するときの安定性を主張した。対応する安定性は、線型多段法の場合に G-安定性 （G-stability）で、ルンゲ＝クッタ法の場合に B-安定性 (B-stability){{sfn|Butcher|1975}} と呼ぶ。 一般的なテスト方程式 <math>y' = f(x,y)</math> と単調性条件 <math>\langle f(x,y) - f(x,z), y-z \rangle \le 0</math> を満足する {{Mvar|f}} を考える。もしすべての <math>h \ge 0</math> に対し、 不等式 <math>\|y_{n+1} - z_{n+1}\| \le \|y_n - z_n\|</math> が成立するとき、そのルンゲ＝クッタ法は '''B-安定''' という。 ここで、{{Mvar|''y{{sub|n}}''}} と {{Mvar|''z{{sub|n}}''}} はそれぞれの初期値に対する数値解である。 <math>f(x,y) = \lambda y</math> に方法を適用することで、B-安定性はA-安定性より強い条件であることがわかる{{sfn|Hairer|Wanner|1996|p=181}}。

== 脚注 ==
{{reflist|30em}}

== 参考文献 ==
* {{Citation | last1=Atkinson | first1=Kendall A. | title=An Introduction to Numerical Analysis | publisher=[[John Wiley & Sons]] | location=New York | edition=2nd | isbn=978-0-471-50023-0 | year=1989}}.
* {{Citation| last1 = Butcher | first1 = John C. | author1-link=John C. Butcher | title=Coefficients for the study of Runge-Kutta integration processes | url=http://journals.cambridge.org/action/displayAbstract?fromPage=online&aid=4922056 | doi=10.1017/S1446788700027932 |journal = Journal of the Australian Mathematical Society | volume=3 | issue=2 |date=May 1963 | pages=185–201}}.
* {{Citation | last1=Butcher | first1=John C. | author1-link=John C. Butcher | title= A stability property of implicit Runge-Kutta methods| year=1975 |journal= BIT| volume= 15 |pages= 358–361 | doi=10.1007/bf01931672}}.
* {{Citation | last1=Butcher | first1=John C. | author1-link=John C. Butcher | title=Numerical Methods for Ordinary Differential Equations | publisher=[[John Wiley & Sons]] | location=New York | isbn=978-0-470-72335-7 | year=2008}}.
* {{Citation | last1=Dahlquist | first1=Germund | author1-link=Germund Dahlquist | title=A special stability problem for linear multistep methods | doi=10.1007/BF01963532 | year=1963 | journal=BIT | issn=0006-3835 | volume=3 | pages=27–43}}.
* {{Citation | last1=Hairer | first1=Ernst | last2=Nørsett | first2=Syvert Paul | last3=Wanner | first3=Gerhard | title=Solving ordinary differential equations I: Nonstiff problems | publisher=[[Springer-Verlag]] | location=Berlin, New York | isbn=978-3-540-56670-0 | year=1993}}.
* {{Citation | last1=Hairer | first1=Ernst | last2=Wanner | first2=Gerhard | title=Solving ordinary differential equations II: Stiff and differential-algebraic problems | publisher=[[Springer-Verlag]] | location=Berlin, New York | edition=2nd | isbn=978-3-540-60452-5 | year=1996}}.
* {{Citation|title=A First Course in the Numerical Analysis of Differential Equations (Second Edition)|year=2008|last1=Iserles|first1=Arieh|publisher=[[Cambridge University Press]]|isbn=978-0-521-73490-5}}.
* {{Citation | last1=Kaw | first1=Autar | last2=Kalu | first2=Egwu | title=Numerical Methods with Applications | url=http://numericalmethods.eng.usf.edu/topics/textbook_index.html | publisher=autarkaw.com | edition=1st | year = 2008}}.
* {{Citation | last1=Kutta | first1=Martin Wilhelm | author1-link=Martin Wilhelm Kutta | title=Beitrag zur näherungsweisen Integration totaler Differentialgleichungen | journal=Zeitschrift für Mathematik und Physik | volume=46 | pages=435–453 | year=1901}}.
* {{Citation | last1=Press | first1=William H. | last2=Flannery | first2=Brian P. | last3=Teukolsky | first3=Saul A. | author3-link=Saul Teukolsky | last4=Vetterling | first4=William T. | title=[[Numerical Recipes|Numerical Recipes: The Art of Scientific Computing]] | chapter=Section 17.1 Runge-Kutta Method | chapter-url=http://apps.nrbook.com/empanel/index.html#pg=907 | publisher=[[Cambridge University Press]] | edition=3rd | isbn=978-0-521-88068-8 | year=2007}}. Also, [http://apps.nrbook.com/empanel/index.html#pg=910 Section 17.2. Adaptive Stepsize Control for Runge-Kutta].
* {{Citation | last1=Stoer | first1=Josef | last2=Bulirsch | first2=Roland | title=Introduction to Numerical Analysis | publisher=[[Springer-Verlag]] | location=Berlin, New York | edition=3rd | isbn=978-0-387-95452-3 | year=2002}}.
* {{Citation | last1=Süli | first1=Endre | last2=Mayers | first2=David | title=An Introduction to Numerical Analysis | publisher=[[Cambridge University Press]] | isbn=0-521-00794-1 | year=2003}}.

* John C. Butcher: "B-Series : Algebraic Analysis of Numerical Methods", [[スプリンガー#出版|Springer]]（SSCM, volume 55), ISBN 978-3030709556 (April, 2021).

* [https://ocw.kyoto-u.ac.jp/wp-content/uploads/2000/04/1999_introduction-to-simulation_11.pdf Runge-Kutta 法についてのノート (早川尚男)] # 公式を導出する例を示している。

==関連項目==
{{Commonscat|Runge-Kutta methods}}
* [[ルンゲ＝クッタ法のリスト]]
* [[常微分方程式の数値解法]]
* [[線型多段法]]

{{DEFAULTSORT:るんけくつたほう}}
[[Category:解析学]]
[[Category:数値解析]]
[[Category:数値微分方程式]]
[[Category:数学のエポニム]]
[[Category:数学に関する記事]]