数值积分算法与MATLAB实现毕业论文设计

毕业设计（论文）设计（论文）题目：数值积分算法与MATLAB实现重庆邮电大学本科毕业设计（论文）-PAGE60-摘要在求一些函数的定积分时，由于原函数十分复杂难以求出或用初等函数表达，导致积分很难精确求出，只能设法求其近似值，因此能够直接借助牛顿-莱布尼兹公式计算定积分的情形是不多的。数值积分就是解决此类问题的一种行之有效的方法。积分的数值计算是数值分析的一个重要分支；因此，探讨近似计算的数值积分方法是有着明显的实际意义的。本文从数值积分问题的产生出发，详细介绍了一些数值积分的重要方法。本文较详细地介绍了牛顿-科特斯求积公式，以及为了提高积分计算精度的高精度数值积分公式，即龙贝格求积公式和高斯-勒让德求积公式。除了研究这些数值积分算法的理论外，本文还将这些数值积分算法在计算机上通过MATLAB软件编程实现，并通过实例用各种求积公式进行运算，分析比较了各种求积公式的计算误差。【关键词】数值积分牛顿-科特斯求积公式高精度求积公式MATLAB软件

ABSTRACTWhenthesolutionofthedefiniteintegralofsomefunctionvalues，becausetheoriginalfunctionisverycomplexanddifficulttofindtheelementaryfunctionexpression,theintegralisdifficulttoaccuratelycalculate,onlymanagedtofindtheapproximatevalue,andthecaseissmallthatallowstodirectinterfacewiththeNewton-Leibnizformulatocalculatethedefiniteintegral.Numericalintegrationisaneffectivemethodtosolvesuchproblems.Thenumericalintegrationisanimportantbranchofnumericalanalysis;therefore,exploringtheapproximatecalculationofthenumericalintegrationmethodhasobviouspracticalsignificance.Thisarticledeparturefromthenumericalintegrationproblem,describedindetailsomeimportantnumericalintegrationmethods.ThispaperhasintroduceddetailtheNewton-Coatesquadratureformula,andinordertoimprovethecalculationaccuracyofnumericalintegrationformulas,MorepreciseformulashaveRombergquadratureformulasandtheGauss-Legendrequadratureformula.Inadditiontothestudyofthesenumericalintegrationalgorithmtheory,thearticlealsoinvolvewhatthesenumericalintegrationalgorithmbeprogrammedbymatlabsoftwareonthecomputer,andanexampleiscalculatedwithavarietyofquadratureformulas,finallyanalysisandcomparisontovariousquadratureformulascalculationerror.【Keywords】NumericalintegrationNewton-CotesquadratureformulaHigh-precisionquadratureformulaMatlabsoftware

目录前言 1第一章牛顿-科特斯求积公式 2第一节数值求积公式的构造 2第二节复化求积公式 9第三节本章小结 12第二章高精度数值积分算法 13第一节梯形法的递推 13第二节龙贝格求积公式 14第三节高斯求积公式 17第四节高斯-勒让德求积公式 19第五节复化两点高斯-勒让德求积公式 22第六节本章小结 23第三章各种求积公式的MATLAB编程实现与应用 24第一节几个低次牛顿-科特斯求积公式的MATLAB实现 24第二节复化求积公式的MATLAB实现 28第三节龙贝格求积公式的MATLAB实现 33第三节高斯-勒让德求积公式的MATLAB实现 34第五节各种求积算法的分析比较 36第六节本章小结 38结论 39致谢 40参考文献 41附录 43一、英文原文 43二、英文翻译 52前言对于定积分，在求某函数的定积分时，在一定条件下，虽然有牛顿-莱布里茨公式可以计算定积分的值，但在很多情况下的原函数不易求出或非常复杂。被积函数的原函数很难用初等函数表达出来，例如等；有的函数的原函数存在，但其表达式太复杂，计算量太大，有的甚至无法有解析表达式。因此能够借助牛顿-莱布尼兹公式计算定积分的情形是不多的。另外，许多实际问题中的被积函数往往是列表函数或其他形式的非连续函数，对这类函数的定积分，也不能用不定积分方法求解，只能设法求其近似值。因此，探讨近似计算的数值积分方法是有明显的实际意义的，即有必要研究定积分的数值计算方法，以解决定积分的近似计算。而数值积分就是解决此类问题的一种有效的方法，它的特点是利用被积函数在一些节点上的信息求出定积分的近似值。微积分的发明是人类科学史上一项伟大的成就，在科学技术中，积分是经常遇到的一个重要计算环节。数值积分是数学上重要的课题之一，是数值分析中重要的内容之一，也是应用数学研究的重点。随着计算机的出现，近几十年来，对于数值积分问题的研究已经成为一个很活跃的研究领域。现在，数值积分在计算机图形学，积分方程，工程计算，金融数学等应用科学领域都有着相当重要的应用，所以研究数值积分问题有着很重要的意义。国内外众多学者在数值积分应用领域也提出了许多新方法。在很多实际应用中，只能知道积分函数在某些特定点的取值，比如天气测量中的气温、湿度、气压等，医学测量中的血压、浓度等等。通过这个课题的研究，我们将会更好地掌握运用数值积分算法求特殊积分函数的定积分的一些基本方法、理论基础；并且通过matlab软件编程的实现，应用于实际生活中。

第一章牛顿-科特斯求积公式第一节数值求积公式的构造大多数实际问题的积分是需要用数值积分方法求出近似结果的。数值积分原则上可以用于计算各种被积函数的定积分，无论被积函数是解析解形式还是数表形式，其基本原理都是用多项式函数近似代替被积函数，用多项式的积分结果近似代替对被积函数的积分。由于所选多项式形式的不同，可以有许多种数值积分方法。而利用插值多项式来构造数值求积公式是最常用的一种方法。对于积分，用一个容易积分的函数去代替被积函数，这样的自然以多项式为最佳，因为多项式能很好的逼近任何连续函数，而且容易求出其原函数。求积公式的推导在积分区间上取有限个点，作的次插值多项式，其中，为次插值基函数。用近似代替被积函数，则得若记则得数值求积公式其中称为求积系数，称为求积节点。则称该求积公式为插值型求积公式。知道了插值型求积公式以及其构造方法。为了便于计算与应用，常将积分区间的等分点作为求积节点，这样构造出来的插值型求积公式就称为牛顿-科特斯（Newton-Cotes）求积公式。在积分区间上取个等距节点，其中，做次拉格朗日插值多项式，因为，所以记截去第二项得显然与无关，只与节点有关。令，则当时，，于是而从而得记则故求积公式可写成这就是牛顿-科特斯求积公式，其中称为科特斯系数。部分科特斯系数取值如下表1.1科特斯系数具有以下特点[1]（1）（2）（3）当³8时，出现负数，稳定性得不到保证。而且当较大时，由于Runge现象，收敛性也无法保证。故一般不采用高阶的牛顿-科特斯求积公式。（4）当£7时，牛顿-科特斯公式是稳定的。表1.1部分科特斯系数表知道了什么是牛顿-科特斯求积公式，下面我们来看它的误差估计，首先来看看牛顿-科特斯求积公式的截断误差。我们知道牛顿-科特斯求积公式是一个插值型数值求积公式，当用插值多项式代替进行积分时，其截断误差即积分真值和近似值之差，推导如下，由插值多项式的误差估计可知，用次拉格朗日多项式逼近函数时产生的误差为其中。对上式两边从到作定积分，便可得出它的截断误差二、几个低次牛顿-科特斯求积公式从上面的讨论可知，用多项式近似代替被积函数进行数值积分时，虽然最高次数可是8，但是8次多项式的计算式非常繁杂的。常用的是下面介绍的几种低次多项式。矩形求积公式定义1.1在牛顿-科特斯求积公式中，如果取QUOTEn=0，用零次多项式（即常数）代替被积函数，即用矩形面积代替曲边梯形的面积，则有称式为矩形求积公式根据牛顿-科特斯求积公式的误差理论式，矩形求积公式的误差估计为2、梯形求积公式定义1.2[1]在牛顿-科特斯求积公式中，如果取QUOTEn=1，用一次多项式代替被积函数，即用梯形面积代替曲边梯形的面积，则有其中，QUOTEfx0=fa,fx1=fb,查表可得QUOTEc01称式为梯形求积公式由于用一次多项式QUOTEp1(x）近似代替被积函数，QUOTEfx，所以它的精度是1。也就是说，只有当被积函数是一次多项式时，梯形求积公式才是准确的。根据牛顿-科特斯求积公式的误差理论式，梯形求积公式的误差估计为QUOTEf''ξ是被积函数QUOTEfx二阶导数在QUOTEx=ξ点的取值，QUOTEξϵ[a,b]3、辛浦生求积公式定义1.3[2]QUOTE313在牛顿-科特斯求积公式中，如果取QUOTEn=2，用二次多项式代替被积函数，即曲边用抛物线代替，则有其中，，，QUOTEx0=a，x1=a+b2,x2=b,查表可得QUOTEc02=称式为辛浦生求积公式，也称抛物线求积公式。它的几何意义是：用过3个点，，QUOTEa,fa,a+b2,fa+b2,b,fb的抛物线和，QUOTEx=a,x=b构成的曲边梯形面积，近似地代替了被积函数QUOTEfx形成的曲边和，QUOTEx=a,x=b构成的曲线梯形面积。下面对辛浦生求积公式的误差进行估计。由于辛浦生求积公式是用二次多项式逼近被积函数推得的，原则上它的代数精度为2.但因多项式次数是偶数，根据定理1.1可知，它的代数精度为3过，和QUOTEb,fb3个点，构造一个QUOTEfx的三次Lagrange插值多项式QUOTEp3(x),且使QUOTEp3‘a+b2=f'a+b2。根据Lagrange插值余项定理得对上式两边从QUOTEa到QUOTEb进行积分，即可得到根据定积分中值定理可知，在QUOTE[a,b]上总有一点满足下述关系：通过变量代换QUOTEt=x-a,dt=dx,,，很容易求得把这个结果代入式QUOTE，便得出辛浦生求积公式的误差估计式4、科特斯求积公式由于和时具有相同的迭代精度，但是时计算量小，故的Newton-Cotes积分公式用的很少。定义1.4QUOTE1.42[18]在牛顿-科特斯求积公式中，如果取QUOTEn=4时，牛顿—科特斯公式为称式为科特斯求积公式。同理根据式可求得其误差估计式三、求积公式的代数精度如果被积函数为任意一个次数不高于次的多项式时，数值求积公式一般形式的截断误差；而当它是次多项式时，，则说明数值求积公式具有次代数精度。一个数值求积公式的代数精度越高，表示用它求数值积分时所需逼近被积函数的多项式次数越高。定理1.QUOTEQUOTE111[3]牛顿-科特斯求积公式的代数精度等于，当为偶数时，牛顿-科特斯求积公式的代数精度等于证明如果被积函数是一个不大于次的多项式，则，即；而当时任意一个次多项式时，，故所以，按照代数精度的定义可知，一般情况下，牛顿-科特斯求积公式的代数精度等于。当为次多项式时，牛顿-科特斯求积公式的代数精度至少等于。若设是一个次多项式，这时为一常数，而因此，只要证明在是偶数时，，，定理就可得证。为此，设，令，于是由于为偶数，不妨设，为正整数，则，于是再引进变换，则，，，代入上式右侧，得出最后的积分中被积函数是奇函数，所以积分结果等于零，定理1.1得证。第二节复化求积公式前面导出的误差估计式表明，用牛顿-科特斯公式计算积分近似值时，步长越小，截断误差越小。但缩小步长等于增加节点，亦即提高插值多项式的次数。龙格现象表明，这样做并不一定能提高精度。理论上已经证明，当QUOTEn→∞时，牛顿-科特斯公式所求得的近似值不一定收敛于积分的准确值，而且随着的增大，牛顿-科特斯公式是不稳定的。因此，实际中不常用高阶牛顿-科特斯公式。为了提高计算精度，可考虑对被积函数用分段低次多项式插值，由此导出复化求积公式。一、复化梯形求积公式在实际应用中，若将积分区间分成若干个小区间，在各个小区间上采用低次的求积分式（梯形公式或辛浦生公式），然后再利用积分的区间可加性，把各区间上的积分加起来，便得到新的求积公式，这就是复化积分公式的基本思想。以梯形面积近似曲边梯形面积，即用梯形公式求小区间上积分的近似值。这样求得近似值显然比用梯形公式计算精度高。定积分存在定理表明，只要被积函数连续，当小区间的长度趋于零时，小梯形面积之和即就趋于曲边梯形面积的准确值，即定积分的准确值。定义1.5[4]QUOTE1.52将积分区间QUOTEa,b进行N进行等分，记为QUOTEh=b-aN,QUOTExk=a+kh(k=0,1,⋯,N-1)在每个小区间QUOTE[xk,xk+1](k=0,1,⋯,N-1)若将所得的近似值记为QUOTETN，整理得称式为复化梯形公式。记为QUOTETN当时，即收敛于如果QUOTEfx∈C(2)[a,b]，则在小区间QUOTE[xk,xk+1]因此因为QUOTEf''x区间QUOTE[a,b]上连续，由介值定理知存在QUOTEξkϵ[a,b]，使得从而有这就是复化梯形公式的截断误差。下面讨论复化梯形公式的数值稳定性。设计算函数值QUOTEfxk时产生的误差为QUOTEδk(k=0,1,⋯,N-1)，则用式计算结果的误差为因此，无论QUOTEN为多大，复化梯形公式都是稳定的。二、复化辛浦生求积公式如果用分段二次插值函数近似被积函数，即在小区间上用Simpson公式计算积分近似值，就可导出复化Simpson公式。定义1.6[5]QUOTE1.615将积分区间QUOTE[a,b]分成QUOTEN=2m等分，分点为，在每个小区间QUOTEx2k-2,x2kk=0,1,⋯,n-1上。用上。用Simpson公式求积分，则有求和得整理后得到式就称为复化辛浦生求积公式。记为QUOTESN如果，QUOTEfxϵC(4)a,b,则由因为QUOTEf4ξ为连续，故存在QUOTEξϵ[a,b],使得代入上式得式表明，步长h越小，截断误差越小。与复化梯形公式的分析相类似，可以证明，当QUOTEN=2m⇢∞时，用复化Simpson公式所求得的近似值收敛于积分值，而且算法具有数值稳定性。三、复化科特斯求积公式定义1.7将积分区间等分为个子区间，每个子区间的中点QUOTEh=b-aN,QUOTExk=a+k4(k=0,1,⋯,N-1)，子区间长度，在每个子区间QUOTE[xk,xk+1](k=0,1,⋯,2N-1)上用科特斯公式求和，得式就称为复化科特斯求积公式，记为，式中，，类似地可以推出复化科特斯公式的截断误差为第三节本章小结本章节开篇介绍了数值求积公式的构造，主要是用运用插值多项式。接着介绍了几个低次的牛顿-科特斯求积公式，即梯形公式、辛浦生公式、科特斯公式，以及牛顿-科特斯求积公式的改进复化求积公式，并对各个求积公式进行了相应的误差分析。第二章高精度数值积分算法复化求积公式是提高精确度的一种行之有效的方法，但是在使用复化型求积公式之前必须先给出步长。步长太大精度难以保证，步长太小则又会导致计算量的增加，而事先给出一个合适的步长往往是困难的。在实际计算中常常采用变步长的计算方法，即在步长逐次减半的过程中，反复利用复化求积公式进行计算，并同时查看相继两次计算结果的误差是否达到要求，直到所求得的积分值满足要求为止。下面以梯形公式为例第一节梯形法的递推在变步长的过程中探讨梯形法的计算规律如下：设将积分区间分为等分，则一共有个等分点，，这里用表示复化梯形法求得的积分值，其下标表示等分数。由余项公式可知，积分值为再将各子区间分半，使得区间成等分。此时所得积分近似值记为，则再由余项公式可知，积分值为假定在上变化不大，即有，于是得，左式也可以写成为这说明用作为积分的近似值时，其误差近似为。计算过程中常用是否满足作为控制计算精度的条件。如果满足，则取作为的近似值；如果不满足，则再将区间分半，直到满足要求为止。实际计算中的递推公式为在给定控制参数后，当满足时，则以作为积分的近似值。通过类似的推导，还可以得到下面的结论对于辛浦生公式，假定在上变化不大，则有对于科特斯公式，假定在上变化不大，则有第二节龙贝格求积公式梯形法的算法简单，单精度低，收敛的速度缓慢。如何提高收敛速度以节省计算量，这是人们极为关心的课题。由此引出了龙贝格公式。由梯形的递推法可以看出，将积分区间等分时，用复化梯形公式计算的结果作为积分的近似值，其误差近似值为。可以设想，如果用这个误差作为的一种补偿，即将作为积分的近似值，可望提高其精确度。直接根据复化求积公式，不难验证这说明，将区间对分前后两次复化梯形公式的值，按式作线性组合恰好等于复化辛浦生公式的值，它比更接近于近似值。同样，根据式用于作线性组合会得到比更精确的值，且通过直接验证可得再由式用与作线性组合，又可得到比更精确的值，通常记为,即式就称为龙贝格求积公式[2]。上述用若干个积分近似值推算出更为精确的积分近似值的方法，称为外推方法。我们将序列,和分别称为梯形序列、辛浦生序列、科特斯序列和龙贝格序列。由龙贝格序列当然还可以继续进行外推，得到新的求积序列。但由于在新的求积序列中，其线性组合的系数分别为与。因此，新的求积序列与前一个序列结果相差不大。故通常外推到龙贝格序列为止。利用龙贝格序列求积的算法称为龙贝格算法。这种算法具有占用内存少、精确度高的优点。因此，成为实际中常用的求积方法。下面给出龙贝格求积算法的计算步骤：第1步：算出和的值，根据公式计算；第2步：将分半，算出后，根据公式计算,再根据公式计算；第3步：再将区间分半，算出及的值，并根据公式和计算出及，再由公式计算出；第4步：将区间再次分半，计算，，，并由公式计算；第5步：将区间再次分半，类似上述过程计算，，，。重复上述过程可计算得到，，，一直算到龙贝格序列中前后两项之差的绝对值不超过给定的误差限为止。定义2.1[6]上述用若干个积分近似值算出更精确的积分近似值的方法，称之为外推法。上述计算步骤也可用表2.1表示表2.1龙贝格计算步骤表QUOTE2.12K梯形序列辛浦生序列科特斯序列龙贝格序列012345﹕﹕﹕﹕﹕第三节高斯求积公式前面介绍的个节点的Newton-Cotes求积公式，其特征是节点是等距的。这种特点使得求积公式便于构造，复化求积公式易于形成。但同时也限制了公式的精度。是偶数时，代数精度为，是奇数时，代数精度为；我们知道个节点的插值型求积公式的代数精确度不低于。设想：能不能在区间上适当选择个节点使插值求积公式的代数精度高于？答案是肯定的，适当选择节点，可使公式的精度最高达到，这就是本节所要介绍的高斯型求积公式。例：其中，固定在，可以适当选取，只有两个自由度，得到的是梯形公式，其代数精度只有1。如果对求积节点也进行适当的选取，将有四个自由度，得到如下公式：这个积分公式的代数精确度为3，这就是高斯型求积公式，上面的求积节点称为高斯点。一、高斯型求积公式和高斯点对于含有个参数，的求积公式：，适当选取这个参数，可以使得数值积分公式的代数精确度达到，我们称这一类求积公式为高斯型求积公式，称这类求积公式的积分节点为高斯点，系数称为高斯系数。定义2.2[7]QUOTE2.1[4]如果个求积节点的求积公式的代数精度为，则这个求积节点称为高斯点。因为高斯求积公式也是插值型求积公式，故有结论：个节点的插值型求积公式的代数精度，满足二、高斯点的特征定理2.1[7]QUOTE2.1[8]设是个相异点，以这个点为零点的次多项式为则是高斯点的充要条件是对于任意不超过次的多项式，成立证明1）必要性。设为高斯点，则对任意不超过次的多项式，均有，则对于任意次数不超过次多项式，是次数不超过的多项式，且注意到，，故2）充分性。设对于一切次数不超过次的多项式，成立，又设是次数不超过次的多项式，用去除，商，余，即，可知和均是不超过次的多项式，从而，又因求积公式是插值多项式的构造导出的，由的选取，其代数精确度可以达到，而是次数不超过次的多项式，因此成立，因，所以，，故而也即由于是次数不超过次的多项式，因此该积分公式的代数精确度至少为，因而由定义2.2知节点是高斯点。定义2.3对于关系，我们称之为正交性，即与任意次多项式正交，而这样的多项式类称为正交多项式。三、高斯型求积公式的余项QUOTE项[4]其中，，。积分，可以计算得到第四节高斯-勒让德求积公式常用的高斯型求积公式有高斯-勒让德公式、高斯-切比雪夫公式、高斯-拉盖尔公式、高斯-埃尔米特公式，下面着重介绍高斯-勒让德公式一、Legendre多项式，称式为勒让德（Legendre）多项式，其具有前面提到的正交性性质，即对于任意次数不超过的多项式，成立因此，多项式的零点就是相应的高斯-勒让德求积公式的高斯点。勒让德多项式的前几项如下：，，，勒让德多项式的性质：性质1勒让德多项式的首项系数为性质2当为奇数时，为奇函数；当为偶数，为偶函数性质3对一切次数不高于次的多项式，有二、高斯-勒让德求积公式1）当时，，其零点为，易得，的高斯-勒让德求积公式是，其代数精确度为12）当时，，其零点为，，此时，的高斯-勒让德求积公式是：，其代数精确度为33）当时，，其零点为和0，设高斯-勒让德求积公式为：令，得到方程eq\o\ac(○,1)令，得到方程eq\o\ac(○,2)令，得到方程eq\o\ac(○,3)联立方程eq\o\ac(○,1)、eq\o\ac(○,2)、eq\o\ac(○,3)可解得，所以，其代数精确度为54）当时，，由，得，故四个零点为，即，，，相应的，可以解出，，高斯-勒让德求积公式为其代数精确度为7综上，高斯-勒让德求积公式的积分节点和积分系数表如表2.2表2.2积分节点和积分系数表三、高斯-勒让德求积公式的余项对阶求积公式（共有个求积节点）其中，，为多项式的零点。积分，可以计算得到：四、一般区间上的高斯-勒让德求积公式前面讲解的高斯-勒让德求积公式是计算标准区间上的定积分，对于一般的有界区间上的定积分，可以通过变量代换转化为区间的定积分。即在积分中，令，当时，；当时，；且。这个变换为线性变换，其反变换为于是成立然后，我们可以用区间上的高斯-勒让德求积公式。第五节复化两点高斯-勒让德求积公式本节对两点公式进行推广和复化，得到了新的数值积分公式同时分析了它的积分误差和收敛阶定义2.4[8]QUOTE2.23设，是一种复化求积公式，若当时成立则称求积公式是阶收敛的。例如，复化的Trapezoid公式和复化的Simpson公式分别具有二阶和四阶收敛性。定理2.2[8]QUOTE2.23设，，则复化两点高斯-勒让德求积公式为的余项表达式为该方法具有四阶收敛性。考虑积分，准确解为0.74684204…用复化两点高斯-勒让德公式可求得=0.74680332，其绝对误差为=0.00001960，与复化梯形公式，复化Simpson公式以及原两点Gauss-Legendre公式相比，这里构造的复化两点高斯-勒让德求积公式方法在运算量和精度方面都有显著的改善。第六节本章小结本章主要介绍了精确度比较高的两种数值求积公式，即龙贝格求积公式和高斯型求积公式，对其进行了相应的误差分析。其中高斯型求积公式主要介绍了高斯-勒让德求积公式，并对两点高斯-勒让德求积公式进行了复化。第三章各种求积公式的MATLAB编程实现与应用MATLAB是由MathWorks公式开发的一种主要用于数值计算及可视化图形处理的工程语言，是当今最优秀的科技应用软件之一。它将数值计算、矩阵运算、图形图像处理、信号处理和仿真等诸多强大的功能集成在较易使用的交互计算机环境中，为科学研究、工程应用提供了一种功能强、效率高的编程工具。下面我们将各种求积算法通过MATLAB软件编程实现，以下程序均用MATLAB7.0编写，运行坏境：1、硬件环境CPU（intelCorei3-2310M,2.1GHz），内存（2GB昱联）,2、软件环境windows7（32位）操作系统。以下总共编写了六个算法程序，部分代码参考文献[10-13],为了体现程序的正确性，以下程序都以为例进行运算。原积分的精确值为几个低次牛顿-科特斯求积公式的MATLAB实现先用M文件定义一个名为f1.m的函数：%i是要调用第几个被积函数g(i)，x是自变量functionf=f1(i,x)g(1)=sqrt(x);ifx==0g(2)=1;elseg(2)=sin(x)/x;endg(3)=4/(1+x^2);f=g(i);程序一QUOTE一[12]：function[C,g]=NCotes(a,b,n,m)%a，b分别为积分的上下限；%n是子区间的个数；%m是调用上面第几个被积函数；%当n=1时计算梯形公式；当n=2时计算辛浦生公式，以此类推；i=n;h=(b-a)/i;z=0;forj=0:ix(j+1)=a+j*h;s=1;ifj==0s=s;elsefork=1:js=s*k;endendr=1;ifi-j==0r=r;elsefork=1:(i-j)r=r*k;endendifmod((i-j),2)==1q=-(i*s*r);elseq=i*s*r;endy=1;fork=0:iifk~=jy=y*(sym('t')-k);endendl=int(y,0,i);C(j+1)=l/q;z=z+C(j+1)*f1(m,x(j+1));endg=(b-a)*z1）当输入，时，即在MATLAB命令窗口输入>>NCotes(0,1,1,2)即可得用梯形公式的积分值和相应科特斯系数如图3.12）当输入，时，即在MATLAB命令窗口输入>>NCotes(0,1,2,2)即可得用辛浦生公式的积分值和相应科特斯系数如图3.23）当输入，时，即在MATLAB命令窗口输入>>NCotes(0,1,4,2)即可得用科特斯公式的积分值和相应科特斯系数如图3.3图3.1图3.2图3.3复化求积公式的MATLAB实现一、复化梯形求积公式的MATLAB实现通过QUOTEfx的QUOTEn+1个等步长节点逼近积分其中，，QUOTExk=a+kh,x0=a,x程序二：functions=trapr1(f,a,b,n)%f是被积函数；%a,b分别为积分的上下限；%n是子区间的个数；%s是梯形总面积；h=(b-a)/n;s=0;fork=1:(n-1)x=a+h*k;s=s+feval('f',x);endformatlongs=h*(feval('f',a)+feval('f',b))/2+h*s;先用M文件定义一个名为f.m的函数：functiony=f(x)ifx==0y=1;elsey=sin(x)/x;end在MATLAB命令窗口中输入>>trapr1('f',0,1,4)回车得到如图3.4图3.4若取子区间的个数在MATLAB命令窗口中输入>>trapr1('f',0,1,8)回车得到如图3.5图3.5复化辛浦生求积公式的MATLAB实现程序三QUOTE三[13]：functions=simpr1(f,a,b,n)%f是被积函数；%a,b分别为积分的上下限；%n是子区间的个数；%s是梯形总面积,即所求积分数值；h=(b-a)/(2*n);s1=0;s2=0;fork=1:nx=a+h*(2*k-1);s1=s1+feval('f',x);endfork=1:(n-1)x=a+h*2*k;s2=s2+feval('f',x);ends=h*(feval('f',a)+feval('f',b)+4*s1+2*s2)/3;先用M文件定义一个名为f.m的函数：functiony=f(x)ifx==0y=1;elsey=sin(x)/x;end在MATLAB命令窗口中输入>>simpr1('f',0,1,4)回车得到如图3.6图3.6若取子区间个数时在MATLAB命令窗口中输入>>simpr1('f',0,1,8)回车得到如图3.7图3.7复化科特斯求积公式的MATLAB实现程序四：functions=cotespr1(f,a,b,n)%f是被积函数；%a,b分别为积分的上下限；%n是子区间的个数；%s是梯形总面积,即所求积分数值；h=(b-a)/n;s1=0;s2=0;s3=0;s4=0;fork=1:nx=a+(4*k-3)*h/4;s1=s1+feval('f',x);endfork=1:nx=a+(4*k-2)*h/4;s2=s2+feval('f',x);endfork=1:nx=a+(4*k-1)*h/4;s3=s3+feval('f',x);endfork=1:(n-1)x=a+4*k*h/4;s4=s4+feval('f',x);ends=h*(7*feval('f',a)+7*feval('f',b)+32*s1+12*s2+32*s3+14*s4)/90;同样先用M文件定义一个名为f.m的函数：functiony=f(x)ifx==0y=1;elsey=sin(x)/x;end在MATLAB命令窗口中输入>>cotespr1('f',0,1,4)回车得到如图3.8图3.8第三节龙贝格求积公式的MATLAB实现构造数表来逼近积分其中。表示数表的最后一行，最后一列的值。程序五QUOTE五[16]：function[R,quad,err,h]=romber(f,a,b,n,delta)%f是被积函数%a,b分别是积分的上下限%n+1是T数表的列数%delta是允许误差%R是T数表%quad是所求积分值M=1;h=b-a;err=1J=0;R=zeros(4,4);R(1,1)=h*(feval('f',a)+feval('f',b))/2while((err>delta)&(J<n))|(J<4)J=J+1;h=h/2;s=0;forp=1:Mx=a+h*(2*p-1);s=s+feval('f',x);endR(J+1,1)=R(J,1)/2+h*s;M=2*M;forK=1:JR(J+1,K+1)=R(J+1,K)+(R(J+1,K)-R(J,K))/(4^K-1);enderr=abs(R(J,J)-R(J+1,K+1));endquad=R(J+1,J+1)先用M文件定义一个名为f.m的函数：functiony=f(x)ifx==0y=1;elsey=sin(x)/x;end在MATLAB命令窗口中输入>>romber('f',0,1,5,0.5*(10^(-8)))回车得到如图3.9图3.9高斯-勒让德求积公式的MATLAB实现程序六QUOTE六7[8]：function[A,x]=Guass1(N)i=N+1;f=((sym('t'))^2-1)^i;f=diff(f,i);t=solve(f);forj=1:ifork=1:iX(j,k)=t(k)^(j-1);endifmod(j,2)==0B(j)=0;elseB(j)=2/j;endendX=inv(X);forj=1:iA(j)=0;x(j)=0;fork=1:iA(j)=A(j)+X(j,k)*B(k);x(j)=x(j)+t(j);endx(j)=x(j)/k;endfunctiong=GuassLegendre(a,b,n,m)%a,b分别是积分的上下限；%n+1为节点个数；%m是调用f1.m中第几个被积函数；[A,x]=Guass1(n);g=0;fori=1:n+1y(i)=(b-a)/2*x(i)+(a+b)/2;f(i)=f1(m,y(i));g=g+(b-a)/2*f(i)*A(i);end用M文件分别把上面两个自定义函数定义为名为Guass1.m函数和GuassLegendre.m函数用M文件定义一个名为f1.m的函数functionf=f1(i,x)g(1)=sqrt(x);ifx==0g(2)=1;elseg(2)=sin(x)/x;endg(3)=4/(1+x^2);f=g(i);在MATLAB命令窗口中输入>>GuassLegendre(0,1,2,2)>>GuassLegendre(0,1,3,2)回车得到如图3.10图3.10各种求积算法的分析比较例1分别用不同的方法计算积分，并作比较。下面用各种求积公式分别计算积分，并给出了相应的计算误差，进行比较，结果如下：1、用Newton-Cotes公式当n=1时，即用梯形公式，用程序一在MATLAB命令窗口中输入>>NCotes(0,1,1,2)得0.927035492403950.01904757796323当n=2时,即用Simpson公式,用程序一在MATLAB命令窗口中输入>>NCotes(0,1,2,2)得0.946145882273590.000062811906407当n=4时,即用科特斯公式，用程序一在MATLAB命令窗口中输入>>NCotes(0,1,4,2)得0.946083004063670.0000000663035132、用复化梯形公式令h=1/8=0.125，用程序二在MATLAB命令窗口中输入>>trapr1('f',0,1,8),得=0.94569086352700.0003922068401823、用复化辛浦生公式令h=1/8=0.125，用程序三在MATLAB命令窗口中输入>>simpr1('f',0,1,8),得=0.946083085384950.0000000150177674、用Romberg公式，用程序五在MATLAB命令窗口中输入>>romber('f',0,1,5,0.5*(10^(-8)))，得0.946083070367180.0000000000000025、用高斯-勒让德求积公式，令，用2个节点的Gauss公式0.94604115827633用3个节点的Gauss公式，用程序六在MATLAB命令窗口中输入>>GuassLegendre(0,1,2,2)，得0.9460831340784730.000000063711290算法比较此例题的精确值为0.946083070367183由例题的各种求积算法可知：对Newton-cotes公式，当n=1时只有1位有效数字，当n=2时有3位有效数字，当n=4时有7位有效数字。对复化梯形公式有2位有效数字，对复化辛浦生公式有7位有效数字。用复化梯形公式，对积分区间二分了11次用2049个函数值，才可得到7位准确数字。用Romberg公式对区间二分3次，用了9个函数值，能得到同样的结果；二分4次，用了14个函数值，却能得到14位有效数字。用高斯-勒让德求积公式仅用了3个函数值，就能得到同样比较精确的6位有效数字。第六节本章小结本章主要对各种求积算法，即牛顿-科特斯、复化求积、龙贝格、高斯-勒让德求积公式的算法在计算机上运用MATLAB软件进行编程实现，从而简化了运算。并且通过实例计算，误差分析对各种求积公式进行了分析比较。

结论本文主要讨论了数值积分的计算方法并通过MATLAB软件编程实现，通过前面的研究我们知道求数值积分近似值的计算方法很多，包括Newton-Cotes求积公式、复化求积公式、Romberg求积公式、高斯求积公式等等。其中Newton-Cotes方法是一种利用插值多项式来构造数值积分的常用方法，这其中梯形积分方法的误差最大，近似效果最差，Simpson方法的精度比梯形积分高了一个数量级，它的代数精度比梯形积分的代数精度高，能更好地近似积分值；Cotes积分方法的误差比Simpson积分精度高两个数量级。因此，一般情况下，代数精度越高，积分公式计算精度也越高。但是高阶的Newton-Cotes方法的收敛性没有保证，因此，在实际计算中很少使用高阶的Newton-Cotes公式。复化梯形积分方法比单独的梯形积分精度高，它的积分精度和被积函数有关，还和复化积分时的步长有关。复化Simpson积分公式比单独的Simpson积分公式高近7个数量级，效果明显。Romberg方法收敛速度快、计算精度较高，但是计算量较大。Gauss求积方法积分精度高、数值稳定、收敛速度较快，但是节点与系数的计算较麻烦、而且要求已知积分函数。一般来说，Newton-Cotes方法的代数精度越高，数值积分的效果越好、越精确。当积分区间比较大的时候，可以采用复化积分方法可以得到更好的效果；变步长积分方法不仅可以很好地控制计算误差，并且可以寻找到适当的积分步长；Romberg积分方法可以更好地利用变步长复化积分公式得到的积分序列从而得到更为精确的数值结果，是一个较好的数值积分方法。高斯求积方法精确度高，收敛性快也是一种很优秀的数值积分方法。

致谢行文至此，我的这篇论文已接近尾声；岁月如梭，我四年的大学时光也即将敲响结束的钟声。离别在即，站在人生的又一个转折点上，心中难免思绪万千，一种感恩之情油然而生。首先感谢重庆邮电大学四年来对我的培养，是博学的老师们教会了我学习的方法、锻炼了我思考的能力、指明了我未来奋斗的方向，从而使我进一步明确了人生的目标。其次，我要感谢我的指导老师—郑继明教授，他的严谨细致、一丝不苟的作风一直是我工作、学习中的榜样；他的循循善诱的教导和不拘一格的思路给予我无尽的启迪。在撰写整个毕业论文的过程当中，他为我们考虑到了每一个细节，从开题报告到毕业论文的拟定修改上，郑老师更是不厌其烦的为我们做好每一步的细心指导。对此，我表示衷心地感谢。没有郑老师，我的论文也不可能这么顺利的完成。同时，我也要感谢每一位给过我帮助的老师和同学，在我撰写论文的过程当中同样给了我大量有益的建议，在此一并向他们表示真诚的感谢，感谢他们对我的支持和帮助。最后感谢这篇论文所涉及到的各位学者，本文引用了数位学者的研究文献，如果没有各位学者的研究成果带给我的的帮助和启发，我将很难完成本篇论文的写作。由于我的学术水平有限，所写论文难免有不足之处，恳请各位老师和学友批评指正。最后，衷心感谢评阅论文及参加答辩的各位老师！

参考文献张德丰等．MATLAB数值计算方法［M］．北京：机械工业出版社，2010．郑继明等．数值计算方法［M］．重庆：重庆大学出版社，2005．樊守芳．Newton-Cotes数值求积公式的注记［J］．枣庄学院学报，2011，28(2)：186-190．刘小伟．基于MATLAB的复合梯形数值积分法的研究与实验［J］．甘肃联合大学学报(自然科学版)，2010，24(4)：20-23．余丹．用Simpson公式进行数值积分［J］．华北电力大学数理系学报，2010，(23)：189-190．Davis，P．J．andP．Rahinowitz．MethodsofNumericalintegration(secondedition)[M]．AcademicPress.NewYork，1984．曹丽华．一类广义Gauss型求积公式［J］．深圳大学数学物理学报，2007，27(3)：524-534．郭晓斌．复化两点Gauss-Legendre公式及其误差分析［J］．西北师范大学学报，2010，4(3)：49-51．王建强．多种数值积分方法比较分析［J］．武汉大学测绘轩辕学报，2010，2(1)：104-106．刘玉娟，陈应祖．龙贝格积分法及其应用编程［J］．重庆科技学院学报，2007，9(1)：97-99．朱叶志．MATLAB数值分析与应用［M］．北京：机械工业出版社，2009．LiF，LiX．Theneighbor-scatteringnumbercanbecomputedinpolynomialtimeforintervalgraphs[J]．ComputersandMathematicswithApplications，2007，54：679-686．伍丽华，周玲丽．数学软件教程［M］．广州：中山大学出版社，2008．熊华，杨国孝．一类振荡函数的数值积分方法［J］．北京理工大学学报，1999，19(3)：280-284．Varma．AK．Onoptimalquadratureformulae［J］．StudiaSciMathHungar，1984，19：437-446．戴立辉．一些数值积分公式中间点的渐近性［J］．闽江学院学报，2011，32(2)：5-7．WANGXH，HEYG，ZENGZz．Numericalintegrationstudybasedontrianglebasisneuralnetworkalgorithm[J]．JournalofElectronics&InformationTechnology．2004．26(3)：394-399．GUORu．ANoteonNewton-Cote’SNumericalIntegralFormula［J］.JournalofXinjiangUniversity(NaturalScienceEdition)，2010，27(2)：186-190．刘小伟．基于MATLAB的变步长梯形数值积分法的研究与实验［J］．廊坊师范学院学报(自然科学版)，2010，10(1)：39-42.

附录一、英文原文：AnOptimalAdaptativeNumericalIntegrationMethodA.Nicolet,A.Genon,W.Legros,M.Ume,F.Delince*ABSTRACTThispaperpresentsanadaptativeextensionoftheGaussianintegrationmethod.ItiswellknownthattheGaussianintegrationmethodisoptimalforsufficientlysmoothfunctions(i.e.whichmaybeapproximatedbyapolynomial)inthesensethatitgivesthemaximumaccuracyforagivennumberofnodes.Unfortunatelyitisnotalwayspossibletochooseapriorithenumberofnodesfortheintegration.OnealternativeistotrysuccessiveGaussianformulaewithanincreasingnumberofpointsuntiltheyagreewiththerequiredaccuracy.Inthiscase,mostoftheadvantagesofthemethodarelost.AlessaccuratebutnaturallyadaptativemethodsuchastheRombergmethodmaybecomeabettersolution.Theideaoftheoptimaladaptativemethodistofindaseriesofintegrationformulatewithanincreasingnumberofnodesinorderthatthesetofabscissaeoflowerorderformulaeisasubsetofabscissaeofhigherorderformulae.Then,thesequentialevaluationofformulaeofincreasingorderonlyrequirestheadditionofnewpoints.Underthisconstraint,theremainingdegreesoffreedom(thenewabscissaandalltheweightfactors)areusedtoobtainformulaeofmaximumorder.INTRODUCTIONSomewell-knownandclassicalmethods[1]arereviewedinordertosituatethenewmethod.Thosemethodsarevalidforsufficientlysmoothintegrands.TrapezoidalmethodAmongstthesimplest,thismethodconsistsinchoosingequally-spacedpointsbetweentheendpointsandtoapproximatethefunctionbypiecewiselinearfunctions.Thetrapezoidalruleis:With:andRombergmethodThebasicideaistousetheresultsfromsuccessiverefinementsofthetrapezoidalrule:TheRichardsonextrapolationisappliedtothissequenceinordertoeliminatehighordererrorterms[2].Thismethodisnaturallyadaptative.Foreachrefinement,anewtrapezoidalapproximationiscomputedwiththenumberofpointsmultipliedbytwoandreusingthepreviouslycomputedvaluesofthefunction.ThentheRichardsonextrapolationisappliedtothenewsequence.Thisprocessisrepeateduntiltherequiredaccuracyisreached.Thenumberoffunctionevaluationsisnotknownapriorianddependsontheintegrand.GeneralruleMostoftheintegrationrulestointegratethefollowingexpression:havethefrom:Theapproximationisalinearcombination,withweightfactors,ofvaluesofthefunction,fornabscissae.Intherestofthepaper,theintegral(3),i.e.theintegralontheinterval[a,b]oftheproductofthefunctionwiththekernel,willbereferredas“theintegralof”.Formulae(4)aretabulatedintheliteratureforsomekernelsandanassociatedinterval[a,b].Orthogonalpolynomials[3]associatedtotheset‘kernel-interval[a,b]’playanimportantroleinthetheory(table1).KernelaIntervalbAssociatedorthogonalpolynomials1-11Legendre0Laguerre-Hermite-11Tschebychew(firstkind)01Orthogonalpolynomialsassociatedtoon[0,1](Berthod-zaborowskiformulae)Table1.Orthogonalpolynomials[4]Theorthogonalityofpolynomials,isexpressedby[3]:Anumericalintegrationruleischaracterizedbyanintegerpsuchthatallthepolynomialsoforderklessthanorequaltopareexactlyintegrated:GaussianmethodGaussianintegrationrulesareoptimalinthesensethatpismaximal.Annpointrulethathas2ndegreesoffreedom(nabscissaeandnweightfactors)canintegrateexactlyallthepolynomialsuptotheorder2n-1(i.e.polynomialsthathaveupto2ncoefficients).Thenecessaryandsufficientconditionisthatallthepowersofxuptotheorder2n-1arecorrectlyintegratedontheconsideredinterval:Therelations(7)constituteasystemof2nequationswith2nunknowns(theandthe)whosesolutiongivestheparametersofintegrationrule.Classicalapproachistointroduceanauxiliarypolynomialwhoserootsaretheabscissae:Sumsofequations(7)weightedbythecoefficientsareconstructed,inordertodoappear:forRelations(9)constituteasystemofnequationswithaunknownsthatgivesthecoefficients:forThesolutionofthesystem(10)determinesthepolynomialwhoserootsaretheabscissaeoftheintegrationrule.Arelationshipbetweenthispolynomialandtheorthogonalpolynomialsassociatedtotheproblemmaybefound.istheorthogonalpolynomialofdegreenassociatedtotheproblem.Anypolynomialofdegree2n-1maybeexpressedas:Withaquotientpolynomialandarestpolynomial,bothofdegreeatmostequalton-1.Thepolynomialisintegratedexactly:Thepolynomialbeingorthogonaltoallthepolynomialsofdegreelessorequalton,thefirsttermofthelefthandmemberisequaltozero.Therelation(12)isalwaystrueonlyifthefirsttermoftherighthandmemberiscancelled,whatisassuredifisequaltozerofor.Thusisequaltoapartfromaconstantfactor.Thetheoryoftheorthogonalpolynomialsassuresthattheirrootsaresimple,realandsituatedintheinterval[a,b].Theabscissaebeingdetermined,themustbecomputedinorderthatanypolynomialoforderlessthanore