數值積分--自适應辛普森法

定積分的定義

簡單來說，函數 f ( x ) f(x) f(x) 在區間 [ l , r ] [l,r] [l,r] 上的定積分 ∫ l r f ( x ) d x \int_{l}^{r}f(x)\mathrm{d}x ∫lr​f(x)dx 指的是 f ( x ) f(x) f(x) 在區間 [ l , r ] [l,r] [l,r] 中與 x x x 軸圍成的區域的面積（其中 x x x 軸上方的部分為正值， x x x 軸下方的部分為負值）。

很多情況下，我們需要高效，準确地求出一個積分的近似值。下面介紹的 辛普森法，就是這樣一種求數值積分的方法。

辛普森法

這個方法的思想是将被積區間分為若幹小段，每段套用二次函數的積分公式進行計算。

對于一個二次函數 f ( x ) = A x 2 + B x + C f(x)=Ax^2+Bx+C f(x)=Ax2+Bx+C，有：

∫ l r f ( x ) d x = ( r − l ) ( f ( l ) + f ( r ) + 4 f ( l + r 2 ) ) 6 \int_l^r f(x) {\mathrm d}x = \frac{(r-l)(f(l)+f(r)+4 f(\frac{l+r}{2}))}{6} ∫lr​f(x)dx=6(r−l)(f(l)+f(r)+4f(2l+r​))​

推導過程：

對于一個二次函數 f ( x ) = A x 2 + B x + C f(x)=Ax^2+Bx+C f(x)=Ax2+Bx+C；

求積分可得 F ( x ) = ∫ 0 x f ( x ) d x = a 3 x 3 + b 2 x 2 + c x + D F(x)=\int_0^x f(x) {\mathrm d}x = \frac{a}{3}x^3+\frac{b}{2}x^2+cx+D F(x)=∫0x​f(x)dx=3a​x3+2b​x2+cx+D 在這裡 D 是一個常數，那麼

∫ l r f ( x ) d x = F ( r ) − F ( l ) = a 3 ( r 3 − l 3 ) + b 2 ( r 2 − l 2 ) + c ( r − l ) = ( r − l ) ( a 3 ( l 2 + r 2 + l r ) + b 2 ( l + r ) + c ) = r − l 6 ( 2 a l 2 + 2 a r 2 + 2 a l r + 3 b l + 3 b r + 6 c ) = r − l 6 ( ( a l 2 + b l + c ) + ( a r 2 + b r + c ) + 4 ( a ( l + r 2 ) 2 + b ( l + r 2 ) + c ) ) = r − l 6 ( f ( l ) + f ( r ) + 4 f ( l + r 2 ) ) \begin{aligned} \int_l^r f(x) {\mathrm d}x &= F(r)-F(l) \\ &= \frac{a}{3}(r^3-l^3)+\frac{b}{2}(r^2-l^2)+c(r-l) \\ &=(r-l)(\frac{a}{3}(l^2+r^2+lr)+\frac{b}{2}(l+r)+c) \\ &=\frac{r-l}{6}(2al^2+2ar^2+2alr+3bl+3br+6c)\\ &=\frac{r-l}{6}((al^2+bl+c)+(ar^2+br+c)+4(a(\frac{l+r}{2})^2+b(\frac{l+r}{2})+c)) \\ &=\frac{r-l}{6}(f(l)+f(r)+4f(\frac{l+r}{2})) \end{aligned} ∫lr​f(x)dx​=F(r)−F(l)=3a​(r3−l3)+2b​(r2−l2)+c(r−l)=(r−l)(3a​(l2+r2+lr)+2b​(l+r)+c)=6r−l​(2al2+2ar2+2alr+3bl+3br+6c)=6r−l​((al2+bl+c)+(ar2+br+c)+4(a(2l+r​)2+b(2l+r​)+c))=6r−l​(f(l)+f(r)+4f(2l+r​))​

我們這樣考慮：假如有一段圖像已經很接近二次函數的話，直接帶入公式求積分，得到的值精度就很高了，不需要再繼續分割這一段了。

于是我們有了這樣一種分割方法：每次判斷目前段和二次函數的相似程度，如果足夠相似的話就直接代入公式計算，否則将目前段分割成左右兩段遞歸求解。

現在就剩下一個問題了：如果判斷每一段和二次函數是否相似？

我們把目前段直接代入公式求積分，再将目前段從中點分割成兩段，把這兩段再直接代入公式求積分。如果目前段的積分和分割成兩段後的積分之和相差很小的話，就可以認為目前段和二次函數很相似了，不用再遞歸分割了。

上面就是自适應辛普森法的思想。在分治判斷的時候，除了判斷精度是否正确，一般還要強制執行最少的疊代次數

代碼解釋：

網上有很多闆子，有的在第一個asr函數中有*15的操作，有的卻沒有；有的在第一個sar有eps/2操作的，有的沒有。

這裡給出一種比較正确的代碼模闆

const double eps = 10e-12;//控制精度，按題目來給
double f(double x) {
    double y = /*函數*/;
    return y;
}
double simpson(double l, double r) {//Simpson公式
    double mid = (l + r) / 2;
    return (f(l) + 4 * f(mid) + f(r)) * (r - l) / 6;
}
double asr(double l, double r, double eps, double ans) {
    double mid = (l + r) / 2;
    double L = simpson(l, mid), R = simpson(mid, r);
    if (fabs(L + R - ans) <= eps * 15) return L + R + (L + R - ans) / 15;     //确認精度
    return asr(l, mid, eps / 2, L) + asr(mid, r, eps/ 2, R);     //精度不夠則遞歸調用
}
double asr(double l, double r, double eps) {
    return asr(l, r, eps, simpson(l, r));
}
           

參考文獻：A Review of Error Estimation in Adaptive Quadrature