NEERC 2014 D题 Damage Assessment

罗博士

于 2020-11-02 22:26:07 发布

阅读量339

点赞数 2

分类专栏： ACM/ICPC ACM高等数学

本文链接：https://blog.csdn.net/u012061345/article/details/109456172

版权

ACM/ICPC 同时被 2 个专栏收录

55 篇文章 8 订阅

订阅专栏

ACM高等数学

4 篇文章 0 订阅

订阅专栏

NEERC 2014 D题 Damage Assessment

题目描述
直接计算定积分
辛普森方法
Romberg方法

题目描述

这里是NEERC2014的相关资料，包括榜、题目、标程，但是好像没有解题报告。
这是CF上的重现Gym100553。
当然在ICPC Live Archive上也有，这里是链接，但是好像不能提交了。所以去CF上比较好。
牛客网上也有这道题，有的代码两边都能过，但有的代码只能过其中一边。数据应该是一样的，也许是SPJ与时限写的不同。
如图，一个油罐，两端是球缺，中间是一个圆柱体。倾斜以后，再给定一个液面高度，求液面的体积。
df23
按如下示意建立一个坐标系，然后求一个定积分即可（定积分的计算可以看这里）
$\int_{-h_{球缺}}^{l+h_{球缺}}f(x)dx$
其中，积分区域跟球缺的高度以及 $l$ 有关，而 $f (x)$ 则是 $x$ 处液体的横截面积，该面积要么为 $0$ ，要么必然为一个弓形（当然整圆也可）。
在这里插入图片描述

直接计算定积分

这是CF上可以AC的代码，但是在牛客网T了两个点。

#include <bits/stdc++.h>
using namespace std;

double const EPS = 3E-3;
double const PI = acos(-1.0);
double const INF = 1E100;

inline int sgn(double x){return x>EPS?1:(x<-EPS?-1:0);}
inline bool is0(double x){return 0==sgn(x);}
inline bool isEq(double x,double y){return is0(x-y);}
inline double myasin(double x){
    if(x<-1.0)x=-1.0;
    else if(x>1.0)x=1.0;
    return asin(x);
}
inline double myacos(double x){
    if(x<-1.0)x=-1.0;
    else if(x>1.0)x=1.0;
    return acos(x);
}

//计算定积分
double integral(double(*f)(double),double x1,double x2){
    double sum = 0.0;
    for(double x=x1;x<=x2;x+=EPS){
        sum += f(x);
    }
    return sum * EPS;
}

//计算弓形的高度，半径为r，高度为h
inline double area(double r, double h){
    if(is0(r)||sgn(h)<=0) return 0.0;
    if(sgn(h-r-r)>=0) return PI * r * r;
    double a = r - h;
    double theta = myacos(a/r);
    return (theta - 0.5 * sin(theta+theta)) * r * r;
}

double D;//球冠底面的直径
double L;//圆柱体的长度
double R,R2;//球冠的球的半径
double T;//倾斜以后,圆柱体母线的高度,即倾斜角就是arcsin(T/L)
double H;//剩下的油,其液面高度，是从圆柱体的底端开始计算
double Theta;
double hOfQiuque,xLeft,xRight,rOfYuanzhu;
double K,B;//y=kx+B;当k是INF表示x=B

//计算x处的横截面积，这是一个弓形
double f(double x){
    double r = rOfYuanzhu;
    if(x<0){
        r = sqrt(R2-(x-xLeft)*(x-xLeft));
    }else if(x>L){
        r = sqrt(R2-(x-xRight)*(x-xRight));
    }
    double y0 = K*x + B;
    if(isEq(K,INF)) y0 = (x<=B?INF:-INF);
    return area(r,y0+r-rOfYuanzhu);
}

double proc(){
    if(isEq(T,L)) K = INF, B = H;
    else Theta = asin(T/L), K = -tan(Theta), B = H / cos(Theta);

    rOfYuanzhu = 0.5 * D;
    xLeft = sqrt(R2-rOfYuanzhu*rOfYuanzhu);
    hOfQiuque = R - xLeft;
    xRight = L + hOfQiuque - R;
    return integral(f,-hOfQiuque,L+hOfQiuque);
}

int main(){
    freopen("damage.in","r",stdin);
    freopen("damage.out","w",stdout);
    //freopen("1.txt","r",stdin);
    while(5==scanf("%lf%lf%lf%lf%lf",&D,&L,&R,&T,&H)){
        R2 = R * R;
        printf("%.6f\n",proc()/1E6);
    }
    return 0;
}

这个是牛客网上可以AC的，但是在CF上WA。

#include <bits/stdc++.h>
using namespace std;

double const EPS = 1E-4;
double const PI = acos(-1.0);
double const INF = 1E100;

inline int sgn(double x){return x>EPS?1:(x<-EPS?-1:0);}
inline bool is0(double x){return 0==sgn(x);}
inline bool isEq(double x,double y){return is0(x-y);}
inline double myasin(double x){
    if(x<-1.0)x=-1.0;
    else if(x>1.0)x=1.0;
    return asin(x);
}
inline double myacos(double x){
    if(x<-1.0)x=-1.0;
    else if(x>1.0)x=1.0;
    return acos(x);
}

//计算定积分
double integral(double(*f)(double),double x1,double x2){
    double sum = 0.0;
    double fx;
    for(double x=x1;x<=x2;x+=EPS){
        sum += (fx = f(x));
        if(x>=0&&is0(fx)) break;
    }
    return sum * EPS;
}

//计算弓形的高度，半径为r，高度为h
inline double area(double r, double h){
    if(is0(r)||sgn(h)<=0) return 0.0;
    if(sgn(h-r-r)>=0) return PI * r * r;
    double a = r - h;
    double theta = myacos(a/r);
    return (theta - 0.5 * sin(theta+theta)) * r * r;
}

double D;//球冠底面的直径
double L;//圆柱体的长度
double R,R2;//球冠的球的半径
double T;//倾斜以后,圆柱体母线的高度,即倾斜角就是arcsin(T/L)
double H;//剩下的油,其液面高度，是从圆柱体的底端开始计算
double Theta;
//double A,B,C;//表示直线Ax+By=C;
double hOfQiuque,xLeft,xRight,rOfYuanzhu;
double K,B;//y=kx+B;当k是INF表示x=B

//计算x处的横截面积，这是一个弓形
double f(double x){
    double r = rOfYuanzhu;
    if(x<0){
        r = sqrt(R2-(x-xLeft)*(x-xLeft));
    }else if(x>L){
        r = sqrt(R2-(x-xRight)*(x-xRight));
    }
    double y0 = K*x + B;
    if(isEq(K,INF)) y0 = (x<=B?INF:-INF);
    return area(r,y0+r-rOfYuanzhu);
}

double proc(){
    /*
    if(is0(T)) A=0.0,B=1.0,C=H;
    else if(isEq(T,L))A=1.0, B=0.0, C=H;
    else Theta = myasin(T/L),A=sin(Theta),B=cos(Theta),C=H;
        //*/
    if(isEq(T,L)) K = INF, B = H;
    else Theta = asin(T/L), K = -tan(Theta), B = H / cos(Theta);

    rOfYuanzhu = 0.5 * D;
    xLeft = sqrt(R2-rOfYuanzhu*rOfYuanzhu);
    hOfQiuque = R - xLeft;
    xRight = L + hOfQiuque - R;
    return integral(f,-hOfQiuque,L+hOfQiuque);
}

int main(){
    //freopen("1.txt","r",stdin);
    //freopen("2.txt","w",stdout);
    while(5==scanf("%lf%lf%lf%lf%lf",&D,&L,&R,&T,&H)){
        D /= 1E3, L /= 1E3, R /= 1E3,T /= 1E3, H /= 1E3;
        R2 = R * R;
        printf("%.3f\n",1E3*proc());
    }
    return 0;
}

直接计算的话，精度合适的范围很窄，不合适的话，要么就T，要么就WA。

辛普森方法

要注意写成分段积分的形式，因为如果整个积分的话，有可能存在起点、中点、终点均为 $0$ 的情况，但起点到中点之间却有一小段体积。另外，递归的时候 $e p s$ 不要乘 $0.5$ ，否则会超时。直接设定一个值比较好。这是CF上能够AC的。

#include <bits/stdc++.h>
using namespace std;

double const EPS = 1E-1;
double const PI = acos(-1.0);
double const INF = 1E100;

inline int sgn(double x){return x>EPS?1:(x<-EPS?-1:0);}
inline bool is0(double x){return 0==sgn(x);}
inline bool isEq(double x,double y){return is0(x-y);}
inline double myasin(double x){
    if(x<-1.0)x=-1.0;
    else if(x>1.0)x=1.0;
    return asin(x);
}
inline double myacos(double x){
    if(x<-1.0)x=-1.0;
    else if(x>1.0)x=1.0;
    return acos(x);
}

double simpson(double(*f)(double),double a,double b){
    double mid = ( a + b ) * 0.5;
    return (b-a)*(f(a)+f(b)+4.0*f(mid)) / 6.0;
}
//自适应辛普森积分的递归调用
double asr(double(*f)(double),double a,double b,double eps,double ans){
    double mid = ( a + b ) * 0.5;
    double lans = simpson(f,a,mid);
    double rans = simpson(f,mid,b);
    //精度够直接返回
    if(fabs(lans+rans-ans)<=eps) return ans;
    //精度不够递归调用
    return asr(f,a,mid,eps,lans) + asr(f,mid,b,eps,rans);
}
//自适应辛普森积分，区间[a, b]，精度eps，原函数f
double asr(double (*f)(double),double a,double b,double eps){
    return asr(f,a,b,eps,simpson(f,a,b));
}

//计算弓形的高度，半径为r，高度为h
inline double area(double r, double h){
    if(is0(r)||sgn(h)<=0) return 0.0;
    if(sgn(h-r-r)>=0) return PI * r * r;
    double a = r - h;
    double theta = myacos(a/r);
    return (theta - 0.5 * sin(theta+theta)) * r * r;
}

double D;//球冠底面的直径
double L;//圆柱体的长度
double R,R2;//球冠的球的半径
double T;//倾斜以后,圆柱体母线的高度,即倾斜角就是arcsin(T/L)
double H;//剩下的油,其液面高度，是从圆柱体的底端开始计算
double Theta;
double hOfQiuque,xLeft,xRight,rOfYuanzhu;
double K,B;//y=kx+B;当k是INF表示x=B

//计算x处的横截面积，这是一个弓形
double f(double x){
    double r = rOfYuanzhu;
    if(x<0){
        r = sqrt(R2-(x-xLeft)*(x-xLeft));
    }else if(x>L){
        r = sqrt(R2-(x-xRight)*(x-xRight));
    }
    double y0 = K*x + B;
    if(isEq(K,INF)) y0 = (x<=B?INF:-INF);
    return area(r,y0+r-rOfYuanzhu);
}

double proc(){
    if(isEq(T,L)) K = INF, B = H;
    else Theta = asin(T/L), K = -tan(Theta), B = H / cos(Theta);

    rOfYuanzhu = 0.5 * D;
    xLeft = sqrt(R2-rOfYuanzhu*rOfYuanzhu);
    hOfQiuque = R - xLeft;
    xRight = L + hOfQiuque - R;
    return asr(f,-hOfQiuque,0.0,EPS)
        + asr(f,0.0,L,EPS)
        + asr(f,L,L+hOfQiuque,EPS);
}

int main(){
    /*以下两行注释掉就能在牛客上AC*/
    freopen("damage.in","r",stdin);
    freopen("damage.out","w",stdout);
    //freopen("1.txt","r",stdin);
    while(5==scanf("%lf%lf%lf%lf%lf",&D,&L,&R,&T,&H)){
        R2 = R * R;
        printf("%.6f\n",proc()/1E6);
    }
    return 0;
}

Romberg方法

Romberg方法采用了更多段的分段积分。因为在CF上如下数据很难通过。
10000 10000 5000 10000 0
这一组数据本质上就是计算一个球缺的体积，因此其积分曲线的一阶导数变化比较多样，有的地方比较平缓，而有的地方比较陡峭。使用Romberg方法需要将 $e p s$ 设置的比较小，而这又会T其他点。分段积分就比较好。

#include <bits/stdc++.h>
using namespace std;

double const EPS = 1E-1;
double const PI = acos(-1.0);
double const INF = 1E100;

inline int sgn(double x){return x>EPS?1:(x<-EPS?-1:0);}
inline bool is0(double x){return 0==sgn(x);}
inline bool isEq(double x,double y){return is0(x-y);}
inline double myasin(double x){
    if(x<-1.0)x=-1.0;
    else if(x>1.0)x=1.0;
    return asin(x);
}
inline double myacos(double x){
    if(x<-1.0)x=-1.0;
    else if(x>1.0)x=1.0;
    return acos(x);
}

//Romberg积分,求f(x)在[x1,x2]上的定积分
//r[0]表示T序列，r[1]表示S序列,r[2]表示C序列,r[3]表示R序列
double Romberg(double(*f)(double), double x1, double x2,double r[][4]){
    double h = x2 - x1;
    int m = 1;
    r[0][3] = INF;

    //首先计算T0
    r[0][0] = 0.5 * h * (f(x1)+f(x2));

    int const n = 20;
    for(int i=1;i<n;++i){
        double *now = *(r + i);
        double *prev = *(r + i - 1);
        double h2 = h;
        h *= 0.5, m <<= 1;

        //计算Ti
        double x = x1 + h;
        double& sum = now[0];
        sum = 0.0;
        for(int j=1;j<m;j+=2){
            sum += f(x);
            x += h2;
        }
        sum = 0.5 * prev[0] + h * sum;

        //计算Si
        now[1] = (4.0*now[0]-prev[0]) / 3.0;
        //计算Ci
        now[2] = (16.0*now[1]-prev[1]) / 15.0;
        //计算Ri
        now[3] = (64.0*now[2]-prev[2]) / 63.0;
        //cout<<i<<" "<<r[i][0]<<" "<<r[i][3]<<endl;
        //判断
        if(isEq(prev[3],now[3])) return now[3];
    }
    return r[n-1][3];
    throw runtime_error("XXXXXX");
}//*/

//计算弓形的高度，半径为r，高度为h
inline double area(double r, double h){
    if(is0(r)||sgn(h)<=0) return 0.0;
    if(sgn(h-r-r)>=0) return PI * r * r;
    double a = r - h;
    double theta = myacos(a/r);
    return (theta - 0.5 * sin(theta+theta)) * r * r;
}

double D;//球冠底面的直径
double L;//圆柱体的长度
double R,R2;//球冠的球的半径
double T;//倾斜以后,圆柱体母线的高度,即倾斜角就是arcsin(T/L)
double H;//剩下的油,其液面高度，是从圆柱体的底端开始计算
double Theta;
double hOfQiuque,xLeft,xRight,rOfYuanzhu;
double K,B;//y=kx+B;当k是INF表示x=B
double Tmp[1004][4];

//计算x处的横截面积，这是一个弓形
double f(double x){
    double r = rOfYuanzhu;
    if(x<0){
        r = sqrt(R2-(x-xLeft)*(x-xLeft));
    }else if(x>L){
        r = sqrt(R2-(x-xRight)*(x-xRight));
    }
    double y0 = K*x + B;
    if(isEq(K,INF)) y0 = (x<=B?INF:-INF);
    return area(r,y0+r-rOfYuanzhu);
}

double proc(){
    if(isEq(T,L)) K = INF, B = H;
    else Theta = asin(T/L), K = -tan(Theta), B = H / cos(Theta);

    rOfYuanzhu = 0.5 * D;
    xLeft = sqrt(R2-rOfYuanzhu*rOfYuanzhu);
    hOfQiuque = R - xLeft;
    xRight = L + hOfQiuque - R;
    int n = 10;
    double delta = (L+hOfQiuque+hOfQiuque)/n;
    double x = -hOfQiuque,ans=0.0;
    for(int i=0;i<n;++i){
        ans += Romberg(f,x,x+delta,Tmp);
        x += delta;
    }
    return ans;
}

int main(){
 	/*以下两行注释掉就能在牛客上AC*/
    freopen("damage.in","r",stdin);
    freopen("damage.out","w",stdout);
    //freopen("1.txt","r",stdin);
    while(5==scanf("%lf%lf%lf%lf%lf",&D,&L,&R,&T,&H)){
        R2 = R * R;
        printf("%.6f\n",proc()/1E6);
    }
    return 0;
}

罗博士

关注

2
点赞
踩
1

收藏

觉得还不错? 一键收藏
0
评论
NEERC 2014 D题 Damage Assessment

NEERC 2014 D题 Damage Assessment题目描述直接计算定积分辛普森方法Romberg方法题目描述这里是NEERC2014的相关资料，包括榜、题目、标程，但是好像没有解题报告。这是CF上的重现Gym100553。当然在ICPC Live Archive上也有，这里是链接，但是好像不能提交了。所以去CF上比较好。牛客网上也有这道题，但是数据与CF上的好像不同，有的代码两边都能过，但有的代码只能过其中一边。数据应该是一样的，也许是SPJ与时限写的不同。如图，一个油罐，两端是球缺
复制链接

扫一扫