1、1（20_届）本科毕业设计信息与计算科学浅谈MATLAB在数值积分中的应用2摘要本文归纳总结了求数值积分的几个重要方法，探讨了这些方法在MATLAB上的实现。借助MATLAB中的GUI功能，设计一个简单实用的积分操作界面。在该界面上可以实现数据处理和相关图形的功能，从而实现数值积分的可视化。关键字数值积分MATLABGUI可视化NUMERICALINTERGRATIONOFMATLABAPPLICATIONABSTRACTTHISARTICLESUMMARIZESSEVERALIMPORTANTNUMERICALINTEGRATIONMETHODS,THENDISCUSSESTHESEMETH
2、ODSINMATLABSIMPLEMENTATIONWECANDESIGNASIMPLEANDPRACTICALUSERINTERFACEWITHGUITHEINTERFACECANBEACHIEVEDINDATAPROCESSINGANDRELATEDGRAPHICSFEATURES，THATENABLESTHEVISUALIZATIONOFNUMERICALINTEGRATIONKEYWORDSNUMERICALINTEGRATIONMATLABGUIVISUALIZATION前言在微积分中，我们已知函数FX在闭区间,AB上连续，且存在原函数FX，求FX在该区间上的定积分可用牛顿莱布尼兹公
3、式求解，即BBAAFXDXFXFBFA。用牛顿莱布尼兹公式计算积分的方法在理论上和解决实际问题中起到很大的作用，但在工程技术领域常遇到复杂的情况无法用牛顿莱布尼兹公式求解。如某些被积函数FX的原函数FX不是初等函数；或者被积函数FX非常复杂；或者被积函数的原函数结构非常复杂；再者就是没有具体的FX表达形式，只有一些零散的测试数据、表格或图形。诸如此类情况，很多可采用“数值积分”的方法求出积分近似值。在积分所限定的区间中，被积函数FX的插值式在N个点处求值，这些点称为节点。节点用IX表示，并假设它们有序且各不相同（即，12NAXXXB）。如果1XA且NXB，那么区间就是闭的，用符号,AB表示，如
4、果1XA且NXB，那么区间就是开的，用符号,AB表示。我们已经知道，分段多项式的插值比单个多项式的全局插值更加优秀，这个结论对数值积分方法也适用。在每个小段上对FX进行低阶多项式逼近。再对每个小段上的逼近多项式积分时，就得到基本公式。基本公式只涉及用足够的,XFX来定义分段多项式的某一段，将此公式应用到每个小段并把结果相加就得到了复合公式，或称为扩展公式。31几个重要的数值积分方法数值积分的基本思路用数值方法近似地求一个函数FX在区间,AB上定积分BAFXDX的基本思路，可以归结到定积分的定义1LIMNBKANKBAIFXDXFN其中1,KKKXXAB1,2,KN。当N的取值充分大时，I的数值
5、积分就是1,NBNKAKBAFXDXIFN（1）1,NBKAKBAFXDXFRFN（2）（1）和（2）式都称为数值求积分式，其中RF表示截断误差，称为数值求积公式的余项。如果数值求积公式（1）式对N次多项式精确地成立，则称此公式具有N次代数精度。显然，K取值不同，数值积分NI的结果就不同，当然，定积分I的值是一样的。这种做法相当于用相对简单的阶梯函数,1,2,KFKN代替FX作积分。实际上各种不同的数值积分方法就主要在于研究选择什么样的简单函数代替FX，使得既能保证一定的精度，又能使计算量足够小。11矩形公式我们知道，如果在区间,AB上，则定积分0FX表示曲线I下的面积，不妨先从图形上看看如何
6、近似计算这块面积（参见图1）。将区间,AB进行N等分，积分步长为/HBAN。记01KNAXXXXB。在小区间上用矩形面积近似FX下面曲边梯形的面积，矩形公式中可以取做端点函数值为小矩形的高，或取右端点的函数值为小矩形的高，于是在整个区间,AB内构成台阶形。容易看出，两个台阶围城面积分别为10NNKHLHFX，/HBAN341NNKHRHFX，/HBAN4在图1中，两个台阶形分别小于和大于所求面积。（3）和（4）式就是计算定积分的矩形数值积分，简称矩形公式。图112梯形公式将矩形公式的两个公式平分，则每个小区间上的小矩形变为小梯形（图1中线段1IIMM）。设函数FX在区间,AB上的1N个等距节点
7、0,1,2,IBAXAKHKNHN处有定义，且其函数值为IFX，则在闭区间0,1,2,KN上存在101,2NNKNKHBATHFXFXFXHN（51）或11,2NNKKKHBATFXFXHN（52）使得,BNNAIFXDXTRFT（6）其中,NRFT是NT的截断误差，即余项，（51）和（52）式就是常用的复合梯形公式，它相当于用分段线性插值函数作为FX的近似。设函数FX在区间,AB上的1N个等距节点0,1,2,IBAXAKHKNHN处的函XY5数值为,1,2,KFXKN，且FX在,AB上具有连续的二阶导数，则复合梯形公式（51）或52是2阶收敛的，且存在,AB，使得（51）的截断误差,NRFT
8、为2,12BNNNAHBARFTFXDXTF如果记2MAX,AXBBAMFXHN则在,AB上的绝对误差为22,12NHRFTMBA。13辛普森公式为提高精度可以用分段二次插值函数代替21N。由于每段要用到相邻两个小区间端点的三个函数值，所以区间数目必须是偶数，记2,0,1,NMKM。在第K段的两个区间上用三个节点222121,KKKKXFXXFX，222,KKXX作二次插值函数KSX，得2222212243KKXKKKKXHSKDXFXFXFX，求M段之和就得到整个区间上的近似积分。函数FX在区间,AB上的21M个等距节点0,1,2,2,2KBAXAKHKMHM处有定义，且其函数值为0,1,2
9、,2KFXKM，则在,AB上有12211124,32MMNKKKKHBASFAFBFXFXHM（7）使得,BNNAIFXDXSRFS（8）其中,NRFS是NS的截断误差，（7）式称为辛普森公式。函数FX在区间,AB上的21M个等距节点0,1,2,2,2KBAXAKHKMHM处的函数值为FX0,1,2KM，且FX在,AB上具有连续的四阶导数，则复合辛普森公式是四阶收敛的，且存在,AB，使得（7）式的截断误差为44,180BNNAHBARFSFXDXSF（9）如果记44MAX,AXBMFX则式（7）在,AB上的绝对误差为644,180BNNAHRFSFXDXSMBA（10）14牛顿科茨公式我们用简
10、单的、便于积分的、又逼近于被积函数FX的函数PX代替FX来构造求积公式。由于多项式不但计算方便，而且容易积分，因此，常取PX为一个多项式。函数FX在闭区间,AB上的1N个等距节点0,1,2,KBAXAKHKNHN处有定义，且其函数值为0,1,2,KFXKN，则在闭区间,AB上存在N阶牛顿科茨公式NNNKKKONCFBACFX（11）使得,BNNAFXDXNCFRFNCF（12）其中科茨系数0,1,2,NKCKN为01111NKNNKCTTTKTKTNDTNKNK（13）从N阶科茨系数（13）式可以看出，科茨系数只与N有关，与积分区间和被积函数无关。因此，只要给出N,就可以算出科茨系数，从而写出
11、N阶牛顿科茨公式计算定积分BAFXDX时，一方面由于它是由（12）式去掉余项,NRFNCF得到的，因而产生截断误差,NRFNCF；另一方面，由于计算机的字长是有限的，函数值可能带有误差，并且计算NNCF的过程中还会有舍入误差。如果将区间,AB分成N等份，则步长BAHN,10,1,2,KXAKHKN则有复合科茨求积公式111100003732123214790424NNNNBKKKKAKKKKHHHHFXDXFAFXFXFXFXFB1415龙贝格公式为了提高复合求积公式的精度，可以利用区间逐次减半的方法。我们看到近似求积公式的误差随着N的增加而减小，但是对于给定的FX和精度，在求积过程中，用二分
12、法每次将N增加一倍，直到满足精度要求。7根据复合辛普森公式的截断误差公式推出222441NNNSSC，得到收敛速度更快的科茨序列20KKC1248,CCCC在科茨序列20KKC基础上还可以重复同样的方法，根据复合科茨公式的截断误差公式推出323441NNNCCR得到龙贝格序列20KKR如1248,RRRR为了统一上面的符号，记11,12,MNMNTT，则2,122,MNMNTT。将上面的一般形式表示为1,11,1,14,2,3,41JMJMJMJJTTTJMJ（M）15或,11,1,11,2,3,41MJMJMJMJJTTTTJMJ（M）16从而得到龙贝格求积公式。函数FX在区间,AB上有定义
13、且在1N个节点012NAXXXXB处的函数值为0,1,2,KFXKN则在区间AB上用二分法每次将N增加一倍，可得到龙贝格求积公式。11,11221,11,111111,11,1,1,2111,2,3,222224,2,3,41IIIIIIIIKJMJMJMJJHTFAFBHBAHBATTFAKHHHITTTJMJ,（M）17函数FX在闭区间,AB上具有2J阶连续导数，则在闭区间,AB上用二分法每次将N增加一倍，得到龙贝格求积公式中的式15是2J收敛的，即2,2,3,BJMJAFXDXTOHJM（MJ）,18截断误差为8222,BJJJMJJMJAFXDXTOHBHF19,2,3,JM（MJ）其
14、中1,2MJMBAHB为依赖于J的常数，且,MJAB。在实际计算中，并不需要作出序列,1MT后，再作序列,2MT。而是在算出1,T和2,1T后就计算2,2T，计算得2,2T和3,2T后就计算3,3T，以此类推。1,T，2,2T，3,3T，,JJT，计算MJT的顺序详见龙贝格积分表，龙贝格算法是数值稳定的。当,1,JJJJTT或1,1,JJJJTT时，计算停止。16高斯勒让德积分公式如果函数FX在区间1,1上连续，则有代数精度为3N的两点高斯勒让德积分公21133GFFF20如果函数FX在区间1,1上具有4阶连续导数，则1221FXDXGFRF21其中两点高斯勒让德积分公式的截断误差为24111
15、35FRFN个节点的高斯勒让德积分公式及其误差分析。一般的，N个节点的高斯勒让德积分公式的代数精度为21N，其数值积分公式为,11,22,NNNNNNGFAFXAFXAFX,2211NNFXDXGFRF其中截断误差为2121112NNNFRFXDXN（23）其中12NNXXXXXXX。91,11NNKKNKFXDXAFXRFN横坐标KX系数,NKA截断误差NRF11213，131，14135F333,0,55585,999615750F40861136311603399810436034785484510652145154983472875F50906179845905394693101000
16、0000000002369268851047862867050568888888810237732650F6093246951420661209386502386191861017132449240360761573004679139346134123261213F70949107912307415311856040584515140000000000001294849662027970539150381830050504179591837154143271415F80960289856507966664774052553240990183434642501012285363022238103
17、450313706645903626837834174163281617F2用MATLAB求数值积分21矩形公式命令SUM用矩形公式计算数值积分可用下面函数SUM和CUMSUM实现，格式如下（一）函数SUM的调用格式SUMX，如果输入向量X，则输出向量X的所有元的和，可用于按矩10形公式（3）和（4）计算积分；如果输入矩阵X，则输出矩阵X的每列向量的元的和。SUMX,DIM，输入ND数组X，输出为X的每个DIM向量的元的和。（二）函数CUMSUM的调用格式CUMSUMX，输入向量X，输出向量X的元依次累加和，可用于按矩形（3）和（4）计算积分；输入矩阵X，输出矩阵为X的每列向量的元依次累加。C
18、UMSUM（X,DIM）输出ND数组X，输出为X的每个DIM向量的元依次累加和。算例一分别用MATLAB和矩形公式计算SIN220XEDX，并与精确值比较。解将0,2分成20等份，运行程序所以，用矩形公式计算的结果为31023，精确值为31044，误差为00021。22梯形公式命令TRAPZ用MATLAB软件和梯形公式计算数值积分和估计误差，可以用自己编写的MATLAB程序，也可以用系统提供的两种计算程序一种是梯形数值积分的程序TRAPZM；另外一种是累加梯形数值积分的程序CUMTAPZM。（一）函数TRAPZ的调用格式TRAPZY，输入Y，输出为按梯形公式51计算的Y的积分的近似值。TRAP
19、ZX，Y，输入X,Y为同长度的数组，输出为按梯形公式（51）计算Y对X的积分。TRAPZ（X,Y,DIM），输出为按梯形公式51计算Y的每个DIM的积分，其中X的长度必须和SIZEY,DIM相同。（二）函数CUMTRAPZ的调用格式CUMTRAPZ（Y），输入Y，输出用梯形方法计算的Y的累加积分的近似值。CUMTRAPZ（X,Y），输入X,Y为同长度的数组，输出用梯形方法计算的Y对X的累加积分。CUMTRAPZ（X,Y,DIM），输出用梯形方法计算的Y的每个DIM的累加积分，其中X的长度必须11和SIZEY,DIM相同。算例二用TRAPZ和CUMTRAPZ分别计算20SINXEXDX，精确到1
20、04。解输入程序可见，梯形公式计算结果与精确值FS03961很接近。但是，函数CUMTRAPZ计算的结果展示了每次累加的计算过程。23辛普森公式命令QUAD用MATLAB软件和梯形公式计算数值积分和估计误差，可以用自己编写的MATLAB程序外，MATLAB系统还提供了QUADM。调用格式QUADFUN,A,B计算函数YFUNX在区间,AB上的数值积分，自动选择步长，输出数值积分。其中函数YFUNX是以FUNM文件命名的M文件函数，或者是FINLINEFUN形式，或者数组，函数YFUNX将接受向量的自变量，并且返回结果是向量Y，即在X的每个元素处的积分估计值。QUADFUN,A,B,TOL，同上
21、，但指定绝对误差限为TOL，默认值为610。Q,FCNTQUAD，输出数值积分值Q和函数赋值的编号FCNT。辛普森数值积分MATLAB程序FUNCTIONYCOMSIMPSONFUN2,A,B,NZ1FEVALFUN2,AFEVALFUN2,BMN/2HBA/2MXAZ20Z30X20X30FORK222MX2XKHZ2Z22FEVALFUN2,X212ENDFORK322MX3XKHZ3Z34FEVALFUN2,X3ENDYZ1Z2Z3H/3算例三用COMSIMPSONM和QUADM分别计算定积分211220XIEDX，取精度为104。解建立FUN2M文件命名的M文件函数输入程序44牛顿公式
22、命令RNCN阶牛顿科茨公式的MATLAB程序（一）估计误差的MATLAB程序FUNCTIONRNCNCENSUK1PN/2FIXN/2IFP0FORK1N2SUKSUKKENDSUKSYMSTABFXN2,SUT2FORU1NSUSUTUENDSUINTFINTSU,T,0,NYDOUBLEINTFRNCBA/NN3FXN2ABSY/SUKELSEFORK1N1SUKSUKKENDSUKSYMSTABFXN1,SUTFORU1NSUSUTUENDSUINTFINTSU,T,0,NYDOUBLEINTFRNCBA/NN2FXN1ABSY/SUK13END二计算科茨系数0,1,2,NKCKN的程序
23、FUNCTIONCN,RNCNNEWCOTENSYMSTABM,FZZEROS1,N1CNZEROS1,N1SUTK1M1M01FORU1NSU1SUTUM01M0USUSU1M0M01ENDSUMOF1SU/T0INTF1INTF1,T,0,NYDOUBLEINTFLCN11N0Y/NM0K1M1FORJ1NK1KJM1MNJFSU/TJINTFINTF,T,0,NYDOUBLEINTFCNJ11NJY/NK1M1WARNINGOFFMATLABDIVIDEBYZEROENDFNSU/TNINTFNINTFN,T,0,NYDOUBLEINTFNCNN1Y/NM0CNSUK1PN/2FIXN/
24、2IFP0FORK1N2SUKSUKKENDSUKSYMSTABFXN2,SUT2FORU1NSUSUTUENDSUINTFINTSU,T,0,NYDOUBLEINTFRNCNBA/NN3FXN2ABSY/SUKELSEFORK1N1SUKSUKKENDSUKSYMSTABFXN1,SUTFORU1NSUSUTUENDSUINTFINTSU,T,0,N,YDOUBLEINTFRNCNBA/NN2FXN1ABSY/SUKEND三在MALAB中提供了用自适递归8阶牛顿科茨公式计算数值积分的程序QUAD1M。算例四用8阶牛顿科茨公式计算定积分SIN20XIEDX，取精度104。解输入程序所以，用8阶
25、牛顿科茨求积公式计算定积分的值，迭代33次，计算结果31044。25龙贝格命令ROMBERG14用MATLAB软件和龙贝格公式计算定积分BAFXDX的近似值有两种常用的方法（一）根据具体题目编写MATLAB程序计算；（二）根据（19）式编写程序，然后计算时调用。ROMBERGMFUNCTIONRT,R,WUGU,HROMBERGFUN3,A,B,WUCHA,MN1HBAWUGU1XAK0RTZEROS4,4RT1,1HFEVALFUN3,AFEVALFUN3,B/2WHILEWUGUWUCHAHH/2S0FORJ1NXAH2J1SSFEVALFUN3,XENDRTK1,1RTK,1/2HSN2
26、NFORI1KRTK1,I14IRTK1,IRTK,I/4I1ENDWUGUABSRTK1,KRTK1,K1ENDRRTK1,K1算例五取精度为810，分别用1,11,JJJJTT和1,1,JJJJTT作为计算停止的条件，用ROMBERGM程序计算15011DXX，龙贝格积分表，最小步长H，取精度为810，并用精确值比较，然后取精度为710，用1,1,JJJJTT作为计算停止的条件，观察用龙贝格求积公式计算的结果与精确值的绝对误差是否满足710。解取精度为810，用1,11,JJJJTT作为计算停止的条件，输入程序再用1,1,JJJJTT作为计算停止的条件，输入程序15由此可见，用作为计算停止
27、的条件，递归6层，得到近似值0916290731972162。用作为计算停止条件，递归7层，计算得到近似值为0916290731874269。26高斯命令GAUSS高斯勒让德积分公式的MATLAB程序FUNCTIONGL,Y,RGNGAUSSR1FUN5,X,ANLENGTHXN22NYFEVALFUN5,XGLSUMAYSUN1SU2N1WOME1SYMSXFORK1NWOMEWOMEXXKENDWOME2WOME2FRINTWOME2,X,1,1FORK1N2SU2NSU2NKENDSYMSMRGNFRM/SU2N算例六用高斯勒让德积分公式计算212112XIEDX，取代精度为3和5，再根
28、据截断误差公式写出误差公式，并将结果与精度值进行比较。解建立以FUN5M为名的M文件输入程序16所以，用代数精度为3和5的高斯勒让德积分公式数值计算I的结果1GL和2GL，与精确值IFS的绝对误差分别为1WGL和2WGL。3数值积分可视化MATLAB是一种集计算功能、符号运算、数据可视化等强大功能的数学工具软件。但是如果只有在我们掌握了足够的MATLAB命令和编程技巧，才能进行数值积分，显然感到很麻烦。而MATLAB中的GUI功能，可以借助软件本身的数值计算、图形处理功能及GUIDE编辑器，设计出一个简洁美观的积分操作界面，然后在界面上实现数据处理和相关图形的功能，这就实现了数值积分可视化。设
29、计思路在MATLABGUI中创建图形界面，然后在图形界面上实现处理数据及显示相关图形功能，添加控件，最后添加控件的事件代码实现处理数据及显示相关结果的功能。数值积分模块不仅可以计算常见的一元、二元、三元函数积分，还能绘制出二元函数的三维图形和积分区域。用户在此界面可以进行积分动画演示；从不同角度查看积分函数图像与积分区域。通过数值积分可视化，可以为一些不熟悉数值积分或者编程的用户提供一个友好的人机交互平台。既能让我们从繁琐的计算过程中解脱出来，17又能提高工作效率。4总结子曰“工欲善其事，必先利其器”。如果有一种有效的工具能解决在教学与研究中遇到的问题，那么MATLAB语言正是这样的一种工具。
30、它可以将使用者从繁琐、无谓的底层编程中解放出来，把有限的宝贵时间更多地花在解决问题中，提高工作效率。目前，MATLAB已经成为国际上最流行的科学与工程计算的软件工具。它已经不仅仅是一个“矩阵实验室”了，它已经成为了一种具有广泛应用前景的全新的计算机高级编程语言了,有人称它为“第四代”计算机语言，在国内外高校和研究部门正扮演着重要的角色。MATLAB语言的功能也越来越强大，不断适应新的要求提出新的解决方法。数值积分和微分是解许多工程和科研问题的很重要的慨念，虽然有许多微积分问题，可以用解析的方式找到答案，但其计算过程庞大且复杂，尤其是有更多的微积分式是无解析解存在的。所以这时候必须用数值方法求近
31、似解。而MATLAB提供庞大的数学运算函数，是数值积分在MATLAB中的计算变得简单容易。可以预见，在科学运算、自动控制与科学绘图领域MATLAB语言将长期保持其独一无二的地位。参考文献1于育民，连冬艳MATLABGUI在数值积分中的应用J南阳理工学院学报，2009，611171192林洽武，MATLAB解复杂一重定积分J广东教育学院学报，广东教育学院，2010，33043463张书欣，付冰冰基于MATLAB的积分求解析解方法的研究J吉林农业科技学院学报，2010，2191211244施吉林，刘淑珍，陈桂芝计算机数值方法（第三版）M高等教育出版社，20081211595余承依基于MATLAB数
32、值计算的GUI设计J漳州师范学院学报（自然科学版），2004，41738406郭瑞，阿布都热西提阿布都外力关于NEWTONCOTES数值积分公式的修正J新疆大学学报,2010,2271861907蔡苗苗MATLAB使用技巧几则J电脑学习，2004644458谢彦红，李金娜，赵旋基于MATLAB的GAUSS积分节点系数研究J沈阳农业大学学报，2010,412254256189刘庄，李有道计算方法（四）R清华大学基础知识讲座，第二章数值积分，535810陈佩宁，刘竞用MATLAB求数值积分的方法J石家庄职业技术学院学报，2008,6（20）586011余丹用SIMPSON公式进行数值积分J华北电力大学数理系，2010,2370212任玉杰数值分析及其MATLAB实现M高等教育出版58572913薛定宇，陈阳泉高等应用数学问题的MATLAB求解（第二版）M清华大学出版社15117514CURTISFGERALD,PATRICKOWHEATLEYAPPLIEDNUMERICALANALYSISM机械工业出版社15GERALDRECKTENWALDNUMERICALMETHODSWITHMATLABIMPLEMENTATIONANDAPPLICATONM机械工业出版社48147
