zoj 2369 Two Cylinders 辛普森积分 几何

链接：http://acm.zju.edu.cn/onlinejudge/showProblem.do?problemId=1345

题意：给两个无限高的圆柱的半径，它们垂直放置，坐标轴相交，求相交部分的体积。

思路：数值积分，用辛普森公式来算。把三重积分转化为一重的。

假设半径为r1的圆柱沿y轴，r2的圆柱沿z轴，那么从x轴看过去，相交部分在yoz平面的切面是个矩形。

矩形在y轴方向的长度为2*y=sqrt(r1² -x²)，在z轴方向上的长度为2*sqrt(r2² - x²)，这样积分算出来的体积是yoz平面上半部分的体积，还有下半部分的，总的体积要乘以2.

如果是从z方向看过去，则平行于xOy平面所做的切面都是一个矩形，x方向长2*sqrt(r1*r1-z*z)，y方向长2*sqrt(r2*r2-z*z)，同样总体积需要乘以2.

由高数知识可知： 有下面的三重积分成立

∫∫∫dxdydz =|Ω|  (Ω为积分区域，|Ω|为区域Ω的体积) 。所以关键在于把Ω表示出来。参考下面的图形：

 

  列出两个圆柱的方程：

  x² + y² = r1²          y=sqrt(r1² -x²)

  x² + z² = r2²          z=sqrt(r2² - x²)

  则 ∫∫∫dxdydz = ∫dx ∫∫dydz = ∫dx ∫ 2*sqrt(r1² -x²) * 2*sqrt(r2² - x²) dydz = 2* ∫^{min( r1 , r2 )}  2*2* sqrt((r1² -x²)*(r2² - x²)) dx ,  

  其中 x ∈ [ 0 , min(r1, r2) ]       是积分yz还是xy还是xz都可以。

下面是zsl学长的代码，使用的是变步长的simpson积分

#include <cstdio>
#include <cmath>
#include <algorithm>
using namespace std;

#define db double
#define rt return
#define cs const

cs db eps = 1e-9;
inline int sig(db x){rt (x>eps)-(x<-eps);}

db r1, r2;

inline db f(db x) {
    rt 4. * sqrt(r1*r1-x*x) * sqrt(r2*r2-x*x);
}

inline db simpson(db a, db b) {
    rt (b-a)*(f(a) + 4.*f((a+b)/2.) + f(b)) / 6.;
}

inline db calc(db a, db b, int depth) {
    db total = simpson(a, b), mid = (a+b) / 2.;
    db tmp = simpson(a, mid) + simpson(mid, b);
    if(sig(total-tmp) == 0 && depth > 3) rt total;
    rt calc(a, mid, depth + 1) + calc(mid, b, depth + 1);
}

int main() {
//    freopen("std.in", "r", stdin);
    int T;
    scanf("%d", &T);
    while(T--) {
        scanf("%lf%lf", &r1, &r2);
        if(sig(r1-r2) > 0) swap(r1, r2);
        printf("%.4lf
", 2. * calc(0., r1, 0));
    }
    rt 0;
}

这是我的代码，自适应simpson积分

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cmath>
#include<cstring>
#include<algorithm>
#include<set>
#include<map>
#include<vector>
#define lson l,m,rt<<1
#define rson m+1,r,rt<<1|1
typedef long long LL;
using namespace std;
const double PI=acos(-1);
const double eps=1e-4;
const int maxn=100;
double r1,r2;

double F(double x)
{
    return 8.*sqrt((r1*r1-x*x)*(r2*r2-x*x));
}
double simpson(double a,double b)
{
    double m=a+(b-a)/2;
    return (b-a)*(F(a)+4*F(m)+F(b))/6;
}
double asr(double a,double b,double e,double A)
{
    double m=a+(b-a)/2;
    double L=simpson(a,m),R=simpson(m,b);
    if(fabs(L+R-A)<=15*e) return L+R+(L+R-A)/15.;
    return asr(a,m,e/2,L)+asr(m,b,e/2,R);
}
double asr(double a,double b,double e)
{
    return asr(a,b,e,simpson(a,b));
}
int main()
{
//    freopen("in.cpp","r",stdin);
    int n;
    scanf("%d",&n);
    while(n--)
    {
        scanf("%lf%lf",&r1,&r2);
        if(r1>r2) swap(r1,r2);
        printf("%.4lf

",asr(0,r1,eps));
    }

    return 0;
}

注意：simpson积分不太稳定，所以精度eps最好比要求的精度高几级，但是要注意不要超时。

自适应simpson积分模板

double F(double x)
{
    。。。
}
double simpson(double a,double b)
{
    double m=a+(b-a)/2;
    return (b-a)*(F(a)+4*F(m)+F(b))/6;
}
double asr(double a,double b,double e,double A)
{
    double m=a+(b-a)/2;
    double L=simpson(a,m),R=simpson(m,b);
    if(fabs(L+R-A)<=15*e) return L+R+(L+R-A)/15.;
    return asr(a,m,e/2,L)+asr(m,b,e/2,R);
}
double asr(double a,double b,double e)
{
    return asr(a,b,e,simpson(a,b));
}

变步长的simpson积分模板

double f(double x)
{
    ...
}

double simpson(double a, doubleb) {
    return (b-a)*(f(a) + 4.*f((a+b)/2.) + f(b)) / 6.;
}

double calc(double a, double b, int depth) {
    double total = simpson(a, b), mid = (a+b) / 2.;
    double tmp = simpson(a, mid) + simpson(mid, b);
    if(fabs(total-tmp)<eps && depth > 3) return total;
    return calc(a, mid, depth + 1) + calc(mid, b, depth + 1);
}

同样可以就几分的还有龙贝格求积分，勒让德-高斯求积分法，以后再去看看这两个。