辛普森积分入门讲解

辛普森积分

引用部分

定义

辛普森积分法就是在积分区间[a,b]上去找三个点a、b和m=(a+b)/2,计算其原函数的在此处的值,然后用抛物线来拟合原函数。

正文

Simpson积分公式

用途:来求一个函数的积分的近似值，用于面积计算等精度要求不是特别苛刻的地方。

其实它就是用一个二次函数曲线不断拟合逼近原函数，然后求得原函数的近似值。

公式说明:

前置： $g (x)$ 为一个关于 $x$ 的二次函数(抛物线)，其中 $g (x) = A \times x^{2} + B \times x + C$ ，对于求定积分 $\int_{0}^{x} g (x) d x$ ，通过求积得其等于 $\frac{A \times x^{3}}{3} + \frac{B \times x^{2}}{2} + C \times x + D$ 其中 $D$ 为常数，可以看做 $0$ 。令 $W (x) = \int_{0}^{x} g (x) d x$ ，所以对于求一段定积分则有 $\int_{a}^{b} g (x) = W (b) - W (a)$ 。

在平面直角坐标系里，由 $(x_{1}, y_{1}), (x_{2}, y_{2}), (x_{3}, y_{3})$ (其中 $x_{3} = \frac{x_{1} + x_{2}}{2}$ )确定的抛物线 $f (x)$ 在区间[x1,x2]的定积分为:

\int_{x_{1}}^{x_{2}} f (x) d x = \frac{1}{6} \times (x_{2} - x_{1}) \times (y_{1} + y_{2} + 4 \times y_{3})

下面给出简单的证明:

令 $g (x) = A \times x^{2} + B \times x + C$ 为拟合后的抛物线，则有

\int_{x_{1}}^{x_{2}} f (x) d x \approx \int_{x_{1}}^{x_{2}} g (x) d x

= W (x_{2}) - W (x_{1})

= \frac{A}{3} \times (x_{2})^{3} + \frac{B}{2} \times (x_{2})^{2} + C \times x_{2} - (\frac{A}{3} \times (x_{1})^{3} + \frac{B}{2} \times (x_{1})^{2} + C \times x_{1})

= \frac{A}{3} \times ((x_{2})^{3} - (x_{1})^{3}) + \frac{B}{2} \times ((x_{2})^{2} - (x_{1})^{2}) + C \times (x_{2} - x_{1})

= \frac{x_{2} - x_{1}}{6} \times (2 \times A \times ((x_{2})^{2} + x_{1} \times x_{2} + (x_{1})^{2}) + 3 \times B \times (x_{2} + x 1) + 6 \times C)

展开化简整理得：

= \frac{x_{2} - x_{1}}{6} \times (A \times (x_{1})^{2} + B \times x_{1} + C + A \times (x_{2})^{2} + B \times x_{2} + C + 4 \times A \times {(\frac{x_{2} + x_{1}}{2})}^{2})

将其组合成完全平方式(配方)后

= \frac{x_{2} - x_{1}}{6} \times (g (x_{1}) + g (x_{2}) + 4 \times g (\frac{x_{1} + x_{2}}{2}))

= \frac{x_{2} - x_{1}}{6} \times (g (x_{1}) + g (x_{2}) + 4 \times g (x_{3}))

于是我们就得到了simpson积分公式

\int_{a}^{b} f (x) d x \approx \frac{b - a}{6} \times [g (a) + 4 \times g (\frac{a + b}{2}) + g (b)]

在实际计算中 $g (x)$ 的值可以用原函数 $f (x)$ 的值来代替，于是就是如下公式:

\int_{a}^{b} f (x) d x \approx \frac{b - a}{6} \times [f (a) + 4 \times f (\frac{a + b}{2}) + f (b)]

代码：

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

自适应辛普森积分法

那么实际程序该如何实现辛普森积分求积呢?

我们如果要求 $\int_{a}^{b} f (x) d x$ 的近似值的话，可以用递归二分区间求解来达到要求精度。
用如下公式:

\int_{a}^{b} f (x) d x = \int_{a}^{m i d} f (x) d x + \int_{m i d}^{b} f (x) d x

其中

m i d = \frac{a + b}{2}

，证明:显然式证明 :)

但是因为是浮点数(小数),那么递归多少层，在什么时候返回值结束递归呢?
我们容易知道如果递归到 $b - a < e p s$ 的话精度虽然很高，但是时间复杂度太高了，但是如果递归少了，精度又得不到保证，那该如何是好呢?

自适应法

自适应法，就是让程序根据实际情况决定如何运行执行操作。~~自己随便下的定义而已~~

这里我们就要用自适应法来解决这个问题啦，让程序自己去决定递归层数，而且又保证精度。

说的很高深，其实很简单。~~还是比较难吧~~

自动化控制区间分割的大小。

实际操作:二分递归，当满足精度就计算返回值，结束递归。

伪代码:

function(l,r,eps,ans):
mid=(l+r)/2;
lval=左边的值,rval=右边的值;
if (满足精度) return 答案;
eps/=2;
else return 左边递归+右边递归;

注意，这里的 $a n s$ 表示上一层计算的整个区间的答案，用来和当前这层来判断精度， $e p s$ 在递归时每次除以2，这是为了消除精度误差叠加效应，当小误差多了就成大误差了，所以每次要缩小精度。

代码:

double asr(double l,double r,double eps,double ans){
    double mid=(l+r)/2;
    double lval=simpson(l,mid),rval=simpson(mid,r);
    if(fabs(lval+rval-ans)<=15*eps) return lval+rval+(lval+rval-ans)/15;
    return asr(l,mid,eps/2,lval)+asr(mid,r,eps/2,rval);
}
double asme(double a,double b,double eps){
    return asr(a,b,eps,simpson(a,b));
}

推荐文章

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

代码

#include<cmath>
#include<cstdio>
#include<cstring>
#include<algorithm>
#define db double
using namespace std;
const db eps=1e-7;
db a,b,c,d,L,R;
db f(db x){return (c*x+d)/(a*x+b);}
db simpson(db l,db r){return (f(l)+f(r)+4*f((l+r)/2))*(r-l)/6;}
db asr(db l,db r,db exps,db val){
    db mid=(l+r)/2;
    db lval=simpson(l,mid),rval=simpson(mid,r);
    if(fabs(lval+rval-val)<=15*exps){return lval+rval+(lval+rval-val)/15;}
    return asr(l,mid,exps/2,lval)+asr(mid,r,exps/2,rval);
}
db asme(db l,db r,db exps){return asr(l,r,exps,simpson(l,r));}
int main(){
    scanf("%lf%lf%lf%lf%lf%lf",&a,&b,&c,&d,&L,&R);
    printf("%lf
",asme(L,R,eps));
    return 0;
}

模板2【模板】自适应辛普森法2

代码


#include<cmath>
#include<cstdio>
#include<cstring>
#include<algorithm>
#define db double
using namespace std;
const db inf=30;
const db eps=1e-7,zero=1e-10;
db a;
db f(db x){return pow(x,a/x-x);}
db simpson(db l,db r){return (f(l)+f(r)+4*f((l+r)/2))*(r-l)/6;}
db asr(db l,db r,db exps,db val){
    db mid=(l+r)/2;
    db lval=simpson(l,mid),rval=simpson(mid,r);
    if(fabs(lval+rval-val)<=15*exps) return lval+rval+(lval+rval-val)/15;
    return asr(l,mid,exps/2,lval)+asr(mid,r,exps/2,rval);
}
db asme(db l,db r,db exps){return asr(l,r,exps,simpson(l,r));}
int main(){
    scanf("%lf",&a);
    if(a<0)puts("orz");
    else printf("%.5lf
",asme(zero,inf,eps));
    return 0;
}

这个虽然求的是不定积分但是，不要被吓到了，因为当 $x$ 大于30左右后，函数值趋近于0，所以可以不计。
然后当 $a < 0$ 时函数不收敛，所以无解。

其他题目[NOI2005]月下柠檬树

simpson的其他用途:

和扫描线结合求圆面积并和其他不规则图形面积等。

【2018.9.7】最近发现的好文章IN