自适应Simpson积分学习笔记

本文转载自洛谷日报微信版第35期，点击文末阅读原文查看原文。

问题引入：求曲边梯形面积

Simpson法的用途就在于求不规则多边形的面积。在这里，我们从曲边梯形说起。

举个例子：曲边梯形ABCD就是下图中曲线AB、线段AC、CD、DB围成的图形，我们想要求出它的面积 。

对于这个问题，在数学中，我们可以将这个多边形在竖直方向切成一条条细线，然后将面积累加，这也是微积分的重要思想。

但是在信息学中，我们无法做到这么精确的事情，于是，我们想办法用一些图形的面积累加求近似解。

利用矩形去近似面积

但是这个方法还是太粗糙了，于是我们考虑优化它——对于每一段，我们取端点中点在函数上的对应点，借助这个点来构造矩形。

方法实现：

设$C(a,0)$，$D(b,0)$。那么：$\underset{a}{\overset{b}{\int }}f(x)\mathrm{d}x\approx \Delta x_i\sum _{i=1}^{n-1}f((i+\frac{1}{2})\Delta x_i)$。

为了方便，我们让每一段的长度相等，即对于每一段，有：$\Delta x=\frac{b-a}{n}$。

那么：$\underset{a}{\overset{b}{\int }}f(x)\mathrm{d}x\approx \Delta x\sum _{i=1}^{n-1}f((i+\frac{1}{2})\Delta x)$。

利用梯形去近似面积

可以证明，这个方法和矩形近似的优化是一样的，不过这种方法的实现简单。

方法实现：

设$C(a,0)$，$D(b,0)$。那么：$\underset{a}{\overset{b}{\int }}f(x)\mathrm{d}x\approx \Delta x_i(\sum _{i=1}^{n-1}f(i\Delta x_i)+\frac{f(a)+f(b)}{2})$。

为了方便，我们让每一段的长度相等，即对于每一段，有$\Delta x=\frac{b-a}{n}$。

则有$\underset{a}{\overset{b}{\int }}f(x)\mathrm{d}x\approx \Delta x(\sum _{i=1}^{n-1}f(i\Delta x)+\frac{f(a)+f(b)}{2})$。

利用抛物线去近似面积（Simpson法）

Simpson法是先将原曲线近似成一段段抛物线，然后再用牛顿-莱布尼茨公式求每一段的面积 。

可以看出，这个方法的效果相当好了。

我们来看看如何实现这种方法：设$C(a,0)$，$D(b,0)$。

为了方便，我们让每一段的长度相等，即对于每一段，有$\Delta x=\frac{b-a}{n}$。

对于每一段，我们如下处理：

设区间起点为$x_{2i-2}$，终点为$x_{2i}$，中点为$x_{2i-1}$。

我们要用过点$(x_{2i-2},f(x_{2i-2}))$，$(x_{2i-1},f(x_{2i-1}))$，$(x_{2i},f(x_{2i}))$的抛物线$g(x)=A{x}^2+Bx+C$来取代$f(x)$。

所以有：

$f(x_{2i-2})=g(x_{2i-2})=A{x}_{2i-2}^2+B{x}_{2i-2}+C$。

$f(x_{2i-1})=g(x_{2i-1})=A{x}_{2i-1}^2+B{x}_{2i-1}+C=A(\frac{x_{2i-2}+x_{2i}}{2}{)}^2+B(\frac{x_{2i-2}+x_{2i}}{2})+C$。

$f(x_{2i})=g(x_{2i})=A{x}_{2i}^2+B{x}_{2i}+C$。

于是$\underset{x_{2i-2}}{\overset{x_{2i}}{\int }}f(x)\mathrm{d}x\approx \underset{x_{2i-2}}{\overset{x_{2i}}{\int }}g(x)\mathrm{d}x$

$=(\frac{A}{3}{x}^3+\frac{B}{2}{x}^2+Cx){|}_{x_{2i-2}}^{x_{2i}}$

$=\frac{\Delta x}{3}[f(x_{2i-2})+4f(x_{2i-1})+f(x_{2i})]$

所以：$\underset{a}{\overset{b}{\int }}f(x)\mathrm{d}x\approx \frac{\Delta x}{3}\sum _{i=0}^{2n-2}[f(x_{2i})+4f(x_{2i+1})+f(x_{2i+2})]$。

不过也有人认为Simpson法是$n=1$时的情形，即$\underset{a}{\overset{b}{\int }}f(x)\mathrm{d}x\approx \frac{b-a}{6}[f(a)+4f(\frac{a+b}{2})+f(b)]$。

下面把这种方法叫做三点Simpson法。

自适应Simpson法

自适应Simpson法就是对Simpson法的一个优化。

对一段区间$[a,b]$，我们做如下操作。

1. 取中点$mid=\frac{a+b}{2}$。

2. 分别对区间$[a,b]$、区间$[a,mid]$、区间$[mid,b]$应用三点Simpson法，设得到的面积分别为$S_0$、$S_1$、$S_2$。

3. 若$S_0$与$S_1+S_2$差别不大，就认为区间$[a,b]$面积的近似值已经求得，否则分别对区间$[a,mid]$、区间$[mid,b]$递归应用本操作。

可以看出这个方法在保证了精度的同时保证了效率。

我们注意到，上述操作中有两个地方含糊不清。

一个是如何确定“差别不大”，一个是面积的近似值已经求得后返回的面积是多少。

我们认为当且仅当$|S_1+S_2-S_0|<15ϵ$时$S_0$与$S_1+S_2$差别不大（乘以$15$这里可以按需调整）。

返回的面积则是$S_1+S_2+\frac{S_1+S_2-S_0}{15}$。

参考代码

#include<iostream>
#include<cstdio>
#include<cmath>
using namespace std;

inline double fun(double x) {//函数带入求值 
    return x*x*x-2*x*x-x+3;
}

inline double sps(double l,double r) {//辛普森积分 
    double mid=(l+r)/2;
    return (fun(l)+4*fun(mid)+fun(r))*(r-l)/6.0;
}

double asr(double l,double r,double eps,double s) {//自适应辛普森积分 
    double mid=(l+r)/2;
    double ls=sps(l,mid);
    double rs=sps(mid,r);
    if (fabs(ls+rs-s)<=(eps*15)) {
        return ls+rs+(ls+rs-s)/15.0;
    }
    return asr(l,mid,eps/2,ls)+asr(mid,r,eps/2,rs);
}

int main() {
    double l,r;
    scanf("%lf%lf",&l,&r);
    printf("%lf",asr(l,r,1e-8,sps(l,r)));
    return 0;
}

经典例题

例1.[LG4525]【模板】自适应辛普森法1

计算积分：$\underset{L}{\overset{R}{\int }}\frac{cx+d}{ax+b}\mathrm{d}x$结果保留至小数点后6位。

数据保证计算过程中分母不为0且积分能够收敛。

[Analysis]

直接套模板就好了。

参考代码：

#include<iostream>
#include<cstdio>
#include<cmath>
using namespace std;

double a,b,c,d;

inline double fun(double x) {//函数带入求值 
    return (c*x+d)/(a*x+b);
}

inline double sps(double l,double r) {//辛普森积分 
    double mid=(l+r)/2;
    return (fun(l)+4*fun(mid)+fun(r))*(r-l)/6;
}

double asr(double l,double r,double eps,double s) {//自适应辛普森积分 
    double mid=(l+r)/2;
    double ls=sps(l,mid);
    double rs=sps(mid,r);
    if (fabs(ls+rs-s)<=eps*15) {
        return ls+rs+(ls+rs-s)/15;
    }
    return asr(l,mid,eps/2,ls)+asr(mid,r,eps/2,rs);
}

int main() {
    double l,r;
    scanf("%lf%lf%lf%lf%lf%lf",&a,&b,&c,&d,&l,&r);
    printf("%.6lf",asr(l,r,1e-8,sps(l,r)));
    return 0;
}

例2.[LG4526]【模板】自适应辛普森法2

计算积分$\underset{0}{\overset{\infty }{\int }}{x}^{\frac{a}{x}-x}\mathrm{d}x$。

保留至小数点后5位。若积分发散，请输出orz。

[Analysis]

这道题的式子看上去有些复杂，尤其是那个无穷大让人很头疼。对于这一类题目，我们可以试着画画函数图像，寻找一下规律。

通过图像，我们可以发现一下规律：

• 当$a<0$时，函数是发散的，直接输出orz即可。
• 当$a>0，x>10$时，函数值就趋近于0了，于是我们只用考虑$x\in (0,10]$的情况了__（注意$x!=0$）__。
• 又因为本题对精度要求较高，要讲右边界改成20才能过。

参考代码：

#include<iostream>
#include<cstdio>
#include<cmath>
using namespace std;

double a;

inline double fun(double x) {//函数带入求值 
    return pow(x,a/x-x);
}

inline double sps(double l,double r) {//辛普森积分 
    double mid=(l+r)/2;
    return (fun(l)+4*fun(mid)+fun(r))*(r-l)/6.0;
}

double asr(double l,double r,double eps,double s) {//自适应辛普森积分 
    double mid=(l+r)/2;
    double ls=sps(l,mid);
    double rs=sps(mid,r);
    if (fabs(ls+rs-s)<=(eps*15)) {
        return ls+rs+(ls+rs-s)/15.0;
    }
    return asr(l,mid,eps/2,ls)+asr(mid,r,eps/2,rs);
}

int main() {
    scanf("%lf",&a);
    if (a<0) {
        puts("orz");
    }
    else {
        printf("%.5lf",asr(1e-8,20,1e-8,sps(1e-8,20)));
    }
    return 0;
}

另外以上方法还可以解决多边形与圆的面积并问题，这样就不用写扫描线啦~\(≧▽≦)/~。

阅读原文