自适应辛普森积分与误差证明

(好像学到了什么不得了且没用的东西)

定积分

365体育投注函数$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}$,取365体育投注区间的右端点$\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$,那么这两种方法365bet体育在线们就都无法使用,因此就引入了数值积分,最常用的就是自适应辛普森积分

辛普森公式

 二次函数积分公式:

令$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$$365bet体育在线们就可以计算365体育投注小区间$[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$的限制,365bet体育在线们应该找365体育投注更加合适的方法。

现在唯一的问题就是如何进行分段。如果段数少了计算误差就大,段数多了时间效率又会低。365bet体育在线们需要找到一个准确度和效率的平衡点。

365bet体育在线们这样考虑:假如有一段图像已经很接近二次函数的话,直接带入公式求积分,得到的值精度就很高了,不需要再继续分割这一段了。

于是365bet体育在线们有了这样365体育投注分割方法:每次判断当前段和二次函数的相似程度,如果足够相似的话就直接代入公式计算,否则将当前段分割成左右两段递归求解。

现在就剩下一个问题了:如果判断每一段和二次函数是否相似?

365bet体育在线们把当前段直接代入公式求积分,再将当前段从中点分割成两段,把这两段再直接代入公式求积分。如果当前段的积分和分割成两段后的积分之和相差很小的话,就可以认为当前段和二次函数很相似了,不用再递归分割了。

上面就是自适应辛普森法的思想。在分治判断的时候,除了判断精度是否正确,一般还要强制执行最少的迭代次数。

代码实现

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$$365bet体育在线们进行$Hermite$插值.由于$Hermite$插值是牛顿插值的极限情形,为此365bet体育在线们先进行牛顿插值.365bet体育在线们设立插值点$$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$$

posted @ 2021-06-15 20:59  Tenderfoot  阅读(134)  评论(0编辑  收藏  举报