光、爆発！
光、爆発！

ルンゲクッタ法

2015-02-01

ルンゲ＝クッタ法とは

　物体が物理的法則にしたがってどのように動くのかをシミュレーションできたら面白いし、また物理現象を理解するためにも役立ちそうです。調べると、（4次の）ルンゲ＝クッタ法というものがあり、計算が簡単に実行できるわりに精度が良いようです。これは、初期条件と時刻tでのyの変化率が与えられている、次の微分方程式
\frac{dy}{dt}=f(t,y)、y(t_0)=y_0
において、時間の刻み幅をhとして、ある時刻のy_nから、次の時刻のy_{n+1}を、次々と計算していく方法です。具体的には
y_{n+1}=y_n+\frac{h}{6}(k_1+2k_2+2k_3+k_4)
と計算します。ここでk_1、k_2、k_3、k_4は、以下のように定義されます。
\begin{cases}k_1=f(t_n,y_n)\\k_2=f(t_n+\frac{h}{2},y_n+\frac{h}{2}k_1)\\k_3=f(t_n+\frac{h}{2},y_n+\frac{h}{2}k_2)\\k_4=f(t_n+h,y_n+k_3)\end{cases}

計算式の意味

上の計算式を見ると、y_{n+1}は、y_nから、変化率\frac{k_1+2k_2+2k_3+k_4}{6}で時間hだけ進んだ状態として計算されます。
定義から、k_1は、y_nでの変化率です。k_2は、変化率をk_1とした場合に、時刻\frac{h}{2}後の地点における変化率です。また、k_3は、変化率をk_2とした場合に、時刻\frac{h}{2}後の地点における変化率です。そしてk_4は、変化率をk_3とした場合に、時刻h後の地点における変化率です。
これら4点における変化率k_1、k_2、k_3、k_4の重みつき平均をとることにより、正しい変化率を\frac{k_1+2k_2+2k_3+k_4}{6}と近似計算していることになります。

証明

ルンゲ＝クッタ法を導出することは難しいらしいのですが、一度導出されたものを、正しかどうか確認することはできますので、ここでは確認作業を行います。ただ、確認するだけでもかなり多くの計算をすることになります。

A.厳密に計算されるy_n

簡単のため、t_n=t、t_{n+1}=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のベキ級数で表せました。
ここで、\frac{df}{dt}、\frac{d^2f}{dt^2}、\frac{d^3f}{dt^3}を順番に計算すると、
\begin{align}\frac{df}{dt}&=f_t+f_y\frac{dy}{dt}\\&=f_t+f_yf\end{align}
\begin{align}\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\end{align}
となり、これをさらに微分して、
\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
となります。

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

表記を簡潔にするために、f(t_n,y_n)をfと書くことにすると、
k_1=f
です。次に、k_2をテイラー展開すると、
\begin{align}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)\end{align}
となります。また、k_3についてもテイラー展開をh^3の項まで行うと、
\begin{align}k_3&=f(t_n+\frac{h}{2},y_n+\frac{h}{2}k_2)\\&=f+\frac{h}{2}(f_t+k_2f_y)+\frac{h^2}{8}(f_{tt}+2k_2f_{ty}+k_2^2f_{yy})+\frac{h^3}{48}(f_{ttt}+3k_2f_{tty}+3k_2^2f{tyy}+k_2^3f_{yyy})+O(h^4)\\&=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)\end{align}
となります。最後のイコールでは、上で計算したk_2を代入しました。k_4についてもk_3と同様に、\begin{eqnarray}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)\end{eqnarray}
となります。これらにより\frac{k_1+2k_2+2k_3+k_4}{6}を計算すると、Aで計算したものとh^4の係数まで一致することが確認できます。

多変数の場合

上記は一変数の場合、すなわち時間tと変数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)です。
証明も、一変数の場合と全く同様にできます。というのは、例えば、\frac{d^3\mathbf{y}}{dt^3}を例にとると、
\begin{align}\frac{d^3\mathbf{y}}{dt^3}&=\frac{d}{dt}(\frac{\partial \mathbf{f}}{\partial t}+\mathbf{f}\cdot\frac{\partial \mathbf{f}}{\partial \mathbf{y}})\\&=\frac{\partial^2 \mathbf{f}}{\partial t^2}+2\mathbf{f}\cdot\frac{\partial^2 \mathbf{f}}{\partial \mathbf{y}\partial t}+\frac{\mathbf{f}}{\partial t}\cdot\frac{\partial \mathbf{f}}{\partial \mathbf{y}}+\mathbf{f}\cdot\frac{\mathbf{f}}{\partial \mathbf{y}}\cdot\frac{\partial \mathbf{f}}{\partial \mathbf{y}}+\mathbf{f}\cdot \mathbf{f}\frac{\partial^2 \mathbf{f}}{\partial \mathbf{y}^2}\end{align}
となり、ベクトルと行列の計算ではありますが、式の形は一変数の場合の\frac{d^3y}{dt^3}と全く同じです。なので、一変数のときと全く同様の計算が成り立ち、厳密に計算したものと、ルンゲ＝クッタ法から計算されるy_nが、h^4の係数まで一致することが分かります。
「ルンゲクッタ法」をウィキ内検索
最終更新：2015年02月01日 10:04