数値積分

数値積分（Numerical Integration）は、定積分

の値を数値的に近似計算する手法です。解析的に不定積分を求めることが困難な関数や、離散的なデータ点としてしか与えられていない関数（実験データなど）の積分値を求めるために用いられます。物理学では、運動方程式の積分によるエネルギー計算、電磁場のポテンシャル計算、分配関数の計算など、極めて広範な応用を持ちます。

本節では、最も基本的かつ実用的な手法である台形則 (Trapezoidal Rule) と シンプソン則 (Simpson’s Rule) について解説します。これらの手法は「ニュートン・コーツ (Newton-Cotes) の公式」と呼ばれる手法群の一種で、積分区間を等間隔に分割して関数値を評価します。

複合台形則 (Trapezoidal Rule)

原理

台形則は、関数を各小区間で一次関数（直線）で近似し、その下の面積を台形の面積として計算する手法です。

Trapezoidal Rule

積分区間を等分し、刻み幅をとします。分点を ( ) とすると、各区間の面積は以下のように近似されます。

これらを全区間にわたって足し合わせることで、全体の積分値の近似式（複合台形則）が得られます。

シンプソン則を適用するには、分割数が偶数である必要があります。近似式は以下のようになります。

係数がのパターンで並ぶのが特徴です。

誤差

シンプソン則の誤差は、刻み幅の4乗に比例します（）。つまり、分割数を2倍にすると、誤差は約 1/16 に激減します。これが、シンプソン則が広く使われる理由です。

補足: それぞれの公式の導出

台形則の導出 区間において、関数をを通る1次多項式（ラグランジュ補間多項式）で近似します。

これを区間で積分すると、台形則の公式が得られます。

シンプソン則の導出 区間において、3点を通る2次多項式で近似します。

ここではラグランジュ基底多項式です。これを区間（幅）で積分すると、シンプソン則の公式が得られます。

Rustによる実装

それでは、これらの手法をRustで実装し、精度を比較してみましょう。例として、を区間で積分します。解析解はです。

use std::f64::consts::PI;

/// 台形則
fn trapezoidal_rule<F>(f: F, a: f64, b: f64, n: usize) -> f64
where
    F: Fn(f64) -> f64,
{
    let h = (b - a) / n as f64;
    let sum: f64 = (1..n).map(|i| f(a + i as f64 * h)).sum();

    // 両端の点は重み 1/2
    h * (0.5 * f(a) + sum + 0.5 * f(b))
}

/// シンプソン則 (nは偶数でなければならない)
fn simpsons_rule<F>(f: F, a: f64, b: f64, n: usize) -> f64
where
    F: Fn(f64) -> f64,
{
    assert!(n % 2 == 0, "Simpson's rule requires an even number of intervals.");

    let h = (b - a) / n as f64;

    let mut sum_odd = 0.0;
    let mut sum_even = 0.0;

    for i in 1..n {
        let x = a + i as f64 * h;
        if i % 2 == 0 {
            sum_even += f(x);
        } else {
            sum_odd += f(x);
        }
    }

    h / 3.0 * (f(a) + 4.0 * sum_odd + 2.0 * sum_even + f(b))
}

fn main() {
    let f = |x: f64| x.sin();
    let a = 0.0;
    let b = PI;
    let exact = 2.0;

    println!("{:<10} {:<20} {:<20} {:<20} {:<20}",
        "N", "Trapezoidal", "Error (Trap)", "Simpson", "Error (Simp)");
    println!("{}", "-".repeat(95));

    let n_values = [10, 20, 40, 80, 160];

    for &n in &n_values {
        let trap = trapezoidal_rule(f, a, b, n);
        let simp = simpsons_rule(f, a, b, n);

        println!(
            "{:<10} {:.15}    {:.5e}           {:.15}    {:.5e}",
            n, trap, (trap - exact).abs(), simp, (simp - exact).abs());
    }
}

実行結果：

N          Trapezoidal          Error (Trap)         Simpson              Error (Simp)
-----------------------------------------------------------------------------------------------
10         1.983523537509455    1.64765e-2           2.000109517315004    1.09517e-4
20         1.995885972708715    4.11403e-3           2.000006784441801    6.78444e-6
40         1.998971810497066    1.02819e-3           2.000000423093183    4.23093e-7
80         1.999742972445836    2.57028e-4           2.000000026428759    2.64288e-8
160        1.999935744350136    6.42556e-5           2.000000001651570    1.65157e-9

結果の考察

台形則: が2倍（10 -> 20）になると、誤差は約 1/4 ( ) になっています。これは理論通りの収束です。
シンプソン則: が2倍（10 -> 20）になると、誤差は約 1/16 ( ) になっています。これも理論通りの収束を示しています。
同じ分割数で比較すると、シンプソン則の方が圧倒的に高精度であることがわかります。

どちらを使うべきか？

関数が滑らかな場合: シンプソン則が推奨されます。少ない計算コストで高い精度が得られます。
実験データなどノイズを含む場合: 高次の近似がかえってノイズを増幅させる可能性があるため、台形則の方が安定する場合があり、一般的によく用いられます。
周期関数の全周期積分: 周期関数の場合、台形則が例外的に極めて高い精度（指数関数的収束）を示すことが知られています。

まとめ

台形則は一次近似であり、誤差は。シンプルでロバスト。
シンプソン則は二次近似であり、誤差は。滑らかな関数に対して非常に高精度。
シンプソン則を使うには、分割数が偶数である必要がある。

次節では、より高度な積分手法である「ガウス求積法」について学びます。

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