辛普森积分学习笔记

白书上也有讲，就是一个在精度要求不高的时候计算积分的算法。

我感觉其实功利一点的话，不用完全理解它的原理，大概明白然后会套板子做题即可（蒟蒻口胡）。

洛谷P4525

模板题，直接上模板。

#include<bits/stdc++.h>
using namespace std;
const double eps=1e-7;
double a,b,c,d,L,R;

double f(double x) {  //具体计算公式 
    return (c*x+d)/(a*x+b);
}

double simpson(double l,double r) {  //辛普森积分公式 
    return (f(l)+f(r)+4*f((l+r)/2))*(r-l)/6;
}

double asr(double l,double r,double exps,double val) {  //自适应辛普森积分 
    double mid=(l+r)/2;
    double lval=simpson(l,mid),rval=simpson(mid,r);
    if (fabs(lval+rval-val)<=15*exps) return lval+rval+(lval+rval-val)/15;
    else return asr(l,mid,exps/2,lval)+asr(mid,r,exps/2,rval);
}

int main()
{
    scanf("%lf%lf%lf%lf%lf%lf",&a,&b,&c,&d,&L,&R);
    printf("%.6lf
",asr(L,R,eps,simpson(L,R)));
    return 0;
}

View Code

洛谷P4526

另一道模板题，不过辛普森积分好像不适用于发散函数和反常积分（积分区间无穷大）。所以我们要想办法处理这两种情况。①可以证明a<0函数发散,a>0函数收敛。②x到底大于某个数之后变得非常非常小，于是之后的就可以不积分了。于是这题就可以做了。

#include<bits/stdc++.h>
using namespace std;
const double eps=1e-7;
double a;

double f(double x) {  //具体计算公式 
    return pow(x,a/x-x);
}

double simpson(double l,double r) {  //辛普森积分公式 
    return (f(l)+f(r)+4*f((l+r)/2))*(r-l)/6;
}

double asr(double l,double r,double exps,double val) {  //自适应辛普森积分 
    double mid=(l+r)/2;
    double lval=simpson(l,mid),rval=simpson(mid,r);
    if (fabs(lval+rval-val)<=15*exps) return lval+rval+(lval+rval-val)/15;
    else return asr(l,mid,exps/2,lval)+asr(mid,r,exps/2,rval);
}

int main()
{
    scanf("%lf",&a);
    if (a<0) puts("orz");
    else printf("%.5lf
",asr(1e-10,20,eps,simpson(1e-10,20)));
    return 0;
}

View Code

HDU-1724

求椭圆面积一部分。精度要求很低，用辛普森积分。

#include<bits/stdc++.h>
using namespace std;
const double eps=1e-7;
double a,b,l,r;

double f(double x) {  //具体计算公式 
    return 2*b*sqrt(1-x*x/a/a);
}

double simpson(double l,double r) {  //辛普森积分公式 
    return (f(l)+f(r)+4*f((l+r)/2))*(r-l)/6;
}

double asr(double l,double r,double exps,double val) {  //自适应辛普森积分 
    double mid=(l+r)/2;
    double lval=simpson(l,mid),rval=simpson(mid,r);
    if (fabs(lval+rval-val)<=15*exps) return lval+rval+(lval+rval-val)/15;
    else return asr(l,mid,exps/2,lval)+asr(mid,r,exps/2,rval);
}

int main()
{
    int T; cin>>T;
    while (T--) {
        scanf("%lf%lf%lf%lf",&a,&b,&l,&r);
        printf("%.3lf
",asr(l,r,eps,simpson(l,r)));
    }
    return 0;
}

View Code

BZOJ-2178

求圆的面积并。直接在x轴上积分，然后f函数就是直线经过的线段长度。其他还是套辛普森积分即可。

但是这题有非常多细节：第一，要先把被包含的圆去掉。第二，可以看出每次计算f函数都要耗费大量时间，所以要想办法通过尽量复用f函数值来减少计算次数。第三，这也是蒟蒻没想懂的一点，好像用我上面板子的控制精度方法是会超时的（即每次二分将精度也除以2的方法），这题可以用一个确切的精度1e-13然后一直用这个精度控制即可。

#include<bits/stdc++.h>
using namespace std;
const double eps=1e-13;
const int N=1e3+10;
int n;
double x[N],y[N],r[N];
bool ban[N];

struct dat{
    double x,y;
    bool operator < (const dat &rhs) const {
        return x<rhs.x;
    }
}A[N<<1];
double f(double P) {  //计算直线P经过线段长度 
    int num=0;
    for (int i=1;i<=n;i++) {
        if (fabs(x[i]-P)>=r[i]) continue;
        double d=fabs(x[i]-P);
        double l=sqrt(r[i]*r[i]-d*d);
        A[++num].x=y[i]-l; A[num].y=y[i]+l;
    }
    sort(A+1,A+num+1);
    double ret=0,lst=-0x3f3f3f3f;
    for (int i=1;i<=num;i++)
        if (A[i].x>lst) ret+=A[i].y-A[i].x,lst=A[i].y;
        else if (A[i].y>lst) ret+=A[i].y-lst,lst=A[i].y;
    return ret;
}

double simpson(double l,double r) {  //辛普森积分公式 
    return (f(l)+f(r)+4*f((l+r)/2))*(r-l)/6;
}

double cal(double l,double r,double fl,double fr,double fm){
    return (fl+fr+4*fm)*(r-l)/6;
}

double asr(double l,double r,double val,double L,double R,double M) {  //自适应辛普森积分 
    double mid=(l+r)/2,LM=f((l+mid)/2),RM=f((mid+r)/2);
    double lval=cal(l,mid,L,M,LM),rval=cal(mid,r,M,R,RM);
    if (fabs(lval+rval-val)<=eps) return lval+rval+(lval+rval-val)/15;
    return asr(l,mid,lval,L,M,LM)+asr(mid,r,rval,M,R,RM);
}

int main()
{
    cin>>n;
    double L=0x3f3f3f3f,R=-0x3f3f3f3f;
    for (int i=1;i<=n;i++) {
        scanf("%lf%lf%lf",&x[i],&y[i],&r[i]);
        L=min(L,x[i]-r[i]);
        R=max(R,x[i]+r[i]);
    }
    
    for (int i=1;i<=n;i++)
        for (int j=i+1;j<=n;j++) 
        if (!ban[i] && !ban[j]) {
            double dist=sqrt((x[i]-x[j])*(x[i]-x[j])+(y[i]-y[j])*(y[i]-y[j]));
            if (dist>fabs(r[i]-r[j])) continue;
            if (r[j]<=r[i]) ban[j]=1; else ban[i]=1;
        }
    int t=0;
    for (int i=1;i<=n;i++)
        if (!ban[i]) { t++; x[t]=x[i]; y[t]=y[i]; r[t]=r[i]; }
    n=t;        
    
    double M=(L+R)/2;
    printf("%.3lf
",asr(L,R,simpson(L,R),f(L),f(R),f(M)));
    return 0;
}

View Code

洛谷P4207/BZOJ-1502