這道題顯然考察積分，但有些積分是困難的（
所以要平時多做題，能判斷哪些是自己能夠手算積出來的），
能積出來當然盡量算出來，但對于那些困難的，就可以模擬積分的過程， 這里微元可以選擇橫截面，對于每一個橫截面的面積 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 閱讀(399)
評論(1) 編輯 收藏 引用