hdu5060 (数学题，辛普森公式)

War

Time Limit: 8000/5000 MS (Java/Others)    Memory Limit: 32768/32768 K (Java/Others)
Total Submission(s): 151    Accepted Submission(s): 45
Special Judge

Problem Description

Long long ago there are two countrys in the universe. Each country haves its own manor in 3-dimension space. Country A's manor occupys x^2+y^2+z^2<=R^2. Country B's manor occupys x^2+y^2<=HR^2 && |z|<=HZ. There may be a war between them. The occurrence of a war have a certain probability.
We calculate the probability as follow steps.
1. VC=volume of insection manor of A and B.
2. VU=volume of union manor of A and B.
3. probability=VC/VU

 

Input

Multi test cases(about 1000000). Each case contain one line. The first line contains three integers R,HR,HZ. Process to end of file.

[Technical Specification]
0< R,HR,HZ<=100

 

Output

For each case，output the probability of the war which happens between A and B. The answer should accurate to six decimal places.

 

Sample Input

  
  

   
   1 1 1
2 1 1
  
  

 

Sample Output

  
  

   
   0.666667
0.187500
  
  

题意：

就是计算圆柱体和球的相交体积。数学分类讨论，计算。

这里有官方的解题思路，画的图灰常的好，我就借过来用一下辣，哈哈哈哈。。

1003 War

首先要求出球体和圆柱体之间的相交的体积。
然后总体积减去相交的体积就是并的体积。
求交的体积可以如上图，想像成一个矩阵和一个圆相交部分绕X轴旋转而成的一个东西，就是上图阴影部分绕X轴出来的一个东西。
然后体积公式为 
  
  
   
   ∫HR0 2πy∗2R2−y2−−−−−−√dy=−43π(R2−y2)3/2|HR0
  
  
总共有五种情况，还有四种情况是

这个直接计算中间的球体即可。

直接计算圆柱体。

这个可以直接套用公式，把HR改成R。

分成黄公和灰色两部分计算，灰色部分是圆柱体好算。	
黄色部分可以套用公式：

  
  
   
   ∫HRr0 2πy∗2R2−y2−−−−−−√dy=−43π(R2−y2)3/2|HRr0
  
  
其中
  
  
   
   r0=R2−HZ2−−−−−−−−√
  
  

用到了积分，好久没看过高数了。。。乍一看，还真是不懂，后来还是回想起来辣，这次很高兴滴认识了辛普森~~~

这里有介绍 http://baike.baidu.com/view/2710883.htm?fr=aladdin 其实就是一个公式，直接套用就可可以了。

AC代码：

//设拟柱体的高（两底面α，β间的距离）为H，如果用平行于底面的平面γ去截该图形，
//所得到的截面面积是平面γ与平面α之间距离h的不超过3次的函数，
//那么该拟柱体的体积V为V = H (S_1 + 4S_0 + S_2) /6.
//式中，S_1和S_2是两底面的面积，S_0是中截面的面积（即平面γ与平面α之间距离h=H/2时得到的截面的面积）。 

//开森地认识了辛普森 公式  
#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
using namespace std;
const double pi=3.1415926,eps=1e-7;
double r,hr,hz;

double f(double n)
{
	return pi*(r*r-n*n); //求微元x的面积  
} 

double simpson(double a,double b)// a b 是y的取值范围 
{
	return (b-a)/6.0*(f(a)+4*f((a+b)/2)+f(b));
}

double cal(double a,double b)
{
	double sum=simpson(a,b);
	double t=simpson(a,(a+b)/2)+simpson((a+b)/2,b);
	if(abs(sum-t)<eps)return sum;
	return simpson(a,(a+b)/2)+simpson((a+b)/2,b);
} 

int main()
{
	while(~scanf("%lf%lf%lf",&r,&hr,&hz))
	{
		double v=4.0*pi*r*r*r/3.0;
		double hv=pi*hr*hr*hz*2.0;
		double sum=v+hv;
		double ans=0.0;
		double v0=0.0,vm=0.0;
		
		if(hr*hr+hz*hz<=r*r)
		{
			ans=hv/v;
		}
		else if(hr>=r && hz>=r)
		{
			ans=v/hv;
		}
		else if(hz>=r && hr<=r)
		{
			double d=sqrt(r*r-hr*hr);
			v0=2.0*cal(d,r)+pi*hr*hr*d*2.0;
			ans=v0/(sum-v0);
		}
		else if(hr>=r && hz<r)
		{
		//	double d=sqrt(r*r-hr*hr);
			v0=2.0*cal(hz,r);
			ans=(v-v0)/(hv+v0);
		}
		else 
		{
			double d=sqrt(r*r-hr*hr);
			v0=pi*hr*hr*2.0*d+2.0*cal(d,hz);
			
			ans=v0/(sum-v0);
		}
		
		printf("%.6f\n",ans);
	}
}