自適應辛普森積分與誤差證明

Tenderfoot發表於2021-06-15

(好像學到了什麼不得了且沒用的東西)

定積分

表示函式$f(x)$在區間$[a,b]$上積分和的極限,記作$$\int_{a}^{b}f(x)dx$$通俗地講,就是該段函式與$x$軸圍成的面積

計算方法一:

分割區間$$\int_{a}^{b}f(x)dx=\lim_{n\rightarrow ∞}\sum_{i=1}^{n}\frac{b-a}{n}f(a+\frac{b-a}{n}i)$$

計算方法二:

牛頓-萊布尼茨公式$$\int_{a}^{b}f'(x)dx=f(b)-f(a)$$

例如:計算$\int_{0}^{1}e^xdx$

法一:在$[0,1]$上插入$n-1$個分點$x_k(1 \leqslant k \leqslant n-1)$,即$$0=x_0<x_1<x_2<\dots <x_{n-1}<x_n=1$$將區間$n$等分,則每段區間的長度均為$\Delta x_k=\frac{1}{n}$,第$k$個分點為$x_k=\frac{k}{n}$,取每個區間的右端點$\xi_k=\frac{k}{n}$作積分和,就有$$\int_{0}^{1}e^xdx=\lim_{n\rightarrow ∞}\frac{1}{n}\sum_{k=1}^{n}e^{\frac{k}{n}}\\=\lim_{n\rightarrow ∞}\frac{e^{\frac{1}{n}}(e-1)}{n(e^{\frac{1}{n}}-1)}\\=\lim_{n\rightarrow ∞}\frac{e^{\frac{1}{n}}(e-1)}{n·\frac{1}{n}}\\=e-1$$

法二:令$f(x)=e^x$,則可設$F(x)=\int e^xdx=e^x$,即$$F'(x)=f(x)\\ \int_{0}^{1}e^xdx=F(1)-F(0)\\=e-1$$

如果所給的$f(x)$不容易求出$n\rightarrow ∞$時的積分和或者不容易求出$\int f(x)dx$,那麼這兩種方法我們就都無法使用,因此就引入了數值積分,最常用的就是自適應辛普森積分

辛普森公式

 二次函式積分公式:

令$f(x)=Ax^2+Bx+C$,對其進行積分$F(x)=\int_{0}^{x}f(x)dx=\frac A 3x^3+\frac B 2x^2+Cx+D$,$D$是一個常數,那麼就有$$\int_{l}^{r}f(x)dx=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))$$

給定一個自然數$n$,將區間$[l,r]$,分成$2n$個等長的區間$x$。其中$$x_i=l+ih,i=0\dots 2n,h=\frac {r-l}2n$$我們就可以計算每個小區間$[x_{2i-2},x_{2i}],i=1\dots n$的積分值,將所有區間的積分值相加即為總積分

對於$[x_{2i-2},x_{2i}],i=1\dots n$的一個區間,選取其中的三個點$(x_{2i-2},x_{2i-1},x_{2i})$就可以得到一條經過該三點的唯一的拋物線$P(x)$。於是問題從計算原函式在該區間的積分值變成了計算新的二次函式$P(x)$在該區間積分值,就可以用辛普森公式近似計算之$$\int_{x_{2i-2}}^{x_{2i}}f(x)dx\approx \int_{x_{2i-2}}^{x_{2i}}P(x)dx=\frac h 3(f(x_{2i-2})+f(x_{2i})+4f(x_{2i-1}))$$將其分段求和可得到如下結論$$\int_l^r f(x)dx\approx \frac h 3(f(x_0)+4f(x_1)+2f(x_2)+4f(x_3)+2f(x_4)+\dots +4f(x_{2N-1})+f(x_{2N}))$$誤差為(後文有推導)$$-\frac 1{90}(\frac{r-l}2)^5f^{(4)}(\xi)$$其中$\xi$是位於區間$[l,r]$的某個值

程式碼實現

const int N = 1000 * 1000;

inline double Simpson_Integration(double a, double b) 
{
  double h = (b - a) / N;
  double s = f(a) + f(b);
  for(int i = 1; i <= N - 1; i++) 
  {
    double x = a + h * i;
    s += f(x) * ((i & 1) ? 4 : 2);
  }
  s *= h / 3;
  return s;
}

自適應辛普森積分

普通的方法為保證精度在時間方面無疑會受到 $n$的限制,我們應該找一種更加合適的方法。

現在唯一的問題就是如何進行分段。如果段數少了計算誤差就大,段數多了時間效率又會低。我們需要找到一個準確度和效率的平衡點。

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

於是我們有了這樣一種分割方法:每次判斷當前段和二次函式的相似程度,如果足夠相似的話就直接代入公式計算,否則將當前段分割成左右兩段遞迴求解。

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

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

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

程式碼實現

double Simpson(double l, double r) 
{
    double mid = (l + r) / 2;
    return (r - l) * (f(l) + 4 * f(mid) + f(r)) / 6;
}
double Asr(double a , double b , double Eps , double A)
{
    double Middle = a + (b  -a) / 2;
    double L = Simpson(a  ,Middle) , R = Simpson(Middle , b);
    if(fabs(L + R - A) <= 15 * Eps) 
        return L + R + (L + R - A) / 15;
    else 
        return Asr(a , Middle , Eps / 2 , L) + Asr(Middle , b , Eps  / 2 , R);
}

誤差證明

$f$在$(x_0,x_2)$上四階可導,且在$[x_0,x_2]$上三階導函式連續,求證$$\int_{x_0}^{x_2}f(x)dx=\frac h 3[f(x_0)+4f(x_1)+f(x_2)]-\frac{h^5}{90}f^{(4)}(\xi)$$其中$h=x_1-x_0=x_2-x_1$

引入

差商

給定函式$f(x)$和插值節點$x_0,x_1,\cdots ,x_n$,則函式$f(x)$關於節點$x_0,x_1,\cdots ,x_n$的$k$階差商可遞迴定義為$$f[x_0,x_1,\dots x_k]=\frac{f[x_1,x_2,\cdots,x_k]-f[x_0,x_1,\cdots ,x_{k-1}]}{x_k-x_0}$$,其中$f(x)$關於節點$x_i$的$0$階差商定義為其函式值,即$f[x_i]=f(x_i)$

牛頓插值公式

$$f(x)=f(x_0)+(x-x_0)f[x_0,x_1]+(x-x_0)(x-x_1)f[x_0,x_1,x_2]+\cdots +(x-x_0)(x-x_1)\cdots (x-x_{n-1})f[x_1,x_2,\cdots,x_n]$$可記作$$f(x)=N_n(x)+R_n(x)$$當$n\rightarrow \infty$時$R_n(x)=0$  

埃爾米特插值多項式

不少實際的插值問題不但要求在節點上的函式值相等,而且還要求對應的導數值也相等,甚至要求高階導數也相等,滿足這種要求的插值多項式就是埃爾米特插值多項式。$Hermite$插值是牛頓插值的極限情形

證明

進行牛頓插值,設插值點為$$x_0,x_1,x_2$$則經過此三點的牛頓插值多項式為$$f(x_0)+(x-x_0)f[x_0,x_1]+(x-x_0)(x-x_1)f[x_0,x_1,x_2]$$,插值餘項為$$(x-x_0)(x-x_1)(x-x_2)f[x,x_0,x_1,x_2]$$先來看$$\int_{x_0}^{x_2}\left \{f(x_0)+(x-x_0)f[x_0,x_1]+(x-x_0)(x-x_1)f[x_0,x_1,x_2]\right \}dx$$可以先算$$\int_{x_0}^{x_2}f(x_0)dx=2hf(x_0)\\\int_{x_0}^{x_2}(x-x_0)f[x_0,x_1]dx=\frac 1 2 (x_2-x_0)^2f[x_0,x_1]\\=\frac 1 2(2h)^2\frac{f(x_1)-f(x_0)}{h}\\=2h(f(x_1)-f(x_0))\\\int_{x_0}^{x_2}(x-x_0)(x-x_1)f[x_0,x_1,x_2]dx\\= \frac 1 6 f[x_0,x_1,x_2](x_2-x_0)^2(2x_2+x_0-3x_1)$$其中$$f[x_0,x_1,x_2]=\frac{f(x_0)}{2h^2}+\frac {f(x_1)}{-h^2}+\frac{f(x_2)}{2h^2}$$因此$$\int_{x_0}^{x_2}(x-x_0)(x-x_1)f[x_0,x_1,x_2]dx=\frac{[f(x_0)-2f(x_1)+f(x_2)]h}{3}\\\int_{x_0}^{x_2}\left \{f(x_0)+(x-x_0)f[x_0,x_1]+(x-x_0)(x-x_1)f[x_0,x_1,x_2]\right \}dx\\=\frac{f(x_0)+4f(x_1)+f(x_2)h}{3}$$現在來證明$$\int_{x_0}^{x_2}(x-x_0)(x-x_1)(x-x_2)f[x,x_0,x_1,x_2]dx=\frac{-h^5}{90}f^{(4)}(\xi)$$令$$g(t)=\int_{x_0}^{t}(x-x_0)(x-x_1)(x-x_2)f[x,x_0,x_1,x_2]dx$$則$$g'(t)=(t-x_0)(t-x_1)(t-x_2)f[t,x_0,x_1,x_2]\\g'(x_0)=g'(x_1)=g'(x_2)=0,g(x_0)=0$$我們進行$Hermite$插值.由於$Hermite$插值是牛頓插值的極限情形,為此我們先進行牛頓插值.我們設立插值點$$x_0,x'_0,x_1,x'_1,x_2,x'_2$$經過這幾個插值點的牛頓插值多項式為$$p(x)=g(x_0)+(x-x_0)g[x_0,x'_0]\\+(x-x_0)(x-x'_0)g[x_0,x'_0,x_1]\\+(x-x_0)(x-x'_0)(x-x_1)g[x_0,x'_0,x_1,x'_1]\\+(x-x_0)(x-x'_0)(x-x_1)(x-x'_1)g[x_0,x'_0,x_1,x'_1,x_2]\\+(x-x_0)(x-x'_0)(x-x_1)(x-x'_1)(x-x_2)g[x_0,x'_0,x_1,x'_1,x_2,x'_2]$$令$x'_0\rightarrow x_0,x'_1\rightarrow x_1,x'_2\rightarrow x_2$可得到$Hermite$插值多項式為$$q(x)=g(x_0)+(x-x_0)g[x_0,x'_0]\\+(x-x_0)^2g[x_0,x_0,x_1]\\+(x-x_0)^2(x-x_1)g[x_0,x_0,x_1,x_1]\\+(x-x_0)^2(x-x_1)^2g[x_0,x_0,x_1,x_1,x_2]\\+(x-x_0)^2(x-x_1)^2(x-x_2)g[x_0,x_0,x_1,x_1,x_2,x_2]$$令$$k(x)=g(x)-q(x)$$可得$$k(x_0)=k(x_1)=k(x_2)=0,k'(x_0)=k'(x_1)=k'(x_2)=0$$使用多次羅爾中值定理可以得到$$k^{(5)}(\xi)=0$$即$$g^{(5)}(\xi)=5!g[x_0,x_0,x_1,x_1,x_2,x_2]$$易得$$g^{(5)}(\xi)=f^{(4)}(\xi)$$,因此$$f^{(4)}(\xi)=5!g[x_0,x_0,x_1,x_1,x_2,x_2]$$而$$g[x_0,x_0,x_1,x_1,x_2,x_2]=\frac{g[x_0,x_0,x_1,x_1,x_2]-g[x_0,x_1,x_1,x_2,x_2]}{x_0-x_2}\\=\frac{\frac{g[x_0,x_0,x_1,x_1]-g[x_0,x_1,x_1,x_2]}{x_0-x_2}-\frac{g[x_0,x_1,x_1,x_2]-g[x_1,x_1,x_2,x_2]}{x_0-x_2}}{x_0-x_2}\\=\frac{g[x_0,x_0,x_1,x_1]-2g[x_0,x_1,x_1,x_2]+g[x_1,x_1,x_2,x_2]}{(x_0-x_2)^2}\\=\frac{\frac{g[x_0,x_0,x_1]-g[x_0,x_1,x_1]}{x_0-x_1}-2\frac{g[x_0,x_1,x_1]-g[x_1,x_1,x_2]}{x_0-x_2}+\frac{g[x_1,x_1,x_2]-g[x_1,x_2,x_2]}{x_1-x_2}}{(x_0-x_2)^2}\\=\frac{g[x_0,x_0,x_1]-2g[x_0,x_1,x_1]+2g[x_1,x_1,x_2]-g[x_1,x_2,x_2]}{(x_0-x_2)^2(x_0-x_1)}\\=\frac{3\frac{g(x_1)}{(x_0-x_1)^2}-3\frac{g(x_1)-g(x_2)}{(x_1-x_2)^2}}{(x_0-x_2)^2(x_0-x_1)}\\=\frac{3g(x_2)}{(x_1-x_2)^2(x_0-x_2)^2(x_0-x_1)}\\=\frac{3g(x_2)}{-4h^2}$$可見$$-4h^5\frac{f^{(4)}(\xi)}{5!\times 3}=g(x_2)$$即$$g(x_2)=\frac{-h^5}{90}f^{(4)}(\xi)$$即$$\int_{x_0}^{x_2}(x-x_0)(x-x_1)(x-x_2)f[x,x_0,x_1,x_2]dx=\frac{-h^5}{90}f^{(4)}(\xi)$$成立$$Q.E.D$$

相關文章