自适应Simpson法P4525 【模板】自适应辛普森法1

何为自适应Simpson法

自适应辛普森算法(Adaptive Simpson′s rule)是一类近似算法(Approximation algorithm)，主要用于在信息计算时求较难反导的函数的定积分。其思想是利用二次函数曲线来不断**拟合(Overfitting)所求曲线（显然这传统的直接用矩形或直边梯形作为微元更精准），而所谓的Adapative（自适应）**则是用于优化时间复杂度的方法。 ----wikipedia

在计算几何里，显然可在特殊情况下利用自适应Simpson法求图形的面积。

Simpon公式及推导

Simpson公式

∫abf(x)dx≈(b−a)(f(a)+f(b)+4f(a+b2))6\int_a^bf(x)dx\approx\frac{(b-a)(f(a)+f(b)+4f(\frac{a+b}{2}))}{6}∫abf(x)dx≈6(b−a)(f(a)+f(b)+4f(2a+b))

推导

设f(x)f(x)f(x)为被积函数，g(x)=Ax2+Bx+Cg(x)=Ax^2+Bx+Cg(x)=Ax2+Bx+C为用于拟合f(x)f(x)f(x)的函数，即g(x)≈f(x)g(x)\approx f(x)g(x)≈f(x)，则有：

∫abf(x)dx≈∫abg(x)dx=∫ab(Ax2+Bx+C)dx=A3(b3−a3)+B2(b2−a2)+C(b−a)\int_{a}^{b}f(x)dx\approx \int_{a}^{b}g(x)dx=\int_{a}^{b}(Ax^2+Bx+C)dx=\frac{A}{3}(b^3-a^3)+\frac{B}{2}(b^2-a^2)+C(b-a)∫abf(x)dx≈∫abg(x)dx=∫ab(Ax2+Bx+C)dx=3A(b3−a3)+2B(b2−a2)+C(b−a)

=A3(b−a)(b2+ab+a2)+B2(b−a)(b+a)+C(b−a)=(b−a)[A3(b2+ab+a2)+B2(b+a)+C]=\frac{A}{3}(b-a)(b^2+ab+a^2)+\frac{B}{2}(b-a)(b+a)+C(b-a)=(b-a)[\frac{A}{3}(b^2+ab+a^2)+\frac{B}{2}(b+a)+C]=3A(b−a)(b2+ab+a2)+2B(b−a)(b+a)+C(b−a)=(b−a)[3A(b2+ab+a2)+2B(b+a)+C]

=(b−a)6[2A(b2+ab+a2)+3B(b+a)+6C]=\frac{(b-a)}{6}[2A(b^2+ab+a^2)+3B(b+a)+6C]=6(b−a)[2A(b2+ab+a2)+3B(b+a)+6C]

=(b−a)6[Aa2+Ba+C+Ab2+Bb+C+Aa2+Ab2+2Aab+2Ba+2Bb+4C]=\frac{(b-a)}{6}[Aa^2+Ba+C+Ab^2+Bb+C+Aa^2+Ab^2+2Aab+2Ba+2Bb+4C]=6(b−a)[Aa2+Ba+C+Ab2+Bb+C+Aa2+Ab2+2Aab+2Ba+2Bb+4C]

=(b−a)6[f(a)+f(b)+A(a+b)2+2B(a+b)+4C]=\frac{(b-a)}{6}[f(a)+f(b)+A(a+b)^2+2B(a+b)+4C]=6(b−a)[f(a)+f(b)+A(a+b)2+2B(a+b)+4C]

=(b−a)6[f(a)+f(b)+4A(a+b2)2+4B(a+b2)+4C]=\frac{(b-a)}{6}[f(a)+f(b)+4A(\frac{a+b}{2})^2+4B(\frac{a+b}{2})+4C]=6(b−a)[f(a)+f(b)+4A(2a+b)2+4B(2a+b)+4C]

=(b−a)(f(a)+f(b)+4f(a+b2))6=\frac{(b-a)(f(a)+f(b)+4f(\frac{a+b}{2}))}{6}=6(b−a)(f(a)+f(b)+4f(2a+b))

上述推导过程除了求了个一元二次多项式的积分外，运用的都是初等数学的知识。

除此推导方法之外，此公式还有一种更常规的推导方法，那就是利用牛顿-柯特斯求积公式。有篇相关文章写的不错，链接如下：

https://blog.csdn.net/qq_43069547/article/details/83153673

代码实现—自适应

如果想要利用Simpson公式求定积分，显然我们要把积分区间分成多个足够小的区间，但问题就在于如何去划分才能既精确效率又高。

其实，Simpson法的程序实现很简单，就是二分递归的过程，但关键在于“自适应”。所谓自适应，说的直白点，无非就是需要多分治几次的地方，多分治几次；不需要的则可以少分治几次。

那么，我们可以得出自适应Simpson法的核心实现过程：在二分递归划分区间的过程中，如果满足了精度需要，则可以终止当前递归，而这就是自适应辛普森法能够自动控制区间分割大小的手段。

inline double queryans(double l,double r,double eps,double ans)//ans：利用Simpson法得到的以[l,r]为积分区间的f(x)的定积分； eps:精度
{register double mid=(l+r)/2.0;register double ansl=Simpson(l,mid),ansr=Simpson(mid,r);if(fabs(ansl+ansr-ans)<=15.0*eps)//Problem 1    return ansl+ansr+(ansl+ansr-ans)/15.0; return queryans(l,mid,eps/2.0,ansl)+queryans(mid,r,eps/2.0,ansr);//Problem 2
}

上述代码有两点需要说明的地方：

Problem 1：为什么判断是否还需要对当前区间继续划分拟合时的判断条件是误差≤15\le15≤15倍的精度值？

答：据下面链接的文章说，wekipedia里有相关证明，同时贴了张图：

Problem 2：为什么递归求折半积分区间的积分时精度？

如果要求当前区间UUU计算的返回值误差eps(U)<keps(U)<keps(U)<k的话，显然其折半区间L,RL,RL,R的返回值的误差应有eps(L)≤k2,eps(R)≤k2eps(L)\le \frac{k}{2},eps(R)\le \frac{k}{2}eps(L)≤2k,eps(R)≤2k这些不等关系，只有这样才能保证eps(U)=eps(L)+eps(R)≤k2+k2=keps(U)=eps(L)+eps(R)\le \frac{k}{2}+\frac{k}{2}=keps(U)=eps(L)+eps(R)≤2k+2k=k。

最后，需要说明的是，本文大部分内容学习或借鉴自https://www.cnblogs.com/pks-t/p/10277958.html，这篇文章讲的确实不错，许多其他相关文章都没提到的令人疑惑的地方点的较到位。

Code

#include<cstdio>
#include<cmath>
#include<iostream>
using namespace std;double A,B,C,D,L,R;inline double f(double x) { return (C*x+D)/(A*x+B); }
inline double Simpson(double a,double b) { return (f(a)+4*f((a+b)/2.0)+f(b))*(b-a)/6.0; }inline double queryans(double l,double r,double eps,double ans)
{register double mid=(l+r)/2.0;register double ansl=Simpson(l,mid),ansr=Simpson(mid,r);if(fabs(ansl+ansr-ans)<=15.0*eps)   return ansl+ansr+(ansl+ansr-ans)/15.0;//15.0是科学推算得 return queryans(l,mid,eps/2.0,ansl)+queryans(mid,r,eps/2.0,ansr);//精度别忘除以2（列式等价得）
}int main()
{scanf("%lf%lf%lf%lf%lf%lf",&A,&B,&C,&D,&L,&R);printf("%.6lf",queryans(L,R,1e-6,Simpson(L,R)));return 0;
}