POJ 3289 重积分|二分_poj3289-优快云博客

本文链接：https://blog.youkuaiyun.com/CrySleeper/article/details/123012367

博客作者分享了一道涉及复杂体积积分的数学问题，该问题来源于POJ.org的一个竞赛题目。作者详细解释了如何利用积分计算解决圆锥体被切割后的体积，并针对不同情况讨论了实现算法，包括分步积分和查积分表的方法。虽然实现过程繁琐，但作者强调思路并不复杂，只是需要耐心。文章以C++代码展示了具体实现，适合对数学和编程感兴趣的读者。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

poj.org/problem?id=3289

老年人用手积分完全积不出来，只能翻翻高数课本查积分表了.

主要需要解决的体积积分如下：一个圆锥体被一个水平平面和一个垂直平面切割后剩下的部分

$\int (arccos(d/(k*h))*k^2*h^2-d*\sqrt{k^2*h^2-d^2})) d_h$

其中k是常数，为圆锥体顶角的一半的tan值， d为中心轴与垂直平面的距离，作为二分变量也是常数

上述积分的第一部分先使用一次分步积分，然后查积分表，第二部分则直接使用积分表，最终的积分是

$\int (k^2*h^3*arccos(d/(k*h))/3-d*k/6*((h*\sqrt{h^2-(d/k)^2}+(d/k)^2*ln(\left | h+\sqrt{h^2-(d/k)^2)} \right |) - d*k/2*(h*\sqrt{h^2-(d/k)^2}-(d/k)^2*ln(\left | h+\sqrt{h^2-(d/k)^2)} \right |))$

。。。。

本来积分写起来就挺乱的，题目还要分多种情况讨论

放倒瓶子之后

1. 水位线没有达到瓶口

被一个平面切割的圆柱体+被两个垂直平面切割的圆锥体

2. 水位线完全掩盖瓶口

被一个平面切割的圆柱体+被两个垂直平面切割的圆锥体+一个圆柱体

3. 水位线在瓶口中间

被一个平面切割的圆柱体+被两个垂直平面切割的圆锥体+被一个平面切割的圆柱体

实现：

积分

double iv2(double d, double h1, double h2) {
    if (h2<EPS) return 0;
    double b=d/gk;
    double a, aa, v1, dd;
    if (h1<EPS) v1=0; else {
        a=sqrt(h1*h1-b*b);
        aa=acos(b/h1);
        v1=gk*gk*h1*h1*h1*aa/3-d*gk/3*h1/2*a-d*gk*h1/2*a;
        dd=(h1+a); if (dd<0) dd=-dd;
        v1+=(d*gk-d*gk/3)/2*b*b*log(dd);
    }
    a=sqrt(h2*h2-b*b);
    aa=acos(b/h2);
    double v2=gk*gk*h2*h2*h2*aa/3-d*gk/3*h2/2*a-d*gk*h2/2*a;
    dd=(h2+a); if (dd<0) dd=-dd;
    v2+=(d*gk-d*gk/3)/2*b*b*log(dd);
    v1-=v2; if (v1<0) v1=-v1;
    return v1;
}

主函数分了好多种情况，特殊处理了圆锥体不出现的两种情形，然后每种情况都要二分定位水位线的位置，写得有点乱

        scanf("%lf %lf %lf %lf %lf %lf", &k, &hb, &db, &hn, &dn, &h);
        if (k+hb+db+hn+dn+h<EPS) break;
        rb=db/2; rn=dn/2;
        dh = h-hb-hn;
        dr = rb-rn;
        if (dr<EPS) {
            v=pi*rb*rb*k;
            d=rb; dd=rb; md=-rb+EPS;
            while(dd>EPS) {
                while(1) {
                    nd=d-dd; if (nd<md) break;
                    s=acos(nd/rb)*rb*rb-nd*sqrt(rb*rb-nd*nd);
                    nv=s*h;
                    if (nv>v) break;
                    d=nd;
                }
                dd/=2;
            } rh=rb-d;
        } else if (dh<EPS) {
            if (k<=hb) v=pi*rb*rb*k;
            else v=pi*rb*rb*hb+pi*rn*rn*(k-hb);
            nd=rn; nv=(acos(nd/rb)*rb*rb-nd*sqrt(rb*rb-nd*nd))*hb;
            if (v<=nv) {
                md=rn; d=rb; dd=rb-rn;
                while(dd>EPS) {
                    while(1) {
                        nd=d-dd; if (nd<md) break;
                        s=acos(nd/rb)*rb*rb-nd*sqrt(rb*rb-nd*nd);
                        nv=s*hb; if (nv>v) break;
                        d=nd;
                    }
                    dd/=2;
                } rh=rb-d;
            } else {
                nd=-rn; nv=(acos(nd/rb)*rb*rb-nd*sqrt(rb*rb-nd*nd))*hb;
                sb=pi*rn*rn*hn; nv+=sb;
                if (v>=nv) {
                    md=-rb+EPS; d=-rn; dd=rb-rn;
                    while(dd>EPS) {
                        while(1) {
                            nd=d-dd; if (nd<md) break;
                            s=acos(nd/rb)*rb*rb-nd*sqrt(rb*rb-nd*nd);
                            nv=s*hb+sb; if (nv>v) break;
                            d=nd;
                        }
                        dd/=2;
                    } rh=rb-d;
                } else {
                    md=-rn+EPS; d=rn; dd=rn+rn;
                    while(dd>EPS) {
                        while(1) {
                            nd=d-dd; if (nd<md) break;
                            s=acos(nd/rb)*rb*rb-nd*sqrt(rb*rb-nd*nd);
                            nv=s*hb;
                            s=acos(nd/rn)*rn*rn-nd*sqrt(rn*rn-nd*nd);
                            nv+=s*hn;
                            if (nv>v) break;
                            d=nd;
                        }
                        dd/=2;
                    } rh=rb-d;
                }
            }
        } else {
            gk=dr/dh;
            if (k<=hb) v=pi*rb*rb*k;
            else if (k<=h-hn) {
                v=pi*rb*rb*hb+iv1(rb/gk-(k-hb), rb/gk);
            } else if (k<=h) {
                v=pi*rb*rb*hb+iv1(rn/gk, rb/gk)+pi*rn*rn*(hn-(h-k));
            } else assert(0);
            nd=rn; nv=(acos(nd/rb)*rb*rb-nd*sqrt(rb*rb-nd*nd))*hb;
            nv+=iv2(nd, rn/gk, rb/gk);
            if (v<=nv) {
                md=rn; d=rb; dd=(rb-rn)/2;
                h2=rb/gk;
                while(dd>EPS) {
                    while(1) {
                        nd=d-dd; if (nd<md) break;
                        s=acos(nd/rb)*rb*rb-nd*sqrt(rb*rb-nd*nd);
                        h1=nd/gk;
                        nv=s*hb+iv2(nd, h1, h2);
                        if (nv>v) break;
                        d=nd;
                    }
                    dd/=2;
                } rh=rb-d;
            } else {
                h2=rb/gk; h1=rn/gk;
                sbb = iv1(h1, h2);
                nd=-rn; nv=(acos(nd/rb)*rb*rb-nd*sqrt(rb*rb-nd*nd))*hb+sbb-iv2(-nd, rn/gk, rb/gk);
                sb=pi*rn*rn*hn; nv+=sb;
                if (v>=nv) {
                    md=-rb+EPS; d=-rn; dd=rb-rn;
                    while(dd>EPS) {
                        while(1) {
                            nd=d-dd; if (nd<md) break;
                            h0=(nd+rb)/gk;
                            s=(pi-acos(-nd/rb))*rb*rb-nd*sqrt(rb*rb-nd*nd);
                            nv=s*hb+sb-iv2(-nd, -nd/gk, h2)+iv1(h1, h2); if (nv>v) break;
                            d=nd;
                        }
                        dd/=2;
                    } rh=rb-d;
                } else {
                    md=-rn+EPS; d=rn; dd=rn+rn;
                    h2=rb/gk; h1=rn/gk;
                    while(dd>EPS) {
                        while(1) {
                            nd=d-dd; if (nd<md) break;
                            s=acos(nd/rb)*rb*rb-nd*sqrt(rb*rb-nd*nd);
                            nv=s*hb;
                            s=acos(nd/rn)*rn*rn-nd*sqrt(rn*rn-nd*nd);
                            nv+=s*hn;
                            nv+=iv2(nd, h1, h2);
                            if (nv>v) break;
                            d=nd;
                        }
                        dd/=2;
                    } rh=rb-d;
                }
            }
        }
        printf("%.2f\n", rh);

这道题思路不难，但是实现起来很麻烦，我写了4个多小时吧，差点就放弃了。。。。。

POJ AC的2729道题，继续寻找水题ing~