自适应辛普森法学习笔记

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

参考资料

刘汝佳 & 陈锋《算法竞赛进阶指南·入门经典》
Ebola 《题解 P4525 【【模板】自适应辛普森法1】》
function_of_zero 《题解 P4525 【【模板】自适应辛普森法1】》

引入