辛普森积分 - 爱码网

这两天看了看辛普森积分，觉得这是很好的骗分工具，正好也比较简单，就随便写一写。这个东西对于三次以下的函数是正确的，但是对于三次以上的函数我们可以将其近似为低次函数套用Simpson公式，这个东西学名好像叫自适应Simpson积分。

先附上公式辛普森积分。

算法的大概流程就是不断的递归求对应区间的积分，当当前层和下一层的相对误差比较小时就退出，从而得到比较准确的值，但是这个其实是可以被卡掉的，比如一段波动的区间，我们取到的值有可能都是上顶点，就会直接退出，然而这并不是最后的答案。所以一般这个算法都是用来骗分的，2333。

附上代码

 1 double calc(double l,double r,double fl,double fr,double fm){
 2     return (r-l)/6.0*(fl+fr+4*fm);
 3 }
 4 double simpson(double l,double r,double mid,double fl,double fr,double fm,double s){
 5     double m1=(l+mid)/2.0,m2=(mid+r)/2.0;
 6     double f1=F(m1),f2=F(m2);
 7     double s1=calc(l,mid,fl,fm,f1),s2=calc(mid,r,fm,fr,f2);
 8     if(fabs(s-s1-s2)<eps)return s1+s2;
 9     return simpson(l,mid,m1,fl,fm,f1,s1)+simpson(mid,r,m2,fm,fr,f2,s2);
10 }

我一般在simpson函数里会将三个函数值一并传进去，否则在单次计算f的复杂度比较大时，会重复计算，浪费时间。

hdu1724，纯板子题。

这道题我写的还是不传函数值的版本，可以对照着看看

 1 #include <cstdio>
 2 #include <cstring>
 3 #include <iostream>
 4 #include <algorithm>
 5 #include <cmath>
 6 #define eps 1e-10
 7 using namespace std;
 8 double a,b;
 9 int T;
10 double f(double x){
11     return b*sqrt(1-(x*x)/(a*a));
12 }
13 double calc(double l,double r){
14     double mid=(l+r)/2.0;
15     return (r-l)/6*(f(l)+f(r)+4*f(mid));
16 }
17 double simpson(double l,double r,double now){
18     double mid=(l+r)/2.0;
19     double n1=calc(l,mid),n2=calc(mid,r);
20     if(fabs(now-n1-n2)<eps)return now;
21     return simpson(l,mid,n1)+simpson(mid,r,n2);
22 }
23 int main(){
24     scanf("%d",&T);
25     double l,r,ans;
26     while(T--){
27         scanf("%lf%lf%lf%lf",&a,&b,&l,&r);
28         ans=simpson(l,r,calc(l,r))*2;
29         printf("%0.3f\n",ans);
30     }
31 }

View Code