数值积分算法与MATLAB实现陈悦5133201讲解

东北大学秦皇岛分校数值计算课程设计报告数值积分算法及MATLAB实现学院数学与统计学院专业信息与计算科学学号5133201姓名陈悦指导教师姜玉山张建波成绩教师评语：指导教师签字：2015年07月14日1 绪论数值分析是计算数学的一个主要部分，计算数学是数学科学的一个分支，它研究用计算机求解各种数学问题的数值检索方其理论与软件的实现.而数值分析主要研究数值计算.现科学技术的发展与进步提出了越来越多的复杂的数值计算问题，这些问题的圆满解决已远人工手算所能胜任，必须依靠电子计算机快速准确的数据处理能力.这种用计算机处理数值问题的方法，成为科学计算.今天，科学计算的应用范围非常广泛，天气预报、工程设计、流体计算、经济规划和预测以及国防尖端的一些科研项目，如核武器的研制、导弹和火箭的发射等，始终是科学计算最为活跃的领域.1.1 数值积分介绍数值积分是数值分析的重要环节，实际问题当中常常需要计算积分，有些数值方法，如微分方程和积分方程的求解，也都和积分计算相联系.求某函数的定积分时，在多数情况下，被积函数的原函数很难用初等函数表达出来，因此能够借助微积分学的牛顿-莱布尼兹公式计算定积分的机会是不多的.另外，许多实际问题中的被积函数往往是列表函数或其他形式的非连续函数，对这类函数的定积分，也不能用不定积分方法求解.由于以上原因，数值积分的理论与方法一直是计算数学研究的基本课题.对微积分学做出杰出贡献的数学大师，如I.牛顿、L.欧拉、C.F.高斯、拉格朗日等人都在数值积分这个领域作出了各自的贡献，并奠定了这个分支的理论基础.构造数值积分公式最通常的方法是用积分区间上的n 次插值多项式代替被积函数，由此导出的求积公式称为插值型求积公式.特别在节点分布等距的情形称为牛顿-科特斯公式，例如梯形公式(Trapezoidal Approximations)与抛物线公式(ApproximationsUsing Parabolas)就是最基本的近似公式.但它们的精度较差.龙贝格算法是在区间逐次分半过程中，对梯形公式的近似值进行加权平均获得准确程度较高的积分近似值的一种方法，它具有公式简练、计算结果准确、使用方便、稳定性好等优点，因此在等距情形宜采用龙贝格求积公式(Rhomberg Integration).当用不等距节点进行计算时，常用高斯型求积公式计算，它在节点数目相同情况下，准确程度较高，稳定性好，而且还可以计算无穷积分.数值积分还是微分方程数值解法的重要依据.许多重要公式都可以用数值积分方程导出.现探讨数值积分算法以及运用MATLAB软件的具体实现1.2 MATLAB软件MATLAB是美国MathWorks公司出品的商业数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境，主要包括MATLAB 和Simulink两大部分.MATLAB是matrix&laboratory两个词的组合，意为矩阵工厂（矩阵实验室）.是由美国mathworks公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技计算环境.它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言（如C、Fortran）的编辑模式，代表了当今国际科学计算软件的先进水平.MATLAB和Mathematica、Maple并称为三大数学软件.它在数学类科技应用软件中在数值计算方面首屈一指.MATLAB可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连接其他编程语言的程序等，主要应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、金融建模设计与分析等领域.MATLAB 的基本数据单位是矩阵，它的指令表达式与数学、工程中常用的形式十分相似，故用MATLAB 来解算问题要比用C ，FORTRAN 等语言完成相同的事情简捷得多，并且MATLAB 也吸收了像Maple 等软件的优点，使MATLAB 成为一个强大的数学软件.在新的版本中也加入了对C ，FORTRAN ，C++，JAVA 的支持.2 数值积分的基本概念一般的，我们可以在区间[],a b 上适当选取某些节点k x ，然后用()k f x 的加权平均得到平均高度()f ς的近似值，这样构造出的求积公式具有下列形式：()()0d nakkbk f x x A f x =≈∑⎰，式中kx称为求积节点；k A 称为求积系数，亦称伴随节点k x 的权.权k A 仅仅与节点k x 的选取有关，而不依赖于被积函数()f x 的具体形式. 2.1 代数精度的概念如果某个求积公式对于次数不超过m 的多项式均能准确的成立，但对于1m +次多项式就不准确成立，则称该求积公式具有m 次代数精度（或代数精确值）一般地，欲使求积公式具有m 次代数精度，只要令它对于()1,,,m f x x x =L 都能准确成立，这就要求:()()22+1+1=-,1=-,21=-.+1k k k m m m k k A b a A x b a A x b a m ⎧⎪⎪⎪⎨⎪⎪⎪⎩∑∑∑M 2.2 求积公式的余项令求积公式的余项为[]R f ，其中[]()()b a0=d nk k k R f f x x A f x =-∑⎰.；区间[],a b 可以是有限的或无限的.构造求积公式的问题就是确定j x 和j A 使得[]R f 在某种意义下尽可能地小.3 数值积分方法及MATLAB 实现3.1 复合辛普森公式 3.1.1 插值型求积公式用插值多项式()n L x 替换积分()*d b aI f x x =⎰中的被积函数()f x ，然后计算()d bn aI L x x =⎰作为积分的近似值，这样建立的求积公式称为插值型求积公式.用插值多项式()n L x 的表达式()()0()n n k k k L x l x f x ==∑，代入得()0nk k k I A f x ==∑，其中:011011()()()()d ()()()()bk k n k ak k k k k k n x x x x x x x x A x x x x x x x x x -+-+----=----⎰L L L L .其余项为:[]()(1)*()d (1)!n baf R f I I x x n ξω+=-=+⎰.3.1.2 牛顿-科特斯公式介绍取等距节点，把积分区间[],a b 剖分成n 等分.令步长()b a h n-=，并记0,n x a x b ==，则1n +个节点为0,0,1,k x x kh k n =+=L ，代入得:k A =0(1)(1)(1)(1)()d ()!!n k nh t t t k t k t n t n k k ----+----⎰L L . 这种等距节点的插值型求积公式通常称为牛顿-科特斯公式. 3.1.3 辛普森公式利用牛顿-科特斯公式，取n =2，此时为()4()()62b a a b S f a f f b -+⎡⎤=++⎢⎥⎣⎦，即为辛普森公式，其余项为[]5(4)1()()902s b a R f f η-=-.3.1.4 复合辛普森公式将积分区间[],a b 分成n 等分，分点为(0,1,,)k x a kh k n =+=L ，其中()b a h n-=，记区间[]1,k k x x +的中点为12k x +，在每个小区间[]1,k k x x +上用辛普森公式，则得到所谓的复合辛普森公式：()()1110246n n k k k k h S f x fx f x -++=⎡⎤⎛⎫=++⎢⎥ ⎪⎢⎥⎝⎭⎣⎦∑.()()()111102246n n n k k k k h S f a f b f x f x --+==⎡⎤⎛⎫=+++⎢⎥ ⎪⎢⎥⎝⎭⎣⎦∑∑. 余项为[]()4(4)()1802n b a h R S f ξ-=-，(,)a b ξ∈. 3.1.5 复合辛普森公式的MATLAB 实现 代码如下：function s=xinpusen(fun,a,b,n) h=(b-a)/n s1=0;s2=0; for k=0:(n-1) x=a+h*k;s1=s1+feval(fun,x); endfor k=0:(n-1) x=a+h*(k+1/2); s2=s2+feval(fun,x); ends=h/6*(feval(fun,a)+feval(fun,b)+2*s1+4*s2); 3.2 龙贝格公式 3.2.1 梯形法的递推化将区间[],a b 分成n 等份，共有1n +个分点，如果将求积区间再二分一次，则分点增至21n +个，用复合梯形公式求得该子区间上的积分值为()1122()4k k k h f x fx f x ++⎡⎤⎛⎫++⎢⎥ ⎪⎢⎥⎝⎭⎣⎦，将每个子区间上的积分值相加得2n T =()1111002()()42n n k k k k k h h f x f x f x --++==++⎡⎤⎣⎦∑∑. 得到递推公式：12102122n n n k k h T T fx -+=⎛⎫=+ ⎪⎝⎭∑. 3.2.2 龙贝格算法公式当()f x 在[],a b 上充分光滑时， 可证用1()T h 逼近I 的截断误差是：1()I T h -=21a h 42a h +2k k a h +++L L按理查森外推法：()1()F h F h =，1()()(),(1,2,)1mpm m m m p F qh q F h F h m q+-==-L 其中，q 为满足10(1,2,)m p q m -≠=L 的适当正数 .取序列：()()()142,1,2,41m m n m m h T T h T h m +⎛⎫- ⎪⎝⎭==-L . 用()1m T h +来逼近I 的误差为()()21m O h +，这种算法就是龙贝格算法.3.2.3 龙贝格算法MATLAB 实现 代码如下：function s=longbeige(fun,a,b,tol) if nargin<4，tol=1e-4;end i=1;j=1;h=b-a;T(1,1)=h*(feval(fun,a)+feval(fun,b))/2;T(i+1,j)=T(I,j)/2+sum(feval(fun,a+h/2:h:b-h/2))*h/2; T(i+1,j+1)=(4^j*T(i+1,j)-T(i.j))/(4^j-1); while (abs(T(i+1,i+1)-T(I,i))>tol) i=i+1;h=h/2;T(i+1,1)=T(I,1)/2+sum(feval(fun,a+h/2:h:b-h/2))*h/2; for j=1:iT(i+1,j+1)=(4^j*T(i+1,j)-T(I,j))/(4^j-1); end end Ts=T(i+1,j+1); 3.3 自适应法自适应积分法是一种比较经济而且快速的求积分的方法.他能自动地在被积函数变化剧烈的区域增多节点，而在被积函数变化平缓的地方减少节点.因此它是一种不均匀区间的积分方法.按照子区间上的积分方式它可以分为自适应辛普森积分法和自适应梯形积分法.通常是采用自适应辛普森积分法作为子区间的积分方式.自适应积分法的基本步骤如下：(1) 将积分区间[],a b 分成两个相等的1级子区间1,2a a h ⎡⎤+⎢⎥⎣⎦和1,2a h a h ⎡⎤++⎢⎥⎣⎦，且h b a =-；(2) 在上述两个1级子区间上用辛普森积分得到积分(){}1_,^12I a a h ⎧⎫+⎨⎬⎩⎭和(){}1_,^12I a h a h ⎧⎫++⎨⎬⎩⎭;(3) 将子区间1,2a a h ⎡⎤+⎢⎥⎣⎦分成两个相等的2级子区间21,2a a h ⎡⎤+⎢⎥⎣⎦和211,22a h a h ⎡⎤++⎢⎥⎣⎦;(4) (){}(){}(){}221111_,^2_,^1_,^12222I a a h I a a h I a h a h ⎧⎫⎧⎫⎪⎪⎪⎪⎧⎫+=++++⎨⎬⎨⎬⎨⎬⎩⎭⎪⎪⎪⎪⎩⎭⎩⎭;(5) 比较(){}1_,^22I a a h ⎧⎫+⎨⎬⎩⎭和(){}1_,^12I a a h ⎧⎫+⎨⎬⎩⎭， 如果|(){}1_,^12I a a h ⎧⎫+⎨⎬⎩⎭- (){}1_,^22I a a h ⎧⎫+⎨⎬⎩⎭|<110**2epsi ，其中epsi 为整体积分所需要的精度，则认为子区间1,2a a h ⎡⎤+⎢⎥⎣⎦上的积分(){}1_,^12I a a h ⎧⎫+⎨⎬⎩⎭已达到所需精度，不需要再细分;否则就需要再细分，对每个2级子区间做同样的判断1级子区间1,2a h a h ⎡⎤++⎢⎥⎣⎦的操作过程完全与上面相同. 3.4 高斯法对已知求积公式可以讨论它的代数精度，反之也可以按照代数精度要求导出求积公式.对于求积公式，当求积节点k x （0,1,,k n =L ）固定时，公式有1n +个待定参数，故此时可要求满足对1,,,n x x L “准确”这样1n +个约束条件，从而使之至少具有n 次代数精度.进一步，可考虑将k x ()0,1,,k n =L 也视为待定参数，这样公式的待定参数就有22n +个，从而可望公式的代数精度达到21n +.此类高精度的求积公式称为高斯型公式，而对应的节点k x ()0,1,,k n =L 称为区间[],a b 上的高斯点. 3.4.1 高斯法的MATLAB 实现代码如下：function g=gaosi(fname,a,b,n,m)switch mcase 1t=0;A=1;case 2t=[-1/sqrt(3),1/sqrt(3)];A=[1,1];case 3t=[-sqrt(0.6),0.0,sqrt(0.6)];A=[5/9,8/9,5/9];case 4t=[-0.8612136,-0.339981,0.339981,0.861136];A=[0.347855,0.652145,0.652145,0.347855];case 5t=[-0.906180,-0.538469,0.0,0.538469,0.906180];A=[0.236927,0.478629,0.568889,0.478629,0.236927];case 6t=[-0.932470,-0.661209,-0.238619,0.238619,0.661209,0.932470];A=[0.171325,0.360762,0.467914,0.467914,0.360762,0.171325];otherwiseerrorendx=linspace(a,b,n+1);g=0; for i=1:ng=g+gsint(fname,x(i),x(i+1),A,t); endfunction g=gsint(fname,a,b,A,t)g=(b-a)/2*sum(A.*feval(fname,(b-a)/2*t+(a+b)/2)); 3.5 多重积分法考虑二重积分(,)d Rf x y A ⎰⎰，它是曲面(,)z f x y =与平面区域R 围成的体积，对于矩形区域{}(,),R x y a x b c y d =≤≤≤≤.可将它写成累次积分(,)((,)d )d bdac Rf x y dx f x y y x =⎰⎰⎰⎰若用复合辛普森公式，可分别将[,],[,]a b c d 分成N ，M 等份，步长,b a d ch k N M--==，先对积分(,)d d c f x y y ⎰应用复合辛普森公式，令i y c ik =+，121()2i yc i k +=++.1101012(,)d (,)4(,)2(,)(,)6M M di M ci i i k f x y y f x y f x y f x y f x y --+==⎡⎤=+++⎢⎥⎣⎦∑∑⎰.从而得：1101012(,)d d (,)4(,)2(,)(,)d 6M M b dbb b b i M aca a a a i i i k f x y y x f x y f x y f x y f x y x --+==⎡⎤=+++⎢⎥⎣⎦∑∑⎰⎰⎰⎰⎰⎰.4 数值积分方法比较例：用数值积分方法计算方程1204d 1I x x =+⎰的值. 4.1 复合辛普森公式求解 >> fun=inline('4./(1+x.^2)');>> xinpusen(fun,0,1,10)ans =3.2749259863031184.2 龙贝格算法求解longbeige(inline('4./(1+x.^2)'),0,1,1e-6)T =3.00003.1000 3.13333.1312 3.1416 3.14213.1390 3.1416 3.1416 3.14163.1409 3.1416 3.1416 3.1416 3.14163.1414 3.1416 3.1416 3.1416 3.1416 3.1416 ans =3.1415926536382444.3 高斯算法求解>> gaosi(inline('4./(1+x.^2)'),0,1,2,3)ans =3.141591222382834>> gaosi(inline('4./(1+x.^2)'),0,1,4,4)ans =3.1415956115587354.4 三种方法比较分析结果显示每一个算法都接近真实值，但龙贝格算法相比较复合辛普森算法，高斯算法来说更加接近.对于代数精度来说，复合辛普森的代数精度为11，龙贝格代数精度为11，高斯代数精度为11.可见代数精度相同时，龙贝格的求积精度最小，所以相同条件下龙贝格求积公式最能接近准确值.总结随着数学实验的兴起，对整个数学课程教学改革起到了积极的推动作用，我们要熟悉的运用各种数学软件，解决数学运算中繁琐的问题，实现学习的简单，快捷化.同时意识到用MATLAB编程时，要实现代码的层次性，做到有规有矩，那样才能把MATLAB 运用自如.这次课程设计，用MATLAB实验对数值积分进行了实现，简介了5种不同的数值积分的方法，并且实现了其中的3中方法，实现过程中发现了各种方法之间的区别和联系.并且在实验过程中，使自己对数值积分和MATLAB更加的熟悉.做到了学习和实践相联系.参考文献[1]戈慈水.《数值分析》课程数值积分的MatLab实现问题的教学研究[M].时代教育出版社，2011.7[2]张志涌.精通MATLAB6.5版[M].北京航天航空大学出版社，2003.3.[3]陈杰，孙晓君.MATLAB数学实验[M].高等教育出版社，2006.6.[4]朱叶志.MATLAB数值分析与应用[M].北京：机械工业出版社,2009.[5]王正林.精通MATLAB科学计算[M].北京：电子工业出版社，2007.[6]萧树铁.大学数学实验[M].北京：高等教育出版社，1997.。