• 欢迎光临~

自适应辛普森法 学习笔记

开发技术 开发技术 2022-10-16 次浏览

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

于是

[displaystyleint_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) ]

[= (r - 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} ]

以上即为辛普森公式。

普通辛普森法

为了求 (displaystyleint^r_l f(x) , mathrm{d}x),可以将区间 ([l, r]) 分为(2n)个等长的区间 ([g_i, g_{i + 1}](i in {1, 2, ..., 2n}))。设 (h = 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)) 为过该三点的函数,则

[displaystyleint^r_l f(x) , mathrm{d}x approx displaystylesum_{i = 1}^{n} displaystyleint^{g_{2i + 1}}_{g_{2i - 1}} P(x) , mathrm{d}x ]

[= dfrac{h}{3} (displaystylesum^{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])([l, 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;
}
程序员灯塔
转载请注明原文链接:自适应辛普森法 学习笔记
喜欢 (0)