自适应辛普森积分算法

辛普森积分是数值积分的一种，是中点公式和梯形公式的改进。

——这是我在APIO2017时听某位大神讲课所了解的。

下面就讲讲辛普森积分到底是啥玩意儿。

假定我们要求如下定积分

$$\int_a^bf(x)dx$$
略微懂一点微积分知识的都知道，对于一个黎曼可积的函数，我们要求其在某个闭区间上的定积分，要先求该函数的不定积分，即先求原函数。就是找到一个函数 $F(x)$，使得 $F'(x)=f(x)$，然后根据牛顿—莱布尼茨公式去搞。即 $$\int_a^bf(x)dx=F(b)-F(a)$$
而有时候原函数并不好求，比如要求原的函数长得很复杂，要求的原函数非初等函数。。
这时我们需要“避其锋芒”，采用一种利用“近似”思想的方法。

我们都知道定积分就是曲线$f(x)$下方的面积(假设$f(x)$恒大于0)，我们考虑一种近似求出曲线下方面积的方法——辛普森法。

我们先在区间$[a,b]$中等距地取奇数个点$a=x_0,x_1,x_2,...,x_n=b$，相邻两点间的距离为$\Delta x$，再设$y_i=f(x_i)$，然后套用公式

$$\int_a^bf(x)dx \approx{\Delta x\over3}(y_0+4y_1+y_2)+{\Delta x\over3}(y_2+4y_3+y_4)+...+{\Delta x\over3}(y_{n-2}+4y_{n-1}+y_n)$$
点取得越多，近似程度就越高，但计算量也越大。当只取三个点的时候，有
$$\int_a^bf(x)dx\approx{ b-a\over6}[f(a)+4f({a+b\over2})+f(b)]$$
这就是三点辛普森公式，又简称辛普森公式。

然而如果只取三个点误差肯定大，取太多的点计算量也会上去。那么到底取多少个点合适呢？

还好辛普森积分有一个重要的“变种”，称为自适应辛普森法(adaptive Simpson’s rule)，我们设置一个精度Eps，让算法根据具体情况递归的划分区间，容易近似的地方少划几份，不容易近似的地方多划几份。这就是所谓的“自适应”。

具体的话，就是设区间$[a,b]$的中点为$c$，则当且仅当

$$|S(a,c)+S(c,b)-S(a,b)|<15Eps$$ 时(加$=$亦可)直接返回结果，否则递归调用，再次划分区间。递归调用时精度Eps也要相应地减小一半。
ps：
$$\int_a^bf(x)dx=\int_a^cf(x)dx+\int_c^bf(x)dx\ (a<c<b)$$
这里的$S(a,b)$就是$[a,b]$的三点辛普森公式值，返回的结果是$S(c,b)+S(a,c)+\Delta s$，$Δs=[S(a,c)+S(c,b)-S(a,b)]/15|$。顺便提醒一句，如果直接返回$S(a,b)$的话，可能存在的误差超出可接受范围，就是比较可能会出错。。还有，这个$15$是怎么来的别问我。。。

以上就是我对自适应辛普森积分的初步理解，或者可以说是最初理解。本人对关于微积分的知识比较有兴趣，今天初学了辛普森积分又有了新的收获，然而这样的日子恐怕已不多。

---

Code

double f(double x){
    return x * x + x;//写要求辛普森积分的函数
}

double simpson(double L, double R){//三点辛普森积分法，要求f(x)是全局函数
    double mid = (L + R) / 2.0;
    return (f(L) + 4.0 * f(mid) + f(R)) * (R - L) / 6.0;
}

double integral(double L, double R, double Eps){//自适应辛普森积分递归过程
    double mid = (L + R) / 2.0;
    double ST = simpson(L, R), SL = simpson(L, mid), SR = simpson(mid, R);
    if(fabs(SL + SR - ST) <= 15.0 * Eps)  return SL + SR + (SL + SR - ST) / 15.0;//直接返回结果
    return integral(L, mid, Eps/2.0) + integral(mid, R, Eps/2.0);//对半划分区间
}

---

我依然只站在你那一边
我明白那份孤独伤痛