台形則・シンプソン則による数値積分｜Python

数値積分は、解析的に積分できない関数や、データ点しか与えられていない場合に積分値を近似する手法だ。台形則とシンプソン則は最も基本的なアルゴリズムであり、原理を理解しておくと SciPy などのライブラリを使う際にも役立つ。

台形則の原理

台形則は、関数のグラフを台形で近似して面積を求める。区間 $[a, b]$ を $n$ 等分したとき、刻み幅は $h = (b - a) / n$ となり、積分は次のように近似される。

$\int_{a}^{b} f (x) d x \approx h [\frac{f ( a ) + f ( b )}{2} + i = 1 \sum n - 1 f (a + ih)]$

def trapezoidal(f, a, b, n):
    """台形則による数値積分"""
    h = (b - a) / n
    result = (f(a) + f(b)) / 2
    for i in range(1, n):
        result += f(a + i * h)
    return result * h

import math

# ∫₀^π sin(x) dx = 2
f = lambda x: math.sin(x)
print(trapezoidal(f, 0, math.pi, 10))   # 1.9835235375094546
print(trapezoidal(f, 0, math.pi, 100))  # 1.9998355038874436
print(trapezoidal(f, 0, math.pi, 1000)) # 1.9999983550656624

$n$ を増やすと真の値 $2$ に近づく。台形則の誤差は $O (h^{2})$ で、二次精度を持つ。

シンプソン則の原理

シンプソン則は、関数を二次多項式（放物線）で近似する。台形則より高精度で、誤差は $O (h^{4})$ となる。 $n$ は偶数でなければならない。

$\int_{a}^{b} f (x) d x \approx \frac{h}{3} [f (a) + f (b) + 4 奇数 i \sum f (x_{i}) + 2 偶数 i \sum f (x_{i})]$

def simpson(f, a, b, n):
    """シンプソン則による数値積分（n は偶数）"""
    if n % 2 != 0:
        raise ValueError("n must be even")
    h = (b - a) / n
    result = f(a) + f(b)
    for i in range(1, n):
        x = a + i * h
        if i % 2 == 0:
            result += 2 * f(x)
        else:
            result += 4 * f(x)
    return result * h / 3

import math

f = lambda x: math.sin(x)
print(simpson(f, 0, math.pi, 10))   # 2.0000067844418012
print(simpson(f, 0, math.pi, 100))  # 2.000000000067847

$n = 10$ でも台形則の $n = 1000$ より正確だ。シンプソン則の威力がわかる。

台形則

二次精度 $O (h^{2})$ 。実装が単純。データ点の個数に制約がない。

シンプソン則

四次精度 $O (h^{4})$ 。同じ $n$ でより高精度。 $n$ は偶数でなければならない。

精度の比較

具体的な例で精度を比較してみよう。 $\int_{0}^{1} e^{x} d x = e - 1 \approx 1.71828$ を計算する。

import math

f = lambda x: math.exp(x)
true_value = math.e - 1

for n in [4, 8, 16, 32]:
    trap = trapezoidal(f, 0, 1, n)
    simp = simpson(f, 0, 1, n)
    print(f'n={n:2d}: 台形則 誤差={abs(trap - true_value):.2e}, '
          f'シンプソン則 誤差={abs(simp - true_value):.2e}')

$n$ を 2 倍にすると、台形則の誤差は約 $1/4$ に、シンプソン則の誤差は約 $1/16$ になる。

複合則と適応型積分

上で示した公式は「複合台形則」「複合シンプソン則」と呼ばれ、区間を分割して適用している。さらに高度な手法として、誤差を推定しながら分割を調整する「適応型積分」がある。SciPy の quad はこの適応型積分を使っている。

ロンバーグ積分

ロンバーグ積分は、台形則の結果を外挿して高精度化する手法だ。リチャードソン外挿を繰り返し適用する。

def romberg(f, a, b, max_iter=10, tol=1e-10):
    """ロンバーグ積分"""
    R = [[0] * (max_iter + 1) for _ in range(max_iter + 1)]
    
    # 初期値（台形則、n=1）
    h = b - a
    R[0][0] = h * (f(a) + f(b)) / 2
    
    for i in range(1, max_iter + 1):
        h /= 2
        # 台形則の細分化
        sum_mid = sum(f(a + (2*k - 1) * h) for k in range(1, 2**(i-1) + 1))
        R[i][0] = R[i-1][0] / 2 + h * sum_mid
        
        # リチャードソン外挿
        for j in range(1, i + 1):
            R[i][j] = R[i][j-1] + (R[i][j-1] - R[i-1][j-1]) / (4**j - 1)
        
        if i > 1 and abs(R[i][i] - R[i-1][i-1]) < tol:
            return R[i][i]
    
    return R[max_iter][max_iter]

import math
f = lambda x: math.sin(x)
print(romberg(f, 0, math.pi))  # 2.0000000000000004

ロンバーグ積分は少ない関数評価で高精度な結果を得られる。

ガウス求積法

ガウス求積法は、評価点と重みを最適に選ぶことで、少ない点数で高精度を達成する手法だ。 $n$ 点のガウス・ルジャンドル求積法は、 $2 n - 1$ 次までの多項式を正確に積分できる。

# 2 点ガウス・ルジャンドル求積法（区間 [-1, 1]）
def gauss_legendre_2(f):
    """2 点ガウス求積法"""
    x = [-1/math.sqrt(3), 1/math.sqrt(3)]
    w = [1, 1]
    return sum(w[i] * f(x[i]) for i in range(2))

# 区間 [a, b] への変換
def gauss_quad(f, a, b):
    """区間 [a, b] でのガウス求積法"""
    def g(t):
        x = (b - a) / 2 * t + (a + b) / 2
        return f(x) * (b - a) / 2
    return gauss_legendre_2(g)

import math
f = lambda x: x**3  # 3 次多項式は 2 点で正確に積分可能
print(gauss_quad(f, 0, 1))  # 0.25（真の値と一致）

実用上の選択

簡単な関数

シンプソン則で十分。SciPy の quad を使うのが最も簡単。

データ点のみ

台形則（numpy.trapz）やシンプソン則（scipy.integrate.simpson）を使う。

高精度が必要

ロンバーグ積分やガウス求積法。SciPy の romberg や fixed_quad が利用可能。

原理を理解した上でライブラリを使うと、適切な手法を選べるようになる。