【数值积分】自适应Simpson法
includeCPP
·
·
算法·理论
自适应 \mathrm{Simpson} 法
对于一个在区间 [a,b] 上可积的函数，我们有很多方式求其在区间 [a,b] 上的定积分 \int_a^bf(x)\mathrm{d}x ，本文将简单介绍其中的 \mathrm{Simpson} 法及常用于计算机科学中的自适应 \mathrm{Simpson} 法。
我们设 $g(x)=ax^2+bx+c$ 为拟合后的函数，则有：
$$\int_l^rf(x)\mathrm{d}x$$
$$≈ \frac{1}{3}a(r^3-l^3)+\frac{1}{2}b(r^2-l^2)+c(r-l)$$
$$=\frac{(r-l)}{6}(2a(l^2+lr+r^2)+3b(l+r)+6c)$$
$$=\frac{(r-l)}{6}(al^2+bl+c+ar^2+br+c+4a(\frac{l+r}{2})^2+4b(\frac{l+r}{2})+4c)$$
$$=\frac{(r-l)}{6}(f(l)+4f(\frac{l+r}{2})+f(r))$$
这样，我们就得到了 $\mathrm{Simpson}$ 公式，即：
$$\int_l^rf(x)\mathrm{d}x≈\frac{(r-l)}{6}(f(l)+4f(\frac{l+r}{2})+f(r))$$
这个式子非常的好，因为我们对于每一个小区间 $[l,r]$ ，可以用 $O(1)$ 的复杂度计算出拟合这个区间的抛物线的定积分值，那么假设我们把函数 $f(x)$ 分为 $n$ 个等长的区间，对于每个区间用 $\mathrm{Simpson}$ 法拟合求积分值，我们便可以在 $O(n)$ 的复杂度内计算出 $\int_a^bf(x)\mathrm{d}x$ 的近似值。
但是在实际问题中，我们会出现一个问题：该把函数分为几个区间分别拟合呢？
如果区间数小了，拟合的精度就不够了。区间数大了，那速度就太慢了。因此，我们需要一种可以自动调整区间分割大小的方法，来在保证一定精度的前提下最快地计算出定积分的近似值，这就是**自适应** $\mathrm{Simpson}$ **法**。
我们考虑二分递归，每次把区间 $[l,r]$ 分为左右两个区间，然后对左右两个子区间分别求定积分值，然后继续二分递归下去。当我们在某次二分递归一个区间时，发现再进行二分递归，其左右两个子区间的定积分值加起来与本区间的定积分值相差小于 $\varepsilon$ ，那我们就停止递归这个区间，这样，我们就在保证一定精度的前提下快速求出了 $\int_a^bf(x)\mathrm{d}x$ 的近似值。
在实际问题中，要注意根据题面数据范围，合理设置 $\varepsilon$ 的值。
参考模板代码：
```cpp
const double eps=1e-9;//根据题面合理设置值
double f(double x){
//此处填写函数
}
double simpson(double a,double b){//Simpson公式
double c=(a+b)/2;
return (b-a)*(f(a)+4*f(c)+f(b))/6;
}
double work(double l,double r,int stp){//自适应Simpson法
double mid=(l+r)/2;
if(abs(simpson(l,mid)+simpson(mid,r)-simpson(l,r))<eps||stp>=1e5){
return simpson(l,r);
}
return work(l,mid,stp+1)+work(mid,r,stp+1);
}
```