simpson积分公式

看了Simpson积分公式，惊讶于积分方法在几何问题中竟如此强大！

比如几何中的求面积、体积问题，可以用积分的方法，在某个方向上用扫描线或扫描面切过图形，求被覆盖的长度或面积，然后进行积分。

对于这类问题，如果能直接求出面积、体积，或者能列出覆盖长度（面积）关于扫描线（面）的位置的函数，然后手算积分，那当然是再好不过了。但是，如果直接用解析法、公式法比较费时费力，思考难度大，而题目对精度的要求不高（比如精确到0.01，0.001），数据规模不大，可以考虑用积分的方法。

用程序进行积分，一般的方法是矩形（梯形）切割法，但精度比较差。这里用Simpson积分公式，原理是用二次曲线逼近原函数，精度比矩形切割要好。

引理：在平面直角坐标系中，由任意三点(x1, y1), (x2, y2), (x3, y3)（x1<x2<x3，x₂ = (x₁ + x₃)/2）确定的抛物线y= f(x)在[x1, x3]的定积分为

证明不难，暴力拆开即可。

设f(x) = ax2 + bx + c，于是

左边 = (ax33/3 + bx32/2 + cx3) – (ax13/3 + bx12/2 + cx3)

=1/6 × [ 2a(x33 – x13) + 3b(x32 – x12) + 6c(x3 – x1) ]

=1/6 × (x3 – x1) ×[ 2a(x32 + x3x1 + x12) + 2b(x3 + x1) + 6c ]

=1/6 (x3 – x1)

×[ (ax32 + bx3 + c) + (ax12 + bx1 + c) + a(x3+x1)2 + 2b(x3+x1) + 4c ]

=1/6 (x3 – x1) [ y1 + y3 + a(x3+x1)2 + 2b(x3+x1) + 4c ]

注意到x2 = (x1 + x3)/2，所以

左边 = 1/6 (x3 – x1) [ y1 + y3 + 4ax22 + 4bx2+ 4c ]

= 1/6 (x3 – x1) [ y1 + y3 + 4y2 ] = 右边

那么，如果我们把要积分的区域[L, R]划分成n份（n为偶数），份与份之间的分割线的横坐标为x0~xn，其中x0 = L, xn = R，对应的函数值分别为y0~yn，那么，对第0,1,2条分割线；第2,3,4条分割线；第4,5,6条分割线……之间的区域各自进行二次曲线拟合，即应用上面的公式，把结果累加起来，可以得到下面的Simpson积分公式：

(上面的公式中左边的积分号应该是L到R，之前写错了，抱歉)

有时，如果对划分的份数不太确定，可以使用一种叫做“自适应Simpson”的方法。对要积分的区域[L, R]分成两半[L, mid]和[mid, R]，如果用二次曲线对(L, mid, R)求出的积分值，和对(L, (L+mid)/2, mid)和(mid, (mid+R)/2, R)分别积分在加起来的值，相差不超过某个精确度，那么就可以停止划分。

（要注意误差的叠加效应，当划分份数比较多时，每一份的误差虽然很小，但叠加在一起可能很大，所以精确度尽量设小）

如果我们能确定被积函数（就是覆盖长度或面积关于扫描位置的函数）是个一次多项式或者二次多项式，但是我们不想把函数式求出来，那也没关系，用上面的方法积分，因为用了二次曲线拟合，答案还是准确的。

要注意的一点是，被积函数在我们想积分的一段上应该是有连续的取值的，就是说，我们积分的扫描线范围[L, R]，在[L, R]内不能出现的某条扫描线，它没有穿过任何图形。（这样很有可能L, R, mid都没有穿过图形，求出的答案为0）所以一开始要把有图形的“竖线区”预处理出来。预处理的方法是，每个图形都有一个“被穿过”段，把这些段进行合并就可以了。

Simpson积分方法十分暴力好写。之前某题“多个圆与一个凸多边形的面积并”，用普通的扫描线算法写了250+行，而用积分的方法，只写了130行左右，而且思考难度低，特殊情况少。

附上代码：（此题为SPOJ-CIRU，圆的面积并裸题，自适应Simpson）


#include <cstdio>#include <cstring>#include <cmath>#include <algorithm>#include <utility>usingnamespace std;

#define maxN 1055#define eps (1e-6)

inlinedouble sqr(double x){return x * x;}

int N;struct Tcir{   double Ox, Oy, r;   void read(){scanf("%lf%lf%lf",&Ox,&Oy,&r);}} cir[maxN];

pair<double,double> itv[maxN];

double f(double x){   int tot =0;   for(int i =1; i <= N;++ i)//itv[i] = cir[i].Calc(x);   {       double h = sqr(cir[i].r)- sqr(x - cir[i].Ox);       if(h <0)continue;       h = sqrt( fabs(h));       ++ tot;       itv[tot].first = cir[i].Oy - h;       itv[tot].second = cir[i].Oy + h;   }   sort(itv +1, itv + tot +1);   double L =-(1e5), R =-(1e5), res =0;   for(int i =1; i <= tot;++ i)   {       if(itv[i].first > R)       {           res += R - L;           L = R = itv[i].first;       }       if(itv[i].second > R) R = itv[i].second;   }   res += R - L;   return res;}

double Cnt(double len,double fL,double fmid,double fR){   return len *(fL + fR +4* fmid)/6;}

double Simpson(double L,double mid,double R,double fL,double fmid,double fR,double S){   double mid1 =(L + mid)/2;   double fmid1 = f(mid1);   double S1 = Cnt( mid-L, fL,  fmid1, fmid );   double mid2 =(mid + R)/2;   double fmid2 = f(mid2);   double S2 = Cnt( R-mid, fmid, fmid2, fR );   if( fabs(S1 + S2 - S)< eps )return S1 + S2;       elsereturn Simpson(L, mid1, mid, fL, fmid1, fmid, S1)                  +Simpson(mid, mid2, R, fmid, fmid2, fR, S2);}

double Simpson(double L,double R){   double fL = f(L), fR = f(R), mid =(L + R)/2;   double fmid = f(mid);   double S = Cnt(R-L, fL, fmid, fR);   return Simpson( L, mid, R, fL, fmid, fR, S );}

pair<double,double> evt[maxN];

int main(){   freopen("input.txt","r", stdin);   freopen("output.txt","w", stdout);

   scanf("%d",&N);   for(int i =1; i <= N;++ i)   {       cir[i].read();       evt[i].first = cir[i].Ox - cir[i].r;       evt[i].second = cir[i].Ox + cir[i].r;   }

   sort(evt+1, evt+N+1);

   double L =-(1e5), R =-(1e5), ans =0;   for(int i =1; i <= N;++ i)   {       if(evt[i].first > R)       {           ans += Simpson(L, R);           L = R = evt[i].first;       }       if(evt[i].second > R) R = evt[i].second;   }   ans += Simpson(L, R);

   printf("%.3lf\n", ans);

   return0;}精简的核心代码：template<class T>double simpson(const T &f,double a,double b,int n){    const double h = (b-a)/n;  double ans = f(a) + f(b); for(int i = 1;i < n;i+=2) ans += 4*f(a+i*h);   for(int i = 2;i < n;i += 2) ans += 2*f(a+i*h); return ans*h/3;}

simpson积分公式相关推荐

c语言simpson积分计算方法,数值分析复化Simpson积分公式和复化梯形积分公式计算积分的通用程序...
数值分析复化Simpson积分公式和复化梯形积分公式计算积分的通用程序数值分析第五次程序作业 PB09001057 孙琪 [问题] 分别编写用复化Simpson积分公式和复化梯形积分公式计算积分的通 ...
数值分析18 - 通过直接方法得到函数积分近似数值方法（左、右、中、梯形矩形积分公式、Simpson积分公式）
左矩形积分公式右矩形积分公式中矩形积分公式 (1阶精度) 梯形矩形积分公式 (1阶精度) Simpson积分公式 (3阶精度)
python怎么计算曲面的表面积_利用simpson积分公式计算曲面表面积
利用 simpson 积分公式计算曲面表面积夏军剑 ; 张新巍 ; 李维伟 [期刊名称] <科技资讯> [年 ( 卷 ), 期] 2014(012)008 [摘要] 二重积分的数值算法比 ...
【matlab】数值积分公式的程序实现
(一)专题实验(Newton-Cotes积分公式) 1.编写[a,b]上梯形积分公式.Simpson积分公式. 2.利用自己编写的程序计算定积分,计算一下数值解和精确解之间差的绝对值. 梯形积分: f ...
2018年东北农业大学春季校赛题解
[题目链接] 写在前面:从都到尾做了一下这场比赛,似乎好题都是原题,水题都是他们学校自己出的.原题在抄过来的过程中,很多题目的题面.数据范围都出了问题,还有题目数据很水.建议以后这样的比赛不要挂到外面 ...
洛谷P3779 [SDOI2017]龙与地下城（概率论+Simpson+FFT）
题面传送门题解 orz shadowice 正态分布正态分布是随机变量\(X\)的一种概率分布形式.它用一个期望\(\mu\)和方差\(\sigma^2\)就可以描述,记为\(N(\mu,\si ...
CSU 1806 Toll 自适应simpson积分+最短路
分析:根据这个题学了一发自适应simpson积分(原来积分还可以这么求),然后就是套模板了学习自适应simpson积分:http://blog.csdn.net/greatwall1995/arti ...
2021-01-07 matlab数值分析数值积分与数值微分复合梯形公式复合Simpson公式
matlab数值分析数值积分与数值微分 1 复合梯形公式 function I=ftrapz(f,a,b,n) format long %显示15位双精度 h=(b-a)/n; x=linspace ...
HDU 1724 Ellipse ——Simpson积分
[题目分析] 一看题目,直接把椭圆积分起来就可以了嘛. 然后发现椭圆比较难积分,还是算了吧. 用Simpson积分硬上. 大概就是用二次函数去拟合面积. [代码] #include <cstdi ...
C语言实现数值积分之Simpson 1/3法则(附完整源码)
数值积分之Simpson 1/3法则数值积分之Simpson 1/3法则的完整源码(定义,实现,main函数测试) 数值积分之Simpson 1/3法则的完整源码(定义,实现,main函数测试) # ...

simpson积分公式

simpson积分公式相关推荐

最新文章

热门文章