自适应辛普森积分与误差证明
（好像学到了什么不得了且没用的东西）
定积分
表示函数$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 法二：令$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]$的某个值 普通的方法为保证精度在时间方面无疑会受到 $n$的限制，我们应该找一种更加合适的方法。 现在唯一的问题就是如何进行分段。如果段数少了计算误差就大，段数多了时间效率又会低。我们需要找到一个准确度和效率的平衡点。 我们这样考虑：假如有一段图像已经很接近二次函数的话，直接带入公式求积分，得到的值精度就很高了，不需要再继续分割这一段了。 于是我们有了这样一种分割方法：每次判断当前段和二次函数的相似程度，如果足够相似的话就直接代入公式计算，否则将当前段分割成左右两段递归求解。 现在就剩下一个问题了：如果判断每一段和二次函数是否相似？ 我们把当前段直接代入公式求积分，再将当前段从中点分割成两段，把这两段再直接代入公式求积分。如果当前段的积分和分割成两段后的积分之和相差很小的话，就可以认为当前段和二次函数很相似了，不用再递归分割了。 上面就是自适应辛普森法的思想。在分治判断的时候，除了判断精度是否正确，一般还要强制执行最少的迭代次数。 $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$$辛普森公式
二次函数积分公式：
代码实现
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;
}
自适应辛普森积分
代码实现
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);
}
误差证明
引入
差商
牛顿插值公式
埃尔米特插值多项式
证明