python 辛普森积分 辛普森积分法公式

定积分

函数 \(f(x)\) 在区间 \([l,r]\) 上的定积分 \(\int_{l}^{r} f(x) dx\) 指的是 \(f(x)\) 在区间 \([l,r]\) 内与 \(x\) 轴所围成的区域的面积（\(x\)轴上方为正，下方为负）。

我们需要一种高效的求解这种积分的近似值的方法，于是就有了辛普森积分法。

普通辛普森法

辛普森法的基本思想是将求解区间分成若干段，每一段都使用二次函数的积分公式来进行求解。

二次函数积分公式（辛普森公式）：

对于一个二次函数\(f(x)=Ax^2+Bx+C\)，有

$\[\int_l^r f(x) {\mathrm d}x = \frac{(r-l)(f(l)+f(r)+4 f(\frac{l+r}{2}))}{6} \]$

证（from OI-Wiki）：

求积分可得

$\[F(x)=\int_0^x f(x) {\mathrm d}x = \frac{a}{3}x^3+\frac{b}{2}x^2+cx+D \]$

那么则有

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

然后我们通过套公式就可以写出这样一段代码：

double simpson(double l, double r) {
    const double mid = (l + r) / 2;
    return (r - l) * (f(l) + 4 * f(mid) + f(r)) / 6; //f(x)为待求解的函数
}

自适应辛普森

普通的辛普森积分法为了保证精度，在时间效率上会很大地受到区间的限制。

问题在于：区间少了精度不够，区间多了又太慢。

自适应辛普森做到了自动控制拆分区间的大小。

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

到了这里，你就可以去做 Luogu P4525 【模板】自适应辛普森法1了！