码迷,mamicode.com
首页 > 其他好文 > 详细

并不对劲的辛普森积分

时间:2018-03-06 20:21:29      阅读:231      评论:0      收藏:0      [点我收藏+]

标签:tin   space   desc   inpu   else   lap   max   stream   多个   

并不会计算几何的并不对劲的人想必是非常之菜的。

很对劲的太刀流在这里->技术分享图片

辛普森积分本身是把一段函数f(x)当成二次函数算,用(f(l)+4*f((l+r)/2)+f(r))*(r-l)/6算出这段曲线在l到r之间与x轴所夹面积。

但是这听上去非常不精确,所以有人发明了自适应的辛普森积分。   

每一次分别算出这段左一半、右一半和整个的辛普森积分,如果左+右与整个之差大于eps就继续递归处理左、右半部分。

至于为什么不把一段函数当成一次函数算,是因为会被f(m)恰好在直线上这种数据卡:

技术分享图片

而二次函数只能用这种手形数据卡:

技术分享图片

那为什么不用三次函数呢?想必是因为三次函数复杂的让人手不健康而且手形数据不是那么好构造吧。

————————例题分界线————————
bzoj2178: 圆的面积并
Time Limit: 20 Sec Memory Limit: 259 MB
Submit: 2415 Solved: 617
[Submit][Status][Discuss]
Description
给出N个圆,求其面积并
Input
先给一个数字N ,N< = 1000 接下来是N行是圆的圆心,半径,其绝对值均为小于1000的整数
Output
面积并,保留三位小数
————————题解分界线—————————
可以把直线x=x在圆中的部分的长度看成是关于x1的函数。对这个函数求自适应的辛普森积分,但是复杂度非常之玄学。

经某些人的实验,发现eps设为1e-14会TLE,设为1e-12会WA,所以好像没什么可选的…

注意要删去所有被包含的圆、加一堆常数优化才能过。

 

技术分享图片
#include<iostream>
#include<iomanip>
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include<algorithm>
#define fir first
#define sec second
#define eps (1e-13)
#define maxn 1010
#define four 4.0
#define six 6.0
#define two 2.0
#define re register 
using namespace std;
typedef pair<double,double>Pdd; 
typedef struct point{double x,y;}P;
typedef struct circle
{
    P o;
    double r;
    bool operator <(const circle &z)const{ return o.x<z.o.x;}
    Pdd getd(double x)  
    {  
        if(r<=fabs(o.x-x))   return Pdd(0,0);//相离 
        double dis=sqrt(r*r-(o.x-x)*(o.x-x));  
        return Pdd(o.y-dis,o.y+dis); //下、上交点 
    }  
}C;
int n,no[maxn],cntt;
double minl=10101010.0,maxr=-10101010.0;
C cir[maxn],circ[maxn];
Pdd pd[maxn];
double D(P _1,P _2){return sqrt((_1.x-_2.x)*(_1.x-_2.x)+(_1.y-_2.y)*(_1.y-_2.y));}

double getf(double x)
{//此线在多个圆中的部分的并 
    double ans=0,last=-10101010;  
    int cnt=0;  
    for(re int i=1;i<=cntt;i++)  
    {  
        pd[++cnt]=circ[i].getd(x);  
        cnt=(pd[cnt]==Pdd(0,0))?cnt-1:cnt;  
    }  
    sort(pd+1,pd+cnt+1);  
    for(re int i=1;i<=cnt;i++)  
    {  
        if(pd[i].fir>last)  
            ans+=pd[i].sec-pd[i].fir,last=pd[i].sec;   
        else if(pd[i].sec>last)  
            ans+=pd[i].sec-last,last=pd[i].sec;  
    }  
    //printf("***%.8lf %.8lf***",x,ans);
    return ans;  
} 
inline double getsim(double l,double r,double fl,double fm,double fr){return ((r-l)/six)*(fl+four*fm+fr);}
double simpson(double l,double m,double r,double fl,double fm,double fr)
{//是对f(x)表示直线x在圆内的部分的总长度这个函数求积分 
    double lm=(l+m)/two,rm=(m+r)/two,flm=getf(lm),frm=getf(rm);
    double lans=getsim(l,m,fl,flm,fm),rans=getsim(m,r,fm,frm,fr),totans=getsim(l,r,fl,fm,fr);
    //printf("%.8lf %.8lf %.8lf %.8lf\n",l,r,lans+rans,totans);system("pause");
    if(fabs(lans+rans-totans)<=eps)return totans;
    else return simpson(l,lm,m,fl,flm,fm)+simpson(m,rm,r,fm,frm,fr);
}
int main()
{
    memset(no,0,sizeof(no));
    scanf("%d",&n);
    for(re int i=1;i<=n;i++)
    {
        scanf("%lf%lf%lf",&cir[i].o.x,&cir[i].o.y,&cir[i].r);
        if(minl>cir[i].o.x-cir[i].r)minl=cir[i].o.x-cir[i].r;
        if(maxr<cir[i].o.x+cir[i].r)maxr=cir[i].o.x+cir[i].r;
    }
    sort(cir+1,cir+n+1);
    for(re int i=1;i<=n;i++)
    {
        if(no[i])continue;    
        for(int j=i+1;j<=n;j++)
        {
            if(no[j])continue;
            if(D(cir[i].o,cir[j].o)+cir[j].r<=cir[i].r)  
                no[j]=1;  
        }
    }//去除被包含的圆 
    for(re int i=1;i<=n;i++)
        if(!no[i])circ[++cntt]=cir[i];
    double ans=simpson(minl,(minl+maxr)/two,maxr,0,getf((minl+maxr)/two),0);
    printf("%.3lf",ans);
    return 0;
}
并不对劲的辛普森积分

 

宣传一波广义呕吐之光,欢迎加入。

技术分享图片

并不对劲的辛普森积分

标签:tin   space   desc   inpu   else   lap   max   stream   多个   

原文地址:https://www.cnblogs.com/xzyf/p/8515693.html

(0)
(0)
   
举报
评论 一句话评论(0
登录后才能评论!
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码!