シンプソンの公式

学校のレポート，提出して終わりはもったいないなと思ったんで供養（？）

概要

シンプソンの公式（Simpson's rule） は微小区間の関数値を放物線で近似するものである．

区間$[-h,h]$におけるシンプソンの公式

図1: シンプソンの公式と誤差

図1のように$x$軸上に距離$h$で均等に$x_0, x_1, x_2$をとる．して，関数値は$f(x_0), f(x_1), f(x_2)$となる．

計算の簡単のため，中点$x_1$が$0$になるよう座標を$x_1$だけ移動させると

• 左: $(-h, f(x_0))$
• 中: $(0, f(x_1))$
• 右: $(h, f(x_2))$

さて，近似のための放物線を

P(x)=Ax^2+Bx+C\quad(A\ne0)

とおくと

\begin{cases}
&P(-h)=f(x_0)=Ah^2-Bh+C\\
&P(0)=f(x_1)=C\\
&P(h)=f(x_2)=Ah^2+Bh+C
\end{cases}

これを係数について解き

\begin{cases}
&A=\frac{f(-h)-2f(0)+f(h)}{2h^2}\\
&B=\frac{f(h)-f(-h)}{2h}\\
&C=f(0)
\end{cases}

（式(3)とする）

そしてこの近似曲線$f(x)$を$-h$から$h$まで積分して係数を代入すると

\int^{h}_{-h}(Ax^2+Bx+C)dx

偶関数，奇関数の性質より

=2\int^h_0(Ax^2+C)dx \\
=2\left[\frac{1}{3}Ax^3+Cx\right]^{h}_{0} \\
=2\left(\frac{1}{3}Ah^3+Ch\right) \\
=\frac{1}{3}h(2Ah^2+6C)

式(3)より

=\frac{1}{3}h(f(-h)-2f(0)+f(h)+6f(0)) \\
=\frac{1}{3}h(f(-h)+4f(0)+f(h))

$h=\frac{x_2-x_0}{2}$であるから

\int^{h}_{-h}(Ax^2+Bx+C)dx=\frac{x_2-x_0}{6}\left(f(-h)+4f(0)+f(h)\right)

（式(4)とする）

微小区間$h$で$2n$等分されている積分区間$[a, b]$におけるシンプソンの公式

あとは，このニコイチの区間¹の計算を、積分区間全体にわたって足し合わせる．

• 最初，$[x_0, x_2]$の面積: $h/3 (f(x_0) + 4f(x_1) + f(x_2))$
• 次，$[x_2, x_4]$の面積: $h/3 (f(x_2) + 4f(x_3) + f(x_4))$
• その次以降，$[x_4, x_6]$の面積: $h/3 (f(x_4) + 4f(x_5) + f(x_6))+\cdots$
• 最後，$[x_{2n-2}, x_{2n}]$の面積: $h/3 (f(x_{2n-2}) + 4f(x_{2n-1}) + f(x_{2n}))$

これらをすべて足し合わせる．このとき，$x_2, x_4, x_6 \cdots$ といったニコイチ区間の境界となる点（偶数番目の点）は，隣り合う2区間の両方で1回ずつ足されるため，係数が2となる．

したがって以下の公式を得る．

\int_{a}^{b} f(x) dx \approx \frac{b-a}{6} \left( f(x_0) + 4 \sum_{i=1}^{n} f(x_{2i-1}) + 2 \sum_{i=1}^{n-1} f(x_{2i}) + f(x_n) \right)

書き下すと以下のようになる．

\int_a^bf(x)dx\approx{\frac{b-a}{6}(f(x_0)+4f(x_1)+2f(x_2)+\cdots+2f(x_{2n-2})+4f(x_{2n-1})+f(x_{2n}))}

プログラム

from typing import Callable, TypeAlias

interval: TypeAlias = tuple[float, float]
fn: TypeAlias = Callable[[float], float]


def integral(interval: interval, fn: fn, n: int = 1000) -> float:
    a, b = interval
    if n % 2 != 0:
        raise ValueError("n must be even for Simpson's rule")
    h = (b - a) / n

    total = fn(a) + fn(b)
    for i in range(1, n, 2):
        total += 4 * fn(a + i * h)
    for i in range(2, n - 1, 2):
        total += 2 * fn(a + i * h)

    return total * h / 3


def main():
    print(f"{integral((0, 3), lambda x: x**4 + 2 * x) = }")


if __name__ == "__main__":
    main()

実行結果

integral((0, 3), lambda x: x**4 + 2 * x) = 57.60000000003228

考察

シンプソンの公式は放物線で任意の区間を近似する．ここで式(4)を見ると，右辺は，全体を6で割る形をしている．
これには一つ思い当たる節がある．1/6公式である．同式も放物線と$1/6$に関係している．

おそらく2次補完誤差を計算する際に，$(x-a)(x-b)$のような式が現れるため，結果的に1/6公式と同じ構造が現れるものかと思う．

1. 1区間につき微小区間$h$二つ分であることの意．ニコイチ区間単位で数えれば$n$等分である． ↩