6.2 梯形與辛普森積分#

概念說明#

定積分 ∫ᵃᵇ f(x) dx 的幾何意義是曲線下面積。數值積分（求積，quadrature）的基本策略：把 [a, b] 切成許多小段，在每段上用簡單形狀（直線、拋物線）逼近 f，加總面積。

梯形法：每段用直線連兩端點，面積是梯形。複合梯形公式： T(h) = h·[f(a)/2 + f(x₁) + … + f(x_{n−1}) + f(b)/2]
辛普森法（Simpson's 1/3 rule）：每兩段用通過三點的拋物線。複合辛普森公式： S(h) = h/3·[f(a) + 4Σ奇數點 + 2Σ偶數點 + f(b)]

兩者都是第 7 章插值理論的特例：梯形法是線性插值後積分，辛普森法是二次插值後積分。令人驚喜的是，辛普森法的精度比「二次多項式」應有的更好——因為對稱性使三次項也自動抵消，它其實對三次多項式完全精確。

實務上的標準做法是自動積分：從粗糙的 h 開始，每次加密一倍，直到兩次結果的差小於容許誤差。加密一倍時舊點全部可重用，只需計算新中點，這讓加密幾乎是免費的。

演算法步驟#

自動梯形/辛普森：

取 n = 1，算 T₁ 或 S₁。
加密：n ← 2n，只對新增的中點求 f 值，更新積分值。
估誤差：辛普森法可用 |S₂ₙ − Sₙ|/15 作為誤差估計。
若誤差 < ε 或達最大加密次數，回傳；否則回第 2 步。

Python 程式碼範例#

import math

def trapezoid(f, a, b, n):
    h = (b - a) / n
    s = 0.5 * (f(a) + f(b))
    for i in range(1, n):
        s += f(a + i * h)
    return s * h

def simpson(f, a, b, n):
    """n 必須為偶數"""
    h = (b - a) / n
    s = f(a) + f(b)
    for i in range(1, n):
        s += f(a + i * h) * (4 if i % 2 else 2)
    return s * h / 3

def adaptive_simpson(f, a, b, eps=1e-10, max_doubling=30):
    n = 2
    s_old = simpson(f, a, b, n)
    for i in range(max_doubling):
        n *= 2
        s_new = simpson(f, a, b, n)
        err = abs(s_new - s_old) / 15
        if err < eps:
            return s_new, n, err
        s_old = s_new
    return s_new, n, err

f = math.exp
a, b = 0.0, 1.0
exact = math.e - 1
print(f"∫₀¹ eˣ dx = e − 1 = {exact:.15f}\n")
print(f"{'n':>6} {'梯形誤差':>12} {'辛普森誤差':>12}")
for n in [4, 8, 16, 32, 64]:
    et = abs(trapezoid(f, a, b, n) - exact)
    es = abs(simpson(f, a, b, n) - exact)
    print(f"{n:6d} {et:12.3e} {es:12.3e}")

val, n, err = adaptive_simpson(f, a, b)
print(f"\n自動辛普森: {val:.15f}（n={n}，估計誤差 {err:.1e}）")

# 與 scipy 對照
try:
    from scipy.integrate import quad
    print("scipy.quad:", quad(f, a, b)[0])
except ImportError:
    print("（未安裝 scipy，略過對照）")

誤差表會清楚顯示：n 加倍時梯形誤差除以 4（O(h²)）、辛普森誤差除以 16（O(h⁴)）。

誤差分析#

設 f 光滑、h = (b−a)/n。複合梯形法誤差為 −(b−a)h²f″(ξ)/12；複合辛普森法誤差為 −(b−a)h⁴f⁽⁴⁾(ξ)/180。兩個推論：(1) 若 f″ 有界 M，要達誤差 ε，梯形法需 n ≥ (b−a)³M/(12ε) 開根號；(2) 辛普森對次數 ≤ 3 的多項式零誤差（代數精度 3），梯形法代數精度 1。實作層面還有一個現代必修觀念：對端點不可微或劇烈變化的函數（如 √x 在 0 的奇點），均勻網格的收斂會退化，需要自適應細分（scipy 的 quad 基於 adaptive Gauss-Kronrod）或變數變換。另一個優雅的加速是 Romberg 積分：對梯形序列做 Richardson 外推（見 7.3 節），可把誤差階逐級推到 O(h⁸) 以上。

練習#

用梯形與辛普森積分 f(x) = 1/(1+x²) 在 [0, 1]（真值 = π/4），列誤差表並驗證收斂階分別為 2 與 4。
實作 Romberg 積分：維護梯形值的三角表 T[k][j] = (4ʲT[k][j−1] − T[k−1][j−1])/(4ʲ−1)，對 ∫₀¹ e^(−x²) dx 算到誤差 < 1e-12，數一數總共呼叫 f 幾次。
積分 ∫₀¹ √x dx（在 0 有奇點，真值 2/3）：觀察辛普森法的收斂階退化成多少；再代入 x = t² 把積分改寫成 2∫₀¹ t² dt 後重算，說明變數變換為什麼救回了收斂速度。