首页 > 其他分享 >自适应辛普森法从入门到进门

自适应辛普森法从入门到进门

时间：2024-02-16 14:55:05浏览次数：35  
标签：入门 int double mid 辛普森 cir 法从 inline simpson

前言

学数学学的。

simpson

背景

  • 我们要计算这样一个式子：

\[\int_l^r f(x)\text dx \]

显然计算机是很难把柿子推出来的。

函数的拟合

  • 对于一个奇怪的函数，为了对其求导，我们可以用一个图像近似且容易求导的函数来替代，这个过程叫做拟合。

这里我们用二次函数来替代，那么有：

\[\begin{aligned}\int_l^r f(x)\text dx&\approx\int_l^r(ax^2+bx+c)\text dx\\&=\dfrac{a}{3}(r^3-l^3)+\dfrac{b}{2}(r^2-l^2)+c(r-l)\\&=\dfrac{r-l}{6}\left[2a(r^2+rl+l^2)+3b(r+l)+6c\right]\\&=\dfrac{r-l}{6}(2ar^2+2arl+2al^2+3br+3bl+6c)\\&=\dfrac{r-l}{6}\left[(ar^2+br+c)+(al^2+bl+c)+4a\left(\dfrac{r+l}{2}\right)+4b\left(\dfrac{r+l}{2}\right)+4c\right]\\&=\dfrac{r-l}{6}\left[f(r)+f(l)+4f\left(\dfrac{r+l}{2}\right)\right]\end{aligned} \]

code

inline double simpson(double l,double r){
    return (r-l)*(f(r)+f(l)+4.0*f((l+r)/2))/6.0;
}
  • 显然这样的话精度误差会很大，因此我们需要将函数分段拟合。

为了确定合适的分段，我们需要一个“自适应”的操作。

如果这一段函数与拟合的二次函数相似，那么直接把这一段当做这个二次函数，然后套用公式计算；否则应该把这一段分成两半，并递归进行这一过程。

关于函数相似的判断，我们可以将区间的左半和右半部分分别进行计算，与计算的相差小于设定的 \(\varepsilon\) 即可返回值。

在分治判断的时候，除了判断精度是否正确，一般还要强制执行最少的迭代次数。

code

inline double work(double l,double r,double eps,double ans,int depth){
    double mid=(l+r)/2,x=simpson(l,mid),y=simpson(mid,r);
    if(fabs(x+y-ans)<=15*eps and depth<0)
        return (x+y)+(x+y-ans)/15.0;
    return work(l,mid,eps/2,x,depth-1)+work(mid,r,eps/2,y,depth-1);
}

一些题目

P4525 自适应辛普森法 1

传送门

Description

求：

\[\int_l^r \dfrac{cx+d}{ax+b}\text dx \]

Solution

套板子。

code

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

inline double simpson(double l,double r){
    return (r-l)*(f(r)+f(l)+4.0*f((l+r)/2))/6.0;
}

inline double work(double l,double r,double eps,double ans,int depth){
    double mid=(l+r)/2,x=simpson(l,mid),y=simpson(mid,r);
    if(fabs(x+y-ans)<=15*eps and depth<0)
        return (x+y)+(x+y-ans)/15.0;
    return work(l,mid,eps/2,x,depth-1)+work(mid,r,eps/2,y,depth-1);
}

signed main(){
    scanf("%lf%lf%lf%lf%lf%lf",&a,&b,&c,&d,&l,&r);
    printf("%.6f",work(l,r,1e-7,simpson(l,r),12));
    return 0;
}

P4526 自适应辛普森法 2

传送门

Description

求：

\[\int_0^{+\infty} x^{\frac{a}{x}-x}\text dx \]

Solution

注意到当 \(a<0\) 时积分发散，\(a>0\) 时积分收敛。

特判后选一个适合的区间计算即可。

code

inline double f(double x){
    return pow(x,a/x-x);
}

inline double simpson(double l,double r){
    return (r-l)*(f(r)+f(l)+4.0*f((l+r)/2))/6.0;
}

inline double work(double l,double r,double eps,double ans,int depth){
    double mid=(l+r)/2,x=simpson(l,mid),y=simpson(mid,r);
    if(fabs(x+y-ans)<=15*eps and depth<0)
        return (x+y)+(x+y-ans)/15.0;
    return work(l,mid,eps/2,x,depth-1)+work(mid,r,eps/2,y,depth-1);
}

signed main(){
    scanf("%lf",&a);
    if(a<0)
        puts("orz");
    else{
        double l=1e-7,r=23.3;
        printf("%.5lf",work(l,r,1e-7,simpson(l,r),12));
    }
    return 0;
}

HDU-1724 Ellipse

传送门

Description

给定 \(a,b,l,r\)，求椭圆：\(\dfrac{x^2}{a^2}+\dfrac{y^2}{b^2}=1\) 与直线：\(x=l,x=r\) 所组成的封闭图形的面积。

Solution

还是套板子。

code

inline double f(double x){
    return b*(sqrt(1.0-(x*x)/(a*a)));
}

inline double simpson(double l,double r){
    return (r-l)*(f(r)+f(l)+4.0*f((l+r)/2))/6.0;
}

inline double work(double l,double r,double eps,double ans,int depth){
    double mid=(l+r)/2,x=simpson(l,mid),y=simpson(mid,r);
    if(fabs(x+y-ans)<=15*eps and depth<0)
        return (x+y)+(x+y-ans)/15.0;
    return work(l,mid,eps/2,x,depth-1)+work(mid,r,eps/2,y,depth-1);
}

signed main(){
    int T=read();
    while(T--){
        scanf("%lf%lf%lf%lf",&a,&b,&l,&r);
        printf("%.3f",2.0*work(l,r,1e-7,simpson(l,r),12));puts("");
    }
    return 0;
}

P4207 [NOI2005] 月下柠檬树

传送门

Description

给定一棵由圆台、圆锥组成的树和平行光线的角度，求投影面积。

Solution

首先注意到若点 \(p\) 的坐标为 \((0,x)\)，那么投影后的坐标为 \(p'(x\cot \alpha,0)\)。

圆的平行投影还是圆，圆的半径投影后不会改变，只有高度与圆台投影形成的梯形尺寸会改变。这一部分可以画图用相似解决。

把圆和梯形的数据确定后，直接套板子即可。

code

inline double sq(double x){return x*x;}

inline int sgn(double x){
    if(fabs(x)<1e-8)
        return 0;
    return x<0?-1:1;
}

struct circle{
    double x,y,r;
    circle(){}
    circle(double _x,double _y,double _r):x(_x),y(_y),r(_r){}
}cir[N];

struct Trape{
    double l,r,al,ar;
    Trape(){}
    Trape(double _l,double _r,double _al,double _ar):l(_l),r(_r),al(_al),ar(_ar){}
}trape[N];

inline double f(double x){
    double res=0;
    for(int i=1;i<=n;i++)
        if(sgn(cir[i].r-fabs(x-cir[i].x))>0){
            double tmp=sqrt(sq(cir[i].r)-sq(fabs(x-cir[i].x)));
            res=max(res,tmp*2);
        }
    for(int i=1;i<=n;i++)
        if(sgn(x-trape[i].l)>0 and sgn(trape[i].r-x)>0){
            double d=x-trape[i].l,len=trape[i].r-trape[i].l;
            double Len=len*trape[i].al/(trape[i].al-trape[i].ar);
            double tmp=(Len-d)/Len*trape[i].al;
            res=max(res,tmp*2);
        }
    return res;
}

inline double simpson(double l,double r){
    return (r-l)*(f(r)+f(l)+4.0*f((l+r)/2))/6.0;
}

inline double work(double l,double r,double eps,double ans,int depth){
    double mid=(l+r)/2,x=simpson(l,mid),y=simpson(mid,r);
    if(fabs(x+y-ans)<=15*eps and depth<0)
        return (x+y)+(x+y-ans)/15.0;
    return work(l,mid,eps/2,x,depth-1)+work(mid,r,eps/2,y,depth-1);
}

signed main(){
    n=read();scanf("%lf",&alpha);
    for(int i=0;i<=n;i++)
        scanf("%lf",&h[i]);
    for(int i=1;i<=n;i++)
        scanf("%lf",&r[i]);
    Tan=1.0/tan(alpha);
    for(int i=1;i<=n;i++){
        cir[i]=circle(Tan*H,0,r[i]);
        double len=Tan*h[i],ll=(r[i+1]-r[i])*r[i]/len,rr=(r[i+1]-r[i])*r[i+1]/len;
        trape[i]=Trape(Tan*H-ll,Tan*(H+h[i])-rr,sqrt(sq(r[i])-sq(ll)),sqrt(sq(r[i+1])-sq(rr)));
        H+=h[i];
        L=min(L,cir[i].x-cir[i].r);
        R=max(R,cir[i].x+cir[i].r);
    }
    R=max(R,Tan*H);
    printf("%.2lf",work(L,R,5e-4,simpson(L,R),12));
    return 0;
}

BZOJ2178 圆的面积并

传送门

Description

  • 给出 \(n\) 个圆，求其面积并。

Solution

套板子，但有几个优化：

  1. 将被完全包含的圆删除。

  2. 将圆按 \(x\) 排序，每个区间分别算积分。

code

struct circle{
    double x,y,r;
}cir[N];

inline bool cmp(circle a,circle b){
    return a.x-a.r<b.x-b.r;
}

struct Segment{
    double l,r;
    bool operator <(const Segment yy) const {
        return l<yy.l;
    }
}tmp[1005];
int ts;
map<double,double>F;

inline double f(double x){
    if(F.count(x))
        return F[x];
    ts=0;
    for(int i=1;i<=n;i++){
        if(cir[i].r<=fabs(cir[i].x-x))
            continue;
        double ff=cir[i].r*cir[i].r-(cir[i].x-x)*(cir[i].x-x);
        if(ff<=1e-4)
            continue;
        ff=sqrt(ff);
        tmp[++ts]=(Segment){cir[i].y-ff,cir[i].y+ff};
    }
    sort(tmp+1,tmp+ts+1);
    double nows=0,nowr=-1e9;
    for(int i=1;i<=ts;i++){
        if(tmp[i].r<nowr)
            continue;
        if(tmp[i].l<=nowr){
            nows+=tmp[i].r-nowr;
            nowr=tmp[i].r;
        }else{
            nows+=tmp[i].r-tmp[i].l;
            nowr=tmp[i].r;
        }
    }
    return F[x]=nows;
}

inline double simpson(double l,double r){
    return (r-l)*(f(r)+f(l)+4.0*f((l+r)/2))/6.0;
}

inline double work(double l,double r,double eps,double ans,int depth){
    double mid=(l+r)/2,x=simpson(l,mid),y=simpson(mid,r);
    if(fabs(x+y-ans)<=15*eps and depth<0)
        return (x+y)+(x+y-ans)/15.0;
    return work(l,mid,eps/2,x,depth-1)+work(mid,r,eps/2,y,depth-1);
}

signed main(){
    n=read();
    for(int i=1;i<=n;i++)
        cin>>cir[i].x>>cir[i].y>>cir[i].r;
    sort(cir+1,cir+n+1,cmp);
    double nowr=-inf,ans=0;
    for(int i=1;i<=n+1;i++){
        if(cir[i].x-cir[i].r>nowr){
            ans+=work(cir[i].x-cir[i].r,cir[i].x+cir[i].r,1e-5,simpson(cir[i].x-cir[i].r,cir[i].x+cir[i].r),12);
            nowr=cir[i].x+cir[i].r;
        }else if(cir[i].x+cir[i].r>nowr){
            ans+=work(nowr,cir[i].x+cir[i].r,1e-4,simpson(nowr,cir[i].x+cir[i].r),12);
            nowr=cir[i].r+cir[i].x;
        }
    }
    printf("%.3lf\n",ans);
    return 0;
}

其他面积并的东西都是同理的，懒得粘了。

后记

乱搞的东西你为什么要学？

标签：入门,int,double,mid,辛普森,cir,法从,inline,simpson
From： https://www.cnblogs.com/QcpyWcpyQ/p/18017149

相关文章

  • 数论从入门到进门
    gcd朴素欧几里得（辗转相除法）这一节我们默认：\(a,b\in\mathbb{Z}\)对于求解\(\gcd(a,b)\)需要用朴素欧几里得定理。\[\gcd(a,b)=\gcd(b,a\bmodb)\]gcd\(\gcd\)为\(\text{greatestcommondivisor}\)的缩写，即最大公约数。（或称最大公因数）对于\(\gcd\)函数，有以下......
  • 多项式从入门到进门
    多项式全家福多项式一个以\(x\)为变量的多项式定义在一个代数域\(F\)上，将函数\(A(x)\)表示为形式和：\[A(x)=\sum\limits_{i=0}^{n-1}a_ix^i\]我们称\(a_0,a_1,\cdots,a_n\)为多项式的系数，所有系数都属于数域\(\mathbbF\)。如果一个多项式的最高次的非零系数是\(k......
  • 生成函数从入门到进门
    引入先看下面这个例子：\[(1+a_1x)\times(1+a_2x)\times\cdots\times(1+a_nx)\]拆括号得：\[1+(a_1+a_2+\cdots+a_n)x+(a_1a_2+a_1a_3+\cdots+a_{n-1}a_n)x^2+\cdots+(a_1a_2\cdotsa_n)x^n\]其中\(x^2\)的系数包含了从\(n\)个元素\(\{a_1,a_2,\cdots,a_n\}\)中选取两个的......
  • 概率期望从入门到进门
    概率的线性性定义：\(\mathbb{E}(X)=\sum_i\timesP(X=i)\)。其中\(x\)为变量。线性性\[\begin{aligned}\mathbb{E}(aX+bY)&=\sumi\timesP(aX+bY=i)\\&=\sumi\sum_j\sum_k[j+k=i]P(aX=j)P(bY=k)\\&=\sum_j\sum_k(j+k)P(aX=j)\cdotP(bY=k)\\&......
  • HydroOJ 从入门到入土（13）批量修改题号前缀
    题库的管理，无论是用前缀来分组，还是用域来分组，都有不好管理的地方，尤其是题号。有的时候导入了一堆题，导入完发现题号不是自己想要的，但删起来很麻烦，一个一个改更不现实，真是欲哭无泪。本文主要记录了一次批量修改题号前缀的过程，仅供参考。修改中涉及数据库操作，修改前一定要现在虚......
  • 拓扑排序入门
    目录写在前面一些概念算法步骤字典序最大/最小的拓扑序列？模板例题3704.排队家谱树奖金P1983[NOIP2013普及组]车站分级1639.拓扑顺序写在前面昨晚cfdiv3的F就是一道基本上可以说板子的拓扑排序的题目，没有做出来感觉图论很早之前就看了，但是基本没有刷过什么题，开始补一下图论......
  • 掌握C语言文件操作：从入门到精通的完整指南！
    ✨✨欢迎大家来到贝蒂大讲堂✨✨......
  • pytorch深度学习入门（8）之-Torchaudio使用Tacotron2 文本转语音
    https://blog.csdn.net/ajunbin859/article/details/134380417?ops_request_misc=&request_id=&biz_id=102&utm_term=pytorch%E7%89%88%E6%9C%AC%E7%9A%84tacotron%E8%AF%A6%E7%BB%86%E5%AE%89%E8%A3%85%E6%95%99%E7%A8%8B&utm_medium=distribute.pc_search_r......
  • day04_操作系统入门
    今日笔记学操作系统基础概念linux系统linux系统(centos)+vmware安装起来（网络配置，磁盘分区）ubuntu安装xshell服务器的远程连接服务器网站的前后端，数据库app的前后端，数据库微信、腾讯微信的服务器移动端设备上，安装的微信客户端在线笔记笔记对运维来说，就是一个宝藏，mar......
  • day21_乌班图入门
    .请解释yum缓存，如何理解、如何管理去网络源下载软件rpm包，会涉及网络延时，网络资源消耗1.解决，关于yum缓存包的理解（自己搭建yum仓库）11.当你拿到一个初始化的机器，默认安装的软件（centos上的rpm格式的软件）数量可能很少导致你后期使用各种工具，会报错，比如python调用gzip解压缩功能s......

赞助商

阅读排行

网站主页关于我们联系我们网站地图

本网站内容转载自其他媒体，侵权联系[admin##ips99.com]。

Copyright © 2020-2023 IPS99 版权所有 IPS99