zoj 2369 Two Cylinders

最新推荐文章于 2019-01-06 16:09:43 发布

dengyun9338

最新推荐文章于 2019-01-06 16:09:43 发布

阅读量86

点赞数

原文链接：http://www.cnblogs.com/Duahanlang/p/3219159.html

版权

zoj 2369 Two Cylinders

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

题意：已知两个无限长的圆柱的半径，它们轴线相交且垂直，求相交部分的体积。

思路：以相交点为坐标原点，建立三维坐标系，列出两个圆柱的公式，把三维转化为一维积分公式

　　　　V₁ : x² + y² = r₁²　　　　==> y = sqrt ( r₁² - x² );

　　　　V₂ : x² + z² = r₂²　　　　==> z = sqrt ( r₂² - x² );

　　　　V = ∫∫∫ dv = ∫∫∫ dxdydz = ∫(∫∫ dydz)dx　　　　积分区间：( 0 <= x <= r₁ , 0 <= x <= r₂ )

　　　　所以： V = ∫₀^{min(r₁, r₂)} sqrt(( r₁²-x²) * (r₂²-x²));

计算这个积分，我用了好几种方法，但是只有变步长simpon 积分达到精度要求，我也不清楚什么原因。下面是几种求积分的算法：

测试数据：

正确答案：

5.3333
12.1603
70.9361
123.0975
444.2909
778.2605
2488.0144
1869197.4256
498415.2860
164517.4621

变步长 simpon 积分：（得到正确结果）

 1 #include <iostream>
 2 #include <stdio.h>
 3 #include <algorithm>
 4 #include <queue>
 5 #include <map>
 6 #include <set>
 7 #include <vector>
 8 #include <cstring>
 9 #include <math.h>
10 using namespace std;
11 
12 const double eps = 1e-4;
13 double r1, r2;
14 
15 double f(double x)
16 {
17     return 8. * sqrt((r1*r1 - x*x) * (r2*r2 - x*x));
18 }
19 
20 double simpson(double a, double b)
21 {
22     int n, k;
23     double h, t1, t2, s1, s2, ep, p, x;
24     n = 1, h = b - a;
25 
26     t1 = h * ( f(a) + f(b) ) / 2.;        //计算T1 = (b-a)/2*[f(a) + f(b)]
27     s1 = t1;        //T1代替S1
28     ep = eps + 1.;
29     while(ep >= eps)
30     {
31         p = 0;
32         for(k = 0; k <= n-1; k++)
33         {
34             x = a + (k + 0.5) * h;
35             p += f(x);
36         }
37         t2 = (t1 + h * p) / 2.;    //计算T2n = Tn/2 + h/2 * [(累加)f(xk + 0.5)];
38         s2 = (4. * t2 - t1) / 3.; //计算Sn = (4*Tn - Tn) / 3;
39         ep = fabs(s2 - s1);  //计算精度
40         t1 = t2, s1 = s2, n += n, h /= 2.;
41     }
42     return s2;
43 }
44 
45 int main()
46 {
47     int CASE;
48     freopen("zoj2369.txt", "r", stdin);
49     cin >> CASE;
50     while(CASE--)
51     {
52         scanf("%lf%lf", &r1, &r2);
53         if(r2 < r1)    swap(r1, r2);
54         double t = simpson(0., r1);
55         printf("%.4lf\n", t);
56 
57     }
58     return 0;
59 }
60 
61 
62 /* 结果：
63 5.3333
64 12.1603
65 70.9361
66 123.0974
67 444.2908
68 778.2604
69 2488.0144
70 1869197.4256
71 498415.2860
72 164517.4621
73 */

勒让德-高斯求积（精度不够）：

 1 #include <iostream>
 2 #include <stdio.h>
 3 #include <algorithm>
 4 #include <queue>
 5 #include <map>
 6 #include <set>
 7 #include <vector>
 8 #include <cstring>
 9 #include <math.h>
10 using namespace std;
11 
12 const double eps = 1e-4;
13 double r1, r2;
14 
15 double f(double x)
16 {
17     return 8. * sqrt((r1*r1 - x*x) * (r2*r2 - x*x));
18 }
19 
20 double lrgs(double a, double b)
21 {
22     int m, i, j;
23     double s, p, ep, aa, bb, w, x, g, h;
24     static double t[5] = {-0.9061798459, -0.5384693101, 0.0, 0.5384693101, 0.9061798459};
25     static double c[5] = {0.2369268851, 0.4786286705, 0.5688888889, 0.4786286705, 0.2369268851};
26     m = 1;
27     h = b - a; s = fabs(0.001 * h);
28     p = 1.0e+35, ep = eps + 1.0;
29     while(ep >= eps && fabs(h) > s)
30     {
31         g = 0.;
32         for(i = 1; i <= m; i++)
33         {
34             aa = a + (i-1.) * h,  bb = a + i * h;
35             w = 0.;
36             for(j = 0; j < 5; j++)
37             {
38                 x = ((bb - aa) * t[j] + (bb + aa)) / 2.;
39                 w += f(x) * c[j];
40             }
41             g += w;
42         }
43         g = g * h / 2.;
44         ep = fabs(g-p) / (1. + fabs(g));
45         p = g, m += 1, h = (b-a) / m;
46     }
47     return g;
48 }
49 
50 int main()
51 {
52     int t;
53     freopen("zoj2369.txt", "r", stdin);
54     scanf("%d", &t);
55     while(t--)
56     {
57         scanf("%lf%lf", &r1, &r2);
58         if(r1 > r2)    swap(r1, r2);
59         printf("%.4lf\n", lrgs(0, r1));
60     }
61     return 0;
62 }
63 
64 /*结果
65 5.3333
66 12.1619
67 70.9440
68 123.1138
69 444.3504
70 778.3592
71 2488.3679
72 1869425.2091
73 498482.1853
74 164540.5195
75 */

龙贝格求积法（精度不够）：

 1 #include <iostream>
 2 #include <stdio.h>
 3 #include <algorithm>
 4 #include <vector>
 5 #include <cstring>
 6 #include <math.h>
 7 using namespace std;
 8 
 9 const double eps = 1e-4;
10 double r1, r2;
11 
12 double f(double x)
13 {
14     return 8. * sqrt((r1*r1 - x*x) * (r2*r2 - x*x));
15 }
16 
17 double romb(double a, double b)
18 {
19     int n, m, i, k;
20     double y[10], h, ep, p, x, s, q;
21     h = b - a;
22     y[0] = h * ((f(a) + f(b))) / 2.;
23     m = 1, n = 1, ep = eps + 1.0;
24     while(ep >= eps )
25     {
26         p = 0.0;
27         for(i = 0; i < n; i++)
28         {
29             x = a + (i+0.5) * h;
30             p += f(x);
31         }
32         p = (y[0] + h * p) /2.;
33         s = 1.0;
34         for(k = 1; k <= m; k++)
35         {
36             s = 4. * s;
37             q = (s*p - y[k-1]) / (s - 1.);
38             y[k - 1] = p, p = q;
39         }
40         ep = fabs(q - y[m-1]);
41         m = m + 1, y[m - 1] = q, n += n, h /= 2.;
42     }
43     return q;
44 }
45 
46 int main()
47 {
48     int t;
49     freopen("zoj2369.txt", "r", stdin);
50     freopen("zzin1.txt", "w", stdout);
51     scanf("%d", &t);
52     while(t--)
53     {
54         scanf("%lf%lf", &r1, &r2);
55         if(r1 > r2)    swap(r1, r2);
56         printf("%.4lf\n", romb(0, r1));
57     }
58     return 0;
59 }
60 
61 /*结果
62 5.3333
63 12.1529
64 70.8978
65 123.0697
66 444.1897
67 778.0925
68 2487.8023
69 1869191.3654
70 498410.2627
71 164515.7323
72 */

也不知最后两个代码那里写搓了，真心无力。

下面给出标程：

 1 #include <cmath>
 2 #include <cstdio>
 3 #include <algorithm>
 4 
 5 using namespace std;
 6 
 7 const long double H = 5e-5;
 8 
 9 struct F {
10     long double rr1, rr2;
11     F(long double r1, long double r2) : rr1(r1 * r1), rr2(r2 * r2) {
12     }
13     long double operator()(long double x) const {
14         x = x * x;
15         return sqrt((rr2 - x) * (rr1 - x));
16     }
17 };
18 
19 int main() {
20     int re;
21 
22     scanf("%d", &re);
23     for (int ri = 1; ri <= re; ++ri) {
24         long double r1, r2;
25         scanf("%Lf%Lf", &r1, &r2);
26         if (r1 > r2) {
27             swap(r1, r2);
28         }
29         int n = int(r1 / H);
30         // int n = 1000000;
31         long double h = r1 / n / 2;
32         F f(r1, r2);
33         long double s = 0;
34         s += f(0) + f(r1);
35         for (int i = 1; i < 2 * n; ++i) {
36             if (i & 1) {
37                 s += 4 * f(i * h);
38             } else {
39                 s += 2 * f(i * h);
40             }
41         }
42         s /= 6 * n;
43         s *= 8 * r1;
44         printf("%.6Lf\n", s);
45     }
46 
47     return 0;
48 }

网上看的大牛的代码，0ms秒过，我的用了290ms，标程用了310ms

 1 #include <cstdio>
 2 #include <cmath>
 3 #include <algorithm>
 4 using namespace std;
 5 
 6 #define db double
 7 #define rt return
 8 #define cs const
 9 
10 cs db eps = 1e-9;
11 inline int sig(db x){rt (x>eps)-(x<-eps);}
12 
13 db r1, r2;
14 
15 inline db f(db x) {
16     rt 4. * sqrt(r1*r1-x*x) * sqrt(r2*r2-x*x);
17 }
18 
19 inline db simpson(db a, db b) {
20     rt (b-a)*(f(a) + 4.*f((a+b)/2.) + f(b)) / 6.;
21 }
22 
23 inline db calc(db a, db b, int depth) {
24     db total = simpson(a, b), mid = (a+b) / 2.;
25     db tmp = simpson(a, mid) + simpson(mid, b);
26     if(sig(total-tmp) == 0 && depth > 3) rt total;
27     rt calc(a, mid, depth + 1) + calc(mid, b, depth + 1);
28 }
29 
30 int main() {
31 //    freopen("std.in", "r", stdin);
32     int T;
33     freopen("zoj2369.txt", "r", stdin);
34     scanf("%d", &T);
35     while(T--) {
36         scanf("%lf%lf", &r1, &r2);
37         if(sig(r1-r2) > 0) swap(r1, r2);
38         printf("%.4lf\n", 2. * calc(0., r1, 0));
39     }
40     rt 0;
41 
42 }

转载于:https://www.cnblogs.com/Duahanlang/p/3219159.html

dengyun9338

关注

0
点赞
踩
0

收藏

觉得还不错? 一键收藏
0
评论
zoj 2369 Two Cylinders

zoj 2369 Two Cylinders链接：http://acm.zju.edu.cn/onlinejudge/showProblem.do?problemCode=2369题意：已知两个无限长的圆柱的半径，它们轴线相交且垂直，求相交部分的体积。思路：以相交点为坐标原点，建立三维坐标系，列出两个圆柱的公式，把三维转化为一维积分公式　　　　V1 : x2 + y2 ...
复制链接

扫一扫