数値シミュレーション －ルンゲ＝クッタ法－
数値シミュレーション解説 第5回
Published: 8/17/2018
Updated: 12/23/2025
約6分
更新履歴
- 2025/12/23: プログラムを C# 14.0 に修正
こんにちは、Rafkaです。
今回は 4 次精度のルンゲ＝クッタ（Runge-Kutta）法のアルゴリズムの紹介と実装を行います。
アルゴリズム
ルンゲ＝クッタ法では以下の更新式で値を更新し、数値解を求めていきます。
導出はホイン法の時と同様に、
次の時間ステップである を４次の項まで展開したものと、
に
任意定数を掛けて総和を取ったものとで、
同じになるように定数を決定すればできると思います。
かなり大変な計算になりそうですが。
そもそも は
なぜこの形になるのかという話もあると思います。
以下のサイトに、4 次に限らない 次のルンゲ＝クッタ法の定義や導出が書かれていたので、
気になる方は見に行ってみてください。
実装
上で示した更新式を以下のように実装しました。
void RungeKutta(ModelFunction f, Span<double> x, double t, double dt) {
var n = x.Length;
Span<double> k1 = stackalloc double[n];
Span<double> k2 = stackalloc double[n];
Span<double> k3 = stackalloc double[n];
Span<double> k4 = stackalloc double[n];
Span<double> tempX = stackalloc double[n];
var halfDt = dt * 0.5;
f(x, t, k1);
ScaleAdd(x, k1, halfDt, tempX);
f(tempX, t + halfDt, k2);
ScaleAdd(x, k2, halfDt, tempX);
f(tempX, t + halfDt, k3);
ScaleAdd(x, k3, dt, tempX);
f(tempX, t + dt, k4);
for (var i = 0; i < n; i++) {
var slope = k1[i] + k2[i] * 2.0 + k3[i] * 2.0 + k4[i];
x[i] += slope * (dt / 6.0);
}
}12 行目で を計算しています。
14 行目で、 の計算に使用するための、
を計算して
tempX に保存しています。
15 行目で を計算しています。
17 行目で、 の計算に使用するための、
を計算して tempX に保存しています。
18 行目で を計算しています。
20 行目で、 を計算して tempX に保存しています。
21 行目で を計算しています。
23 行目からのループで、最後の更新式を計算して x に足し込むことで
を計算しています。まとめ
これで常微分方程式の数値解法としてよく用いられる 3 つの手法の紹介が終わりました。
これらの手法は、ホイン法はオイラー法よりも誤差が小さくなるように、
ルンゲ＝クッタ法はホイン法よりも誤差が小さくなるように、
それぞれ開発されたので、次回の記事ではこれらの手法の誤差の比較を行う予定です。
