【学习笔记】辛普森法

辛普森法

背景：计算定积分

\[\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) \]

\(c\) 为 \([l,r]\) 中的某个值。如果 \(f(x)\) 在区间上的四阶导数尽可能小，或者区间长度尽可能短，那么这个结果也就越精确（四次以下的多项式函数完全没有误差）。

关于这个误差分析我还没找到 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)\) ，所如果我们一开始就将区间划得足够小，那么这个假设就可以成立。