浅谈自适应Simpson法

Tiphereth_A

2018-08-12 15:33:58

Algo. & Theory

点击前往我的blog获取更好阅读体验

另一篇讲解

0 前置知识&&写在前面

阅读本篇文章需要掌握定积分

1 先来谈谈求面积

有的读者可能会有疑问：不是要讲自适应Simpson法吗，怎么变成求面积了？

因为Simpson法本身就是用来求不规则图形面积的

说到不规则图形，我们往往都是先从曲边梯形开始

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

对于这个问题，我们的解决方案是：把曲边梯形分成n段，每一段用一些规则的几何图形去近似，然后把每一段的面积加起来，这样我们就得到近似值了

可以看出，上述过程的关键就是选择什么样的几何图形去近似

当然，用不同的几何图形近似，效果是不同的

1.1 用矩形去近似

当然，我们可以看出这种近似方法太粗糙了，针对用矩形近似的方案，我们可以做一些优化：

对于每一段，我们取端点中点在函数上的对应点，借助这个点来构造矩形：

这样看起来就舒服多了，但感觉还是有些粗糙，有没有更好的方法呢？

当然有了！

不过在继续之前，我们先来看看如何实现这种方法

设C(a,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)}

为了方便，我们让每一段的长度相等，即对于每一段，有

\Delta x=\frac{b-a}{n}

那么

\displaystyle\intop_a^bf(x)\mathrm{d}x\thickapprox\Delta x\sum_{i=1}^{n-1}{f((i+\frac{1}{2})\Delta x)}

1.2 用梯形去近似

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

有一些部分看起来已经足够精确了，但感觉还是有些粗糙，有没有更好的方法呢？

不过在继续之前，我们还是先来看看如何实现这种方法

设C(a,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})

为了方便，我们让每一段的长度相等，即对于每一段，有

\Delta x=\frac{b-a}{n}

则有

\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})

1.3 Simpson法

好了，终于到正题了

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

（因为笔者在GeoGebra里没找到根据三点画抛物线的工具，所以这里用圆弧代替抛物线了QwQ）

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

我们来看看如何实现这种方法

设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)=Ax^2+Bx+C来取代f(x)

所以有

f(x_{2i-2})=g(x_{2i-2})=Ax_{2i-2}^2+Bx_{2i-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(x_{2i})=g(x_{2i})=Ax_{2i}^2+Bx_{2i}+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 \qquad\qquad\qquad\qquad\qquad\quad=(\frac{A}{3}x^3+\frac{B}{2}x^2+Cx)\Big|_{x_{2i-2}}^{x_{2i}} \qquad\qquad\qquad\qquad\qquad\qquad\qquad\enspace=\frac{\Delta x}{3}[f(x_{2i-2})+4f(x_{2i-1})+f(x_{2i})]

所以

\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})]

不过笔者在网上查到的资料有一部分认为Simpson法是n=1时的情形，即

\displaystyle\intop_a^bf(x)\mathrm{d}x\thickapprox\frac{b-a}{6}[f(a)+4f(\frac{a+b}{2})+f(b)]

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

2 自适应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\epsilon时S_0与S_1+S_2差别不大（乘以15这里可以按需调整）

（2018.09.14 upd: 乘以15这里是经过一系列误差分析后得出的，具体笔者可能会另写一篇文章咕咕咕，感谢@Marser和@_rqy两位dalao的补充）

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

附程序：

inline double F(double num) {
  //自己按需补充
}
inline double simpson(double a,double b) {
  double c=a+(b-a)/2;//防溢出
  return (F(a)+4*F(c)+F(b))*(b-a)/6
}
double asr(double a,double b,double eps,double S) {
  double c=a+(b-a)/2;
  double lS=simpson(a,c),rS=simpson(c,b);
  if(fabs(lS+rS-S)<=15*eps) return lS+rS+(lS+rS-S)/15.0;
  return asr(a,c,eps/2,lS)+asr(c,b,eps/2,rS);//注意这里eps要除以2
}

3 后记

这篇文章笔者写了4h吧，内容还算简单，希望各位能够愉快地享用~(￣▽￣)~*

btw，洛谷P4525、P4526是模板题ヾ(≧▽≦*)连切两道紫题真开心

4 主要参考书目