自适应辛普森法是一种高效求出已知函数在给定区间的定积分的方法。
Part 1 啥是定积分?
函数 \(f(x)\) 在区间 \([l,r]\) 上的定积分 \(\int_{l}^{r}f(x)\mathrm{d}x\) 指的是 \(f(x)\) 在区间 \([l,r]\) 中与 \(x\) 轴围成的区域的面积(其中 \(x\) 轴上方的部分为正值,\(x\) 轴下方的部分为负值)。
画成图就是这样的:
(这里图先咕了)
不少情况下我们都可以通过推式子+查积分表求解定积分,但计算机并不会推式子,而且我们大多数情况下都只需要求出积分的近似数值,这时候就需要一个通用的数值积分方法来解决问题了。
Part 2 辛普森公式
对于一个二次函数 \(f(x)=Ax^2+Bx+C\),我们有:
$$\int_l^r f(x) {\mathrm d}x = \frac{(r-l)(f(l)+f(r)+4 f(\frac{l+r}{2}))}{6}$$
推导过程较为繁琐,这里暂时略去。
事实上只需稍加改动,辛普森公式也可以求出一般平面图形的近似面积(别忘了积分的实质还是面积)。
(如果你用辛普森公式去求不少常见平面图形面积的话,会发现结果和课堂上学到的面积公式完全一致,这当然不是巧合 :-))
Part 3 自适应辛普森法
拿到了一个函数,如果我们直接带入上面的公式求积分,因为函数图像和二次函数并不相似,误差肯定是很大的。
我们考虑将这个函数分成很多段,每段分别求积分,当一个段足够小,小到和二次函数高度接近的时候,我们带入上面的公式求积分,得到的值就很精确了。
问题在于这里的很多段到底该怎么分呢?分少了计算误差就大,分多了时间效率又会低。我们需要找到一个准确度和效率的平衡点。
我们这样考虑:假如有一段图像已经很接近二次函数的话,直接带入公式求积分,得到的值精度就很高了,不需要再继续分割这一段了。
于是我们有了这样一种分割方法:每次判断当前段和二次函数的相似程度,如果足够相似的话就直接代入公式计算,否则将当前段分割成左右两段递归求解。
现在就剩下一个问题了:如果判断每一段和二次函数是否相似?
我们把当前段直接代入公式求积分,再将当前段从中点分割成两段,把这两段再直接代入公式求积分。如果当前段的积分和分割成两段后的积分之和相差很小的话,就可以认为当前段和二次函数很相似了,不用再递归分割了。
上面描述的过程就被称为自适应辛普森法。代码实现如下:
double simpson(double l,double r) { double mid=(l+r)/2; return (r-l)*(f(l)+4*f(mid)+f(r))/6;//辛普森公式 } double asr(double l,double r,double eqs,double ans) { double mid=(l+r)/2; double fl=simpson(l,mid),fr=simpson(mid,r); if(abs(fl+fr+ans)<=15*eqs)return fl+fr+(fl+fr-ans)/15;//足够相似的话就直接返回 return asr(l,mid,eqs/2,fl)+asr(mid,r,eqs/2,fr);//否则分割成两段递归求解 }
这个代码里事实上有很多重复计算,我们可以将部分结果先存储好,需要调用时直接调用即可。具体实现这里略去。
Part 4 Reference
- xht37 《自适应辛普森法 学习笔记》