自适应Simpson积分

最新推荐文章于 2022-10-16 16:44:23 发布

weixin_30437481

最新推荐文章于 2022-10-16 16:44:23 发布

阅读量66

点赞数

CC 4.0 BY-SA版权

原文链接：http://www.cnblogs.com/TSHugh/p/8543710.html

本文介绍了Simpson公式及其在数值积分中的应用，包括如何使用该公式计算不规则函数的定积分以及求解特定问题的实例代码。文中还讨论了Simpson公式的修正方法，并给出了自适应Simpson积分的实现方式。

Simpson公式的话是一个对于三次及以下函数成立的连续区间求定积分的公式.(好像还有许多其他的公式……)
虽然他是一个普通的定积分公式,但是他不仅可以用来求不规则函数的定积分还可以用来瞎搞求面积,你可以把用这个式子求面积看成是对于定积分本质的运用,也可以直接把图像看成函数,我习惯于按照后者来理解.
http://blog.youkuaiyun.com/greatwall1995/article/details/8639135
上面那个博客已经把Simpson积分说的已经比较明白了.
http://blog.youkuaiyun.com/xl2015190026/article/details/53518077
上面这个博客里提到了对于Simpson积分的修正——I.乘15 II.把精度劈半
但是据我的经验,不修正才稳啊.
对于自适应Simpson积分的代码实现,在递归的过程中,除了左右端点以外,其他的参量可以根据个人喜好与具体情况选择性地传入.
hdu1724:Ellipse *真尼玛裸

#include <cmath>
#include <cstdio>
#include <cstring>
#include <algorithm>
typedef double db;
const db eps=1e-4;
db K,B;
inline db f(db x){
  return std::sqrt(K*x*x+B);
}
inline db calc(db l,db r){
  return (f(l)+4*f((l+r)*0.5)+f(r))*(r-l)/6;
}
inline db simpson(db l,db r){
  db mid=(l+r)*0.5;
  if(std::fabs(calc(l,mid)+calc(mid,r)-calc(l,r))<eps)return calc(l,r);
  return simpson(l,mid)+simpson(mid,r);
}
int main(){
  int T,a,b,l,r;scanf("%d",&T);
  while(T--){
    scanf("%d%d%d%d",&a,&b,&l,&r);
    K=-(db)(b*b)/(a*a),B=b*b;
    printf("%.3f\n",simpson(l,r)*2);
  }
}

这是一种简单清真但是慢的打法

bzoj2178:圆的面积并 *也真尼玛裸(如果你不怕被hack的话)

#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
#define ft first
#define sd second
#define mmp(a,b) (std::make_pair(a,b))
typedef double db;
typedef std::pair<db,db> pdd;
const int N=1010;
const db eps=1e-13;
pdd seg[N];
struct Cir{
  int x,y,r;
  inline db dis(Cir a){
    return std::sqrt((a.x-x)*(a.x-x)+(a.y-y)*(a.y-y));
  }
  inline void read(){scanf("%d%d%d",&x,&y,&r);}
  inline bool in(db pos){return x-r<pos&&x+r>pos;}
  inline bool in(Cir a){
    return r>=dis(a)+a.r;
  }
  inline pdd get(db pos){
    db len=std::sqrt(r*r-(pos-x)*(pos-x));
    return mmp(y-len,y+len);
  }
}cir[N];
bool die[N];
int n;
inline db f(db pos){
  int i,len=0;db ret=0,l,r;
  for(i=1;i<=n;++i)
    if(cir[i].in(pos))
      seg[++len]=cir[i].get(pos);
  std::sort(seg+1,seg+(len+1));
  for(i=1;i<=len;++i){
    l=seg[i].ft,r=seg[i].sd;
    if(i!=1&&seg[i-1].sd>l)l=seg[i-1].sd;
    if(l<r)ret+=r-l;
    else seg[i].sd=l;
  }
  return ret;
}
inline db calc(db l,db r,db fl,db fr,db fm){
  return (r-l)/6*(fl+4*fm+fr);
}
inline db Simpson(db l,db r,db mid,db fl,db fr,db fm,db A){
  db m1=(l+mid)/2,m2=(mid+r)/2;
  db fm1=f(m1),fm2=f(m2);
  db L=calc(l,mid,fl,fm,fm1),R=calc(mid,r,fm,fr,fm2);
  if(std::fabs(L+R-A)<=eps)return A;
  return Simpson(l,mid,m1,fl,fm,fm1,L)+Simpson(mid,r,m2,fm,fr,fm2,R);
}
int main(){
  scanf("%d",&n);
  int i,j,l=2000,r=-2000;
  for(i=1;i<=n;++i){
    cir[i].read();
    l=std::min(l,cir[i].x-cir[i].r);
    r=std::max(r,cir[i].x+cir[i].r);
  }
  for(i=1;i<=n;++i)
    if(!die[i])
      for(j=1;j<=n;++j)
        if(j!=i&&!die[j])
          if(cir[i].in(cir[j]))
            die[j]=true;
  j=0;
  for(i=1;i<=n;++i)
    if(!die[i])
      cir[++j]=cir[i];
  n=j;
  db fl=f(l),fr=f(r),fm=f((l+r)/2.);
  printf("%.3f\n",Simpson(l,r,(l+r)/2.,fl,fr,fm,calc(l,r,fl,fr,fm)));
  return 0;
}

这是一种看起来傻逼但是比较快的打法

bzoj1502:[NOI2005]月下柠檬树 ***
首先,你需要高超的空间想象能力(一层一层、一片一片地去想)和大胆的猜想能力(猜是切线).
其次,你还需要一定的高中几何基础(三角函数).
然后,去码吧(板子).

#include <cmath>
#include <cstdio>
#include <cstring>
#include <algorithm>
typedef double db;
const db eps=1e-5;
const db Pi=std::acos(-1.);
const int N=510;
int n,m;
db alpha;
struct Cir{
  db x,r,l;
  bool in(Cir a){return r>a.r+std::fabs(a.x-x);}
  bool in(db pos){return x-r<pos&&pos<x+r;}
  db h(db pos){return std::sqrt(r*r-(x-pos)*(x-pos));}
}cir[N];
inline bool comp1(Cir a,Cir b){return a.l<b.l;}
struct Line{
  db l,r,k,b;
  bool in(db pos){return l<pos&&pos<r;}
  db h(db pos){return k*pos+b;}
}line[N];
inline bool comp2(Line a,Line b){return a.l<b.l;}
inline Line get(db aks1,db w1,db aks2,db w2){
  Line ret;
  ret.l=aks1,ret.r=aks2;
  ret.k=(w1-w2)/(aks1-aks2),ret.b=w1-aks1*ret.k;
  return ret;
}
inline db f(db pos){
  db ret=0.;int i;
  for(i=1;i<=n;++i){
    if(cir[i].in(pos))
      ret=std::max(ret,cir[i].h(pos));
    if(pos<cir[i].l)break;;
  }
  for(i=1;i<=m;++i){
    if(line[i].in(pos))
      ret=std::max(ret,line[i].h(pos));
    if(pos<line[i].l)break;
  }
  return ret*2.;
}
inline db calc(db l,db r,db fl,db fr,db fm){
  return (r-l)*(fl+4*fm+fr)/6;
}
inline db Simpson(db l,db r,db mid,db fl,db fr,db fm,db A){
  db m1=(l+mid)*0.5,m2=(mid+r)*0.5;
  db fm1=f(m1),fm2=f(m2);
  db L=calc(l,mid,fl,fm,fm1),R=calc(mid,r,fm,fr,fm2);
  if(std::fabs(L+R-A)<=eps)return L+R;
  return Simpson(l,mid,m1,fl,fm,fm1,L)+Simpson(mid,r,m2,fm,fr,fm2,R);
}
int main(){
  scanf("%d%lf",&n,&alpha);
  ++n,alpha=std::tan(alpha);
  int i;db h,s=0.,l=1e18,r=-1e18;
  for(i=1;i<=n;++i)
    scanf("%lf",&h),cir[i].x=(s+=h)/alpha;
  for(i=1;i<=n;++i){
    if(i==n)h=0.;
    else scanf("%lf",&h);
    l=std::min(l,cir[i].x-h);
    r=std::max(r,cir[i].x+h);
    cir[i].r=h;
    cir[i].l=(cir[i].x-h);
  }
  for(i=1;i<n;++i){
    if(cir[i].in(cir[i+1])||cir[i+1].in(cir[i]))continue;
    if(cir[i].r==cir[i+1].r)line[++m]=get(cir[i].x,cir[i+1].x,cir[i].r,cir[i].r);
    else{
      h=std::fabs(cir[i].x-cir[i+1].x);
      h=h+h*(std::min(cir[i].r,cir[i+1].r))/(std::fabs(cir[i].r-cir[i+1].r));
      if(cir[i].r>cir[i+1].r)s=std::acos(cir[i].r/h);
      else s=std::acos(-cir[i+1].r/h);
      line[++m]=get(cir[i].x+cir[i].r*std::cos(s),cir[i].r*std::sin(s),cir[i+1].x+cir[i+1].r*std::cos(s),cir[i+1].r*std::sin(s));
    }
  }
  std::sort(cir+1,cir+(n+1),comp1);
  std::sort(line+1,line+(m+1),comp2);
  db mid=(l+r)*0.5,fl=f(l),fr=f(r),fm=f(m);
  printf("%.2f\n",Simpson(l,r,mid,fl,fr,fm,calc(l,r,fl,fr,fm)));
  return 0;
}