自适应辛普森法从入门到进门

前言

学数学学的。

simpson

背景

我们要计算这样一个式子：

\[\int_l^r f(x)\text dx \]

显然计算机是很难把柿子推出来的。

函数的拟合

对于一个奇怪的函数，为了对其求导，我们可以用一个图像近似且容易求导的函数来替代，这个过程叫做拟合。

这里我们用二次函数来替代，那么有：

\[\begin{aligned}\int_l^r f(x)\text dx&\approx\int_l^r(ax^2+bx+c)\text dx\\&=\dfrac{a}{3}(r^3-l^3)+\dfrac{b}{2}(r^2-l^2)+c(r-l)\\&=\dfrac{r-l}{6}\left[2a(r^2+rl+l^2)+3b(r+l)+6c\right]\\&=\dfrac{r-l}{6}(2ar^2+2arl+2al^2+3br+3bl+6c)\\&=\dfrac{r-l}{6}\left[(ar^2+br+c)+(al^2+bl+c)+4a\left(\dfrac{r+l}{2}\right)+4b\left(\dfrac{r+l}{2}\right)+4c\right]\\&=\dfrac{r-l}{6}\left[f(r)+f(l)+4f\left(\dfrac{r+l}{2}\right)\right]\end{aligned} \]

code

inline double simpson(double l,double r){
    return (r-l)*(f(r)+f(l)+4.0*f((l+r)/2))/6.0;
}

显然这样的话精度误差会很大，因此我们需要将函数分段拟合。

为了确定合适的分段，我们需要一个“自适应”的操作。

如果这一段函数与拟合的二次函数相似，那么直接把这一段当做这个二次函数，然后套用公式计算；否则应该把这一段分成两半，并递归进行这一过程。

关于函数相似的判断，我们可以将区间的左半和右半部分分别进行计算，与计算的相差小于设定的 \(\varepsilon\) 即可返回值。

在分治判断的时候，除了判断精度是否正确，一般还要强制执行最少的迭代次数。

code

inline double work(double l,double r,double eps,double ans,int depth){
    double mid=(l+r)/2,x=simpson(l,mid),y=simpson(mid,r);
    if(fabs(x+y-ans)<=15*eps and depth<0)
        return (x+y)+(x+y-ans)/15.0;
    return work(l,mid,eps/2,x,depth-1)+work(mid,r,eps/2,y,depth-1);
}

一些题目

P4525 自适应辛普森法 1

传送门

Description

求：

\[\int_l^r \dfrac{cx+d}{ax+b}\text dx \]

Solution

套板子。

code

inline double f(double x){
    return (c*x+d)/(a*x+b);
}

inline double simpson(double l,double r){
    return (r-l)*(f(r)+f(l)+4.0*f((l+r)/2))/6.0;
}

inline double work(double l,double r,double eps,double ans,int depth){
    double mid=(l+r)/2,x=simpson(l,mid),y=simpson(mid,r);
    if(fabs(x+y-ans)<=15*eps and depth<0)
        return (x+y)+(x+y-ans)/15.0;
    return work(l,mid,eps/2,x,depth-1)+work(mid,r,eps/2,y,depth-1);
}

signed main(){
    scanf("%lf%lf%lf%lf%lf%lf",&a,&b,&c,&d,&l,&r);
    printf("%.6f",work(l,r,1e-7,simpson(l,r),12));
    return 0;
}

P4526 自适应辛普森法 2

传送门

Description

求：

\[\int_0^{+\infty} x^{\frac{a}{x}-x}\text dx \]

Solution

注意到当 \(a<0\) 时积分发散，\(a>0\) 时积分收敛。

特判后选一个适合的区间计算即可。

code

inline double f(double x){
    return pow(x,a/x-x);
}

inline double simpson(double l,double r){
    return (r-l)*(f(r)+f(l)+4.0*f((l+r)/2))/6.0;
}

inline double work(double l,double r,double eps,double ans,int depth){
    double mid=(l+r)/2,x=simpson(l,mid),y=simpson(mid,r);
    if(fabs(x+y-ans)<=15*eps and depth<0)
        return (x+y)+(x+y-ans)/15.0;
    return work(l,mid,eps/2,x,depth-1)+work(mid,r,eps/2,y,depth-1);
}

signed main(){
    scanf("%lf",&a);
    if(a<0)
        puts("orz");
    else{
        double l=1e-7,r=23.3;
        printf("%.5lf",work(l,r,1e-7,simpson(l,r),12));
    }
    return 0;
}

HDU-1724 Ellipse

传送门

Description

给定 \(a,b,l,r\)，求椭圆：\(\dfrac{x^2}{a^2}+\dfrac{y^2}{b^2}=1\) 与直线：\(x=l,x=r\) 所组成的封闭图形的面积。

Solution

还是套板子。

code

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

inline double simpson(double l,double r){
    return (r-l)*(f(r)+f(l)+4.0*f((l+r)/2))/6.0;
}

inline double work(double l,double r,double eps,double ans,int depth){
    double mid=(l+r)/2,x=simpson(l,mid),y=simpson(mid,r);
    if(fabs(x+y-ans)<=15*eps and depth<0)
        return (x+y)+(x+y-ans)/15.0;
    return work(l,mid,eps/2,x,depth-1)+work(mid,r,eps/2,y,depth-1);
}

signed main(){
    int T=read();
    while(T--){
        scanf("%lf%lf%lf%lf",&a,&b,&l,&r);
        printf("%.3f",2.0*work(l,r,1e-7,simpson(l,r),12));puts("");
    }
    return 0;
}

P4207 [NOI2005] 月下柠檬树

传送门

Description

给定一棵由圆台、圆锥组成的树和平行光线的角度，求投影面积。

Solution

首先注意到若点 \(p\) 的坐标为 \((0,x)\)，那么投影后的坐标为 \(p'(x\cot \alpha,0)\)。

圆的平行投影还是圆，圆的半径投影后不会改变，只有高度与圆台投影形成的梯形尺寸会改变。这一部分可以画图用相似解决。

把圆和梯形的数据确定后，直接套板子即可。

code

inline double sq(double x){return x*x;}

inline int sgn(double x){
    if(fabs(x)<1e-8)
        return 0;
    return x<0?-1:1;
}

struct circle{
    double x,y,r;
    circle(){}
    circle(double _x,double _y,double _r):x(_x),y(_y),r(_r){}
}cir[N];

struct Trape{
    double l,r,al,ar;
    Trape(){}
    Trape(double _l,double _r,double _al,double _ar):l(_l),r(_r),al(_al),ar(_ar){}
}trape[N];

inline double f(double x){
    double res=0;
    for(int i=1;i<=n;i++)
        if(sgn(cir[i].r-fabs(x-cir[i].x))>0){
            double tmp=sqrt(sq(cir[i].r)-sq(fabs(x-cir[i].x)));
            res=max(res,tmp*2);
        }
    for(int i=1;i<=n;i++)
        if(sgn(x-trape[i].l)>0 and sgn(trape[i].r-x)>0){
            double d=x-trape[i].l,len=trape[i].r-trape[i].l;
            double Len=len*trape[i].al/(trape[i].al-trape[i].ar);
            double tmp=(Len-d)/Len*trape[i].al;
            res=max(res,tmp*2);
        }
    return res;
}

inline double simpson(double l,double r){
    return (r-l)*(f(r)+f(l)+4.0*f((l+r)/2))/6.0;
}

inline double work(double l,double r,double eps,double ans,int depth){
    double mid=(l+r)/2,x=simpson(l,mid),y=simpson(mid,r);
    if(fabs(x+y-ans)<=15*eps and depth<0)
        return (x+y)+(x+y-ans)/15.0;
    return work(l,mid,eps/2,x,depth-1)+work(mid,r,eps/2,y,depth-1);
}

signed main(){
    n=read();scanf("%lf",&alpha);
    for(int i=0;i<=n;i++)
        scanf("%lf",&h[i]);
    for(int i=1;i<=n;i++)
        scanf("%lf",&r[i]);
    Tan=1.0/tan(alpha);
    for(int i=1;i<=n;i++){
        cir[i]=circle(Tan*H,0,r[i]);
        double len=Tan*h[i],ll=(r[i+1]-r[i])*r[i]/len,rr=(r[i+1]-r[i])*r[i+1]/len;
        trape[i]=Trape(Tan*H-ll,Tan*(H+h[i])-rr,sqrt(sq(r[i])-sq(ll)),sqrt(sq(r[i+1])-sq(rr)));
        H+=h[i];
        L=min(L,cir[i].x-cir[i].r);
        R=max(R,cir[i].x+cir[i].r);
    }
    R=max(R,Tan*H);
    printf("%.2lf",work(L,R,5e-4,simpson(L,R),12));
    return 0;
}

BZOJ2178 圆的面积并

传送门

Description

给出 \(n\) 个圆，求其面积并。

Solution

套板子，但有几个优化：

将被完全包含的圆删除。
将圆按 \(x\) 排序，每个区间分别算积分。

code

struct circle{
    double x,y,r;
}cir[N];

inline bool cmp(circle a,circle b){
    return a.x-a.r<b.x-b.r;
}

struct Segment{
    double l,r;
    bool operator <(const Segment yy) const {
        return l<yy.l;
    }
}tmp[1005];
int ts;
map<double,double>F;

inline double f(double x){
    if(F.count(x))
        return F[x];
    ts=0;
    for(int i=1;i<=n;i++){
        if(cir[i].r<=fabs(cir[i].x-x))
            continue;
        double ff=cir[i].r*cir[i].r-(cir[i].x-x)*(cir[i].x-x);
        if(ff<=1e-4)
            continue;
        ff=sqrt(ff);
        tmp[++ts]=(Segment){cir[i].y-ff,cir[i].y+ff};
    }
    sort(tmp+1,tmp+ts+1);
    double nows=0,nowr=-1e9;
    for(int i=1;i<=ts;i++){
        if(tmp[i].r<nowr)
            continue;
        if(tmp[i].l<=nowr){
            nows+=tmp[i].r-nowr;
            nowr=tmp[i].r;
        }else{
            nows+=tmp[i].r-tmp[i].l;
            nowr=tmp[i].r;
        }
    }
    return F[x]=nows;
}

inline double simpson(double l,double r){
    return (r-l)*(f(r)+f(l)+4.0*f((l+r)/2))/6.0;
}

inline double work(double l,double r,double eps,double ans,int depth){
    double mid=(l+r)/2,x=simpson(l,mid),y=simpson(mid,r);
    if(fabs(x+y-ans)<=15*eps and depth<0)
        return (x+y)+(x+y-ans)/15.0;
    return work(l,mid,eps/2,x,depth-1)+work(mid,r,eps/2,y,depth-1);
}

signed main(){
    n=read();
    for(int i=1;i<=n;i++)
        cin>>cir[i].x>>cir[i].y>>cir[i].r;
    sort(cir+1,cir+n+1,cmp);
    double nowr=-inf,ans=0;
    for(int i=1;i<=n+1;i++){
        if(cir[i].x-cir[i].r>nowr){
            ans+=work(cir[i].x-cir[i].r,cir[i].x+cir[i].r,1e-5,simpson(cir[i].x-cir[i].r,cir[i].x+cir[i].r),12);
            nowr=cir[i].x+cir[i].r;
        }else if(cir[i].x+cir[i].r>nowr){
            ans+=work(nowr,cir[i].x+cir[i].r,1e-4,simpson(nowr,cir[i].x+cir[i].r),12);
            nowr=cir[i].r+cir[i].x;
        }
    }
    printf("%.3lf\n",ans);
    return 0;
}

其他面积并的东西都是同理的，懒得粘了。

后记

乱搞的东西你为什么要学？

posted @ 2024-02-16 14:51 QcpyWcpyQ 阅读(35) 评论(0) 编辑收藏举报

刷新页面返回顶部

QcpyWcpyQ

自适应辛普森法从入门到进门

前言

simpson

背景

函数的拟合

code

code

一些题目

P4525 自适应辛普森法 1

Description

Solution

code

P4526 自适应辛普森法 2

Description

Solution

code

HDU-1724 Ellipse

Description

Solution

code

P4207 [NOI2005] 月下柠檬树

Description

Solution

code

BZOJ2178 圆的面积并

Description

Solution

code

后记

公告