こんにちは、Rafkaです。
前回 → 数値シミュレーション －ホイン法－
今回は４次精度のルンゲ＝クッタ（Runge-Kutta）法のアルゴリズムの紹介と実装を行います。
アルゴリズム
ルンゲ＝クッタ法では以下の更新式で値を更新し、数値解を求めていきます。
$$\begin{align}\boldsymbol{k}_1 &= \boldsymbol{f}\bigl(t, \boldsymbol{x}(t)\bigr)\\ \boldsymbol{k}_2 &= \boldsymbol{f}\Biggl(t + \frac{\Delta t}{2}, \boldsymbol{x}(t) + \frac{\Delta t}{2}\boldsymbol{k}_1\Biggr)\\ \boldsymbol{k}_3 &= \boldsymbol{f}\Biggl(t + \frac{\Delta t}{2}, \boldsymbol{x}(t) + \frac{\Delta t}{2}\boldsymbol{k}_2\Biggr)\\ \boldsymbol{k}_4 &= \boldsymbol{f}\bigl(t + \Delta t, \boldsymbol{x}(t) + \Delta t\boldsymbol{k}_3\bigr)\\ \boldsymbol{x}(t + \Delta t) &= \boldsymbol{x}(t) + \frac{\Delta t}{6}(\boldsymbol{k}_1 + 2\boldsymbol{k}_2 + 2\boldsymbol{k}_3 + \boldsymbol{k}_4)\end{align}$$
導出はホイン法の時と同様に、次の時間ステップである \(\boldsymbol{x}(t + \Delta t)\) を４次の項まで展開したものと、\(\boldsymbol{k}_1,\;\boldsymbol{k}_2,\;\boldsymbol{k}_3,\;\boldsymbol{k}_4\) に任意定数を掛けて総和を取ったものとで同じになるように定数を決定すればできると思います。かなり大変な計算になりそうですが。そもそも \(\boldsymbol{k}_1,\;\boldsymbol{k}_2,\;\boldsymbol{k}_3,\;\boldsymbol{k}_4\) はなぜこの形になるのかという話もあると思います。以下のサイトに、４次に限らない \(n\) 次のルンゲ＝クッタ法の定義や導出が書かれていたので、気になる方は見に行ってみてください。
ルンゲ=クッタ法を導く (1) ： 準備としての高次テイラー法
実装
上で示した更新式を以下のように実装しました。
ホイン法の時と同様に、変更点は更新値を計算するメソッドの部分だけです。
まとめ
これで常微分方程式の数値解法としてよく用いられる３つの手法の紹介が終わりました。
これらの手法は、ホイン法はオイラー法よりも誤差が小さくなるように、ルンゲ＝クッタ法はホイン法よりも誤差が小さくなるように、それぞれ開発されたので、次回の記事ではこれらの手法の誤差の比較を行う予定です。
ついでに、これら３つの手法は異なる点が更新の部分だけなので、ソースコードのリファクタリングも同時に行うかもしれません。