辛普森积分入门讲解

辛普森积分

引用部分

• 百度一下

• 文库

定义

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

正文

• Simpson积分公式

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

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

---

• 公式说明:

---

前置：$g(x)$为一个关于$x$的二次函数(抛物线)，其中$g(x)=A\times x^2+B\times x+C$，对于求定积分$\int_0^xg(x)dx$，通过求积得其等于$\frac{A\times x^3}{3}+\frac{B\times x^2}{2}+C\times x+D$ 其中$D$为常数，可以看做$0$。令$W(x)=\int_0^xg(x)dx$，所以对于求一段定积分则有$\int_a^bg(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)dx=\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)dx≈\int_{x_1}^{x_2}g(x)dx$$
$$=W(x_2)-W(x_1)$$
$$=\frac{A}{3}\times (x_2)^3+\frac{B}{2}\times (x_2)^2+C\times x_2-\left(\frac{A}{3}\times (x_1)^3+\frac{B}{2}\times (x_1)^2+C\times x_1\right)$$
$$=\frac{A}{3}\times \left((x_2)^3-(x_1)^3\right)+\frac{B}{2}\times \left((x_2)^2-(x_1)^2\right)+C\times (x_2-x_1)$$
$$=\frac{x_2-x_1}{6}\times \left(2\times A\times \left((x_2)^2+x_1\times x_2+(x_1)^2\right)+3\times B\times (x_2+x1)+6\times C\right)$$
展开化简整理得：
$$=\frac{x_2-x_1}{6}\times \left(A\times (x_1)^2+B\times x_1+C+A\times (x_2)^2+B\times x_2+C+4\times A\times \left(\frac{x_2+x_1}{2}\right)^2\right)$$
将其组合成完全平方式(配方)后
$$=\frac{x_2-x_1}{6}\times \left(g(x_1)+g(x_2)+4\times g\left(\frac{x_1+x_2}{2}\right)\right)$$
$$=\frac{x_2-x_1}{6}\times \left(g(x_1)+g(x_2)+4\times g(x_3)\right)$$
于是我们就得到了simpson积分公式
$$\int_a^bf(x)dx≈\frac{b-a}{6}\times \left[g(a)+4\times g\left(\frac{a+b}{2}\right)+g(b)\right]$$

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

$$\int_a^bf(x)dx≈\frac{b-a}{6}\times \left[f(a)+4\times f\left(\frac{a+b}{2}\right)+f(b)\right]$$

代码：

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

---

自适应辛普森积分法

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

我们如果要求$\int_a^bf(x)dx$的近似值的话，可以用递归二分区间求解来达到要求精度。
用如下公式:

$$\int_a^bf(x)dx=\int_a^{mid}f(x)dx+\int_{mid}^bf(x)dx$$
其中 $mid=\frac{a+b}{2}$，证明:显然式证明 :)

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

• 自适应法

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

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

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

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

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

伪代码:

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

注意，这里的$ans$表示上一层计算的整个区间的答案，用来和当前这层来判断精度，$eps$在递归时每次除以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\n",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\n",asme(zero,inf,eps));
    return 0;
}

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

---

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

• simpson的其他用途:

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

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