数値微分

解析的に導関数（微分）を求めることが困難、あるいは不可能な関数に対し、計算機を用いて近似的に微分係数を求める手法を数値微分 (Numerical Differentiation) と呼びます。

物理シミュレーションにおいては、物体の運動方程式 を解く際や、ポテンシャルエネルギー から力 を求める際など、至る所で微分計算が必要となります。

本節では、数値微分の最も基本的な手法である差分法 (Finite Difference Method) について、その原理とRustによる実装、そして数値計算特有の誤差の振る舞いについて解説します。

微分の定義と差分近似

関数 の における微分係数 の数学的な定義は以下の通りです。

コンピュータでは無限小の極限 を直接扱うことはできません。その代わり、十分に小さな有限の値 （刻み幅、ステップサイズ）を用いて、微分係数を近似します。これを差分近似と呼びます。

近似の方法によって、いくつかの種類があります。

1. 前進差分 (Forward Difference)

定義式をそのまま用いて、点 とその少し先の点 の値を使う方法です。

2. 後退差分 (Backward Difference)

点 とその少し手前の点 の値を使う方法です。

3. 中心差分 (Central Difference)

点 を中心とした前後 と の値を使う方法です。

幾何学的には、前進差分と後退差分が「2点間の直線の傾き」で接線の傾きを近似するのに対し、中心差分は「前後2点の中点における傾き」を用いていることになります。一般に、中心差分の方が前進・後退差分よりも高精度になります（理由は後述します）。

Rustによる実装

これらの差分法をRustで実装してみましょう。ここでは、関数 を例にとり、 における微分係数を計算します。解析解は なので、これと比較して精度を確認します。

fn forward_difference<F>(f: F, x: f64, h: f64) -> f64
where
    F: Fn(f64) -> f64,
{
    (f(x + h) - f(x)) / h
}

fn backward_difference<F>(f: F, x: f64, h: f64) -> f64
where
    F: Fn(f64) -> f64,
{
    (f(x) - f(x - h)) / h
}

fn central_difference<F>(f: F, x: f64, h: f64) -> f64
where
    F: Fn(f64) -> f64,
{
    (f(x + h) - f(x - h)) / (2.0 * h)
}

fn main() {
    // 微分対象の関数 f(x) = sin(x)
    let f = |x: f64| x.sin();

    // 微分を計算する点
    let x: f64 = 1.0;

    // 刻み幅 h
    let h: f64 = 0.01;

    // 解析解 (真値): cos(1.0)
    let exact = x.cos();

    // 1. 前進差分
    let diff_forward = forward_difference(f, x, h);

    // 2. 後退差分
    let diff_backward = backward_difference(f, x, h);

    // 3. 中心差分
    let diff_central = central_difference(f, x, h);

    println!("解析解:    {:.10}", exact);
    println!("前進差分:  {:.10} (誤差: {:.2e})", diff_forward, (diff_forward - exact).abs());
    println!("後退差分:  {:.10} (誤差: {:.2e})", diff_backward, (diff_backward - exact).abs());
    println!("中心差分:  {:.10} (誤差: {:.2e})", diff_central, (diff_central - exact).abs());
}

実行結果の例：

解析解:    0.5403023059
前進差分:  0.5360859810 (誤差: 4.22e-3)
後退差分:  0.5445006207 (誤差: 4.20e-3)
中心差分:  0.5402933009 (誤差: 9.00e-6)

この結果から、同じ刻み幅 でも、中心差分の方が圧倒的に誤差が小さいことがわかります。

誤差の解析：なぜ中心差分が良いのか？

テイラー展開を用いて誤差（打ち切り誤差）を評価することで、この理由を数学的に説明できます。

関数 を の周りでテイラー展開すると以下のようになります。

前進差分の誤差

前進差分の式にテイラー展開を代入して整理すると：

つまり、近似値と真値 との差（誤差）は、 の1乗に比例します（オーダー ）。

中心差分の誤差

一方、中心差分の式 に代入すると、偶数次の項が相殺されます。

こちらの誤差は、 の2乗に比例します（オーダー ）。

が小さい場合（例： ）、 となり、 に比べてはるかに小さくなります。これが、中心差分が高精度である理由です。

最適なステップサイズと丸め誤差

「 を小さくすればするほど精度が良くなる」と考えがちですが、数値計算ではそう単純ではありません。コンピュータの浮動小数点数には丸め誤差が存在するためです。

を極端に小さくすると、分子の において、非常に近い値同士の引き算（桁落ち）が発生し、有効桁数が失われてしまいます。

実際に を変化させながら誤差をプロットしてみると、以下のようになります。

1. 打ち切り誤差支配領域: が大きい領域。テイラー展開の近似誤差が支配的で、 を小さくすると精度が向上する。
2. 丸め誤差支配領域: が小さすぎる領域。桁落ちによる丸め誤差が支配的で、 を小さくすると逆に精度が悪化する。

最適な は、これら2つの誤差のバランスが取れる点に存在します。

実験：最適なhを探る

以下のコードで、 を から まで変化させたときの誤差を確認してみましょう。

fn central_difference<F>(f: F, x: f64, h: f64) -> f64
where
    F: Fn(f64) -> f64,
{
    (f(x + h) - f(x - h)) / (2.0 * h)
}

fn central_difference_error<F, G>(f: F, df_exact: G, x: f64, h: f64) -> f64
where
    F: Fn(f64) -> f64,
    G: Fn(f64) -> f64,
{
    (central_difference(f, x, h) - df_exact(x)).abs()
}

fn main() {
    let f = |x: f64| x.sin();
    let df_exact = |x: f64| x.cos();
    let x: f64 = 1.0;

    println!("h,      Error (Central)");

    let mut h = 1.0;
    for _ in 0..16 {
        h /= 10.0;
        let error = central_difference_error(f, df_exact, x, h);
        println!("{:.1e}, {:.2e}", h, error);
    }
}

実行結果：

h,      Error (Central)
1.0e-1, 9.00e-4
1.0e-2, 9.00e-6
1.0e-3, 9.01e-8
1.0e-4, 9.00e-10
1.0e-5, 1.11e-11
1.0e-6, 2.77e-11
1.0e-7, 1.94e-10
1.0e-8, 2.97e-9
...
1.0e-16, 1.48e-2

倍精度浮動小数点数（f64）を用いる場合、中心差分における最適な はおおよそ 〜 程度になります（関数のスケールにもよりますが、マシンイプシロン に対して 程度が目安と言われます）。

Note

実用的な数値計算ライブラリでは、このような最適な を自動的に推定するアルゴリズムが組み込まれている場合もあります。

2階微分の数値計算

運動方程式 などで登場する2階微分 も、中心差分法を応用して求めることができます。

この近似式の誤差も です。

fn second_derivative_central<F>(f: F, x: f64, h: f64) -> f64
where
    F: Fn(f64) -> f64,
{
    (f(x + h) - 2.0 * f(x) + f(x - h)) / (h * h)
}

まとめ

• 数値微分には前進差分、後退差分、中心差分などの手法がある。
• 中心差分は の精度を持ち、他の2つ（ ）より高精度であるため、通常はこれが推奨される。
• ステップサイズ は小さければ良いわけではなく、丸め誤差（桁落ち） に注意する必要がある。
• f64 を用いる場合、中心差分での最適な は 程度が目安となる。

次節では、数値積分について学びます。