数値積分は、解析的に積分できない関数や、データ点しか与えられていない場合に積分値を近似する手法だ。台形則とシンプソン則は最も基本的なアルゴリズムであり、原理を理解しておくと SciPy などのライブラリを使う際にも役立つ。
台形則の原理
台形則は、関数のグラフを台形で近似して面積を求める。区間 を 等分したとき、刻み幅は となり、積分は次のように近似される。
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
を増やすと真の値 に近づく。台形則の誤差は で、二次精度を持つ。
シンプソン則の原理
シンプソン則は、関数を二次多項式（放物線）で近似する。台形則より高精度で、誤差は となる。 は偶数でなければならない。
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
でも台形則の より正確だ。シンプソン則の威力がわかる。
二次精度 。実装が単純。データ点の個数に制約がない。
四次精度 。同じ でより高精度。 は偶数でなければならない。
精度の比較
具体的な例で精度を比較してみよう。 を計算する。
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}')
を 2 倍にすると、台形則の誤差は約 に、シンプソン則の誤差は約 になる。
複合則と適応型積分
上で示した公式は「複合台形則」「複合シンプソン則」と呼ばれ、区間を分割して適用している。さらに高度な手法として、誤差を推定しながら分割を調整する「適応型積分」がある。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
ロンバーグ積分は少ない関数評価で高精度な結果を得られる。
ガウス求積法
ガウス求積法は、評価点と重みを最適に選ぶことで、少ない点数で高精度を達成する手法だ。 点のガウス・ルジャンドル求積法は、 次までの多項式を正確に積分できる。
# 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 が利用可能。
原理を理解した上でライブラリを使うと、適切な手法を選べるようになる。
