ルンゲ＝クッタ法

Note

本節のポイント

• 科学技術計算の標準手法である4次のルンゲ＝クッタ法(RK4)のアルゴリズムを理解する。
• ndarrayを用いて、連立微分方程式をベクトル形式で効率的に解く実装方法を学ぶ。
• 4次精度による圧倒的な誤差減少の効果を体感する。

ルンゲ＝クッタ法(Runge-Kutta Methods) は、オイラー法よりも高精度に微分方程式を解くための一連の手法です。その中でも特に、4次のルンゲ＝クッタ法(RK4) が精度と計算コストのバランスが非常に良く、実用上の標準として広く使われています。

原理：4次のルンゲ＝クッタ法(RK4)

オイラー法が区間の始点での傾きだけを使ったのに対し、RK4では区間内の4点の傾きを計算し、それらを加重平均して次の値を決定します。 ベクトル形式の微分方程式 に対して、時間刻み での更新式は以下の通りです。

幾何学的意味

1. : 始点での傾き（オイラー法と同じ）。
2. : を使って中点へ進み、そこでの傾き。
3. : を使って中点へ進み（修正）、そこでの傾き。
4. : を使って終点へ進み、そこでの傾き。

最後にこれらを の重みで平均します。これにより、テイラー展開の4次の項までが一致し、大域誤差は （4次精度） となります。これは、刻み幅 を半分にすると、誤差が に激減することを意味します。

Rustによる実装（ndarrayの活用）

連立微分方程式（多変数）を扱う場合、ndarray入門で学んだndarrayを使うと、数学的なベクトル演算をそのままコードに落とし込めるため、非常に見通しが良くなります。 例として、単振動（調和振動子）の方程式を解いてみましょう。

これを連立化すると、状態ベクトル に対して以下のようになります。

use ndarray::{Array1, arr1};

/// 4次のルンゲ＝クッタ法による1ステップの更新
fn rk4_step<F>(state: &Array1<f64>, t: f64, h: f64, f: &F) -> Array1<f64>
where
    F: Fn(f64, &Array1<f64>) -> Array1<f64>,
{
    let k1 = f(t, state);
    let k2 = f(t + h * 0.5, &(state + &k1 * (h * 0.5)));
    let k3 = f(t + h * 0.5, &(state + &k2 * (h * 0.5)));
    let k4 = f(t + h, &(state + &k3 * h));

    state + (&k1 + &k2 * 2.0 + &k3 * 2.0 + &k4) * (h / 6.0)
}

fn integrate_rk4<F>(
    state0: &Array1<f64>,
    t0: f64,
    t_max: f64,
    h: f64,
    f: &F,
) -> Array1<f64>
where
    F: Fn(f64, &Array1<f64>) -> Array1<f64>,
{
    let mut t = t0;
    let mut state = state0.clone();

    while t < t_max {
        // ステップ幅が余る場合の調整
        let step_h = if t + h > t_max { t_max - t } else { h };
        state = rk4_step(&state, t, step_h, f);
        t += step_h;
    }

    state
}

fn harmonic_oscillator(_t: f64, state: &Array1<f64>) -> Array1<f64> {
    let x = state[0];
    let v = state[1];
    arr1(&[v, -x])
}

fn main() {
    let x0 = arr1(&[1.0, 0.0]); // 初期条件: x=1, v=0
    let t_max = 2.0 * std::f64::consts::PI; // 1周期

    println!("{:<5} {:<15} {:<15}", "h", "Final x", "Error");
    println!("{}", "-".repeat(40));

    for &h in &[0.5, 0.25, 0.125, 0.0625] {
        let state = integrate_rk4(&x0, 0.0, t_max, h, &harmonic_oscillator);
        let exact = 1.0; // cos(2pi) = 1
        println!(
            "{:<5.3} {:<15.10} {:<15.2e}",
            h,
            state[0],
            (state[0] - exact).abs()
        );
    }
}

実行結果：

h     Final x         Error
----------------------------------------
0.500 0.9987316280    1.27e-3
0.250 0.9999579266    4.21e-5
0.125 0.9999986780    1.32e-6
0.062 0.9999999586    4.14e-8

を半分（1/2）にするごとに、誤差がおよそ 1/16 ずつ減少していることがわかります。オイラー法と比較して、少ない計算回数で驚異的な精度が得られることがRK4の強みです。

注意点と限界

RK4は非常に優秀ですが、万能ではありません。

1. 硬い方程式 (Stiff Equation): 時間スケールが大きく異なる現象（例：非常に速い振動とゆっくりの運動）が混在する系では、RK4は不安定になりやすく、ステップ幅 を極端に小さくする必要があります。このような場合は、陰的解法 (Implicit Method) などが用いられます。
2. エネルギー保存: RK4はエネルギーを厳密には保存しません。長時間計算すると、軌道が徐々にズレていくことがあります。

次節では、計算精度を保ちながら効率よく計算するために、ステップ幅 を自動調整する方法を学びます。