随笔-21  评论-10  文章-21  trackbacks-0
   这道题显然考察积分,但有些积分是困难的(所以要平时多做题,能判断哪些是自己能够手算积出来的),
能积出来当然尽量算出来,但对于那些困难的,就可以模拟积分的过程, 这里微元可以选择横截面,对于每一个横截面的面积 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 阅读(343) 评论(1)  编辑 收藏 引用

评论:
# re: pku 3289 Moonshine 2013-06-02 00:38 | ys
这道题有积分公式,具体见poj 3929  回复  更多评论
  

只有注册用户登录后才能发表评论。
【推荐】超50万行VC++源码: 大型组态工控、电力仿真CAD与GIS源码库
网站导航: 博客园   IT新闻   BlogJava   知识库   博问   管理