自适应辛普森法

模板题：P4525 【模板】自适应辛普森法 1

题意：

试计算积分

$$\displaystyle{\int_L^R\frac{cx+d}{ax+b}\,\mathrm{d}x}$$

结果保留至小数点后 $6$ 位。

数据保证计算过程中分母不为 $0$ 且积分能够收敛。

输入格式：

一行，包含 $6$ 个实数 $a,b,c,d,L,R$。

输出格式：

一行，积分值，保留至小数点后 $6$ 位。

数据范围：

$a,b,c,d\in[-10,10]$，$-100\le L<R\le 100$ 且 $R-L\ge1$。

先给出 Simpson 公式

$$\int_a^b f(x) \,\mathrm{d} x \,\approx\,(b-a)\cdot \frac{f(a)+4 f\left(\frac{a+b}{2}\right)+f(b)}{6}$$

下面给出推导和介绍

Simpson 公式的原理在于用二次函数去拟合区间，并计算定积分。

显然，区间划分的越多，精确度也就相对越高，不过相应的时间复杂度也会增大。

不过我们先来看看 Simpson 公式是怎么推导的。

$$\begin{aligned} \int_a^b f(x) \,\mathrm{d} x &\approx \int_a^b (A x^2+B x+C) \,\mathrm{d} x \\[6pt]&=\frac{A}{3}\left(b^3-a^3\right)+\frac{B}{2}\left(b^2-a^2\right)+C(a-b) \\[6pt]&=\frac{(b-a)}{6}\left(2 A\left(b^2+a b+a^2\right)+3 B(b+a)+6 C\right) \\[6pt]&=\frac{(b-a)}{6}\left(2 A b^2+2 A a b+2 A a^2+3 B b+3 B a+6 C\right) \\[6pt]&=\frac{(b-a)}{6}\left(A a^2+B a+C+A b^2+B b+C+4 A\left(\frac{a+b}{2}\right)^2+4 B\left(\frac{a+b}{2}\right)+4 C\right) \\[6pt]&=\frac{(b-a)}{6}\left(f(a)+f(b)+4 f\left(\frac{a+b}{2}\right)\right) \end{aligned}$$

这里第一行用到了积分的小知识，可以看我之前写的这篇不定积分基本公式

接着我们来讲本题用到的自适应辛普森法。

刚刚我们提到了 Simpson 公式需要根据精度进行分段拟合，而自适应辛普森法就是实现这个操作的算法

算法原理也很简单，就是把原来的区间二分递归，如果精度不够就继续递归，否则直接退出。

时间复杂度取决于原函数 $f(x)$ 的复杂程度。

代码：

#include <bits/stdc++.h>
using namespace std;
#define int long long
#define INF 0x3f3f3f3f3f3f3f3f
void up(int &x,int y) { x < y ? x = y : 0; }
void down(int &x,int y) { x > y ? x = y : 0; }
#define N ((int)())

const double eps = 1e-12;
double a,b,c,d,l,r;
double f(double x) { return (c * x + d) / (a * x + b); }
double simpson(double l, double r)
{
    double mid = (l + r) / 2;
    return (f(l) + 4 * f(mid) + f(r)) * (r - l) / 6;
}
double calc(double l, double r,double ans)
{
    double mid = (l + r) / 2;
    double L = simpson(l, mid), R = simpson(mid, r);
    if(fabs(L + R - ans) <= eps) return L + R + (L + R - ans);
    return calc(l, mid, L) + calc(mid, r, R);
}
signed main()
{
    ios::sync_with_stdio(0);
    cin.tie(0);cout.tie(0);
    // freopen("check.in","r",stdin);
    // freopen("check.out","w",stdout);
    cin >> a >> b >> c >> d >> l >> r;
    cout << fixed << setprecision(6) << calc(l, r, simpson(l, r)) << '\n';
    return 0;
}

---

参考文献：

[1] https://www.luogu.com.cn/blog/beili233/solution-p4525