算法筆記--辛普森積分公式

一圖勝千言。

自適應辛普森模板：

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;
}