数値積分

概要

シンプソンの公式により $\int_{x_a}^{x_b} f(x) \mathrm{d}x$ の近似計算を行う．

まず，ある区間 $[x_l, x_r]$ において $f(x)$ を区間の両端と真ん中 $x_m = (x_l + x_r)/2$ の $3$ 点を通る放物線で近似し，定積分を求めてみる．いわゆる 1/6 公式の形に持ち込めば簡単に計算できる．

\[\int\_{x_l}^{x_r} f(x) \mathrm{d}x = \dfrac{x_r-x_l}{6}[ f(x_l)+4f(x_m)+f(x_r) ]\]

この式をシンプソンの公式と言う．関数を $x_m$ の重みを $4$ 倍して求めた平均をとり，定数関数とみたときの積分 (つまり長方形の面積) という見方ができる．

次に，これだけでは精度が不十分なので $n$ を十分大きい整数とし全体を等間隔の数列 $x_0, \cdots , x_{2n}$ で切った微小区間に分割する．最初と最後は区間の両端に合わせ $x_a=x_0, x_b=x_{2n}$ とする．左から順に隣り合う $2$ つの区間を合わせた区間に対してシンプソンの公式を使って足し合わせる．

\[\begin{eqnarray} \int_{x_a}^{x_b} f(x) \mathrm{d}x &=& \sum_{i=0}^{n-1} \int_{x_{2i}}^{x_{2i+2}} f(x) \mathrm{d}x \nonumber \\ &\simeq& \sum_{i=0}^{n-1} \dfrac{2h}{6} [ f(x_{2i})+4f(x_{2i+1})+f(x_{2i+2}) ] \nonumber \\ &=& \dfrac{h}{3} \sum_{i=0}^{2n} k_i f(x_i) \nonumber \end{eqnarray}\]

ここで $h$ は区間の幅 $(x_i - x_{i-1}) / 2n$ ，$k_i$ は $i$ が $0$ または $2n$ のとき $1$，それ以外の奇数のとき $4$，偶数のとき $2$ である．最後の式を合成シンプソンの公式といい，かなり誤差が小さいらしい．

使い方

引数に与えた f, x1, x2 に対し，$\int_{x_0}^{x_1} f(x) \mathrm{d}x$ の数値積分を計算する．n は許容誤差と制限時間を見て調整．両端で発散する場合は微小量ずらしてごまかす．

実装

ld simpson(ld x0, ld x1, std::function<ld(ld)> f) {
    const int N = 1000;
    const ld dx = (x1 - x0) / N;
    ld s = f(x0) + f(x0 + dx * N);
    for (int i = 1; i < N; i += 2) s += ld(4.0) * f(x0 + dx * i);
    for (int i = 2; i < N; i += 2) s += ld(2.0) * f(x0 + dx * i);
    return s * dx / 3.0;
}

検証

AOJ 2497 http://judge.u-aizu.ac.jp/onlinejudge/review.jsp?rid=1578466#1

参考文献

シンプソンの公式 - Wikipedia https://ja.wikipedia.org/wiki/%E3%82%B7%E3%83%B3%E3%83%97%E3%82%BD%E3%83%B3%E3%81%AE%E5%85%AC%E5%BC%8F