自適應辛普森求積法 C++ 實現：變步長遞歸精度控制完整解析

數值積分是科學計算的基石，辛普森法則以其 O(h⁴) 的高精度廣受青睐。但固定步長難以兼顧效率與精度。本文深度解析一個變步長自適應辛普森求積的 C++ 實現，展示如何通過步長遞減與誤差估計實現精度自動控制。

核心演算法原理

自適應辛普森法基於龍格法則：

1. 計算步長 h 與 h/2 下的辛普森近似值 S(h)、S(h/2)
2. 誤差估計：E ≈ |S(h/2) - S(h)| / 15
3. 若 E < ε 則接受結果；否則步長減半遞歸

本實現採用迭代加倍步數策略：m, 2m, 4m, 8m...，複用前一次計算的函數值，大幅降低計算量。

程式碼結構解析

被積函數定義

double cal0(double x) { return 6 * exp(-x*x); }      // I₁ = ∫₀¹ 6e⁻ˣ² dx
double cal1(double x) { return 3 * log(x); }         // I₂ = ∫₁² 3ln(x) dx
double cal2(double x) { return sin(x) / x; }         // I₃ = ∫₁⁴ sin(x)/x dx

核心求積函數 simpson()

void simpson(double a, double b, int m, double (*g)(double)) {
    double h[3], q[3], p[3], s[3], t;
    int i, k;

    // 初始步長 h₀ = (b-a)/(2m)
    h[0] = (b - a) / (2 * m);
    // 計算 S₀ (m 區間)
    // ... 計算 p[0], s[0] ...

    // 步長減半 h₁ = h₀/2，計算 S₁ (2m 區間)
    // 利用 q[0] 存儲新增中點函數值，複用 p[0]
    // s[1] = h₁ * (p[0] + 4*q[0]) / 3

    // 步長再減半 h₂ = h₁/2，計算 S₂ (4m 區間)
    // q[1] 存儲新增中點，p[1] = p[0] + 2*q[0]
    // s[2] = h₂ * (p[1] + 4*q[1]) / 3

    // 迭代循環：k = 2, 3, ...
    for (k = 2; ; k++) {
        h[(k+1)%3] = (b - a) / (pow(2, k+2) * m);
        // 計算新中點函數值 q[k%3]
        // 更新 p[k%3] = p[(k-1)%3] + 2*q[(k-1)%3]
        // 計算 s[(k+1)%3]
        // 收斂判斷：|s[new] - s[old]| < 1e-7
    }
}

關鍵優化：循環陣列複用

使用大小為 3 的循環陣列 h[3]、q[3]、p[3]、s[3]，透過模運算 k%3 實現滾動儲存，空間複雜度 O(1)，避免動態記憶體分配開銷。

輸入輸出範例

// in.dat
0 1 2      // 積分區間 [0,1], 初始 m=2/2=1
1 2 2      // 積分區間 [1,2], 初始 m=1
1 4 2      // 積分區間 [1,4], 初始 m=1

// out.out
I(1)=4.48094
積分區間數目為:16
I(2)=1.15888
積分區間數目為:32
I(3)=0.81212
積分區間數目為:32

精度與效能分析

現代化改進建議

• 使用 std::function<double(double)> 替代函數指標，支援 Lambda 與仿函數
• 引入 constexpr 誤差容限，支援編譯期配置
• 加入最大遞歸深度防護，避免病態函數無限循環
• 使用 long double 提升中間計算精度
• 考慮 std::vector 替代固定陣列，支援動態精度需求

總結

這個實現展示了經典數值方法的工程化智慧：誤差驅動自適應、中點複用減少函數評估、循環緩衝區恆定記憶體。對於科學計算工程師，理解其背後的龍格外推原理與複用策略，比單純複製程式碼更有價值。下次面對積分問題，不妨嘗試自己實現一個自適應版本——這是通往高性能計算的必經之路。