自适应Simpson积分学习笔记

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

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

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

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

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

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

利用矩形去近似面积

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

方法实现：

设$C(a,0)$C(a,0)，$D(b,0)$D(b,0)。那么：$\displaystyle\intop_a^bf(x)\mathrm{d}x\thickapprox\Delta x_i\sum_{i=1}^{n-1}{f((i+\frac{1}{2})\Delta x_i)}$a∫b​f(x)dx≈Δxi​i=1∑n−1​f((i+21​)Δxi​)。

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

那么：$\color{red}\displaystyle\intop_a^bf(x)\mathrm{d}x\thickapprox\Delta x\sum_{i=1}^{n-1}{f((i+\frac{1}{2})\Delta x)}$a∫b​f(x)dx≈Δxi=1∑n−1​f((i+21​)Δx)。

利用梯形去近似面积

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

方法实现：

设$C(a,0)$C(a,0)，$D(b,0)$D(b,0)。那么：$\displaystyle\intop_a^bf(x)\mathrm{d}x\thickapprox\Delta x_i(\sum_{i=1}^{n-1}{f(i\Delta x_i)}+\frac{f(a)+f(b)}{2})$a∫b​f(x)dx≈Δxi​(i=1∑n−1​f(iΔxi​)+2f(a)+f(b)​)。

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

则有$\color{red}\displaystyle\intop_a^bf(x)\mathrm{d}x\thickapprox\Delta x(\sum_{i=1}^{n-1}{f(i\Delta x)}+\frac{f(a)+f(b)}{2})$a∫b​f(x)dx≈Δx(i=1∑n−1​f(iΔx)+2f(a)+f(b)​)。

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

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

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

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

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

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

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

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

所以有：

$f(x_{2i-2})=g(x_{2i-2})=Ax_{2i-2}^2+Bx_{2i-2}+C$f(x2i−2​)=g(x2i−2​)=Ax2i−22​+Bx2i−2​+C。

$f(x_{2i-1})=g(x_{2i-1})=Ax_{2i-1}^2+Bx_{2i-1}+C=A(\frac{x_{2i-2}+x_{2i}}{2})^2+B(\frac{x_{2i-2}+x_{2i}}{2})+C$f(x2i−1​)=g(x2i−1​)=Ax2i−12​+Bx2i−1​+C=A(2x2i−2​+x2i​​)2+B(2x2i−2​+x2i​​)+C。

$f(x_{2i})=g(x_{2i})=Ax_{2i}^2+Bx_{2i}+C$f(x2i​)=g(x2i​)=Ax2i2​+Bx2i​+C。

于是$\displaystyle\intop_{x_{2i-2}}^{x_{2i}}f(x)\mathrm{d}x\thickapprox\intop_{x_{2i-2}}^{x_{2i}}g(x)\mathrm{d}x$x2i−2​∫x2i​​f(x)dx≈x2i−2​∫x2i​​g(x)dx

$\qquad\qquad\qquad\enspace=(\frac{A}{3}x^3+\frac{B}{2}x^2+Cx)\Big|_{x_{2i-2}}^{x_{2i}}$=(3A​x3+2B​x2+Cx)∣∣∣​x2i−2​x2i​​

$\qquad\qquad\qquad\enspace=\frac{\Delta x}{3}[f(x_{2i-2})+4f(x_{2i-1})+f(x_{2i})]$=3Δx​[f(x2i−2​)+4f(x2i−1​)+f(x2i​)]

所以：$\color{red}\displaystyle\intop_a^bf(x)\mathrm{d}x\thickapprox\frac{\Delta x}{3}\sum_{i=0}^{2n-2}[f(x_{2i})+4f(x_{2i+1})+f(x_{2i+2})]$a∫b​f(x)dx≈3Δx​i=0∑2n−2​[f(x2i​)+4f(x2i+1​)+f(x2i+2​)]。

不过也有人认为Simpson法是$n=1$n=1时的情形，即$\color{blue}\displaystyle\intop_a^bf(x)\mathrm{d}x\thickapprox\frac{b-a}{6}[f(a)+4f(\frac{a+b}{2})+f(b)]$a∫b​f(x)dx≈6b−a​[f(a)+4f(2a+b​)+f(b)]。

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

自适应Simpson法

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

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

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

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

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

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

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

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

我们认为当且仅当$|S_1+S_2-S_0|&lt;15\epsilon$∣S1​+S2​−S0​∣<15ϵ时$S_0$S0​与$S_1+S_2$S1​+S2​差别不大（乘以$15$15这里可以按需调整）。

返回的面积则是$\color{blue}S_1+S_2+\frac{S_1+S_2-S_0}{15}$S1​+S2​+15S1​+S2​−S0​​。

参考代码

#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

计算积分：$\displaystyle\intop_L^R\frac{cx+d}{ax+b}\mathrm{d}x$L∫R​ax+bcx+d​dx结果保留至小数点后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

计算积分$\displaystyle\intop_0^\infty x^{\frac{a}{x}-x}\mathrm{d}x$0∫∞​xxa​−xdx。

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

[Analysis]

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

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

• 当$a&lt;0$a<0时，函数是发散的，直接输出orz即可。
• 当$a&gt;0，x&gt;10$a>0，x>10时，函数值就趋近于0了，于是我们只用考虑$\color{blue}x\in (0,10]$x∈(0,10]的情况了__（注意$\color{red}x!=0$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;
}

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

阅读原文