Simpson公式的应用（HDU 1724/ HDU 1071）

辛普森积分法 - 维基百科，自由的百科全书

Simpson's rule - Wikipedia, the free encyclopedia

　　利用这个公式，用二分的方法来计算积分。

1071 ( The area )

 1 #include <iostream>
 2 #include <algorithm>
 3 #include <cstdio>
 4 #include <cstring>
 5 #include <cmath>
 6 
 7 using namespace std;
 8 
 9 const double EPS = 1e-8;
10 double A, B, C, P, Q;
11 
12 template<class T> T sqr(T x) { return x * x;}
13 inline double cal(double x) { return A * sqr(x) + (B - P) * x + C - Q;}
14 inline double sps(double l, double r) { return (cal(l) + cal(r) + 4 * cal((l + r) / 2)) / 6 * (r - l);}
15 
16 double work(double l, double r) {
17     //cout << l << ' ' << r << endl;
18     double ans = sps(l, r), m = (l + r) / 2;
19     if (fabs(ans - sps(l, m) - sps(m, r)) < EPS) return ans;
20     else return work(l, m) + work(m, r);
21 }
22 
23 
24 int main() {
25     int T;
26     double l, r;
27     double x[3], y[3];
28     cin >> T;
29     while (T--) {
30         for (int i = 0; i < 3; i++) cin >> x[i] >> y[i];
31         double p[2], q[2], d[2];
32         for (int i = 0; i < 2; i++) p[i] = sqr(x[i]) - sqr(x[i + 1]), q[i] = x[i] - x[i + 1], d[i] = y[i] - y[i + 1];
33         A = (q[1] * d[0] - q[0] * d[1]) / (p[0] * q[1] - p[1] * q[0]);
34         B = (p[1] * d[0] - p[0] * d[1]) / (p[1] * q[0] - p[0] * q[1]);
35         C = y[0] - B * x[0] - A * sqr(x[0]);
36         //cout << A << ' ' << B << ' ' << C << endl;
37         P = (y[1] - y[2]) / (x[1] - x[2]);
38         Q = y[1] - P * x[1];
39         //cout << P << ' ' << Q << endl;
40         printf("%.2f\n", work(x[1], x[2]));
41     }
42     return 0;
43 }

View Code

 

1724 ( Ellipse )

 1 #include <iostream>
 2 #include <algorithm>
 3 #include <cstdio>
 4 #include <cstring>
 5 #include <cmath>
 6 
 7 using namespace std;
 8 
 9 const double EPS = 1e-8;
10 double A, B;
11 
12 template<class T> T sqr(T x) { return x * x;}
13 inline double cal(double x) { return 2 * B * sqrt(1 - sqr(x) / sqr(A));}
14 inline double sps(double l, double r) { return (cal(l) + cal(r) + 4 * cal((l + r) / 2)) / 6 * (r - l);}
15 
16 double work(double l, double r) {
17     //cout << l << ' ' << r << endl;
18     double ans = sps(l, r), m = (l + r) / 2;
19     if (fabs(ans - sps(l, m) - sps(m, r)) < EPS) return ans;
20     else return work(l, m) + work(m, r);
21 }
22 
23 int main() {
24     int T;
25     double l, r;
26     cin >> T;
27     while (T-- && cin >> A >> B >> l >> r) printf("%.3f\n", work(l, r));
28     return 0;
29 }

View Code

 

　　之后还有题会继续更新。

UPD：

　　就是因为见过这题，所以才学这个公式的。1y~

 ACM-ICPC Live Archive

 1 #include <cstdio>
 2 #include <cstring>
 3 #include <iostream>
 4 #include <algorithm>
 5 #include <cmath>
 6 
 7 using namespace std;
 8 
 9 double coe[4][11];
10 const double EPS = 1e-8;
11 
12 int k;
13 double cal(double x, double *c) {
14     double ret = c[0];
15     for (int i = 1; i <= k; i++) ret *= x, ret += c[i];
16     return ret;
17 }
18 
19 inline double cal(double x, double *p, double *q) { return cal(x, p) / cal(x, q);}
20 inline double cal(double x, double y, double *p, double *q) { return max(cal(x, p, q) - y, 0.0);}
21 inline double simpson(double y, double l, double r, double *p, double *q) { return (cal(l, y, p, q) + cal(r, y, p, q) + 4 * cal((l + r) / 2, y, p, q)) * (r - l) / 6;}
22 
23 inline double getpart(double y, double l, double r, double *p, double *q) {
24     double sum = simpson(y, l, r, p, q);
25     //cout << l << ' ' << r << ' ' << sum << endl;
26     if (fabs(sum - simpson(y, l, (l + r) / 2, p, q) - simpson(y, (l + r) / 2, r, p, q)) < EPS) return sum;
27     return getpart(y, l, (l + r) / 2, p, q) + getpart(y, (l + r) / 2, r, p, q);
28 }
29 
30 inline double getarea(double y, double l, double r, double *p, double *q) {
31     double ret = 0, d = (r - l) / 100;
32     for (int i = 0; i < 100; i++) {
33         ret += getpart(y, l + d * i, l + d * (i + 1), p, q);
34     }
35     return ret;
36 }
37 
38 double dc2(double l, double r, double a, double w) {
39     double m;
40     while (r - l > EPS) {
41         m = (l + r) / 2.0;
42         //cout << m << ' ' << getarea(m, 0, w, coe[0], coe[1]) - getarea(m, 0, w, coe[2], coe[3]) << endl;
43         if (getarea(m, 0, w, coe[0], coe[1]) - getarea(m, 0, w, coe[2], coe[3]) > a) l = m;
44         else r = m;
45     }
46     return l;
47 }
48 
49 int main() {
50     //freopen("in", "r", stdin);
51     //freopen("out", "w", stdout);
52     double w, d, a;
53     while (cin >> w >> d >> a >> k) {
54         for (int i = 0; i < 4; i++) for (int j = 0; j <= k; j++) cin >> coe[i][j];
55         for (int i = 0; i < 4; i++) reverse(coe[i], coe[i] + k + 1);
56         //cout << getarea(-5.51389, 0, w, coe[0], coe[1]) - getarea(-5.51389, 0, w, coe[2], coe[3]) << endl;
57         //cout << cal(3, coe[0], coe[1]) << endl;
58         printf("%.5f\n", -dc2(-d, 0, a, w));
59     }
60     return 0;
61 }

View Code

 

 

——written by Lyon

转载于:https://www.cnblogs.com/LyonLys/p/simpson_Lyon.html