辛普森法求数值积分

部分内容来自 OI-Wiki。

首先，我们引入二次函数积分公式。

对于二次函数

$$ f(x)=ax^2+bx+c, $$

有

$$ \int_{L}^{R} f(x)\text{d}x=\dfrac{R-L}{6}\left[f(L)+f(R)+4f\left(\dfrac{L+R}{2}\right)\right] $$

这里给出计算过程：

先求出 $f(x)$ 的原函数：

$$ F(x)=\int f(x)\text{d}x=\dfrac{a}{3}x^3+\dfrac{b}{2}x^2+cx+D $$

由牛顿-莱布尼茨公式，

$$ \begin{aligned} \int_{L}^{R} f(x)\text{d}x &=F(R)-F(L)\\ &=\dfrac{a}{3}(R^3-L^3)+\dfrac{b}{2}(R^2-L^2)+c(r-l)\\ &=(R-L)\left(\dfrac{a}{3}(L^2+LR+R^2)+\dfrac{b}{2}(L+R)+c\right)\\ &=\dfrac{R-L}{6}(2L^2+2LR+2R^2+3L+3R+6c)\\ &=\dfrac{R-L}{6}\left\{(aL^2+bL+c)+(aR^2+bR+c)+4\left[a\left(\dfrac{l+r}{2}\right)^2+b\left(\dfrac{l+r}{2}\right)+c\right]\right\}\\ &=\dfrac{R-L}{6}\left[f(L)+f(R)+4f\left(\dfrac{L+R}{2}\right)\right] \end{aligned} $$

需要注意这里其实不要求 $a\ne 0$。

考虑使用二次函数拟合被积函数，将被积函数分成若干段，每段内使用辛普森公式求近似值。

但是，如果我们不断细分，无疑会增大算量；如果细分的程度不够，又会影响精度。

一般，我们先将函数分成几段，接下来对每段检查二次函数是否能良好拟合。

具体地，对于函数的某一段 $[l,r]$，如果拟合 $[l,r]$ 的结果跟 $[l,\frac{l+r}{2}]$ 与 $[\frac{l+r}{2},r]$ 的和接近，那么认为拟合良好。

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

double asr(double l, double r, double eps, double ans, int step) {
  double mid = (l + r) / 2;
  double fl = simpson(l, mid), fr = simpson(mid, r);
  if (abs(fl + fr - ans) <= 15 * eps && step < 0)
    return fl + fr + (fl + fr - ans) / 15;
  return asr(l, mid, eps / 2, fl, step - 1) +
         asr(mid, r, eps / 2, fr, step - 1);
}

double calc(double l, double r, double eps) {
  return asr(l, r, eps, simpson(l, r), 12);
}

其中 15 是一个 Magic Number，在 https://doi.org/10.1145/321526.321537 给出了论述。

需要指出，也存在一些病态函数，直接应用自适应辛普森法可能会出错。