自适应Simpson积分学习笔记
发表更新ICPC/CCPC / Notes15 分钟读完 (大约2196个字)

自适应Simpson积分学习笔记

一种精度较高的积分数值计算方法。

原理简述

（随便捞一本数分的书应该都会讲）Simpson 积分是积分数值计算的一种方法，用抛物线而非传统的矩形来近似图形，精度较高，公式为：

然而实际解题的过程中不好控制划分的数量 $n$，所以提出了自适应 Simpson 积分。所谓自适应，即积分数值在可接受范围内时直接返回积分，否则对半分开递归求解，这里每层在积分时只划分成 $1$ 条即可，也就是说：

即可。

模板

1
2
3
4
5
6
7
8
9
10
double simpson(double l, double r){
    double mid = (l + r) / 2;
    return (f(l) + 4 * f(mid) + f(r)) * (r - l) / 6;
}
double solve(double l, double r, double _eps, double I){
    double mid = (l + r) / 2;
    double Il = simpson(l, mid), Ir = simpson(mid, r);
    if(fabs(Il + Ir - I) <= 15 * _eps)   return Il + Ir + (Il + Ir - I) / 15;
    return solve(l, mid, _eps / 2, Il) + solve(mid, r, _eps / 2, Ir);
}

练习

luoguP4525 【模板】自适应辛普森法1

题目链接

>folded
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
#include<algorithm>
#include<cstdio>
#include<cmath>

using namespace std;

double a, b, c, d;
double f(double x){
    return (c * x + d) / (a * x + b);
}
double simpson(double l, double r){
    double mid = (l + r) / 2;
    return (f(l) + 4 * f(mid) + f(r)) * (r - l) / 6;
}
double solve(double l, double r, double eps){
    double mid = (l + r) / 2;
    double Il = simpson(l, mid), Ir = simpson(mid, r), I = simpson(l, r);
    if(fabs(Il + Ir - I) <= 15 * eps)   return Il + Ir + (Il + Ir - I) / 15;
    return solve(l, mid, eps / 2) + solve(mid, r, eps / 2);
}

int main(){
    double l, r;
    scanf("%lf%lf%lf%lf%lf%lf", &a, &b, &c, &d, &l, &r);
    printf("%.6f\n", solve(l, r, 1e-7));
    return 0;
}

luoguP4526 【模板】自适应辛普森法2

题目链接

>folded
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
#include<cstdio>
#include<cmath>

using namespace std;

double a;
double f(double x){
    return pow(x, a / x - x);
}
double simpson(double l, double r){
    double mid = (l + r) / 2;
    return (f(l) + 4 * f(mid) + f(r)) * (r - l) / 6;
}
double solve(double l, double r, double eps){
    double mid = (l + r) / 2;
    double Il = simpson(l, mid), Ir = simpson(mid, r), I = simpson(l, r);
    if(fabs(Il + Ir - I) <= 15 * eps)   return Il + Ir + (Il + Ir - I) / 15;
    return solve(l, mid, eps / 2) + solve(mid, r, eps / 2);
}

int main(){
    scanf("%lf", &a);
    if(a < 0)   puts("orz");
    else    printf("%.5f\n", solve(1e-8, 20, 1e-7));
    return 0;
}

hdu1724 Ellipse

题目链接

>folded
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
#include<algorithm>
#include<cstdio>
#include<cmath>

using namespace std;

int T;
double a, b;

double f(double x){ // the function
    return 2 * sqrt(b * b * (1 - x * x / a / a));
}
double simpson(double l, double r){
    double mid = (l + r) / 2;
    return (f(l) + 4 * f(mid) + f(r)) * (r - l) / 6;
}
double solve(double l, double r, double _eps, double I){
    double mid = (l + r) / 2;
    double Il = simpson(l, mid), Ir = simpson(mid, r);
    if(fabs(Il + Ir - I) <= 15 * _eps)   return Il + Ir + (Il + Ir - I) / 15;
    return solve(l, mid, _eps / 2, Il) + solve(mid, r, _eps / 2, Ir);
}

int main(){
    scanf("%d", &T);
    while(T--){
        scanf("%lf%lf", &a, &b);
        double l, r; scanf("%lf%lf", &l, &r);
        printf("%.3f\n", solve(l, r, 1e-5, simpson(l, r)));
    }
    return 0;
}

SP8073 CIRU - The area of the union of circles

题目链接

求圆的面积并。将 $x=x_0$ 处圆覆盖的线段长度作为函数值 $f(x_0)$，这个函数值可以 $O(n\lg n)$ 求得。对这个函数 $f(x)$ 在整个区间内作 Simpson 积分即可。但是当圆的分布较为稀疏时，函数 $f(x)$ 可能有大量的 $0$ 点，在自适应 Simpson 积分时可能不会递归求解下去。所以我们选取一系列有圆覆盖的连续区间分段积分。

一些优化：小圆在大圆内部先删去；单独一个没有与其他圆有交的圆先算出答案并删去。

>folded
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
#include<algorithm>
#include<cstdio>
#include<cmath>

using namespace std;

const double eps = 1e-8;
const double PI = acos(-1);
const double INF = 1e16;
inline int sgn(double x){
    if(fabs(x) < eps)	return 0;
    else if(x > 0)	return 1;
    else	return -1;
}
inline int cmp(double x, double y){
    if(fabs(x-y) < eps)	return 0;
    else if(x > y)	return 1;
    else	return -1;
}

struct Vector{
    double x, y;
    Vector(double x = 0, double y = 0):x(x), y(y){}
    void read(){ scanf("%lf%lf", &x, &y); }
};
typedef Vector Point;
Vector operator + (Vector A, Vector B){ return Vector{A.x + B.x, A.y + B.y}; }
Vector operator - (Vector A, Vector B){ return Vector{A.x - B.x, A.y - B.y}; }
double operator * (Vector A, Vector B){ return A.x * B.x + A.y * B.y; } // dot product
double Length(Vector A){ return sqrt(A * A); }
struct Line{
    Point p; Vector v;
};
struct Circle{
    Point p;
    double r;
    bool operator < (const Circle &C) const{ return r > C.r; }
};

// ------------------------------------------------------------------------------- //

const int N = 1005;

int n;
Circle c[N];
bool b[N];
double ans;

struct Segment{
    double l, r;
    bool operator < (const Segment &A) const{ return l < A.l; }
}a[N], seg[N];

double f(double x){
    int id = 0;
    for(int i = 1; i <= n; i++){
        double d = c[i].r * c[i].r - (c[i].p.x - x) * (c[i].p.x - x);
        if(sgn(d) <= 0) continue;
        d = sqrt(d);
        a[++id] = (Segment){c[i].p.y - d, c[i].p.y + d};
    }
    sort(a+1, a+id+1);
    double res = 0, pre = -1e9;
    for(int i = 1; i <= id; i++){
        if(a[i].l > pre)    res += a[i].r - a[i].l, pre = a[i].r;
        else if(a[i].r > pre)   res += a[i].r - pre, pre = a[i].r;
    }
    return res;
}
double simpson(double l, double r){
    double mid = (l + r) / 2;
    return (f(l) + 4 * f(mid) + f(r)) * (r - l) / 6;
}
double solve(double l, double r, double _eps){
    double mid = (l + r) / 2;
    double Il = simpson(l, mid), Ir = simpson(mid, r), I = simpson(l, r);
    if(fabs(Il + Ir - I) <= 15 * _eps)   return Il + Ir + (Il + Ir - I) / 15;
    return solve(l, mid, _eps / 2) + solve(mid, r, _eps / 2);
}

int main(){
    scanf("%d", &n);
    for(int i = 1; i <= n; i++)
        scanf("%lf%lf%lf", &c[i].p.x, &c[i].p.y, &c[i].r);

    sort(c+1, c+n+1);
    int cid = 0;
    for(int i = 1; i <= n; i++){
        if(sgn(c[i].r) == 0)    continue;
        bool in = false;
        for(int j = 1; j <= cid; j++){
            if(cmp(Length(c[j].p - c[i].p), c[j].r - c[i].r) <= 0){
                in = true;
                break;
            }
        }
        if(!in) c[++cid] = c[i];
    }
    n = cid; cid = 0;

    for(int i = 1; i <= n; i++)
        for(int j = 1; j < i; j++)
            if(cmp(Length(c[i].p - c[j].p), c[i].r + c[j].r) < 0)
                b[i] = b[j] = 1;
    for(int i = 1; i <= n; i++)
        if(!b[i])   ans += PI * c[i].r * c[i].r;
        else    c[++cid] = c[i];
    n = cid;

    int id = 0;
    for(int i = 1; i <= n; i++)
        seg[++id] = (Segment){c[i].p.x - c[i].r, c[i].p.x + c[i].r};
    sort(seg+1, seg+id+1);
    double pre = -1e9;
    for(int i = 1; i <= id; i++){
        if(seg[i].l > pre)  ans += solve(seg[i].l, seg[i].r, 1e-5), pre = seg[i].r;
        else if(seg[i].r > pre) ans += solve(pre, seg[i].r, 1e-5), pre = seg[i].r;
    }
    printf("%.3f\n", ans);
    return 0;
}

[NOI2005]月下柠檬树

题目链接

圆的投影仍然是圆，圆台的投影是两圆及其外公切线。套用 Simpson 积分即可。

>folded
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
#include<algorithm>
#include<cstdio>
#include<cmath>

using namespace std;

const double eps = 1e-8;
const double PI = acos(-1);
const double INF = 1e16;
inline int sgn(double x){
    if(fabs(x) < eps)	return 0;
    else if(x > 0)	return 1;
    else	return -1;
}
inline int cmp(double x, double y){
    if(fabs(x-y) < eps)	return 0;
    else if(x > y)	return 1;
    else	return -1;
}

struct Vector{
    double x, y;
    Vector(double x = 0, double y = 0):x(x), y(y){}
    void read(){ scanf("%lf%lf", &x, &y); }
};
typedef Vector Point;
Vector operator + (Vector A, Vector B){ return Vector{A.x + B.x, A.y + B.y}; }
Vector operator - (Vector A, Vector B){ return Vector{A.x - B.x, A.y - B.y}; }
Vector operator * (Vector A, double k){ return Vector{A.x * k, A.y * k}; }
double operator * (Vector A, Vector B){ return A.x * B.x + A.y * B.y; } // dot product
double Length(Vector A){ return sqrt(A * A); }
double Angle(Vector A){ return atan2(A.y, A.x); } // polar angle of vector A
struct Line{
    Point p; Vector v;
};
struct Circle{
    Point p; double r;
    Point getPoint(double alpha){ return Point(p.x + cos(alpha) * r, p.y + sin(alpha) * r); }
};
void getTangents(Circle C1, Circle C2, Point c1[], Point c2[], int &resn){
    // resn is the number of tangent lines
    // c1[] and c2[] are relevant points on C1 and C2
    resn = 0;
    if(cmp(C1.r, C2.r) < 0)   swap(C1, C2), swap(c1, c2);
    double d = Length(C1.p - C2.p);
    if(sgn(d) == 0 && cmp(C1.r, C2.r) == 0){ resn = -1; return; }  // two circles are the same
    if(cmp(C1.r - C2.r, d) > 0) return;                         // contained
    double base = Angle(C2.p - C1.p);
    if(cmp(C1.r - C2.r, d) == 0){                               // internally tangent
        c1[++resn] = C1.getPoint(base), c2[resn] = C2.getPoint(base);
        return;
    }
    double ang = acos((C1.r - C2.r) / d);
    c1[++resn] = C1.getPoint(base - ang), c2[resn] = C2.getPoint(base - ang);
    c1[++resn] = C1.getPoint(base + ang), c2[resn] = C2.getPoint(base + ang);
//    if(cmp(C1.r + C2.r, d) == 0)                                // externally tangent
//        c1[++resn] = C1.getPoint(base), c2[resn] = C2.getPoint(base + PI);
//    else if(cmp(C1.r + C2.r, d) < 0){							// separated
//        ang = acos((C1.r + C2.r) / d);
//        c1[++resn] = C1.getPoint(base - ang), c2[resn] = C2.getPoint(base - ang + PI);
//        c1[++resn] = C1.getPoint(base + ang), c2[resn] = C2.getPoint(base + ang + PI);
//    }
}

// ------------------------------------------------------------------------------- //

const int N = 1005;

int n;
Circle c[N];
double alpha, H[N], R[N], mxR, mnL = 1e9;
pair<Point, Point> seg[N]; int sid;

double f(double x){
    double mx = 0, res = 0;
    for(int i = 1; i <= sid; i++){
        if(cmp(seg[i].first.x, x) <= 0 && cmp(seg[i].second.x, x) >= 0){
            Point p = seg[i].first + (seg[i].second - seg[i].first) * ((x - seg[i].first.x) / (seg[i].second.x - seg[i].first.x));
            mx = max(mx, p.y);
        }
    }
    for(int i = 1; i <= n; i++)
        if(cmp(c[i].r, fabs(x - c[i].p.x)) > 0)
            mx = max(mx, sqrt(c[i].r * c[i].r - fabs(x - c[i].p.x) * fabs(x - c[i].p.x)));
    res += mx;
    return res;
}
double simpson(double l, double r){
    double mid = (l + r) / 2;
    return (f(l) + 4 * f(mid) + f(r)) * (r - l) / 6;
}
double solve(double l, double r, double _eps, double I){
    double mid = (l + r) / 2;
    double Il = simpson(l, mid), Ir = simpson(mid, r);
    if(fabs(Il + Ir - I) <= 15 * _eps)   return Il + Ir + (Il + Ir - I) / 15;
    return solve(l, mid, _eps / 2, Il) + solve(mid, r, _eps / 2, Ir);
}

int main(){
    scanf("%d%lf", &n, &alpha);
    double Tan = 1 / tan(alpha);
    for(int i = 1; i <= n + 1; i++){
        scanf("%lf", &H[i]);
        H[i] += H[i-1];
    }
    for(int i = 1; i <= n + 1; i++){
        if(i <= n)  scanf("%lf", &R[i]);
        c[i].p = Point(H[i] * Tan, 0);
        c[i].r = R[i];
        mnL = min(mnL, c[i].p.x - c[i].r);
        mxR = max(mxR, c[i].p.x + c[i].r);
        if(i > 1){
            Point p1[10], p2[10]; int pid;
            getTangents(c[i-1], c[i], p1, p2, pid);
            for(int j = 1; j <= pid; j++){
                if(sgn(p1[j].y) <= 0)   continue;
                if(p1[j].x > p2[j].x)   swap(p1[j], p2[j]);
                seg[++sid] = make_pair(p1[j], p2[j]);
            }
        }
    }
    printf("%.2f\n", 2 * solve(mnL, mxR, 1e-6, simpson(mnL, mxR)));
    return 0;
}

自适应Simpson积分学习笔记

http://xyfjason.github.io/blog-xcpc/2020/03/18/自适应Simpson积分学习笔记/

作者

xyfJASON

发布于

2020-03-18

更新于

2021-02-25

许可协议

#

评论

关注我

分类

最新文章

标签

AtCoder2
CCPC2
Codeforces57
EOJ1
GMCPC1
HDU11
ICPC3
Nowcoder12
POI1
hihoCoder1
学习笔记17
随笔1
×