浅谈Simpson积分
NaCly_Fish
2018-10-02 22:26:36
### 注意：在看这篇文章之前，请确保你已经有微积分的基础知识
Simpson 积分有什么用？
别急，先来看一下这题：[P4525](https://www.luogu.org/problemnew/show/P4525)
像这么一坨式子：
$$\int\limits_{L}^{R}\frac{cx+d}{ax+b} \text dx$$
直接积分好像很难，~~反正我是不会~~。
ps：求积分就可以理解为 `在一段区间内，函数图像与x轴围成的面积，不过在x轴上方为正，下方为负`。
按照一般方法来求的话：
$$\int\limits_{L}^{R}f(x) \text dx = F(R)-F(L)$$
算出来 $f(x)$ 的不定积分很可能就不是初等函数，计算难度也很大。
等等，你说什么？用矩形来近似？像这样？
那我们来试试：
```cpp
#include<cstdio>
#include<iostream>
using namespace std;
double a,b,c,d,l,r;
double f(double x){
return (c*x+d)/(a*x+b);
}
double simpson(double a,double b){
double c = a+(b-a)/2;
return (f(a)+4*f(c)+f(b))*(b-a)/6;
}
int main(){
scanf("%lf%lf%lf%lf%lf%lf",&a,&b,&c,&d,&l,&r);
printf("%.6lf",simpson(l,r));
return 0;
}
```
这代码两分钟就能打出来，样例都过了（`flag*1`），我交一发，你看。。。
[评测记录1](https://www.luogu.org/recordnew/show/11403373) `AC*4 WA*5 TLE*1`
~~为什么会变成这样呢？~~
看来用矩形来近似，不仅精度低，而且慢。
这个时候我们就要用 Simpson 积分来近似啦！
Simpson 积分，可以看成是用矩形近似的升级版。
这回我们从曲线上取 $3$ 个点，用二次函数来拟合曲线。
后面推导的式子较多，建议拿出纸和笔，自己跟着算一算。
设原函数为 $f(x)$，用来拟合的函数为 $g(x)$。
$g(x)$ 是二次函数，设为 $g(x)=ax^2+bx+c$
因为 $g(x)$ 是用来拟合 $f(x)$ 的，所以有：
$$\int\limits_L^Rf(x) \text dx\approx \int\limits_L^Rax^2+bx+c\space \text dx$$
$$\int f(x) \text dx \approx \frac13ax^3+\frac12bx^2+cx+C$$
然后再带入 $R$ 和 $L$：
$$\int\limits_L^Rf(x) \text dx \approx \frac13a(R^3-L^3)+\frac12b(R^2-L^2)+c(R-L)$$
然后提公因式，原式为：
$$\frac{R-L}{6}(2a(R^2+LR+L^2)+3b(R+L)+6c)$$
把里面展开来：
$$\frac{R-L}{6}(2aR^2+2aLR+2aL^2+3bR+3bL+6c)$$
重新整理一下式子：
$$\frac{R-L}{6}((aR^2+bR+c)+(aL^2+bL+c)+aR^2+aL^2+2aLR+2bR+2bL+4c)$$
再整理：
$$\frac{R-L}{6}\left((aR^2+bR+c)+(aL^2+bL+c)+4\left(a\left(\frac{R+L}{2}\right)^2+b\left(\frac{R+L}{2}\right)+c\right)\right)$$
为啥我要加一些看似多余的括号呢？
看一下括号里的式子，是不是可以用 $f$ 来替代？
于是就能得到这三个式子：
$$aR^2+bR+c\approx f(R)$$
$$aL^2+bL+c\approx f(L)$$
$$4\left(a\left(\frac{R+L}{2}\right)^2+b\left(\frac{R+L}{2}\right)+c\right)\approx 4f\left(\frac{R+L}{2}\right)$$
把这三个式子带回去
最后我们就得到了 $\text{Simpson}$ 积分公式：
$$\int\limits_L^Rf(x) \text dx \approx \frac{R-L}{6}\left(f(L)+4f\left(\frac{L+R}{2}\right)+f(R)\right)$$
这里给出形象化的图解$\text{qwq}$：
图中红色的曲线就是我们要求积分的函数：
$$f(x) = x^{\frac{2}{x}-x}$$
在这里，我们求它在 $[1,2]$ 上的定积分
怎么办呢？按照我们之前的步骤，直接取 $L=1,R=2$ 套公式就能算出来了。
那近似的精确度怎么样呢？经过计算，误差在 $10^{-2}$ 范围内 。
于是就有了这段代码：
```cpp
#include<cstdio>
#include<iostream>
using namespace std;
double a,b,c,d,l,r;
double f(double x){
return (c*x+d)/(a*x+b);
}
double simpson(double a,double b){
double c = a+(b-a)/2;
return (f(a)+4*f(c)+f(b))*(b-a)/6;
}
int main(){
scanf("%lf%lf%lf%lf%lf%lf",&a,&b,&c,&d,&l,&r);
printf("%.6lf",simpson(l,r));
return 0;
}
```
这次仍然是轻松水过样例，肯定没问题（`flag*2`）
[评测记录2](https://www.luogu.org/record/show?rid=11358070) `AC*5 WA*5`
看来精度还是不够，怎么办呢？
之前我们看到，Simpson 只用了三个点来拟合函数 $f(x)$，显然精度还是不够的，可以考虑分治解决。
对于要积分的区间，再划分成两半，这样就解决精度问题啦！
于是我们就可以这样：
double simpson(double L, double R,double eps){
double mid = (L+R)/2;
if(R-L<eps) return (f(L)+4*f(mid)+f(R))*(R-L)/6;
return simpson(L,mid,eps)+simpson(mid,R,eps);
}
(这里 eps 表示的就是误差范围)
结果呢？
[评测记录3](https://www.luogu.org/recordnew/show/11462902) `AC*9 WA*1`
看来还是可以继续优化：
我们不难注意到，对于函数 $f(x)$，在它比较平缓，变化幅度小的地方，完全可以少分几份，精度也不会差的太多；而对于变化幅度大的地方，分固定的那么多份精度也不够。
所以，如果计算积分时，分治或不分治的计算结果差别不大，那就没必要继续往下分了；否则就继续分。
最后就得到了近似求定积分比较准确的方法：
**自适应** Simpson
double integral(double L, double R, double eps){
//计算[L,R]上的定积分
double mid = (L+R)/2;
double ST = simpson(L,R); //整体
double SL = simpson(L,mid); //左半边
double SR = simpson(mid,R); //右半边
if(fabs(SL+SR-ST) <= 15*eps) return SL+SR+(SL+SR-ST)/15; //误差足够小，直接返回结果
return integral(L, mid, eps/2) + integral(mid, R, eps/2); //继续分
}
有人问 $15$ 那个数是怎么回事，我也是说不太清到底怎么来的。。。
总之就是可以比较准确地调节误差。
这段代码交上去，成功AC
[评测记录4](https://www.luogu.org/record/show?rid=11360443)