オイラー法

Note

本節のポイント

常微分方程式の最も基本的な解法であるオイラー法の原理を理解する。

1次精度手法の誤差の蓄積（大域誤差）について学ぶ。

2階微分方程式を連立1階微分方程式に帰着させる手法を習得する。

オイラー法 (Euler’s Method) は、常微分方程式の数値解法の中で最も基本的かつ単純な手法です。実用的な精度は低いものの、微分方程式を数値的に解く（時間発展させる）ための基礎的な概念を理解するために不可欠です。

原理

1階の常微分方程式の初期値問題を考えます。

微分の定義により、時間刻み幅が十分に小さければ、導関数は前方差分で近似できます。

これをについて解くと、現在の状態から次の時刻の状態を求める更新式（漸化式）が得られます。

Euler Method

幾何学的には、「現在の点における接線の傾きに沿って、時間だけ直線的に進む」ことを意味します。このため、関数の曲率が大きい場合、ステップを進めるごとに真の解から大きく外れていくことになります。

誤差の性質

オイラー法の誤差には、1ステップごとに発生する 局所誤差 と、計算終了までに蓄積される 大域誤差 があります。

局所誤差: テイラー展開をの項で打ち切っているため、となります。
大域誤差: 計算ステップ数はに比例するため、誤差の蓄積によりとなります。

このように、全体の誤差がの1乗に比例する手法を 1次精度(First-order accuracy) と呼びます。刻み幅を半分にしても、誤差は半分にしかなりません。実用的な科学技術計算では、より高次（通常は4次以上）の手法が求められます。

Rustによる実装

まずは、1変数の微分方程式を解くプログラムを実装してみましょう。例として、以下の放射性崩壊などを記述する方程式を考えます。

解析解はです。

/// オイラー法による1ステップの更新
fn euler_step<F>(x: f64, t: f64, h: f64, f: F) -> f64
where
    F: Fn(f64, f64) -> f64,
{
    x + h * f(t, x)
}

fn euler_integrate<F>(x0: f64, t0: f64, t_max: f64, h: f64, f: F) -> f64
where
    F: Fn(f64, f64) -> f64 + Copy,
{
    let mut t = t0;
    let mut x = x0;

    while t < t_max {
        let step_h = if t + h > t_max { t_max - t } else { h };
        x = euler_step(x, t, step_h, f);
        t += step_h;
    }

    x
}

fn main() {
    let f = |_t: f64, x: f64| -x;
    let x0 = 1.0;
    let t_max = 2.0;

    println!(
        "{:<5} {:<15} {:<15} {:<15}",
        "h", "Numerical", "Exact", "Error"
    );
    println!("{}", "-".repeat(55));

    for &h in &[0.4, 0.2, 0.1, 0.05] {
        let x = euler_integrate(x0, 0.0, t_max, h, f);
        let exact = (-t_max).exp();
        println!(
            "{:<5.2} {:<15.8} {:<15.8} {:<15.2e}",
            h,
            x,
            exact,
            (x - exact).abs()
        );
    }
}

実行結果：

h     Numerical       Exact           Error
-------------------------------------------------------
0.40  0.07776000      0.13533528      5.76e-2
0.20  0.08589935      0.13533528      4.94e-2
0.10  0.12157665      0.13533528      1.38e-2
0.05  0.12851216      0.13533528      6.82e-3

を半分にするごとに、誤差もおよそ半分になっていることが確認できます。これが1次精度の特徴です。

2階微分方程式への適用（連立化）

ニュートンの運動方程式のような2階微分方程式は、速度を導入することで、2つの連立1階微分方程式に書き換えて解くことができます。

これをオイラー法で書くと：

一般に、多変数あるいは高階の微分方程式は、状態をベクトルとして扱うことで、単一のベクトル微分方程式として統一的に記述できます。これについては次節のルンゲ＝クッタ法で詳しく扱います。

まとめ

オイラー法は、微分を前方差分で近似する最も単純な解法。
1次精度 ( )であり、実用的な精度を得るには非常に小さな刻み幅が必要。
2階以上の微分方程式は、変数を増やすことで1階の連立方程式に変換できる。

次節では、オイラー法の精度を大幅に改善したルンゲ＝クッタ法を学びます。

Keyboard shortcuts

Rust計算物理学チュートリアル

オイラー法

原理

誤差の性質

Rustによる実装

2階微分方程式への適用（連立化）

まとめ