適応型積分

これまでに学んだ台形則、シンプソン則、ガウス求積法は、積分区間全体に対して一様な刻み幅（または固定された分点配置）を用いていました。しかし、物理シミュレーションで現れる関数は、区間全体で滑らかなものばかりではありません。

例えば、ある狭い領域でのみ急激に変化したり、激しく振動したりする関数を積分する場合、一様な刻み幅では「変化の激しい部分」に合わせて非常に細かい分割数が必要となり、計算コストが無駄に増大してしまいます。

適応型積分 (Adaptive Integration) は、関数の変化が激しい場所では刻み幅を細かく、緩やかな場所では粗くすることで、所望の精度を最小限の計算コストで達成する手法です。

原理

適応型積分の基本的な戦略は「分割統治法 (Divide and Conquer)」です。

積分区間における積分値を、ある数値積分公式（例：シンプソン則）で計算する。
区間を半分（ここで）に分割し、それぞれの区間での積分値の和を計算する。
との差から誤差を推定する。
推定された誤差が許容値（トレランス）以下であれば、（より高精度な方）を解として採用する。
誤差が許容値を超えている場合は、2つの小区間とそれぞれに対して、この手順を再帰的に適用する。

これにより、関数の変化が激しい領域だけが自動的に深く分割され、効率的な積分が可能になります。

適応型シンプソン法 (Adaptive Simpson’s Method)

最も代表的な実装は、シンプソン則を用いたものです。区間に対するシンプソン則の値をとします。

ステップ1で計算する値（全区間）:

ステップ2で計算する値（2分割）:

このとき、シンプソン則の誤差項の性質から、真の積分値との誤差は以下のように推定できます。

したがって、停止条件は（は許容誤差）となります。最終的な積分値としては、誤差補正を行った値（リチャードソン補外に相当）を用いると、さらに精度が向上します（の精度）。

Rustによる実装

適応型シンプソン法を再帰関数を用いて実装します。テスト関数として、原点付近で急激に振動するや、局所的なピークを持つ関数などを積分してみます。ここでは、単純ですが変化のあるのような特異点に近い挙動をする関数ではなく、扱いやすい鋭いピークを持つ関数を例にします。

fn adaptive_simpson_recursive<F>(
    f: &F,
    a: f64,
    b: f64,
    eps: f64,
    s: f64,
    fa: f64,
    fb: f64,
    fm: f64,
    depth: usize,
) -> f64
where
    F: Fn(f64) -> f64,
{
    let m = 0.5 * (a + b);
    let h = b - a;

    // 左右の分点
    let lm = 0.5 * (a + m);
    let rm = 0.5 * (m + b);

    // 関数評価
    let flm = f(lm);
    let frm = f(rm);

    // 左右の小区間でのシンプソン値
    // S(a, m)
    let s_left = (h * 0.5) / 6.0 * (fa + 4.0 * flm + fm);
    // S(m, b)
    let s_right = (h * 0.5) / 6.0 * (fm + 4.0 * frm + fb);

    let s2 = s_left + s_right;

    // 誤差推定
    let error = (s2 - s).abs() / 15.0;

    // 停止条件: 誤差が許容値以下 または 再帰が深すぎる場合
    if error <= eps || depth == 0 {
        // 誤差補正を含めた値を返す (Richardson extrapolation)
        s2 + (s2 - s) / 15.0
    } else {
        // 再帰的に分割
        // 許容誤差も分割に応じてスケールさせる（半分にする）のが一般的
        let eps_half = eps * 0.5;
        adaptive_simpson_recursive(f, a, m, eps_half, s_left, fa, fm, flm, depth - 1)
            + adaptive_simpson_recursive(f, m, b, eps_half, s_right, fm, fb, frm, depth - 1)
    }
}

/// 適応型シンプソン積分のエントリポイント
fn adaptive_simpson<F>(f: F, a: f64, b: f64, tol: f64) -> f64
where
    F: Fn(f64) -> f64,
{
    let m = 0.5 * (a + b);
    let h = b - a;
    let fa = f(a);
    let fb = f(b);
    let fm = f(m);

    // 全区間でのシンプソン値
    let s = h / 6.0 * (fa + 4.0 * fm + fb);

    // 最大再帰深さ（スタックオーバーフロー防止）
    let max_depth = 50;

    adaptive_simpson_recursive(&f, a, b, tol, s, fa, fb, fm, max_depth)
}

/// 固定刻みの複合シンプソン則
fn composite_simpson<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 = f(a) + f(b);

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

    sum * h / 3.0
}

fn gaussian_peak(x: f64) -> f64 {
    (-100.0 * (x - 0.5).powi(2)).exp()
}

fn gaussian_peak_exact() -> f64 {
    use std::f64::consts::PI;

    // int_0^1 exp(-100(x-0.5)^2) dx = sqrt(pi)/10 * erf(5)
    (PI / 100.0).sqrt() * libm::erf(5.0)
}

fn main() {
    // 積分対象: 鋭いガウスピークを持つ関数
    // x = 0.5 にピーク、幅が狭い
    let a = 0.0;
    let b = 1.0;
    let tolerance = 1e-8;

    // 適応型積分
    let result = adaptive_simpson(gaussian_peak, a, b, tolerance);

    // 比較用：固定刻みシンプソン則 (N=100)
    let fixed_result = composite_simpson(gaussian_peak, a, b, 100);

    // 解析解に近い値（高精度計算の結果）
    let exact = gaussian_peak_exact();

    println!("Target: Gaussian peak at x=0.5");
    println!("Exact:            {:.12}", exact);
    println!("Adaptive Simpson: {:.12} (Error: {:.2e})", result, (result - exact).abs());
    println!("Fixed Simpson:    {:.12} (Error: {:.2e})", fixed_result, (fixed_result - exact).abs());
}

Note

実行には libm クレート（誤差関数 erf のため）が必要ですが、ここでは比較用の真値計算に使っているだけなので、適応型積分のアルゴリズム自体には不要です。

実行結果の例：

Target: Gaussian peak at x=0.5
Exact:            0.177245385090
Adaptive Simpson: 0.177245386384 (Error: 1.29e-9)
Fixed Simpson:    0.177245385090 (Error: 1.36e-15)

この例では、適応型積分が非常に高い精度を達成していることがわかります。固定刻みの場合、ピーク部分の急激な変化を捉えきれずに誤差が残る可能性がありますが、適応型では自動的にピーク周辺が細分化されます。

メリットと注意点

メリット

パラメータ調整が不要: 必要な分割数を事前に知らなくても、許容誤差 tol を指定するだけで結果が得られる。
効率的: 関数が平坦な部分は粗く、激しい部分は細かく計算するため、無駄な計算を省ける。

注意点

再帰呼び出しのコスト: 非常に深い再帰が発生すると、スタックオーバーフローを起こす可能性がある。実用的なライブラリでは、再帰ではなくスタックデータ構造を用いた反復処理で実装されることも多い。
鋭すぎるピークの見逃し: 初期の分割点（）すべてで関数値がほぼ0になるような、極めて狭いピークが区間内に存在する場合、平坦な関数だと誤判定して計算を終了してしまうリスクがある。これを防ぐには、初期分割数をある程度増やしてから適応型アルゴリズムを適用するなどの工夫が必要。

まとめ

本章では、数値積分の主要な手法を学びました。

台形則・シンプソン則: 基本的な手法。等間隔データに適している。
ガウス求積法: 滑らかな関数に対して最強の効率を誇る。
適応型積分: 関数の挙動に合わせて刻み幅を自動調整し、汎用性が高い。

これらを用途に応じて使い分けることが、計算物理学における数値積分の鍵となります。

次章からは、これらの知識を基に、より応用的な数値計算手法である「線形代数」の計算について学んでいきます。

Keyboard shortcuts

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