这道题显然考察积分，但有些积分是困难的（
所以要平时多做题，能判断哪些是自己能够手算积出来的），
能积出来当然尽量算出来，但对于那些困难的，就可以模拟积分的过程， 这里微元可以选择横截面，对于每一个横截面的面积 S 是 一个弓形，这是好求的，再沿轴线方向积分却困难了，边看代码边解释：
waterloo的标程
1 #include <stdio.h>
2 #include <math.h>
3 #include <assert.h>
4
5 #define N 1000
6
7 double s,delta,k,hb,db,hn,dn,h,V,hs,v;
8
9 //求每个横截面微元弓形的面积
10 double area(double r, double s) {
11 double t = db/2-s;
12 double theta = t<r?acos(t/r):0;
13 double wedge = theta * r * r;
14 double triangles = r * cos(theta) * r * sin(theta);
15 return wedge - triangles;
16 }
17 //求横放 液面高度为 s时 的体积
18 double volume(double s) {
19 if (s*2 > db) return V - volume(db-s);
20 double vb = area(db/2, s) * hb;
21 double vn = area(dn/2, s) * hn;
22 double vc = area(db/2, s) + area(dn/2, s);
23 int i;
24 /*================*/
25 /*精彩的地方
26
27 因为求不出那个定积分，所以只能分成 N 段逼近，当然 N 越大越准确，越大时间越多
28 如果对被积函数的特性，如单调性， 斜率的变化情况快慢，大致曲线的形状，越清楚，
29 就越准确，像在这里，微元 area 以 r 作为变量，却不是正比于 r^2，所以直接用
30 for (i=1;i<N; i++) {
31 vc += area((db + (dn-db)*i/N)/2, s);
32 }
33 vc *= (h-hn-hb)/N;
34 去模拟那个积分过程效果就要差的多，N==5000时才能得到正确答案，而下面得方法 N==20
35 就能得到正确答案，差距不是一点点！！
36 */
37 for (i=1;i<N;i+=2) {
38 vc += 4 * area((db + (dn-db)*i/N)/2, s);
39 }
40 for (i=2;i<N;i+=2) {
41 vc += 2 * area((db + (dn-db)*i/N)/2, s);
42 }
43 vc *= (h-hn-hb)/3/N;
44 /*====================*/
45 return vb + vn + vc;
46 }
47
48 doit(){
49 int hnx = hn, hbx = hb, hx = h, dnx = dn;
50 V = 2 * volume(db/2);//总的容积
51 /*=======================*/
52 //这段处理时修改瓶子，使得“瓶子”总是装满液体的，这样就可以用同一个公式求出液体的体积
53 if (k <= hb) {
54 hn = 0;
55 hb = k;
56 h = k;
57 } else if (k < h-hn) {
58 dn = db + (dn-db) * (k-hb) / (h-hn-hb);
59 hn = 0;
60 h = k;
61 } else {
62 hn = hn - h + k;
63 h = k;
64 }
65 v = 2 * volume(db/2);//液体的体积
66 /*========================*/
67 hn = hnx; hb = hbx; h = hx; dn = dnx;
68 //将瓶子还原回来
69
70 /*==========================*/
71 //这个二分看得不怎么习惯
72 s = db/2;
73 for (delta=s/2;delta > .00001; delta /=2) {
74 if (volume(s) > v) s -= delta;
75 else s += delta;
76 }
77 /*======================*/
78 printf("%0.2lf\n",s);
79 }
80
81 main(){
82 while (6 == scanf("%lf%lf%lf%lf%lf%lf",&k,&hb,&db,&hn,&dn,&h)) {
83 if (!(k||hb||db||hn||dn||h)) return;
84 doit();
85 }
86 }
87
my code：
1 #include<iostream>
2 #include<cmath>
3 using namespace std;
4 #define pi 3.1415926535897932384626433832795
5 #define integral_Num 20//数据特殊才取20，照理还是取大点好
6 #define eps 1e-8
7 int sig(double x){return x < -eps ? -1 : x > eps ; }
8 double k ,hb ,db ,hn ,dn ,h, V;
9 double area(double height, double r)
10 {
11 if( sig(height - db/2 + r) < 0 ) return 0;
12 if( sig(height - db/2 - r) >= 0 ) return pi*r*r;
13 /*
14 这里要确保acos的值为-1.0 ~1.0之间， 甚至为1.0+1e-10都不行,
15 此外 r 不能为 0，r == 0 的情况必须在上面就处理掉
16 */
17 double theta = acos( ( db/2 - height) / r );
18 return r * r * theta - r * r * sin(theta) * cos(theta) ;
19 }
20
21 double volumn(double height)
22 {
23 int i;
24 double
25 vb = area(height, db/2) * hb,
26 vn = area(height, dn/2) * hn,
27 vc = area(height, db/2) + area(height, dn/2);
28 for(i = 1; i < integral_Num; i+=2)
29 vc+=4 * area(height , 0.5 * db - 0.5 * (db - dn) * i / integral_Num );
30 for(i = 2; i < integral_Num; i+=2)
31 vc+=2 * area(height , 0.5 * db - 0.5 * (db - dn) * i / integral_Num );
32 vc*=(h - hb - hn) / 3 / integral_Num;
33 return vc+vn+vb;
34 }
35 void doit()
36 {
37 double tk = k, thb = hb, tdb = db, thn = hn, tdn = dn, th = h;
38 if(k <= hb){ hn=0; hb=k; h=hb; }
39 else if(k > h - hn) { hn=k+hn-h; h=k; }
40 //要防止 h - hn - hb == 0 的情况运行到这一步
41 else {dn = db + (hb - k)*(db - dn)/(h - hn - hb); h=k; hn=0; }
42
43 double v, l = 0 , r=db, mid;
44 v=2*volumn(db/2);
45 k = tk, hb = thb, db = tdb, hn = thn, dn = tdn, h = th;
46 while(r - l >= 1e-3)
47 {
48 mid=(r+l)/2;
49 if(volumn(mid) > v) r = mid;
50 else l = mid;
51 }
52 printf("%.2lf\n",mid);
53 }
54 int main()
55 {
56 /*freopen("C.1.dat","r",stdin);
57 freopen("ans.txt","w",stdout);*/
58 while(scanf("%lf%lf%lf%lf%lf%lf",&k, &hb, &db, &hn, &dn, &h) !=EOF)
59 {
60 if(!(k || hb || db || hn || dn || h) )break;
61 doit();
62 }
63 }
64
posted on 2009-02-27 22:45
wangzhihao 阅读(390)
评论(1) 编辑 收藏 引用