Chlience

【算法】辛普森积分法
怎么求图形面积并呢?比如说在 三角形面积并 这道题中,可以通过扫描线的方法算出交点,然后直接用梯形面积公式来算但是...
扫描右侧二维码阅读全文
18
2019/01

【算法】辛普森积分法

怎么求图形面积并呢?

比如说在 三角形面积并 这道题中,可以通过扫描线的方法算出交点,然后直接用梯形面积公式来算
但是如果是 圆的面积并 这样的题,就没有一个很快速的方法去求出每一段内的面积

一般来说,按照积分的思想,可以用很多东西来拟合图像面积,比如说矩形或者是梯形
在将原图形划分为无穷多块时,显然是可以计算出精确的面积的
但是计算机没有无限微分的能力,若区间划分太小便会导致无法在规定时间内出解

所以说我们需要一个在相对而言比较大的区间划分方案中能够尽可能精确的拟合图形面积的方法
辛普森积分法就是这样一种工具

辛普森法则(Simpson's rule)是一种数值积分方法,是牛顿-寇次公式的特殊形式,以二次曲线近的方式取代矩形或梯形积分公式,以求得定积分的数值近似解。其近似值如下:

$$ \int_{a}^bf(x)dx\approx\frac{b-a}{6}\left[f(a)+4f(\frac{a+b}{2})+f(b)\right] $$

这个公式非常的妙啊,在一定的精度范围下它能够帮我们很好的拟合图像的面积
有多妙呢?它甚至能够做到精确的求某些图形的面积/体积,具体可以去 辛普森积分法 查看

如何使用这个公式?

要计算 $[a,b]$ 区间内的面积,只需要三个值:$f(a),f(b),f(c)$,即用 $x=a,b,c$ 分别去切这个图像,直线在图形内部的长度
当然,如果是求体积就是切出 $x=a,b,c$ 时图形的面积
至于如何计算 $f(x)$,只需要求出每个规则的图形和直线/平面的相交位置,然后排序暴力合并

一般来说,在 $OI$ 中会使用自适应 $Simpson$ 法

假设当前区间 $[l,r]$

若 $getSimpson(l,r)$ 和 $getSimpson(l,mid)+getSinpson(mid,r)$ 差距比较小,或者说在 $eps$ 范围内时,直接返回答案
否则递归向下寻找答案

并且值得注意的是,如果递归层数太小,区间太大,请强制向下递归几层,防止被卡

然后就可以愉快的做题了

辛普森积分法虽然是个很有用的工具,但是毕竟求得的仅仅是一个近似解
还是要注意精度问题啊!

BZOJ 1845 [CQOI2005]三角形面积并

#include <bits/stdc++.h>
using namespace std;
#define mid (l + r) / 2
const int N = 110;
const double eps = 1e-8;
typedef pair <double , double> Pair;
struct Vector {
    double x, y;
    Vector operator + (const Vector a) const {return Vector{x + a.x, y + a.y};}
    Vector operator - (const Vector a) const {return Vector{x - a.x, y - a.y};}
    double operator * (const Vector a) const {return x * a.y - y * a.x;}
    Vector operator * (double a) {return Vector{x * a, y * a};}
};
struct Triangle {
    Vector p[3];
    Vector& operator [] (int x) {return p[x];}
}t[N];
Pair sta[N];
int n, cnt, top;

int read();
bool crossCheck(Vector, Vector, Vector, Vector);
bool parallelCheck(Vector, Vector, Vector, Vector);
Vector crossPos(Vector, Vector, Vector, Vector);
double f(double);
double calSimpson(double, double);
double simpson(double, double);

int main() {
    scanf("%d", &n);
    double min_x = 1000000, max_x = - 1000000;
    for(int i = 1; i <= n; ++ i)
        for(int j = 0; j < 3; ++ j) {
            scanf("%lf %lf", &t[i][j].x, &t[i][j].y);
            max_x = max(max_x, t[i][j].x);
            min_x = min(min_x, t[i][j].x);
        }
    printf("%.2f\n", simpson(min_x, max_x) - eps);
}
bool crossCheck(Vector s1, Vector e1, Vector s2, Vector e2) {
    if(max(s1.x, e1.x) < min(s2.x, e2.x) || max(s2.x, e2.x) < min(s1.x, e1.x)) return 0;
    if(max(s1.y, e1.y) < min(s2.y, e2.y) || max(s2.y, e2.y) < min(s1.y, e1.y)) return 0;
    if(((e1 - s1) * (s2 - s1)) * ((e1 - s1) * (e2 - s1)) > 0) return 0;
    if(((e2 - s2) * (s1 - s2)) * ((e2 - s2) * (e1 - s2)) > 0) return 0;
    return 1;
}
bool parallelCheck(Vector s1, Vector e1, Vector s2, Vector e2) {
    return fabs((e1 - s1) * (e2 - s2)) < eps;
}
Vector crossPos(Vector p1, Vector v1, Vector p2, Vector v2) {
    return v1 * (((p2 - p1) * v2) / (v1 * v2)) + p1;
}
double f(double x) {
    top = 0;
    for(int i = 1; i <= n; ++ i) {
        double U = - 1000000, D = 1000000;
        Vector a = Vector{x, - 1000000};
        Vector b = Vector{x, 1000000};
        for(int j = 0; j < 3; ++ j) {
            if(crossCheck(t[i][j], t[i][(j + 1) % 3], a, b) && !parallelCheck(t[i][j], t[i][(j + 1) % 3], a, b)) {
                Vector p = crossPos(t[i][j], t[i][(j + 1) % 3] - t[i][j], a, b - a);
                U = max(U, p.y); D = min(D, p.y);
            }
        }
        if(U >= D)
            sta[++ top] = {D , U};
    }
    sort(sta + 1, sta + top + 1);
    double ans = 0 , last = - 1000000;
    for(int i = 1; i <= top; ++ i) {
        if(last >= sta[i].second) continue;
        ans += sta[i].second - max(sta[i].first, last);
        last = sta[i].second;
    }
    return ans;
}
double calSimpson(double l, double r) {return (f(l) + f(r) + 4 * f(mid)) * (r - l) / 6;}
double simpson(double l, double r) {
    if(r - l > 1.0)
        return simpson(l, mid) + simpson(mid, r);
    double L = calSimpson(l, mid), R = calSimpson(mid, r), ANS = calSimpson(l, r);
    if(fabs(L + R - ANS) <= eps) return ANS;
    else return simpson(l, mid) + simpson(mid, r);
}
Last modification:January 19th, 2019 at 09:57 pm

Leave a Comment