【学习笔记】辛普森法

· · 个人记录

背景：计算定积分

\int_{l}^{r}f(x)dx

应用场景：数值积分，计算特殊平面图形面积

注：本文主要讲解的是自适应辛普森法的各种边界条件推导。

前置知识

辛普森公式

对于二次函数 f(x)=ax^2+bx+c 有公式

\int_{l}^{r}f(x)dx=\frac{(r-l)(f(l)+f(r)+4f(\frac{l+r}{2}))}{6}=S[l,r]

这个等式证明非常简单，可以自己试试看（提示：拉格朗日插值法）。

而对于一个非二次的区间可积函数我们可以用一个抛物线去拟合它。

记 h=r-l，对于一般的 [l,r] 可积函数 f(x) 有

\int_{l}^rf(x)dx=S[l,r]-O(h^5)

其中余项（也就是误差）为

-\frac{h^5}{2880}f^{(4)}(c) 关于这个误差分析我还没找到 reference，先 mark 住。 ### 普通辛普森法 我们可以将区间划分为 $n$ 段，对每一段分别进行辛普森。 由于对于三次以下多项式函数辛普森公式都能完美拟合，所以显然当 $n$ 越大误差越小。 普通辛普森法可以有效提高数值积分的精度，但有时间上的短板。 ### 自适应辛普森法 自适应辛普森法可以根据实际情况选择是否将区间继续拆分细化。 那要怎么判断呢？ 我们每次先将区间一分为二，如果两个小区间的和与真实值的误差小于 eps，那么就返回两个区间的和，否则继续递归拆分。 不过我们并不知道当前的误差大小，所以我们要通过数学方法估计。 对于区间 $[l,r]$，我们对其使用辛普森公式： $$ S=\int_{l}^{r}f(x)dx=S[l,r]-\frac{h^5}{2880}f^{(4)}(c) $$ 我们拆分成左右区间分别套用辛普森公式： $$ S_l=\int_{l}^{\frac{l+r}{2}}f(x)dx=S[l,\frac{l+r}{2}]-\frac{1}{32}\cdot\frac{h^5}{2880}f^{(4)}(c_1)\\ S_r=\int_{\frac{l+r}{2}}^{r}f(x)dx=S[\frac{l+r}{2},r]-\frac{1}{32}\cdot\frac{h^5}{2880}f^{(4)}(c_2)\\ $$ 此时我们假定 $f^{(4)}(c)\approx f^{(4)}(c_1) \approx f^{(4)}(c_2)$ 并记 $\rho=\frac{h^5}{2880}f^{(4)}(c)

我们作差得

\begin{aligned} S-(S_l+S_r)&\approx\int_{l}^{r}f(x)dx-\int_{l}^{\frac{l+r}{2}}f(x)dx-\int_{\frac{l+r}{2}}^{r}f(x)dx\\ &+\rho-\frac{\rho}{16}\\ &=\frac{15}{16}\rho \end{aligned}

也就是说整个区间计算的结果与拆分后两段之和的差，大约为两端之和与真实值误差的 15 倍。

因此如果有

|S-(S_l+S_r)|<15\times EPS

就证明了 S_l+S_r 与真实值误差不超过 eps，无需继续拆分。

对于返回的结果，我们需要作以下修正：

\begin{aligned} \int_{l}^{r}f(x)dx&\approx S_l+S_r-\frac{\rho}{16}\\ &=S_l+S_r+(S_l+S_r-S)/15 \end{aligned}

这样可以更接近真实值。

代码：

double simpson(double l,double r){
    double mid=(l+r)/2.0;
    return (r-l)*(f(l)+f(r)+4*f(mid))/6.0;
}

double asr(double l,double r,double eps,double now,int step){
    double mid=(l+r)/2;
    double fl=simpson(l,mid),fr=simpson(mid,r);
    if(fabs(fl+fr-now)<=eps*15&&step<0){
        return fl+fr+(fl+fr-now)/15.0;
    }
    return asr(l,mid,eps/2.0,fl,step-1)+asr(mid,r,eps/2.0,fr,step-1);
}

double calc(double l,double r,double eps){
    return asr(l,r,eps,simpson(l,r),12);
}

注意到我们还维护了一个 step 表示最小拆分次数。这是因为上述推导中有一个关键的假设 f^{(4)}(c)\approx f^{(4)}(c_1) \approx f^{(4)}(c_2) ，所如果我们一开始就将区间划得足够小，那么这个假设就可以成立。