自适应辛普森法

pufanyi

发布于：Feb 4, 2019

并没想打算开篇小小结，只是随便写写（感觉开小小结的成本太高了啊……）

感觉还是挺妙的。以前听 zyy 大佬讲课的时候说过，函数的积分用二次函数拟合一下。

加入我们要算 $\int_l^r F(x)\mathrm{d}x$ ，我们用二次函数 $f(x)=ax^2+bx+c$ 来拟合。

我们发现：

\begin{aligned} \int_l^r f(x)\mathrm{d}x&=\frac{a}{3}r^3+\frac{b}{2}r^2+cr-\frac{a}{3}l^3-\frac{b}{2}l^2-cl\\ &=\frac{a}{3}\left(r-l\right)\left(r^2+l^2+lr\right)+\frac{b}{2}\left(r-l\right)\left(r+l\right)+c\left(r-l\right)\\ &=\frac{r-l}{6}\left(2ar^2+2al^2+2alr+3br+3bl+6c\right)\\ &=\frac{r-l}{6}\left(ar^2+br+c+al^2+bl+c+ar^2+al^2+2alr+2br+2bl+4c\right)\\ &=\frac{r-l}{6}\left(f(l)+f(r)+a(l+r)^2+2b(l+r)+4c\right)\\ &=\frac{r-l}{6}\left(f(l)+f(r)+4\left(a\left(\frac{l+r}{2}\right)^2+b\left(\frac{l+r}{2}\right)+c\right)\right)\\ &=\frac{r-l}{6}\left(f(l)+f(r)+4f\left(\frac{l+r}{2}\right)\right) \end{aligned}

于是我们有：

\int_l^r F(x)\mathrm{d}x\approx\frac{r-l}{6}\left(F(l)+F(r)+4F\left(\frac{l+r}{2}\right)\right)

然后就可以拟合辣！

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

然而就是要么 WA（炸精度）要么 T……

有一种叫自适应辛普森的方法解决这个问题，我们用辛普森计算 $\int_l^rF(x)\mathrm{d}x$ ， $\int_l^{(l+r)/2}F(x)\mathrm{d}x$ ， $\int_{(l+r)/2}^rF(x)\mathrm{d}x$ 的值，如果 $\left|\int_l^rF(x)\mathrm{d}x-\left(\int_l^{(l+r)/2}F(x)\mathrm{d}x+\int_{(l+r)/2}^rF(x)\mathrm{d}x\right)\right|\le \mathrm{eps}$ ，那么直接返回值，否则递归下去算。

upd：忽然发现我学了个假的 simpson。其实eps 应该是相对误差，所以每次递归下去时 $\mathrm{eps}$ 都应该除 $2$ 。听说最后判定时应该是 $\le 15\cdot \mathrm{eps}$ ，虽然我并不知道是为什么。话说这和 eps 乘上 15 有什么区别吗

这样复杂度就是 $\mathcal{O}\left(\text{玄学}\right)$ ，交之前洗把脸就行了。

inline double asr(double l, double r, double ans, double eps)
{
  double mid = (l + r) * .5;
  double ll = simpson(l, mid), rr = simpson(mid, r);
  if(fabs(ll + rr - ans) < 15 * eps)
    return ll + rr;
  return asr(l, mid, ll, eps * .5) + asr(mid, r, rr, eps * .5);
}