自适应Simpson积分笔记 Adaptive Simpson

发表于2020-08-11更新于2023-09-15字数统计4k阅读时长42分

Adaptive Simpson

辛普森公式

SimpsonSimpson 积分是积分数值计算的一种方法，用抛物线而非传统的矩形来近似图形，精度较高，公式为：

∫abf(x) dx≈b−a6n∑i=1n[f(xi)+4f(xi+xi−12)+f(xi−1)]\int_{a}^{b}f(x)\,dx \approx \frac{b - a}{6n} \sum_{i = 1}^{n} [f(x_i) + 4f(\frac{x_i + x_{i - 1}}{2}) + f(x_{i - 1})]

具体证明见 学习笔记 自适应辛普森(Simpson)积分

实际计算中不好判断划分数 nn，因此提出了自适应辛普森积分，即在积分数值在可接受范围内直接返回积分值，否则对区间二分递归求解，即：

∫abf(x) dx≈b−a6[f(a)+4f(a+b2)+f(b)]\int_{a}^{b}f(x)\,dx \approx \frac{b - a}{6} [f(a) + 4f(\frac{a + b}{2}) + f(b)]

模板

不能 O(1)O(1) 求出 f(x)f(x) 时，保存已有结果能降低时间复杂度，进而可提高精度。可用unordered_map实现优化。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
unordered_map<int, int> mp;
//f根据要求不同实现不同
inline long double simpson(long double l, long double r) {
    double mid = (l + r) / 2;
    return (f(l) + f(mid) * 4 + f(r)) * (r - l) / 6;
}

long double integral(long double l, long double r, long double eps， double i) {
    long double mid = (l + r) / 2;
    long double iL = simpson(l, mid), iR = simpson(mid, r);
    if (fabs(iL + iR - i) <= eps * 15) return iL + iR + (iL + iR - i) / 15;
    return integral(l, mid, eps / 2, iL) + integral(mid, r, eps / 2, iR);
}

integral(l, r, 1e-5, simpson(l, r));

面积并

自适应辛普森积分可以用来解决一些面积并问题，即平面上几何图形覆盖面积的总和。

求面积并：对于 xx，求出该处覆盖的线长 (xx 与图形相交得到一系列线段，按左端点排序计算可得)，那么总线长就是面积并了 (联想二重积分)。

注意：保证线长是连续的，如果不连续需要分段求。此外还需注意区间范围，过大无法保证精度，可能需要其它方法 (如三角形面积并一般用扫描线)。

积分例题

A 自适应辛普森法1

Description

题目链接

Code

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
#include <bits/stdc++.h>

using namespace std;

long double a, b, c, d, L, R;

unordered_map<long double, long double> mp;
inline long double f(long double x) {
    return (c * x + d) / (a * x + b);
}

inline long double simpson(long double l, long double r) {
    double mid = (l + r) / 2;
    return (f(l) + f(mid) * 4 + f(r)) * (r - l) / 6;
}

long double integral(long double l, long double r, long double eps) {
    long double mid = (l + r) / 2;
    long double iL = simpson(l, mid), iR = simpson(mid, r);
    long double i = simpson(l ,r);
    if (fabs(iL + iR - i) <= eps * 15) return iL + iR + (iL + iR - i) / 15;
    return integral(l, mid, eps / 2) + integral(mid, r, eps / 2);
}

int main () {
    std::ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
    cout << setprecision(6) << fixed;
    cin >> a >> b >> c >> d >> L >> R;
    cout << integral(L, R, 1e-7) << endl;
    return 0;
}

B 自适应辛普森法2

Description

题目链接

Solution

推导可得 a<0a < 0 时发散。a≥0a \ge 0 时， x<20x < 20 时就已收敛

Code

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 <bits/stdc++.h>

using namespace std;

long double a;

unordered_map<long double, long double> mp;
inline long double f(long double x) {
    return pow(x, a / x - x);
}

inline long double simpson(long double l, long double r) {
    double mid = (l + r) / 2;
    return (f(l) + f(mid) * 4 + f(r)) * (r - l) / 6;
}

long double integral(long double l, long double r, long double eps) {
    long double mid = (l + r) / 2;
    long double iL = simpson(l, mid), iR = simpson(mid, r);
    long double i = simpson(l ,r);
    if (fabs(iL + iR - i) <= eps * 15) return iL + iR + (iL + iR - i) / 15;
    return integral(l, mid, eps / 2) + integral(mid, r, eps / 2);
}

int main () {
    std::ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
    cout << setprecision(5) << fixed;
    cin >> a;
    if (a < -1e-8) cout << "orz" << endl;
    else cout << integral(1e-8, 20, 1e-7) << endl;
    return 0;
}

面积并例题

A SP8073 CIRU - The area of the union of circles

Description

题目链接
求圆的面积并

Solution

注意到目标函数值不连续，可能存在大量0，故分段进行辛普森积分。

可以通过判断圆的相交求出联通块，对联通块上的圆求面积并。或者将圆的左右端点坐标存入 vector 排序，每次取相邻且不同的坐标 l,rl, r，对 [l,r][l, r] 进行积分。

本题 f(x)f(x) 不能在常数时间获得，可用 unordered_map 优化。

Code

The area of the union of circles
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
#include <bits/stdc++.h>

using namespace std;

constexpr int MAXN = 1e3 + 3;
constexpr long double EPS = 1e-10;

int n, vis[MAXN];
vector<int> chosen, G[MAXN];

struct Circle {
    /*@param x x of center O @param y y of center O @param R radius*/
    long double x, y, R;
    /*@param l 最左端 @param r 最右端*/
    long double l, r;
    /*按左端点排序 相同则按右端点排序*/
    friend inline bool operator<(const Circle& a, const Circle& b) {
        return a.l < b.l - EPS;
    }
} circle[MAXN];

struct Seg {
    /*@param l 左端点 @param r 右端点*/
    long double l, r;
    /*按左端点排序 相同则按右端点排序*/
    friend inline bool operator<(const Seg& a, const Seg& b) {
         return a.l < b.l - EPS;
    }
} seg[MAXN];

unordered_map<double, double> mp;
inline long double f(long double x) { //算x处覆盖的线段总长
    int tot = 0;
    for (auto& i : chosen) {
        long double tmp = pow(circle[i].R, 2) - pow(x - circle[i].x, 2);
        if (tmp < EPS) continue;
        seg[++tot].l = circle[i].y - sqrt(tmp);
        seg[tot].r = circle[i].y + sqrt(tmp);
    }
    sort(seg + 1, seg + tot + 1);
    long double l = seg[1].l, len = .0;
    for (int i = 1; i <= tot; i++) {
        l = max(l, seg[i].l);
        if (l < seg[i].r - EPS) len += seg[i].r - l;
        l = max(l, seg[i].r);
    }
    return len;
}

inline long double simpson(long double l, long double r) {
    long double mid = (l + r) / 2;
    return (f(l) + f(mid) * 4 + f(r)) * (r - l) / 6;
}

long double integral(long double l, long double r, long double eps) {
    long double mid = (l + r) / 2;
    long double iL = simpson(l, mid), iR = simpson(mid, r);
    long double i = simpson(l ,r);
    if (fabs(iL + iR - i) < eps * 15) return iL + iR + (iL + iR - i) / 15;
    return integral(l, mid, eps / 1.5) + integral(mid, r, eps / 1.5);
}

inline void init() {
    for (int i = 1; i <= n; i++) { //circle已按l排序
        if (vis[i] || circle[i].R < EPS) { //半径为0跳过
            vis[i] = true; continue;
        }
        for (int j = i + 1; j <= n; j++) {
            if (vis[j] || circle[j].R < EPS) { //半径为0跳过
                vis[j] = true; continue;
            }
            long double tmp = pow(circle[i].x - circle[j].x, 2)
                            + pow(circle[i].y - circle[j].y, 2);
            if (tmp < pow(circle[i].R - circle[j].R, 2) + EPS) {//被包含
                vis[j] = true; continue; //包含关系不加边
            }
            if (tmp < pow(circle[i].R + circle[j].R, 2) + EPS)
                G[i].push_back(j), G[j].push_back(i); //i与j联通
        }
    }
}

inline void dfs(int cur, long double& r) { //处理连通块
    vis[cur] = true, chosen.push_back(cur);
    for (auto& to : G[cur]) {
        if (vis[to]) continue;
        dfs(to, r = max(r, circle[to].r)); //求最大右端点
    }
}

inline long double calc() { //计算面积并
    long double l, r, ans = .0;
    for (int i = 1; i <= n; i++) {
        if (vis[i]) continue;
        chosen.clear();
        l = circle[i].l, r = circle[i].r;
        dfs(i, r); //处理连通块
        if (chosen.size() == 1) ans += acos(-1) * circle[i].R * circle[i].R;
        else ans += integral(l, r, 1e-5);
    }
    return ans;
}

int main () {
    std::ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
    cout << setprecision(3) << fixed;
    cin >> n;
    for (int i = 1; i <= n; i++) {
        cin >> circle[i].x >> circle[i].y >> circle[i].R;
        circle[i].l = circle[i].x - circle[i].R;
        circle[i].r = circle[i].x + circle[i].R;
    }
    sort(circle + 1, circle + n + 1);
    init();
    cout << calc();
    return 0;
}

B NOI2005 月下柠檬树

Prolbem

题目链接

Solution

投影为一系列圆以及梯形构成，其中梯形的斜边为相邻圆的公切线 (可能不存在)。

切点求法：两圆心分别与对应的切点连线，圆心向连线做垂线，切点向轴做垂线，通过相似三角形求出。

注意：由于投影具有对称性，只需对一半进行积分，否则可能超时。

Code

NOI2005 月下柠檬树
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
#include <bits/stdc++.h>

using namespace std;

constexpr int MAXN = 3e3 + 3;
constexpr double EPS = DBL_EPSILON;

int n;
double h[MAXN], sumH, alpha;

struct Circle {
    /*@param x x of center O @param y y of center O @param R radius*/
    double x, y, R;
    /*@param l 最左端 @param r 最右端*/
    double l, r;
} circle[MAXN];

struct Seg {
    double l, r;
    friend inline bool operator<(const Seg& a, const Seg& b) {
        return a.l < b.l - EPS;
    }
} seg[MAXN];

struct Trapezoid {      //梯形类
    double x1, y1; //右上顶点
    double x2, y2; //右下顶点
    double deltaX = 0, height = 0;
} trape[MAXN];

unordered_map<double, double> mp;
inline double f(double x) {  //算x处覆盖的线段总长
    if (mp.count(x)) return mp[x];
    int tot = 0;
    for (int i = 1; i <= n; i++) { //圆形
        double tmp = pow(circle[i].R, 2) - pow(x - circle[i].x, 2);
        if (tmp < EPS) continue;
        seg[++tot].l = circle[i].y - sqrt(tmp);
        seg[tot].r = circle[i].y + sqrt(tmp);
    }
    for (int i = 1; i <= n; i++) {  //梯形
        double l = min(trape[i].x1, trape[i].x2);
        double r = max(trape[i].x1, trape[i].x2);
        if (r < x + EPS) continue;
        if (x < l + EPS) {
            seg[++tot].l = trape[i].y2, seg[tot].r = trape[i].y1;
            continue;
        }
        double delta = trape[i].x2 - x;
        double ratio = abs(delta / trape[i].deltaX);
        if (trape[i].deltaX > EPS) {
            seg[++tot].r = trape[i].y2 + trape[i].height * ratio;
            seg[tot].l = trape[i].y2;
        } else if (trape[i].deltaX < EPS) {
            seg[++tot].l = trape[i].y2 + trape[i].height * ratio;
            seg[tot].r = trape[i].y1;
        }
    }
    sort(seg + 1, seg + tot + 1);
    double len = .0, l = seg[1].l;
    for (int i = 1; i <= tot; i++) {
        l = max(l, seg[i].l);
        if (l < seg[i].r - EPS) len += seg[i].r - l;
        l = max(l, seg[i].r);
    }
    return mp[x] = len;
}

inline double simpson(double l, double r) {
    double mid = (l + r) / 2;
    return (f(l) + f(mid) * 4 + f(r)) * (r - l) / 6;
}

double integral(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) <= eps * 15) return iL + iR + (iL + iR - i) / 15;
    return integral(l, mid, eps / 1.6, iL) + integral(mid, r, eps / 1.6, iR);
}

inline double calc() {  //计算面积并
    for (int i = 1; i <= n; i++) {
        double dis = circle[i + 1].y - circle[i].y;  //圆心距离
        if (dis < abs(circle[i].R - circle[i + 1].R) + EPS) continue;
        double theta = acos((circle[i].R - circle[i + 1].R) / dis);
        trape[i].x1 = circle[i + 1].R * sin(theta); //右上顶点
        trape[i].y1 = circle[i + 1].y + circle[i + 1].R * cos(theta);
        trape[i].x2 = circle[i].R * sin(theta); //右下顶点
        trape[i].y2 = circle[i].y + circle[i].R * cos(theta);
        trape[i].deltaX = trape[i].x2 - trape[i].x1; //
        trape[i].height = trape[i].y1 - trape[i].y2; //高度
    }
    double l = circle[1].l, r = circle[1].r;
    for (int i = 2; i <= n; i++)
        l = min(l, circle[i].l), r = max(r, circle[i].r);
    return integral(0, r, 5e-4, simpson(l, r)) * 2; //对称性 对1边积分
}

int main() {
    std::ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
    cout << setprecision(2) << fixed;
    cin >> n >> alpha;
    for (int i = 0; i <= n; i++) cin >> h[i];
    sumH = h[0];
    for (int i = 1; i <= n; i++) {
        cin >> circle[i].R;
        circle[i].x = 0, circle[i].y = sumH / tan(alpha);
        circle[i].l = -circle[i].R, circle[i].r = circle[i].R;
        sumH += h[i];
    }
    circle[n + 1].x = 0, circle[n + 1].y = sumH / tan(alpha);
    circle[n + 1].l = 0, circle[n + 1].r = 0;  //顶点
    cout << calc() << endl;
    return 0;
}

C 三角形

Description

题目链接
求等腰直角三角形的面积并

Solution

注意到即使三角形有交集，目标函数也会不连续。故不能求连通块，而是改为将三角形的左右端坐标存入 vector 排序，分段积分。

Code

三角形
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
#include <bits/stdc++.h>

using namespace std;

constexpr int MAXN = 2e3 + 3;
constexpr double EPS = DBL_EPSILON;

struct Triangle {
    /*@param x1 上顶x坐标 @param y1 上顶y坐标*/
    int x1, y1;
    /*@param x2 下顶x坐标 @param y2 下顶y坐标*/
    int x2, y2;
    friend inline bool operator<(const Triangle& a, const Triangle& b) {
        return a.x1 < b.x1;
    }
} tri[MAXN];

int n, vis[MAXN], st, ed;
vector<int> G[MAXN], posX, chosen;
pair<double, double> seg[MAXN];

unordered_map<double, double> mp;
inline double f(double x) {
    if (mp.count(x)) return mp[x];
    int tot = 0;
    for (auto& i : chosen) {
        if (x < -EPS + tri[i].x1 || tri[i].x2 < -EPS + x) continue;
        seg[++tot].first = 1.0 * tri[i].y2;
        seg[tot].second = 1.0 * tri[i].y2 + tri[i].x2 - x;
    }
    sort(seg + 1, seg + tot + 1);
    double len = 0, l = seg[1].first;
    for (int i = 1; i <= tot; i++) {
        l = max(l, seg[i].first);
        if (l < seg[i].second - EPS) len += seg[i].second - l;
        l = max(l, seg[i].second);
    }
    return mp[x] = len;
}

inline double simpson(double l, double r) {
    double mid = (l + r) / 2;
    return (f(l) + f(mid) * 4 + f(r)) * (r - l) / 6;
}

double integral(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) <= eps * 15) return iL + iR + (iL + iR - i) / 15;
    return integral(l, mid, eps / 2, iL) + integral(mid, r, eps / 2, iR);
}

inline void init(double res = 0) {
    sort(tri + 1, tri + n + 1);
    sort(posX.begin(), posX.end());
    auto edIt = unique(posX.begin(), posX.end());
    for (auto it = posX.begin() + 1; it != edIt; it++) {
        double l = *(it - 1), r = *(it);
        chosen.clear();
        for (int i = 1; i <= n; i++) {
            if ((tri[i].x1 <= l && l <= tri[i].x2) ||
                (tri[i].x1 <= r && r <= tri[i].x2)) {
                chosen.push_back(i);
            }
        }
        l += 2e-9, r -= 2e-9;
        res += integral(l, r, 1e-4, simpson(l, r));
    }
    cout << setprecision(1) << fixed << res << endl;
}

int main () {
    std::ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
    cin >> n;
    for (int i = 1, x, y, m; i <= n; i++) {
        cin >> x >> y >> m;
        tri[i].x1 = x, tri[i].y1 = y + m;
        tri[i].x2 = x + m, tri[i].y2 = y;
        posX.push_back(x), posX.push_back(x + m);
    }
    init();
    return 0;
}

D HNOI2012 三角形覆盖问题

Description

题目链接

Solution

正解为扫面线，但依然可以用自适应 simpson 积分，只不过需要直接设置 eps，而不是每次递归的时候折半。此外还需判断一个三角形是不是被另一个三角形完全覆盖，是则可以将其跳过。

Code

三角形覆盖问题
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
#include <bits/stdc++.h>

using namespace std;

constexpr long double EPS = 1e-13;
constexpr int MAXN = 1e4 + 5;

int n, x, y, m;
vector<int> chosen, pos;
bool vis[MAXN];

struct Triangle {
    int x, y, x2, y2;
    friend inline bool operator<(const Triangle& a, const Triangle& b) {
        return a.x < b.x;
    }
} tri[MAXN];

struct Segment {
    long double l, r;
    friend bool operator<(const Segment& a, const Segment& b) {
        return a.l < b.l;
    }
} seg[MAXN];

unordered_map<long double, long double> mp;
inline long double f(long double x) {
    if (mp.count(x)) return mp[x];
    int tot = 0;
    for (auto& i : chosen) {
        if (x < (long double)tri[i].x - EPS || (long double)tri[i].x2 < x - EPS) continue;
        seg[++tot].l = (long double)tri[i].y2;
        seg[tot].r = (long double)tri[i].y - x + (long double)tri[i].x;
    }
    sort(seg + 1, seg + tot + 1);
    long double l = seg[1].l, len = .0;
    for (int i = 1; i <= tot; i++) {
        l = max(l, seg[i].l);
        if (seg[i].r - l > EPS) len += seg[i].r - l;
        l = max(l, seg[i].r);
    }
    return mp[x] = len;
}

long double simpson(long double l, long double r) {
    long double mid = (l + r) / 2;
    return (f(l) + f(mid) * 4 + f(r)) * (r - l) / 6;
}

long double integral(long double l, long double r, long double eps, long double i) {
    long double mid = (l + r) / 2;
    long double iL = simpson(l, mid), iR = simpson(mid, r);
    if (fabs(iR + iL - i) <= eps) return iL + iR + fabs(iR + iL - i);
    return integral(l, mid, eps, iL) + integral(mid, r, eps, iR);
}

void init() {
    for (int i = 1; i <= n; i++) {
        if (vis[i]) continue;
        if (tri[i].x2 - tri[i].x == 0) {
            vis[i] = 1;
            continue;
        }
        for (int j = i + 1; j <= n; j++) {
            if (vis[j]) continue;
            if (tri[j].x2 - tri[j].x == 0) {
                vis[j] = 1;
                continue;
            }
            if (tri[j].y <= tri[i].y && tri[j].x >= tri[i].x &&
                tri[j].x2 <= tri[i].x2 && tri[j].y2 >= tri[i].y2 &&
                tri[j].y <= tri[i].y - (tri[j].x - tri[i].x)) {
                vis[j] = 1;
            }
        }
    }
}

void calc() {
    sort(pos.begin(), pos.end());
    sort(tri + 1, tri + n + 1);
    auto edIt = unique(pos.begin(), pos.end());
    init();
    long double ans = .0;
    for (auto it = pos.begin() + 1; it != edIt; it++) {
        int l = *(it - 1), r = *it;
        chosen.clear();
        for (int i = 1; i <= n; i++) {
            if (vis[i]) continue;
            if (tri[i].x >= r) break;
            if ((tri[i].x <= l && l <= tri[i].x2) ||
                (tri[i].x <= r && r <= tri[i].x2))
                chosen.emplace_back(i);
        }
        if (!chosen.size()) continue;
        long double L = 2 * EPS + l, R = - 2 * EPS + r;
        ans += integral(L, R, 1e-14, simpson(L, R));
    }
    cout << setprecision(1) << fixed << ans;
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(0), cout.tie(0);
    cin >> n;
    for (int i = 1; i <= n; i++) {
        cin >> x >> y >> m;
        tri[i].x = x, tri[i].y = y + m;
        tri[i].x2 = x + m, tri[i].y2 = y;
        pos.emplace_back(x), pos.emplace_back(x + m);
    }
    calc();
    return 0;
}