自适应辛普森法学习笔记

ChenSh715

已于 2022-10-16 23:43:38 修改

阅读量183

点赞数 1

分类专栏：学习笔记文章标签：学习算法

于 2022-08-07 21:53:31 首次发布

本文链接：https://blog.csdn.net/weixin_44279267/article/details/126217857

版权

学习笔记专栏收录该内容

3 篇文章 0 订阅

订阅专栏

对于一个二次函数 $f(x) = ax^2 + bx + c$ ，积分得 $\displaystyle\int_0^x f(t) \, \mathrm{d}t = \dfrac{a}{3}x^3 + \dfrac{b}{2}x^2 + cx + C$ 。

于是

$\displaystyle\int_l^r f(x) \, \mathrm{d}x = F(r) - F(l)$

$\dfrac{a(r^3 - l^3)}{3} + \dfrac{b(r^2 - l^2)}{2} + c(r - l) + (C - C)$

$l)[\dfrac{a(r^2 + rl + l^2)}{3} + \dfrac{b(r + l)}{2} + c]$

$\dfrac{r - l}{6}(2ar^2 + 2arl + 2al^2 + 3br + 3bl + 6c)$

$\dfrac{r - l}{6}\{(ar^2 + br + c) + (al^2 + bl + c) + 4[a(\dfrac{l + r}{2})^2 + b(\dfrac{l + r}{2}) + c]\}$

$\dfrac{(r - l)(f(l) + 4f(\dfrac{l + r}{2}) + f(r))}{6}$

以上即为辛普森公式。

普通辛普森法

为了求 $\displaystyle\int^r_l f(x) \, \mathrm{d}x$ ，可以将区间 $[l, r]$ 分为 $2 n$ 个等长的区间 $[g_i, g_{i + 1}](i \in \{1, 2, ..., 2n\})$ 。设 $\dfrac{r - l}{2n}$ ，近似地计算过 $g_{i - 1}, f(g_{i - 1}))$ 、 $g_i, f(g_i))$ 、 $g_{i + 1}, f(g_{i + 1}))$ 三点的二次函数在区间 $g_{i - 1}, g_{i + 1}]$ 的积分。当 $n$ 足够大时，计算结果便越接近原函数的积分。设 $P (x)$ 为过该三点的函数，则

$\displaystyle\int^r_l f(x) \, \mathrm{d}x \approx \displaystyle\sum_{i = 1}^{n} \displaystyle\int^{g_{2i + 1}}_{g_{2i - 1}} P(x) \, \mathrm{d}x$

$\dfrac{h}{3} (\displaystyle\sum^{n}_{i = 1} f(g_{2i - 1}) + 4f(g_{2i}) + f(g_{2i + 1}))$

$\dfrac{h}{3}(f(g_1) + 4f(g_2) + 2f(g_3) + 4f(g_4) + ... + 2f(g_{2n - 1}) + 4f(g_{2n}) + f(g_{2n+1}))$

Code:

int n = 1e7;

double calc(double l, double r) {
    double h = (r - l) / (2 * n), res = f(l) + f(r);
    
    for (int i = 2; i < 2 * n; ++i)
        res += f(l + i * h) * (i & 1? 2: 4);
    
    return res * h / 3;
}

自适应辛普森法

考虑到当分割出的二次函数在该区间内已经足够接近原函数时，可以直接使用该积分值；而当二次函数仍不够接近时还需再分割，于是可以考虑分别用辛普森公式计算区间 $[l, r]$ 、 $\dfrac{l + r}{2}]$ 、 $[\dfrac{l + r}{2}, r]$ 的积分值，若大区间的积分值与小区间的积分值和的差足够小，则可以认为二次函数的图像在该区间内足够接近原函数的图像，并直接返回大区间的积分值。否则，递归计算两个小区间。一般还需设定一个次数，当已经递归超过该次数时才开始判断是否可以直接返回。

P4525 【模板】自适应辛普森法 1 AC Code:

#include <cstdio>
#define simpson(l, r) ((r - l) * (f(l) + f(r) + 4 * f((l + r) / 2)) / 6)
#define f(x) (double(c * (x) + d) / (a * (x) + b))
#define abs(x) ((x) >= 0? (x): -(x))

using namespace std;

double a, b, c, d, L, R;

double calc(double l, double r, int step, double sps) {
    double mid = (l + r) / 2;
    double fl = simpson(l, mid), fr = simpson(mid, r);
    
    if (abs(fl + fr - sps) <= 1e-6 && step < 0)
        return sps;
    else
        return calc(l, mid, step - 1, fl) + calc(mid, r, step - 1, fr);
}

int main() {
    scanf("%lf %lf %lf %lf %lf %lf", &a, &b, &c, &d, &L, &R);
    printf("%lf", calc(L, R, 12, simpson(L, R)));
    
    return 0;
}