Simpson自适应公式[with.Willem]

向Willem爷学习了Simpson自适应公式。

先是做 [NOI2004]降雨量 这题，看到题解里面说可以使用Simpson自适应公式，于是向Willem爷了解了一下。膜膜Willem大爷。

先说积分：积分（integral）的几何意义是函数的曲线上 $x​$ 的一段区间与 $x​$ 轴围成的曲边梯形的面积。

例如该函数的积分记作$$\int_a^b {f(x)} \,{\rm d}x​$$。

如何求这个函数的积分呢？Simpson自适应公式就是其中的一种方法。

Simpson公式就是将复杂的函数 $f(x)​$ 近似拟合成二次函数，再通过二次函数算积分的方法对于$f(x)​$的积分进行估算。下面这个就是Simpson公式。

$$\begin{align} \int_a^b {f(x)} \,{\rm d}x\approx\frac {(b−a)\times(f(a)+f(b)+4 f(\frac{a+b}{2}))} {6} \end{align}$$

但是这个公式算出来的只是一个近似值。题目要求精度怎么办？

自适应就是求积分时能够自动控制切割的区间大小和长度，使精度能满足要求。

来看这行代码：

double asr(double l,double r,double v){
	double mid=(l+r)/2,lv=simpson(l,mid),rv=simpson(mid,r);
	if (fabs(v-(lv+rv))<eps) return lv+rv;
	else return asr(l,mid,lv)+asr(mid,r,rv);
}

通过二分的方法来进行精度的保证——在可控的时间内计算出正确的近似值。

来看洛谷的模板：[模板]自适应辛普森法

代码也很短：

#include<bits/stdc++.h>
#define eps 1e-12

using namespace std;

double a,b,c,d,L,R;

double f(double x){
	return (c*x+d)/(a*x+b);
}

double simpson(double l,double r){
	return (r-l)*(f(r)+f(l)+4*f((l+r)/2))/6;
}

double asr(double l,double r,double v){
	double mid=(l+r)/2,lv=simpson(l,mid),rv=simpson(mid,r);
	if (fabs(v-(lv+rv))<eps) return lv+rv;
	else return asr(l,mid,lv)+asr(mid,r,rv);
}

int main(void){
	scanf("%lf%lf%lf%lf%lf%lf",&a,&b,&c,&d,&L,&R);
	printf("%.6lf\n",asr(L,R,simpson(L,R)));
}

洛谷还有另外一个模板题：[模板]自适应辛普森法2

$\int_0^\infty {f(x)}\,{\rm d}x​$ ($f(x)=x^{\frac{a}{x}-x}​$)

这题没有上界——上界是无穷大。怎么办？

我们来看一下这个$f(x)​$：

不难发现，$f(x)$如果将$x$取大，其函数值均趋向于$0$。

那么这道题的策略就是取一个较大的上界，再用Simpson自适应公式求解。

#include<bits/stdc++.h>
#define eps 1e-10

using namespace std;

double a;

double f(double x){
	return pow(x,a/x-x);
}

double simpson(double l,double r){
	return (r-l)*(f(l)+f(r)+4*f((l+r)/2))/6;
}

double asr(double l,double r,double v){
	double mid=(l+r)/2,lv=simpson(l,mid),rv=simpson(mid,r);
	if (fabs(lv+rv-v)<eps) return lv+rv;
	return asr(l,mid,lv)+asr(mid,r,rv);
}

int main(void){
	scanf("%lf",&a);
	if (a<0){
		puts("orz");
	return 0;}
	printf("%.5f",asr(eps,30,simpson(eps,30)));
}

再来看一道利用Simpson自适应公式的基础题：[HDU]Ellipse

大致题目如下：

Math is important!! Many students failed in 2+2’s mathematical test, so let’s AC this problem to mourn for our lost youth..
Look this sample picture:

A ellipses in the plane and center in point O. the L,R lines will be vertical through the X-axis. The problem is calculating the blue intersection area. But calculating the intersection area is dull, so I have turn to you, a talent of programmer. Your task is tell me the result of calculations.(defined PI=3.14159265 , The area of an ellipse A=PIab )

利用椭圆公式，就可以写出要积分的$f(x)$。

#include<bits/stdc++.h>
#define eps 1e-8

using namespace std;

int N;
double a,b,L,R;

double f(double x){
	return sqrt((a*a*b*b-b*b*x*x)/(a*a))*2;
}

double simpson(double l,double r){
	return (r-l)*(f(l)+f(r)+4*f((l+r)/2))/6;
}

double asr(double l,double r,double v){
	double mid=(l+r)/2,lv=simpson(l,mid),rv=simpson(mid,r);
	if (fabs(lv+rv-v)<eps) return lv+rv;
	else return asr(l,mid,lv)+asr(mid,r,rv);
}

int main(void){
	scanf("%d",&N);
	while (N--){
		scanf("%lf%lf%lf%lf",&a,&b,&L,&R);
		printf("%.3f\n",asr(L,R,simpson(L,R)));
	}
}

这些大致就是Simpson自适应公式的基础内容，感谢Willem大爷的谆谆教诲。

感谢！