數值積分

什麼是數值積分

　　數值積分可以用來求定積分的近似值。對于很多函數來說，我們是可以使用初等函數來表示出其積分的，對于這種函數，隻需要求出不定積分然後代入值就能得到定積分了。

　　可是除此之外還有許多難求的函數和沒法使用初等函數表示的函數。當我們想要求出它們的定積分的時候，需要使用數值積分來求解。

　　在ACM中一些題目需要使用數值積分來求解，以下列出一些求數值積分的方法，由簡單到難，而對ACMer來說最重要的是複合Simpson，其精度較高，且可調精度，是亂搞積分幾何的利器。

　　我從這學的：網易公開課 MIT 數值積分

　　學習的契機是想要A掉這題：ZOJ 3898 我的題解

法一·黎曼和

黎曼和是用将區間等長分為 n 段，然後用矩形去逼近函數，每段的長為Δx。

可以選擇每段左側的函數值作為矩形的高，也可以選擇每段右側的函數值作為矩形的高。

若設 n+1 個函數值從左至右為 x0,x1⋯xn ，可得如下公式：

Left Hand Riemann=Right Hand Riemann=Δxf(x0)+Δxf(x1)+⋯+Δxf(xn−1)=Δx∑i=0n−1f(xi)Δxf(x1)+Δxf(x2)+⋯+Δxf(xn)=Δx∑i=1nf(xi)

這種逼近比較粗糙，不過能比較好地傳達數值積分的概念。

法二·梯形法

黎曼和雖然簡單，但是精度堪憂，沒法很好地模拟逼近函數，接下來介紹第二種方法。

可以看出用矩形逼近的時候有很多空缺，使用梯形去逼近，就能大大提高精度了，看上去很像了。

沿用之前的 xi ，由梯形公式我們可以得到如下公式：

Trapezoid=Δx(f(x0)+f(x1))2+Δx(f(x1)+f(x2))2+⋯+Δx(f(xn−1)+f(xn))2=Δx(f(x0)2+f(x1)+⋯+f(xn−1)+f(xn)2)=Left Hand Riemann+Right Hand Riemann2

可以看出梯形法求出的就是左右黎曼和的平均值。

公式中第一個和最後一個需要除以 2 其他都能合出來1個。

法三·Simpson公式

前兩種都比較簡單，但是精度比較差，第三種方法是用抛物線去逼近。

其實我并不知道是怎麼計算的，是從公開課上學來的。

使用 Simpson 公式首先需要 n 為偶數。

将整個區間分為n2段，每段的底為 2Δx ，高比較難搞，但是是有公式的： Simpson height=f(xl)+4f(xm)+f(xr)6

于是有公式：

Simpson=Δx3(f(x0)+4f(x1)+2f(x2)+⋯+2f(xn−2)+4f(xn−1)+f(xn))

法四·複合Simpson

Simpson 公式的精度其實已經相當不錯了，可是相對于ACM中所需的精度仍然有差距。

為了提高精度，我們需要多次重複利用 Simpson 公式。

我們先定義被積函數為 f(x) ，定義 Simpson(l,r)=(r−l)f(l)+4f(l+r2)+f(r)6

這也就是正常的 Simpson 公式，再定義函數 RSimpson 為複合 Simpson 。

當我們要求 RSimpson(l,r) ，先令 m=l+r2

當 Simpson(l,r)≈Simpson(l,m)+Simpson(m,r) 時我們就認為精度夠了傳回其中一個。

當不滿足的時候我們就再分段去求 RSimpson(l,m)+RSimpson(m,r)

這樣得到的精度就比較高了，而且通過定義 ≈ 的範圍可以調整精度。

總結公式如下：

RSimpson(l,r)={Simpson(l,r)                                              approximateRSimpson(l,m)+RSimpson(m,r)       else

複合Simpson的實作

inline double getAppr(double fl, double fm, double fr, double l, double r) {
    return (fl+*fm+fr)*(r-l)/;
}

double Simpson(double l, double r, double fl, double fr) {
    double m = (l+r)/, lm = (l+m)/, rm = (r+m)/;
    double fm = f(m), flm = f(lm), frm = f(rm);
    double vlr = getAppr(fl, fm, fr, l, r);
    double vlm = getAppr(fl, flm, fm, l, m);
    double vrm = getAppr(fm, frm, fr, m, r);
    return fabs(vlr-vlm-vrm) < EPS ? vlr : Simpson(l, m, fl, fm)+Simpson(m, r, fm, fr);
}
           

複合Simpson的實作2(Natureal的代碼)

inline double getAppr(double l,double r){
    return (f(l) + *f((l+r)/) + f(r)) * (r - l) / ;
}

double Simpson(double l,double r){
    double sum = getAppr(l,r);
    double mid = (l+r)/;
    double suml = getAppr(l,mid);
    double sumr = getAppr(mid,r);
    return (fabs(sum - suml - sumr) < EPS) ? sum : Simpson(l, mid) + Simpson(mid, r);
}