ルンゲ＝クッタ法２ - カールの曲がった地平線

　古典的な４次のルンゲ＝クッタ法によって計算される物体の位置 $y_n$ は、厳密解と比べてどの程度近いのでしょうか。ルンゲ＝クッタ法１ - タケウチCarlの地平線で述べたように、両者をテイラー展開して比較すると（これは $h$ のベキ級数ですが）、４次の項（ $h^4$ の項）まで一致することが確かめられます。ここではその確認を行います。また、利便性を図るために多変数の場合へと拡張します。

厳密解における位置 $y_n$

　記号を簡単にするために、 $t_n=t, y_n=y(t)$ 、 $t_{n+1}=t+h, y_{n+1}=y(t+h)$ とおいて、テイラー展開をすると、
$y_{n+1}-y_n=y(t+h)-y(t)=hy'+\frac{h^2}{2}y^{(2)}+\frac{h^3}{6}y^{(3)}+\frac{h^4}{24}y^{(4)}+\cdots$
$=hf+\frac{h^2}{2}(\frac{df}{dt})+\frac{h^3}{6}(\frac{d^2f}{dt^2})+\frac{h^4}{24}(\frac{d^3f}{dt^3})+\cdots$
となります。微小な時間間隔 $h$ のベキ級数で表せました。しかし、係数は $f$ の高次の全微分係数となっています。後々ルンゲ＝クッタ法の $y_n$ と比較しやすいように、 $\frac{df}{dt}$ 、 $\frac{d^2f}{dt^2}$ 、 $\frac{d^3f}{dt^3}$ を、 $f$ の変微分係数を使って書き直しておきましょう。まず、全微分の公式から、
$\frac{df}{dt}=f_t+f_y\frac{dy}{dt}=f_t+f_yf$
です。これを $t$ で全微分して、
$\frac{d^2f}{dt^2}=\frac{d}{dt}(f_t+f_yf)=(f_{tt}+ff_{ty})+(f_t+ff_y)f_y+f(f_{ty}+ff_{yy})$
$=f_{tt}+2ff_{ty}+f^2f_{yy}+f_tf_y+f{f_y}^2$
です。さらにこれを $t$ で全微分して、
$\frac{d^3f}{dt^3}=f_{ttt}+3ff_{tty}+3f^2f_{tyy}+f^3f_{yyy}+f_{tt}f_y$
$+5ff_{ty}f_y+4f^2f_{yy}f_y+3f_tf_{ty}+3ff_tf_{yy}+f_t{f_y}^3+f{f_y}^3$
となります。単純に全微分を適用していけばいいだけなのですが、計算は結構ハードです。

ルンゲ＝クッタ法から計算される $y_n$

　記号の簡略化のために、 $f(t_n,y_n)$ を $f$ と書くことにします。すると、 $k_1=f$ です。
　次に、 $k_2$ をテイラー展開して、
$k_2=f(t_n+\frac{h}{2},y_n+\frac{h}{2}k_1)=f+\frac{h}{2}(f_t+k_1f_y)+\frac{h^2}{8}(f_{tt}+2k_1f_{ty}+k_1^2f_{yy})$
$+\frac{h^3}{48}(f_{ttt}+3k_1f_{tty}+3k_1^2f{tyy}+k_1^3f_{yyy})+O(h^4)$
$=f+\frac{h}{2}(f_t+ff_y)+\frac{h^2}{8}(f_{tt}+2ff_{ty}+f^2f_{yy})$
$+\frac{h^3}{48}(f_{ttt}+3ff_{tty}+3f^2f{tyy}+f^3f_{yyy})+O(h^4)$
となります。
　また、 $k_3$ についても同様に、テイラー展開を $h^3$ の項まで行うと、
$k_3=f+\frac{h}{2}(f_t+k_2f_y)+\frac{h^2}{8}(f_{tt}+2ff_{ty}+f^2f_{yy}+2f_tf_y+2ff_y^2)$
$+\frac{h^3}{48}(f_{ttt}+3ff_{tty}+3f^2f{tyy}+f^3f_{yyy}$
$+3f_yf_{tt}+12ff_yf_{ty}+9f^2f_yf_{yy}+6f_tf_{ty}+6ff_tf_{yy})+O(h^4)$
となります。
　 $k_4$ についても同様にテイラー展開して、
$k_4=f+h(f_t+k_2f_y)+\frac{h^2}{2}(f_{tt}+2ff_{ty}+f^2f_{yy}+f_tf_y+ff_y^2)$
$+\frac{h^3}{24}(4f_{ttt}+12ff_{tty}+12f^2f{tyy}+4f^3f_{yyy}+3f_yf_{tt}$
$+18ff_yf_{ty}+15f^2f_yf_{yy}+6f_tf_y^2+6ff_y^3+12f_tf_{ty}+12ff_tf_{yy})+O(h^4)$
となります。
　以上によって $\frac{h}{6}(k_1+2k_2+2k_3+k_4)$ を計算すると、これが厳密解のテイラー展開と $h^4$ の係数まで一致することが確認できます。

多変数の場合

　ここまでは、位置 $y$ が一変数の場合（方向が一つしかない場合）でした。多変数の場合に拡張しておくと、多次元空間での運動がシミュレーションできるようになりますし、ニュートンの運動方程式のように、２階以上の微分方程式も数値計算できるようになります。
　実は、見た目は一変数の場合と全く変わらず、ルンゲ＝クッタ法の計算公式は、 $y$ を多次元のベクトル $\mathbf{y}$ で置き換えたものになります。すなわち、微分方程式 $\frac{d\mathbf{y}}{dt}=\mathbf{f}(t,\mathbf{y})$ 、 $\mathbf{y}(t_0)=\mathbf{y}_0$ において、時間の刻み幅を $h$ として、ある時刻の位置 $\mathbf{y}_n$ から、次の時刻の $\mathbf{y}_{n+1}$ は、
$\mathbf{y}_{n+1}=\mathbf{y}_n+\frac{h}{6}(\mathbf{k}_1+2\mathbf{k}_2+2\mathbf{k}_3+\mathbf{k}_4)$
によって計算されます。
　ここで、 $\mathbf{y}$ 、 $\mathbf{f}$ 、 $\mathbf{k_{*}}$ は皆ベクトルです。例えば、 $\mathbf{y}=(y_1,y_2,\cdots,y_n)$ です。
　そして、やはり、厳密に計算した位置 $\mathbf{y}$ と、ルンゲ＝クッタ法から計算される $\mathbf{y}$ を比較すると、両者のテイラー展開の $h^4$ の係数まで一致することが確かめられます。
　確認は一変数の場合と全く同様にできます。というのは、一変数のときにやった微分の計算（テイラー展開）は、多変数になるとベクトルと行列の計算になりますが、式の形は一変数の場合と全く同じだからです。