【自适应辛普森·入门篇】积分问题的巧妙解法

\(Preface\)

自适应辛普森，一个非常有趣的算法。

感觉也不需要多么高深的数学知识，果然算法这种东西实践起来总比想象中容易许多。

辛普森积分

考虑对于任何一个函数图象，在\(\Delta x\)极小的时候，都可以将这一范围内的函数视为一段抛物线。

假设这段抛物线拟合后的函数为\(Ax^2+Bx+C\)，然后推导一下：

\[\begin{align}\int_l^r(Ax^2+Bx+C)&=\frac A3(r^3-l^3)+\frac B3(r^2-l^2)+\frac C3(r-l)\\&=\frac{r-l}6(2A(l^2+lr+r^2)+3B(l+r)+6C)\\&=\frac{r-l}6((Al^2+Bl+C)+4(A(\frac{l+r}2)^2+B(\frac{l+r}2)+C)+(Ar^2+Br+C))\\&=\frac{r-l}6(f(l)+4f(\frac{l+r}2)+f(r))\end{align} \]

这就是所谓辛普森积分的公式了：

\[\int_l^rf(x)dx\approx\frac{r-l}6(f(l)+4f(\frac{l+r}{2})+f(r)) \]

考虑到这一公式的前提是\(\Delta x\)极小，因此我们需要不断二分递归以保证精度，不难发现复杂度显然爆炸。

所以我们需要改进这一算法。

自适应辛普森

递归次数少，难以保持精度；递归次数多，复杂度又难以接受。

而所谓的自适应辛普森算法，就是在二分的过程中根据精度，自动控制区间分割的大小。

可以直接看代码。

/*友情提醒：I表示inline，可忽略；DB表示double；Con DB&表示const double&，可以直接视作double*/
namespace SimpsonInt
{
	I DB f(Con DB& x) {return /*需要积分的函数*/;}
	I DB Simpson(Con DB& l,Con DB& r) {return (r-l)*(f(l)+4*f((l+r)/2)+f(r))/6;}//辛普森积分
	I DB ASR(Con DB& l,Con DB& r,Con DB& ans,Con DB& eps)//自适应辛普森算法
	{
		DB mid=(l+r)/2,L=Simpson(l,mid),R=Simpson(mid,r);
		if(fabs(L+R-ans)<=15*eps) return L+R+(L+R-ans)/15;//控制精度（15据说是研究结果）
		return ASR(l,mid,L,eps/2)+ASR(mid,r,R,eps/2);//二分递归，注意eps除以2
	}
}

模板：自适应辛普森法1

• 求\(\int_L^R\frac{cx+d}{ax+b}dx\)。
• \(-100\le L<R\le 100\)且\(R-L\ge 1\)。

真正的模板题。

只要在上面给出代码中的函数\(f(x)\)里填上(c*x+d)/(a*x+b)即可。

模板：自适应辛普森法2

• 求\(\int_0^{\infty}x^{\frac ax-x}dx\)。
• 若积分发散，输出\(orz\)
• \(a\le50\)。

可证明\(a<0\)时积分发散，\(a>0\)时积分收敛。（具体证明这里略去）

这道题看似是要求无限积分，实际上当\(x\)到\(15\)左右\(x^{\frac ax-x}\)就已经小得可以忽略不计了。

所以说，我们只要求\(\int_0^{15}x^{\frac ax-x}dx\)即可。

最后要做的就是在代码中的函数\(f(x)\)里填上pow(x,a/x-x)。

\(Postscript\)

自适应辛普森算法是一个简单而又有用的算法。

具体来说，就是学了并不会造成多大负担，而万一碰到了这种题目又能发挥出巨大的作用。

所以，余认为，这还是非常值得一学的。