zoukankan      html  css  js  c++  java
  • 自适应Simpson积分

    用$g(x)=A*x^2+B*x+C$来代替原函数$f(x)$,设$g(x)$原函数为$G(x)$,显然$G(x)=frac{1}{3}*A*x^3+frac{1}{2}*B*x^2+C*x$

    $int_a^b f(x) dx approx int_a^b g(x) dx=G(b)-G(a)$

    代入化简得

    $$int_a^b f(x) dx approx frac{a-b}{6}*[g(a)+g(b)+4*g(frac{a+b}{2})]$$

    二分区间,然后计算区间积分,通过期望容差来控制二分

    如果$|S(a,mid)+S(mid,b)-S(a,b)|<15*varepsilon $,则中止二分,并且认为

    $$S(a,b)=S(a,m)+S(m,b)+frac{S(a,mid)+S(mid,b)-S(a,b)}{15}$$

    否则继续二分

    其中$mid=frac{a+b}{2},varepsilon为期望容差$

    参考博客:https://www.cnblogs.com/chy-2003/p/11644112.html

    例题:

    luogu 4525 【模板】自适应辛普森法1

    #include <iostream>
    #include <algorithm>
    #include <cstring>
    #include <cstdio>
    #include <cmath>
    
    using namespace std;
    
    const double INF = 20;
    
    double a, b, c, d, l, r;
    
    inline double sqr(double x)
    {
        return x * x;
    }
    
    double f(double x)
    {
        return (c * x + d) / (a * x + b);
    }
    
    double calc(double l, double r)
    {
        double mid = (l + r) / 2;
        return (f(l) + 4 * f(mid) + f(r)) * (r - l) / 6;
    }
    
    double simpson(double l, double r, double eps)
    {
        double mid = (l + r) / 2;
        double st = calc(l, r), sl = calc(l, mid), sr = calc(mid, r);
        if (fabs(sl + sr - st) <= 15 * eps) return sl + sr + (sl + sr - st) / 15;
        return simpson(l, mid, eps / 2) + simpson(mid, r, eps / 2);
    }
    
    int main()
    {
        // freopen("in.txt", "r", stdin);
        // freopen("out.txt", "w", stdout);
        scanf("%lf%lf%lf%lf%lf%lf", &a, &b, &c, &d, &l, &r);
        printf("%.6lf
    ", simpson(l, r, 4e-7));
        return 0;
    }

    luogu 4526 【模板】自适应辛普森法2

    #include <iostream>
    #include <algorithm>
    #include <cstring>
    #include <cstdio>
    #include <cmath>
    
    using namespace std;
    
    const double INF = 20;
    
    double a;
    
    inline double sqr(double x)
    {
        return x * x;
    }
    
    double f(double x)
    {
        return pow(x, a / x - x);
    }
    
    double calc(double l, double r)
    {
        double mid = (l + r) / 2;
        return (f(l) + 4 * f(mid) + f(r)) * (r - l) / 6;
    }
    
    double simpson(double l, double r, double eps)
    {
        double mid = (l + r) / 2;
        double st = calc(l, r), sl = calc(l, mid), sr = calc(mid, r);
        if (fabs(sl + sr - st) <= 15 * eps) return sl + sr + (sl + sr - st) / 15;
        return simpson(l, mid, eps / 2) + simpson(mid, r, eps / 2);
    }
    
    int main()
    {
        // freopen("in.txt", "r", stdin);
        // freopen("out.txt", "w", stdout);
        scanf("%lf", &a);
        if (a < 0) printf("orz
    ");
        else printf("%.5lf
    ", simpson(4e-7, INF, 4e-7));
        return 0;
    }

    hdu 1724 Ellipse

    #include <iostream>
    #include <algorithm>
    #include <cstring>
    #include <cstdio>
    #include <cmath>
    
    using namespace std;
    
    int T;
    double a, b, l, r;
    
    inline double sqr(double x)
    {
        return x * x;
    }
    
    double f(double x)
    {
        double t = sqr(b) - sqr(b) / sqr(a) * sqr(x);
        return 2 * sqrt(t);
    }
    
    double calc(double l, double r)
    {
        double mid = (l + r) / 2;
        return (f(l) + 4 * f(mid) + f(r)) * (r - l) / 6;
    }
    
    double simpson(double l, double r, double eps)
    {
        double mid = (l + r) / 2;
        double st = calc(l, r), sl = calc(l, mid), sr = calc(mid, r);
        if (fabs(sl + sr - st) <= 15 * eps) return sl + sr + (sl + sr - st) / 15;
        return simpson(l, mid, eps / 2) + simpson(mid, r, eps / 2);
    }
    
    int main()
    {
        // freopen("in.txt", "r", stdin);
        // freopen("out.txt", "w", stdout);
        scanf("%d", &T);
        while (T--) {
            scanf("%lf%lf%lf%lf", &a, &b, &l, &r);
            printf("%.3lf
    ", simpson(l, r, 1e-6));
        }
        return 0;
    }
  • 相关阅读:
    NoHttp开源Android网络框架1.0.0之架构分析
    3种浏览器性能測试
    自己定义控件-画板,橡皮擦,刮刮乐
    android优化 清除无效代码 UCDetector
    iOS推送 (百度推送)
    C#中的协变OUT和逆变
    使用反射构造对象实例并动态调用方法
    用反射获取构造函数带参数的实例对象
    自己实现一个IOC(控制翻转,DI依赖注入)容器
    func 和action 委托的使用
  • 原文地址:https://www.cnblogs.com/zzzzzzy/p/13502580.html
Copyright © 2011-2022 走看看