タケウチCarlの地平線

都内在住31歳の独身サラリーマンが、日々木工や読書、散歩などを楽しみつつ、いつか脱サラして小屋暮らしや旅暮らしをすることを夢見るブログ

ルンゲ=クッタ法2

 古典的な4次のルンゲ=クッタ法によって計算される物体の位置y_nは、厳密解と比べてどの程度近いのでしょうか。ルンゲ=クッタ法1 - タケウチCarlの地平線で述べたように、両者をテイラー展開して比較すると(これはhのベキ級数ですが)、4次の項(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が一変数の場合(方向が一つしかない場合)でした。多変数の場合に拡張しておくと、多次元空間での運動がシミュレーションできるようになりますし、ニュートン運動方程式のように、2階以上の微分方程式数値計算できるようになります。
 実は、見た目は一変数の場合と全く変わらず、ルンゲ=クッタ法の計算公式は、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の係数まで一致することが確かめられます。
 確認は一変数の場合と全く同様にできます。というのは、一変数のときにやった微分の計算(テイラー展開)は、多変数になるとベクトルと行列の計算になりますが、式の形は一変数の場合と全く同じだからです。