zoukankan      html  css  js  c++  java
  • 算法笔记--辛普森积分公式

    一图胜千言。

    自适应辛普森模板:

    const double eps=1e-10;
    double f(double x)
    {
        //你的函数。 
    } 
    double simpson(double l,double r)
    {
        return (f(l)+4.0*f((l+r)/2.0)+f(r))/6.0*(r-l);
    }
    double integral(double l,double r)
    {
        double mid=(l+r)/2.0;
        double res=simpson(l,r);
        if(fabs(res-simpson(l,mid)-simpson(mid,r))<eps)return res;
        else return integral(l,mid)+integral(mid,r);
    }

    复合辛普森模板:

    double f(double x) {
        double l = max(s1, x-w);
        double r = min(s2, x+w);
        if(r-l < eps) return eps;
        else return r-l;
    }
    double simpson(double l, double r) {
        return (f(l) + 4.0*f((l+r)/2.0) + f(r))/6.0*(r-l);
    }
    double integral(double l, double r) {
        double m = (l+r)/2.0;
        double res = simpson(l, r);
        if(fabs(res-simpson(l, m)-simpson(m, r)) < eps) return res;
        else return integral(l, m) + integral(m, r);
    }
    double fuhesimpson(double l, double r, int n) {
        double h = (r-l)/n, ans = 0;
        for (int i = 0; i < n; ++i) {
            ans += integral(l+i*h, l+(i+1)*h);
        }
        return ans;
    }

    参考:https://blog.csdn.net/tengweitw/article/details/43311685 

    例题1:HDU 1724 Ellipse

    代码:

    #include<bits/stdc++.h>
    using namespace std;
    #define ll long long
    #define pb push_back
    #define pi acos(-1.0)
    #define mem(a,b) memset(a,b,sizeof(a))
    
    const double eps=1e-10;
    double a,b;
    double f(double x)
    {
        return b*sqrt(1.0-(x*x)/(a*a)); 
    } 
    double simpson(double l,double r)
    {
        return (f(l)+4.0*f((l+r)/2.0)+f(r))/6.0*(r-l);
    }
    double integral(double l,double r)
    {
        double mid=(l+r)/2.0;
        double res=simpson(l,r);
        if(fabs(res-simpson(l,mid)-simpson(mid,r))<eps)return res;
        else return integral(l,mid)+integral(mid,r);
    }
    
    int main()
    {
        ios::sync_with_stdio(false);
        cin.tie(0);
        int T;
        cin>>T;
        double l,r;
        while(T--)
        {
            cin>>a>>b>>l>>r;
            cout<<fixed<<setprecision(3)<<2*integral(l,r)<<endl;
        }
        return 0;
    }
    View Code
  • 相关阅读:
    627. Swap Salary
    176. Second Highest Salary
    596. Classes More Than 5 Students
    183. Customers Who Never Order
    181. Employees Earning More Than Their Managers
    182. Duplicate Emails
    175. Combine Two Tables
    620. Not Boring Movies
    595. Big Countries
    HDU 6034 Balala Power! (贪心+坑题)
  • 原文地址:https://www.cnblogs.com/widsom/p/7694882.html
Copyright © 2011-2022 走看看