simpson法求积分 专题练习

时间：2022-04-24 20:33:32

【xsy1775】数值积分

题意

多组询问，求\(\int_l^r\sqrt{a(1-{x^2\over b})}dx\)

分析

double f(double x) {
    return sqrt(a*(1-x*x/b));
}

double Get(double l,double r) {
    return (r-l)*(f(l)+f(r)+4*f((l+r)/2))/6;
}

double Calc(double l,double r) {
    double gx=Get(l,r);
    double mid=(l+r)/2,gl=Get(l,mid),gr=Get(mid,r),gs=gl+gr;
    if (!cmp(gx-gs))
        return gs;
    else {
        double rl=Calc(l,mid),rr=Calc(mid,r),rs=rl+rr;
        return rs;
    }
}

【bzoj2178】圆的面积并

//...

int cmp(double x);

int n,flag[N];
struct Cir {
    int x,y,r;
    Cir(int _x=0,int _y=0,int _r=0) {
        x=_x,y=_y,r=_r;
    }
    friend int operator < (Cir a,Cir b) {
        return a.r<b.r;
    }
}cir[N];

int tot;
struct Line {
    double l,r;
    Line(double _l=0,double _r=0) {
        l=_l,r=_r;
    }
    friend int operator < (Line a,Line b) {
        return cmp(a.l-b.l)<0;
    }
}line[N];

double res;

//...

int cmp(double x) {
    if (fabs(x)<EPS) return 0;
    return x<0?-1:1;
}

int Contain(Cir a,Cir b) {
    int dx=a.x-b.x,dy=a.y-b.y,dr=a.r-b.r;
    return dx*dx+dy*dy<=dr*dr;
}

#define ld double
#define get f
#define F rep
#define p line
#define segment Line
#define eps EPS
#define a cir
#define inf 1e9
inline bool cmp2(segment a,segment b) {
    return a.l<b.l;
}
inline ld get(ld x) {
    int cnt=0;
    F(i,1,n) if (cmp(fabs(x-a[i].x)-a[i].r)<=0) {
        ld tmp=sqrt(pow(a[i].r,2)-pow(x-a[i].x,2));
        if (!cmp(tmp)) continue;
        p[++cnt]=(segment) {
            a[i].y-tmp,a[i].y+tmp
        };
    }
    sort(p+1,p+cnt+1,cmp2);

    /*
    ld h=-inf,ans=0;
    F(i,1,cnt) {
        if (h<p[i].l) ans+=p[i].r-p[i].l,h=p[i].r;
        else if (h<p[i].r) ans+=p[i].r-h,h=p[i].r;
    }
    return ans;
    */

    double totLen=0,mxp=MIN;
    rep(i,1,cnt) {
//      if (cmp(line[i].r-mxp)<=0) continue;
        if (mxp<p[i].l/*cmp(line[i].l-mxp)>0*/)
            totLen+=line[i].r-line[i].l,mxp=line[i].r;
        else if (mxp<p[i].r/*cmp(mxp-p[i].r)<0*/)
            totLen+=line[i].r-mxp,mxp=line[i].r;
//      mxp=line[i].r;
    }
    return totLen;
}
#undef ld
#undef get
#undef F
#undef p
#undef segment
#undef eps
#undef a
#undef inf

double Get(double l,double r) {
    return (r-l)*(f(l)/6+f(r)/6+f((l+r)/2)*2/3);
}

double Calc(double l,double r) {
    double gx=Get(l,r);
    double mid=(l+r)/2;
    double gl=Get(l,mid);
    double gr=Get(mid,r);
    double gs=gl+gr;
    if (!cmp(gx-gs))
        return gs;
    else {
        double rl=Calc(l,mid);
        double rr=Calc(mid,r);
        double rs=rl+rr;
        return rs;
    }
}

int main(void) {
    //...

    sort(cir+1,cir+n+1);
    rep(i,1,n) rep(j,i+1,n)
        if (Contain(cir[j],cir[i])) {
            flag[i]=1;
            break;
        }
    int cur=0;
    rep(i,1,n) if (!flag[i])
        cir[++cur]=cir[i];
    n=cur;

    res=Calc(L,R);
    printf("%.3lf\n",res);

    //...
}

小结

对于这类求面积的问题，意味着它满足连续性，可以使用积分来求解。

【hdu1724】Ellipse

double f(double x) {
    double y2=(1-pow(x,2)/pow(a,2))*pow(b,2);
    if (cmp(y2)<0) return 0; else return 2*sqrt(y2);
}

double Get(double l,double r) {
    return (r-l)*(f(l)+f(r)+4*f((l+r)/2))/6;
}

double Calc(double l,double r) {
    double mid=(l+r)/2;
    double gl=Get(l,mid);
    double gr=Get(mid,r);
    double gs=Get(l,r);
    if (!cmp(gs-(gl+gr)))
        return gl+gr;
    else {
        double rl=Calc(l,mid);
        double rr=Calc(mid,r);
        return rl+rr;
    }
}

相关文章