C Sharpと数値解析 - ルンゲ・クッタ法

概要

ルンゲ・クッタ法は、微分方程式 ${\dfrac {dy}{dt}}=f(t,y)$ の数値解法として、カール・ルンゲとマルティン・クッタによって開発された手法である。

この手法は、解曲線に沿って複数の点で関数を評価して、それらを適切に組み合わせることにより高精度な近似解を得ることにある。

ルンゲ・クッタ法の理論は、テイラー級数展開に基づいており、任意の次数nのルンゲ・クッタ法は、テイラー展開のn次の項まで一致するように設計されている。
これにより、刻み幅hに関する誤差を $O(h^{n+1})$ まで抑えることができる。

一般的なルンゲ・クッタ法の形式は次式となる。

$y(n+1)=y(n)+h(b_{1}k_{1}+b_{2}k_{2}+\cdots +b_{n}k_{n})$

ここで、$k1,k2,\cdots ,k_{n}$ は各段階での傾きの評価値であり、$b1,b2,\cdots ,b_{n}$ は重み係数である。
各段階での傾きk_iは、以下に示すように計算される。
${\begin{aligned}k1&=f(t_{n},y_{n})\\k2&=f(t_{n}+c_{2}h,y_{n}+h(a_{21}k_{1}))\\k3&=f(t_{n}+c_{3}h,y_{n}+h(a_{31}k_{1}+a_{32}k_{2}))\\\cdots \end{aligned}}$

• 精度制御

  高次の方法ほど理論的な精度は高くなるが、実際の問題では丸め誤差の影響も考慮する必要がある。

• 安定性

  高次の方法は一般に安定性領域が狭くなる傾向があり、特に硬い方程式では注意が必要である。

• 計算効率

  次数が上がるほど関数評価の回数が増加するため、問題に応じて適切な次数を選択する必要がある。

• 適応的ステップサイズ制御

  より高度な実装では、局所誤差の推定に基づいてステップサイズを動的に調整する機能を組み込むことが一般的である。

N次のルンゲ・クッタ法

1次のルンゲ・クッタ法

オイラー法と同等であり、最も単純な形式である。
1点での傾き評価のみを使用するため、精度は低いものの、計算コストは最小限である。

2次のルンゲ・クッタ法

中点法やホイン法として知られる方法が含まれる。
2点での傾き評価を行い、それらを組み合わせることにより、1次よりも高い精度を得ることができる。

3次のルンゲ・クッタ法

3点での傾き評価を行い、より高精度な近似を実現する。
ただし、4次法と比較すると使用頻度は低い。

4次のルンゲ・クッタ法

最も広く使われている方法であり、精度と計算効率のバランスが優れている。
4点での傾き評価を行い、これらを適切な重みで組み合わせる。

高次のルンゲ・クッタ法

5次以上の高次の方法も存在するが、計算コストの増加に対する精度の向上が比較的小さいため、特殊な用途を除いてあまり使用されていない。

ルンゲ・クッタ法の種類

埋め込み型ルンゲ・クッタ法

異なる次数の方法を組み合わせることにより、効率的な誤差推定を可能にする手法である。

陰的ルンゲ・クッタ法

硬い方程式に対して、より優れた安定性を持つ変種である。

シンプレクティック・ルンゲ・クッタ法

ハミルトン系の数値解法として、エネルギー保存等の幾何学的性質を保持する特殊な形式である。

4次のルンゲ・クッタ法

 using System;
 using System.Collections.Generic;
 
 /// <summary>
 /// 4次のルンゲクッタ法による微分方程式ソルバ
 /// </summary>
 public class RungeKutta
 {
    private readonly Func<double, double, double> _function;
    private readonly double _stepSize;
 
    /// <summary>
    /// コンストラクタ
    /// </summary>
    /// <param name="function">微分方程式の右辺を表す関数 f(t, y)</param>
    /// <param name="stepSize">計算のステップサイズ</param>
    public RungeKutta(Func<double, double, double> function, double stepSize)
    {
       _function = function;
       _stepSize = stepSize;
    }
 
    /// <summary>
    /// 指定された区間で微分方程式を解く
    /// </summary>
    /// <param name="t0">開始時刻</param>
    /// <param name="y0">初期値</param>
    /// <param name="tEnd">終了時刻</param>
    /// <returns>計算結果の時刻と解のペアのリスト</returns>
    public List<(double Time, double Value)> Solve(double t0, double y0, double tEnd)
    {
       var result = new List<(double Time, double Value)>();
       var t      = t0;
       var y      = y0;
 
       // 初期値を結果リストに追加
       result.Add((t, y));
 
       // 終了時刻に達するまで計算を続ける
       while (t < tEnd)
       {
          // 4次のルンゲ・クッタ法による1ステップの計算
          var k1 = _function(t, y);
          var k2 = _function(t + _stepSize / 2, y + _stepSize * k1 / 2);
          var k3 = _function(t + _stepSize / 2, y + _stepSize * k2 / 2);
          var k4 = _function(t + _stepSize, y + _stepSize * k3);
 
          // 次のステップの値を計算
          y += _stepSize * (k1 + 2 * k2 + 2 * k3 + k4) / 6;
          t += _stepSize;
 
          // 結果をリストに追加
          result.Add((t, y));
       }
 
       return result;
    }
 }

 // 例 減衰を表す微分方程式 : y' = -2y
 Func<double, double, double> f = (t, y) => -2 * y;
 
 // 4次のルンゲ・クッタ法による解法
 var rk         = new RungeKutta(f, stepSize: 0.1);
 var rkSolution = rk.Solve(t0: 0, y0: 1, tEnd: 5);