【算法】自适应辛普森(Simpson)算法详解

先来看一道题：

给你三条直线(x=a,x=b,y=0)和一个函数

求该函数图像与这三条直线围成的面积。

【输入】

第一行n,a,b。

第二行n+1个整数A0~An。

【输出】

一行一个整数，即这四条先围成的面积，保留两位小数。

【样例输入1】

2 0 1

1 0 0

【样例输出1】

0.33

【样例输入2】

3 -1 0

1 0 0 0

【样例输出2】

-0.25

【备注】

数据范围：1≤n≤7，-10≤a≤b≤10，-10≤Ai≤10

注意：这里指的“面积”是代数面积，即x轴上方的面积是正的，x轴下方的面积是负的。

来看看样例一

通过一个图，我们了解到，题意就是求一个函数图象在[a,b]的区间里下方的面积。

于是我们想到把这个不规则图形剖分成近似的规则图形

即

把每一段的面积加起来就可以得到答案了。

可是，要划分成多少段才好呢？太多了时间慢，太少了答案不准确，因此，自适应辛普森法便权衡了这两点，在容易近似的地方少划分一点，在不容易近似的地方多划分，这样便既快又准了。

代码实现：

double simpson(double a,double b)//simpson，cal是全局函数，用来计算f(x)的值 {double c=(a+b)/2;return (cal(a)+4*cal(c)+cal(b))*(b-a)/6;}double asr(double a,double b,double eps,double A)//自适应simpson递归过程 {double c=(a+b)/2;double L=simpson(a,c),R=simpson(c,b);if(fabs(L+R-A)<=15*eps)return L+R+(L+R-A)/15.0;return asr(a,c,eps/2,L)+asr(c,b,eps/2,R);}double asr(double a,double b,double eps)//主过程,精度要求为eps,区间为[a,b] {return asr(a,b,eps,simpson(a,b));}

于是我们便得到了本题的代码：

#include <iostream>#include <cstdio>#include <cmath>#include <algorithm>using namespace std;const double eps=1e-5;int xi[10],n,a,b;double ans=0;double cal(double x){double re=0;for(int i=n;i>=0;--i){re=re*x+xi[i];}return re;}double simpson(double a,double b){double c=(a+b)/2;return (cal(a)+4*cal(c)+cal(b))*(b-a)/6;}double asr(double a,double b,double eps,double A){double c=(a+b)/2;double L=simpson(a,c),R=simpson(c,b);if(fabs(L+R-A)<=15*eps)return L+R+(L+R-A)/15.0;return asr(a,c,eps/2,L)+asr(c,b,eps/2,R);}double asr(double a,double b,double eps){return asr(a,b,eps,simpson(a,b));}int main(){scanf("%d%d%d",&n,&a,&b);if(b<a)swap(a,b);for(int i=n;i>=0;--i){scanf("%d",xi+i);}printf("%.2lf\n",asr(a,b,0.01));}

并且通过了这道题。

【完】