自适应辛普森法 学习笔记

作者: xht37 分类: 笔记 发布时间: 2019-08-04 23:12

点击数:585

OI 中的微积分——自适应辛普森法

引入

我们要计算这样一个式子:

$$\int_l^r f(x) {\rm d}x$$

我们显然不可能让计算机去推柿子。

那怎么办呢?

二次函数的定积分

对于一个二次函数

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

我们有

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

函数的拟合

对于一个奇怪的函数,为了对其求导,我们可以用一个图像近似且容易求导的函数(这里我们使用二次函数)来代替它,这个过程称之为“拟合”。

这样的话精度误差可能会很大,因此我们需要将函数分段拟合。

自适应辛普森法

显然我们并不知道分成多少段合适,因为段分得越多,精度误差越小,但时间消耗越大;段分得越少,精度误差越大,但时间消耗越小。

我们需要一个“自适应”的操作。

如果这一段函数与拟合的二次函数“相似”,那么我们直接把这一段当做这个二次函数,然后套用公式计算;如果这一段函数与拟合的二次函数“不相似”,那么我们应该把这一段分成两半,并递归进行这一过程。

接下来的问题就是,如何判断是否“相似”呢?

我们可以对整段、左半部分、右半部分分别套用二次函数定积分公式进行计算,结果分别记作 $A,L,R$,若 $L+R$ 与 $A$ 相差的值不超过我们设定的精度了,那么我们就认为这一段函数与拟合的二次函数是“相似”的。

namespace Simpson {
    #define mid ((l + r) / 2)

    inline ld f(ld x) {
        return ...;
    }

    inline ld simpson(ld l, ld r) {
        return (r - l) * (f(l) + f(r) + 4 * f(mid)) / 6;
    }

    inline ld asr(ld l, ld r, ld eps, ld A) {
        ld L = simpson(l, mid), R = simpson(mid, r);
        if (fabs(L + R - A) <= 15 * eps) return L + R + (L + R - A) / 15;
        return asr(l, mid, eps / 2, L) + asr(mid, r, eps / 2, R);
    }

    inline ld main(ld l, ld r) {
        return asr(l, r, 1e-10, simpson(l, r));
    }
}

【练习】P4525 【模板】自适应辛普森法1

优化

仔细观察可以发现,上面那份代码中有很多函数值被计算了多次,当计算函数值所需要的时间较大时,这份代码的常数将会爆炸。

我们可以进行一些简单的常数优化。

namespace Simpson {
    #define mid ((l + r) / 2)

    inline ld f(ld x) {
        return ...;
    }

    inline ld simpson(ld l, ld r, ld fl, ld fr, ld fm) {
        return (r - l) * (fl + fr + 4 * fm) / 6;
    }

    inline ld asr(ld l, ld r, ld eps, ld fl, ld fr, ld fm) {
        ld flm = f((l + mid) / 2), frm = f((r + mid) / 2);
        ld L = simpson(l, mid, fl, fm, flm), R = simpson(mid, r, fm, fr, frm), A = simpson(l, r, fl, fr, fm);
        if (fabs(L + R - A) <= 15 * eps) return L + R + (L + R - A) / 15;
        return asr(l, mid, eps / 2, fl, fm, flm) + asr(mid, r, eps / 2, fm, fr, frm);
    }

    inline ld main(ld l, ld r) {
        return asr(l, r, 1e-10, f(l), f(r), f(mid));
    }
}

【练习】P4526 【模板】自适应辛普森法2

【练习】UVA1356 Bridge

参考资料

发表评论

电子邮件地址不会被公开。 必填项已用*标注