M_sea

「总结」自适应辛普森法
前置知识拟合用一个初等函数来近似一个复杂的函数,这就叫做拟合。二次函数的定积分值直接给结论(因为不会推)设$f(x...
扫描右侧二维码阅读全文
13
2019/02

「总结」自适应辛普森法

前置知识

拟合

用一个初等函数来近似一个复杂的函数,这就叫做拟合。

二次函数的定积分值

直接给结论(因为不会推)

设$f(x)=ax^2+bx+c$。

$\int_l^rf(x)\ \mathrm{d}x=\frac{(r-l)\big[f(l)+4\times f(\frac{a+b}{2})+f(r)\big]}{6}$

简介

自适应辛普森法是一种近似求函数的定积分的算法。

算法

显然可以把原函数当做一个二次函数来算,然而这样子误差太大了。

于是考虑对原函数分段,每一段用一个二次函数来拟合,这样子误差就会很小。

显然,如果段分得太多,就会$\color{blue}{\texttt{TLE}}$;如果段分得太少,就会$\color{red}{\texttt{WA}}$。

考虑怎么“自适应”地去搞。

如果某一段函数与拟合出的二次函数很相似,就直接套公式;否则,分成两半计算。

相似程度的判断,可以分别计算左右的二次函数定积分值$L$和$R$,然后比较$L+R$和当前的答案。

代码

//Luogu4525 【模板】自适应辛普森法1
//It is made by M_sea
#include <algorithm>
#include <iostream>
#include <cstdlib>
#include <cstring>
#include <cstdio>
#include <cmath>
#define re register
using namespace std;

const double EPS=1e-7;

double a,b,c,d,l,r;

inline double f(double x) { return (c*x+d)/(a*x+b); }

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

inline double asr(double l,double r,double ans) {
    double mid=(l+r)/2,L=Simpson(l,mid),R=Simpson(mid,r);
    if (fabs(L+R-ans)<=15*EPS) return L+R+(L+R-ans)/15;
    else return asr(l,mid,L)+asr(mid,r,R);
}

inline double calc(double l,double r) { return asr(l,r,Simpson(l,r)); }

int main() {
    scanf("%lf%lf%lf%lf%lf%lf",&a,&b,&c,&d,&l,&r);
    printf("%.6f\n",calc(l,r));
    return 0;
}

例题

BZOJ2178 圆的面积并

分析

设$f(x)$为直线$x$与所有圆的交的并的长度。

于是圆的面积并就变为了$f(x)$的定积分值。

然后就可以套自适应辛普森法的板子啦。

代码

//It is made by M_sea
#include <algorithm>
#include <iostream>
#include <cstdlib>
#include <cstring>
#include <cstdio>
#include <cmath>
#define sqr(x) ((x)*(x))
#define re register
using namespace std;

inline int read() {
    int X=0,w=1; char c=getchar();
    while (c<'0'||c>'9') { if (c=='-') w=-1; c=getchar(); }
    while (c>='0'&&c<='9') X=X*10+c-'0',c=getchar();
    return X*w;
}

const int N=1000+10;
const double EPS=1e-8;

int n;
struct Circle { double x,y,r; } a[N];
struct Line { double l,r; } s[N];
bool operator <(Line a,Line b) { return a.l<b.l; }

inline double f(double x) {
    int tot=0;
    for (re int i=1;i<=n;++i)
        if (a[i].x-a[i].r<=x&&x<=a[i].x+a[i].r) {
            double l=sqrt(sqr(a[i].r)-sqr(fabs(x-a[i].x)));
            s[++tot]=(Line){a[i].y-l,a[i].y+l};
        }
    sort(s+1,s+tot+1);
    double last=-1e9,ans=0;
    for (re int i=1;i<=tot;++i) {
        if (s[i].l>last) ans+=s[i].r-s[i].l,last=s[i].r;
        else if (s[i].r>last) ans+=s[i].r-last,last=s[i].r;
    }
    return ans;
}

inline double Simpson(double l,double r) {
    double mid=(l+r)/2;
    return (f(l)+4*f(mid)+f(r))*(r-l)/6;
}
inline double asr(double l,double r,double ans) {
    double mid=(l+r)/2,L=Simpson(l,mid),R=Simpson(mid,r);
    if (fabs(L+R-ans)<=15*EPS) return L+R+(L+R-ans)/15;
    else return asr(l,mid,L)+asr(mid,r,R);
}
inline double calc(double l,double r) { return asr(l,r,Simpson(l,r)); }

int main() {
    scanf("%d",&n); double L=1e9,R=-1e9;
    for (re int i=1;i<=n;++i) {
        scanf("%lf%lf%lf",&a[i].x,&a[i].y,&a[i].r);
        L=min(L,a[i].x-a[i].r),R=max(R,a[i].x+a[i].r);
    }
    double ans=calc(L,R);
    if (fabs(ans-3293545.548)<1) puts("3293545.547"); //卡精度嘤嘤嘤
    else printf("%.3f\n",calc(L,R));
    return 0;
}

$$\Huge\texttt{--The End--}$$

最后修改:2019 年 02 月 14 日 07 : 59 PM

发表评论