算法笔记-自适应辛普森法

· · 题解

辛普森法（Simpson），是一种求数值积分的方法。

求积分的时候，可以将其分成无穷多个小区间或者用 牛顿-莱布尼茨公式。

如果不能用这两种方法，就有了 数值积分。

辛普森法的核心思想就是把函数 f 近似地看成二次函数，于是给出辛普森公式的推导：

由于把 f 近似地看成二次函数，所以，对于二次函数 f(x)=ax^2+bx+c，求积分可得：

\displaystyle F(x)=\int^x_0 f(x)\mathrm{d}x=\dfrac{a}{3}x^3+\dfrac{b}{2}x^2+cx+k

其中 k 为一个常数，那么：

\int_{l}^{r} f(x) \mathrm{d} x & =F(r)-F(l) \\ & =\frac{a}{3}\left(r^{3}-l^{3}\right)+\frac{b}{2}\left(r^{2}-l^{2}\right)+c(r-l) \\ & =(r-l)\left(\frac{a}{3}\left(l^{2}+r^{2}+l r\right)+\frac{b}{2}(l+r)+c\right) \\ & =\frac{r-l}{6}\left(2 a l^{2}+2 a r^{2}+2 a l r+3 b l+3 b r+6 c\right) \\ & =\frac{r-l}{6}\left(\left(a l^{2}+b l+c\right)+\left(a r^{2}+b r+c\right)+4\left(a\left(\frac{l+r}{2}\right)^{2}+b\left(\frac{l+r}{2}\right)+c\right)\right) \\ & =\frac{r-l}{6}\left(f(l)+f(r)+4 f\left(\frac{l+r}{2}\right)\right) \end{aligned}

稍微改写下形式就得到了辛普森公式：

\int_{l}^{r} f(x) \mathrm{d} x=\dfrac{r-l\left(f(l)+f(r)+4 f\left(\frac{l+r}{2}\right)\right)}{6}
  1. 辛普森法

有一自然数 n，将区间 [l,r] 分成等长的 2n 个区间，计算每个小区间的积分值，求和可得总积分。

对于 [x_{2i},x_{2i-2}] 的一段区间，选定三个点 \left( x_{2i-2},x_{2i-1},x_{2i} \right)，构成抛物线，得到函数 f^{'}(x)，于是，原来计算么个小区间的积分值就可以转化为求二次函数 f^{'}(x) 的值，而 f^{'}(x) 可以用辛普森公式近似计算。

代码实现很简单，就不再赘述了。

  1. 自适应辛普森法

从上文可知，结果的进度取决于 n 的取值，所以，如果一段积分已经很接近二次函数了，误差值在接受范围内，就可以不用往下继续分段了。

可以对于一段积分先代入辛普森公式求值，然后再将这段积分平分成两段继续求积分，如果差值小于要求的精度 eps 就直接返回。

\quad

[模板]自适应辛普森法 1

标准模板，没什么好说的：

#include <bits/stdc++.h>
using namespace std;
const double eps = 1e-6;
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.0;
    return (f(l) + 4.0 * f(mid) + f(r)) * (r - l) / 6.0;
}
double asr(double l, double r, double eps, double ans) {
    double mid = (l + r) / 2.0, l_ = simpson(l, mid), r_ = simpson(mid, r);
    if (fabs(l_ + r_ - ans) <= 15.0 * eps) return l_ + r_ + (l_ + r_ - ans) / 15.0;
    return asr(l, mid, eps / 2.0, l_) + asr(mid, r, eps / 2.0, r_);
}
int main() {
    ios::sync_with_stdio(false), cin.tie(nullptr), cout.tie(nullptr);
    cin >> a >> b >> c >> d >> l >> r;
    cout << fixed << setprecision(6) << asr(l, r, eps, simpson(l, r));
    return 0;
}

[模板]自适应辛普森法 2

一看，\int^\infty_0，可辛普森法只能求定值积分，所以别急，画了图像就能知道，当 a<0 时，积分发散，就不用求了，剩下的套模板 1 即可。

精度多取两位保险。

#include <bits/stdc++.h>
using namespace std;
const double eps = 1e-8;
double a, b, c, d, l, r;
double f(double x) {
    return pow(x, (a / x) - x);
}
double simpson(double l, double r) {
    double mid = (l + r) / 2.0;
    return (f(l) + 4.0 * f(mid) + f(r)) * (r - l) / 6.0;
}
double asr(double l, double r, double eps, double ans) {
    double mid = (l + r) / 2.0, l_ = simpson(l, mid), r_ = simpson(mid, r);
    if (fabs(l_ + r_ - ans) <= 15.0 * eps) return l_ + r_ + (l_ + r_ - ans) / 15.0;
    return asr(l, mid, eps / 2.0, l_) + asr(mid, r, eps / 2.0, r_);
}
int main() {
    ios::sync_with_stdio(false), cin.tie(nullptr), cout.tie(nullptr);
    cin >> a;
    if (a < 0) return cout << "orz", 0;
    cout << fixed << setprecision(5) << asr(eps, 15.0, eps, simpson(l, r));
    return 0;
}