Simpson积分学习笔记

昨天$ACM$模拟的时候遇到了一道$Simpson$积分相关的题

完全不知道怎么求，我们组$Fishmen$被$Bym$嘲讽了很久

于是今天下午结合各种资料还是看了一下

这个东西不要觉得它看上去讲什么积分很高级 实际上认真推导也不是很难

首先$Simpson$积分的式子是这个东西

积分号的下标代表区间左端点，上标代表区间右端点

这个值的意思就是函数在$[a, b]$区间与$x$轴，$x = a$， $x = b$围成的面积

首先引入一些东西，$Simpson$积分就是用一个二次曲线来无限逼近原曲线

来达到求得近似值的效果

首先引入一个东西 对于一个二次函数$f(x) = ax^2 + bx + c$

求积分可得$F(x) = \displaystyle\int_0^xf(x)dx = \frac{a}{3}x^3 + \frac{b}{2}x^2 + cx + D$

在这里$D$是一个常数 这个东西具体怎么证等我问了物理组的同学再来填坑

那么$\displaystyle\int_L^Rf(x)dx = F(R) - F(L)$
$\displaystyle\int_L^Rf(x)dx = \frac{a}{3}(R^3-L^3) + \frac{b}{2}(R^2-L^2)+c(R-L)$
$\displaystyle\int_L^Rf(x)dx = (R - L)[\frac{a}{3}(L^2+R^2+LR) + \frac{b}{2}(L + R)+c]$
$\displaystyle\int_L^Rf(x)dx = \frac{b - a}{6}(2aL^2+2aR^2+2aLR+3bL+3bR+6c)$
$\displaystyle\int_L^Rf(x)dx = \frac{b - a}{6}\{(aL^2 + bL + c) + (aR^2 + bR + c) + 4[a(\frac{L + R}{2})^2 + b(\frac{L + R}{2}) + c]\}$
$\displaystyle\int_L^Rf(x)dx =\frac{b - a}{6}[f(L) + f(R) + 4f(\frac{L + R}{2})]$

就这样证完了… 然后我们来考虑这个代码怎么实现

double Simpson(double a, double b) {
    int mid = (a + b) / 2;
    return (b - a) / 6 * (f(a) + f(b) + 4 * f(mid));
}

好了写完了
这样子直接计算我们发现误差非常大，毕竟原图像可能不能很好的用二次函数逼近

自适应Simpson积分

那怎么办呢，我们可以考虑模拟微分的过程，把区间分成无穷小份再把得到的答案积分就是最后的答案了

但是一般我们只要求一个近似值，所以我们就递归到满足正常误差即可

还有一个问题 每次递归下去都有精度误差，那么累加起来可能就和正确答案相差甚远了

这个时候我们可以在递归的时候提高精度要求，尽量让小的答案精确这样就可以了

在刘汝佳的代码中判断是这样写的

double _asr(double l, double r, double v, double eps) {
    double mid = (l + r) / 2;
    double v1 = simpson(l, mid), v2 = simpson(mid, r);
    if(fabs(v1+v2-v) <= 15*eps) return v1+v2 + (v1+v2-v)/15;
    else return _asr(l, mid, v1, eps/2) + _asr(mid, r, v2, eps/2);
}

大概意思就是把多出的精度损失也加回去，这样子能够比较好的保证精确度

裸的模板题

这个就不用讲了 套板子即可

Codes

#include<bits/stdc++.h>

using namespace std;
typedef double db;

const db eps = 1e-6;
db a, b, c, d;

db f(db x) {
    return (c * x + d) / (a * x + b);
}

db integral(db x1, db x2) {
    db mid = (x1 + x2) / 2;
    return (x2 - x1) / 6 * (f(x1) + f(x2) + f(mid) * 4);
}

db Simpson(db x1, db x2, db eps) {
    db mid = (x1 + x2) / 2, res1, res2, res;
    res = integral(x1, x2); res1 = integral(x1, mid); res2 = integral(mid, x2);
    if(fabs(res1 + res2 - res) < eps * 15) return res1 + res2 + (res1 + res2 - res) / 15;   
    return Simpson(x1, mid, eps / 2) + Simpson(mid, x2, eps / 2);
}

int main() {
#ifndef ONLINE_JUDGE
    freopen("4525.in", "r", stdin);
    freopen("4525.out", "w", stdout);
#endif
    db L, R;
    cin >> a >> b >> c >> d >> L >> R;
    printf("%.6lf", Simpson(L, R, eps));
    return 0;
}