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

东北大学秦皇岛分校数值计算课程设计报告数值积分算法及MATLAB实现学院数学与统计学院专业信息与计算科学学号5133201姓名陈悦指导教师姜玉山张建波成绩教师评语：指导教师签字：2015年07月14日数学与统计学院课程设计（实习）报告第1页1绪论数值分析是计算数学的一个主要部分，计算数学是数学科学的一个分支，它研究用计算机求解各种数学问题的数值检索方其理论与软件的实现.而数值分析主要研究数值计算.现科学技术的发展与进步提出了越来越多的复杂的数值计算问题，这些问题的圆满解决已远人工手算所能胜任，必须依靠电子计算机快速准确的数据处理能力.这种用计算机处理数值问题的方法，成为科学计算.今天，科学计算的应用范围非常广泛，天气预报、工程设计、流体计算、经济规划和预测以及国防尖端的一些科研项目，如核武器的研制、导弹和火箭的发射等，始终是科学计算最为活跃的领域.1.1数值积分介绍数值积分是数值分析的重要环节，实际问题当中常常需要计算积分，有些数值方法，如微分方程和积分方程的求解，也都和积分计算相联系.求某函数的定积分时，在多数情况下，被积函数的原函数很难用初等函数表达出来，因此能够借助微积分学的牛顿-莱布尼兹公式计算定积分的机会是不多的.另外，许多实际问题中的被积函数往往是列表函数或其他形式的非连续函数，对这类函数的定积分，也不能用不定积分方法求解.由于以上原因，数值积分的理论与方法一直是计算数学研究的基本课题.对微积分学做出杰出贡献的数学大师，如I.牛顿、L.欧拉、C.F.高斯、拉格朗日等人都在数值积分这个领域作出了各自的贡献，并奠定了这个分支的理论基础.构造数值积分公式最通常的方法是用积分区间上的n次插值多项式代替被积函数，由此导出的求积公式称为插值型求积公式.特别在节点分布等距的情形称为牛顿-科特斯公式，例如梯形公式(TrapezoidalApproximations)与抛物线公式(ApproximationsUsingParabolas)就是最基本的近似公式.但它们的精度较差.龙贝格算法是在区间逐次分半过程中，对梯形公式的近似值进行加权平均获得准确程度较高的积分近似值的一种方法，它具有公式简练、计算结果准确、使用方便、稳定性好等优点，因此在等距情形宜采用龙贝格求积公式(RhombergIntegration).当用不等距节点进行计算时，常用高斯型求积公式计算，它在节点数目相同情况下，准确程度较高，稳定性好，而且还可以计算无穷积分.数值积分还是微分方程数值解法的重要依据.许多重要公式都可以用数值积分方程导出.现探讨数值积分算法以及运用MATLAB软件的具体实现数学与统计学院课程设计（实习）报告第2页1.2MATLAB软件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数值积分的基本概念一般的，我们可以在区间,ab上适当选取某些节点kx，然后用kfx的加权平均得到平均高度f的近似值，这样构造出的求积公式具有下列形式：0dnakkbkfxxAfx，式中kx称为求积节点；kA称为求积系数，亦称伴随节点kx的权.权kA仅仅与节点kx的选取有关，而不依赖于被积函数fx的具体形式.2.1代数精度的概念如果某个求积公式对于次数不超过m的多项式均能准确的成立，但对于1m次多项式就不准确成立，则称该求积公式具有m次代数精度（或代数精确值）数学与统计学院课程设计（实习）报告第3页一般地，欲使求积公式具有m次代数精度，只要令它对于1,,,mfxxx都能准确成立，这就要求:22+1+1=-,1=-,21=-.+1kkkmmmkkAbaAxbaAxbam2.2求积公式的余项令求积公式的余项为Rf，其中ba0=dnkkkRffxxAfx.；区间,ab可以是有限的或无限的.构造求积公式的问题就是确定jx和jA使得Rf在某种意义下尽可能地小.3数值积分方法及MATLAB实现3.1复合辛普森公式3.1.1插值型求积公式用插值多项式()nLx替换积分*dbaIfxx中的被积函数fx，然后计算dbnaILxx作为积分的近似值，这样建立的求积公式称为插值型求积公式.用插值多项式()nLx的表达式0()nnkkkLxlxfx，代入得0nkkkIAfx，其中:011011()()()()d()()()()bkknkakkkkkknxxxxxxxxAxxxxxxxxx.其余项为:(1)*()d(1)!nbafRfIIxxn.3.1.2牛顿-科特斯公式介绍取等距节点，把积分区间,ab剖分成n等分.令步长bahn，并记0,nxaxb，则1n个节点为0,0,1,kxxkhkn，代入得:数学与统计学院课程设计（实习）报告第4页kA=0(1)(1)(1)(1)()d()!!nknhtttktktntnkk.这种等距节点的插值型求积公式通常称为牛顿-科特斯公式.3.1.3辛普森公式利用牛顿-科特斯公式，取n=2，此时为()4()()62baabSfaffb，即为辛普森公式，其余项为5(4)1()()902sbaRff.3.1.4复合辛普森公式将积分区间,ab分成n等分，分点为(0,1,,)kxakhkn，其中()bahn，记区间1,kkxx的中点为12kx，在每个小区间1,kkxx上用辛普森公式，则得到所谓的复合辛普森公式：1110246nnkkkkhSfxfxfx.111102246nnnkkkkhSfafbfxfx.余项为4(4)()1802nbahRSf，(,)ab.3.1.5复合辛普森公式的MATLAB实现代码如下：functions=xinpusen(fun,a,b,n)h=(b-a)/ns1=0;s2=0;fork=0:(n-1)x=a+h*k;s1=s1+feval(fun,x);endfork=0:(n-1)x=a+h*(k+1/2);s2=s2+feval(fun,x);数学与统计学院课程设计（实习）报告第5页ends=h/6*(feval(fun,a)+feval(fun,b)+2*s1+4*s2);3.2龙贝格公式3.2.1梯形法的递推化将区间,ab分成n等份，共有1n个分点，如果将求积区间再二分一次，则分点增至21n个，用复合梯形公式求得该子区间上的积分值为1122()4kkkhfxfxfx，将每个子区间上的积分值相加得2nT=1111002()()42nnkkkkkhhfxfxfx.得到递推公式：12102122nnnkkhTTfx.3.2.2龙贝格算法公式当fx在,ab上充分光滑时，可证用1()Th逼近I的截断误差是：1()ITh=21ah42ah2kkah按理查森外推法：1()FhFh，1()()(),(1,2,)1mpmmmmpFqhqFhFhmq其中，q为满足10(1,2,)mpqm的适当正数.取序列：142,1,2,41mmnmmhTThThm.用1mTh来逼近I的误差为21mOh，这种算法就是龙贝格算法.3.2.3龙贝格算法MATLAB实现代码如下：functions=longbeige(fun,a,b,tol)ifnargin4，tol=1e-4;endi=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;数学与统计学院课程设计（实习）报告第6页T(i+1,1)=T(I,1)/2+sum(feval(fun,a+h/2:h:b-h/2))*h/2;forj=1:iT(i+1,j+1)=(4^j*T(i+1,j)-T(I,j))/(4^j-1);endendTs=T(i+1,j+1);3.3自适应法自适应积分法是一种比较经济而且快速的求积分的方法.他能自动地在被积函数变化剧烈的区域增多节点，而在被积函数变化平缓的地方减少节点.因此它是一种不均匀区间的积分方法.按照子区间上的积分方式它可以分为自适应辛普森积分法和自适应梯形积分法.通常是采用自适应辛普森积分法作为子区间的积分方式.自适应积分法的基本步骤如下：(1)将积分区间,ab分成两个相等的1级子区间1,2aah和1,2ahah，且hba；(2)在上述两个1级子区间上用辛普森积分得到积分1_,^12Iaah和1_,^12Iahah;(3)将子区间1,2aah分成两个相等的2级子区间21,2aah和211,22ahah;(4)221111_,^2_,^1_,^12222IaahIaahIahah;(5)比较1_,^22Iaah和1_,^12Iaah，如果|1_,^12Iaah-1_,^22Iaah|110**2epsi，其中epsi为整体积分所需要的精度，则认为子区数学与统计学院课程设计（实习）报告第7页间1,2aah上的积分1_,^12Iaah已达到所需精度，不需要再细分;否则就需要再细分，对每个2级子区间做同样的判断1级子区间1,2ahah的操作过程完全与上面相同.3.4高斯法对已知求积公式可以讨论它的代数精度，反之也可以按照代数精度要求导出求积公式.对于求积公式，当求积节点kx（0,1,,kn）固定时，公式有1n个待定参数，故此时可要求满足对1