自适应辛普森法是一种可以在较快时间内求解定区间数值积分的算法，其思想为将这个定区间内的函数分段并近似为二次函数进行求值。
1. 辛普森公式 辛普森公式可以为一个二次函数求解定积分。
对于一个二次函数 $f(x)$，有：
证明：
设 $f(x)=Ax^2+Bx+C$，则有 $g(x)=\dfrac{A}{3}x^3+\dfrac{B}{2}x^2+Cx$ 使得 $g’(x)=f(x)$。
因此原式可化为：
2. 普通辛普森法 该方法在1743年发表于托马斯·辛普森的一篇论文中。方法如下：
给定一个正整数 $n$，将区间 $[l,r]$ 分成 $2n$ 个子区间，则两个区间的空隙为 $h=\dfrac{r-l}{2n}$，第 $i$ 个区间的结束点为 $p_i=l+h\cdot i$（为方便计算，令 $p_0=l$）。
接着对于 $i\in[1,n]$，依次用二次函数对 $[p_{2i-2)},p_{2i}]$ 区间内的函数值取近似，则此段内的积分值约为 $\dfrac{h(f(p_{2i-2})+4f(p_{2i-1})+f(p_{2i}))}{3}$。
则可算出积分值约为：
查看代码
1 2 3 4 5 6 7 8 9 double simple_simpson (double l,double r,int n = 1e6 ) double h = (r - l) / (n * 2.0 ); double ans = f(l)+f(r); for (int i = 1 ;i < 2 * n;i ++) { double p = l + h * i; ans += f(p) * ((i & 1 ) + 1 ) * 2.0 ; } return ans * h / 3.0 ; }
普通辛普森法的误差约为 $-\frac{1}{90}(\frac{r-l}{2})^5f^4(\xi)$，其中 $l \le \xi \le r$。
3. 自适应辛普森法 很明显，普通辛普森法的时间和精度受限于 $n$ 的大小，这使得我们不能很好控制计算结果的效率和精度。
考虑设计一个可以自我调节精度的算法，当二次函数足够接近原函数值时直接返回，否则将区间分割为左右两部分，递归重复此过程。
注意到，如果区间直接代入公式算出的结果与左右两边代入公式算出结果之和相差不大时，这段函数可以近似认为是二次函数而直接返回。
这就是自适应辛普森法的思想。
查看代码
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 double simpson (double l,double r) return (r - l) * (f(l) + 4 * f((l + r) / 2 ) + f(r)) / 6 ; } double asr (double l,double r,double eps,double A) double mid = (l + r) / 2 ; double 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); } double solve (double l,double r) return asr(l,r,eps,simpson(l,r)); }
4. 例题 1. 洛谷P4525 【模板】自适应辛普森法1 给定 $a,b,c,d,L,R$，计算 $\int_L^R \dfrac{cx+d}{ax+b} \mathrm{d}x$，结果保留至小数点后 $6$ 位。
模板题，直接套上模板即可。
查看代码
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 #include <cstdio> #include <cmath> typedef double D;const D eps = 1e-7 ;D a,b,c,d,l,r; D f (D x) { return (c * x + d) / (a * x + b); } D simpson (D l,D r) { return (r - l) * (f(l) + 4 * f((l + r) / 2 ) + f(r)) / 6 ; } D asr (D l,D r,D A) { D mid = (l + r) / 2 ; D 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,L) + asr(mid,r,R); } D solve (D l,D r) { return asr(l,r,simpson(l,r)); } int main () scanf ("%lf%lf%lf%lf%lf%lf" ,&a,&b,&c,&d,&l,&r); printf ("%.6lf\n" ,solve(l,r)); return 0 ; }
2. 洛谷P45256【模板】自适应辛普森法2 给定 $a(a \le 50)$，计算 $\int_0^\infty x^{\frac{a}{x}-x} \mathrm{d}x$，结果保留至小数点后 $5$ 位。若积分发散，输出 $\text{orz}$。
在 $\text{desmos}$ 上画个图之后，可以注意到原函数当 $a<0$ 时发散，而当 $a \ge 0,x\ge 15$ 时，原函数的值已无限趋近于 $0$，可以忽略不计。因此可以利用此性质将原积分化为 $\int_{eps}^{15} x^{\frac{a}{x}-x} \mathrm{d}x$（$x$ 不能取 $0$，因此从 $eps$ 开始计算积分）。
查看代码
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 #include <cstdio> #include <cmath> typedef double D;const D eps = 1e-7 ;D a; D f (D x) { return pow (x,a / x - x); } D simpson (D l,D r) { return (r - l) * (f(l) + 4 * f((l + r) / 2 ) + f(r)) / 6 ; } D asr (D l,D r,D A) { D mid = (l + r) / 2 ; D 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,L) + asr(mid,r,R); } D solve (D l,D r) { return asr(l,r,simpson(l,r)); } int main () scanf ("%lf" ,&a); if (a < 0 ) { puts ("orz" ); return 0 ; } printf ("%.5lf\n" ,solve(eps,15 )); return 0 ; }
