数值积分的算法比较及其MATLAB实现

哈尔滨工业大学本科毕业设计（论文）编 号： 审定成绩： 重庆邮电大学毕业设计（论文）设计（论文）题目：数值积分算法与 MATLAB 实现学 院 名 称 ： 数理学院学 生 姓 名 ：专 业 ： 数学与应用数学班 级 ：学 号 ： 指 导 教 师 哈尔滨工业大学本科毕业设计（论文）：答辩组 负责人 ：填表时间： 年 月重庆邮电大学教务处制哈尔滨工业大学本科毕业设计（论文）- I -摘 要在求一些函数的定积分时，由于原函数十分复杂难以求出或用初等函数表达，导致积分很难精确求出，只能设法求其近似值，因此能够直接借助牛顿-莱布尼兹公式计算定积分的情形是不多的。数值积分就是解决此类问题的一种行之有效的方法。积分的数值计算是数值分析的一个重要分支；因此，探讨近似计算的数值积分方法是有着明显的实际意义的。本文从数值积分问题的产生出发，详细介绍了一些数值积分的重要方法。本文较详细地介绍了牛顿-科特斯求积公式，以及为了提高积分计算精度的高精度数值积分公式，即龙贝格求积公式和高斯-勒让德求积公式。除了研究这些数值积分算法的理论外，本文还将这些数值积分算法在计算机上通过 MATLAB 软件编程实现，并通过实例用各种求积公式进行运算，分析比较了各种求积公式的计算误差。【关键词】数值积分 牛顿-科特斯求积公式 高精度求积公式 MATLAB 软件哈尔滨工业大学本科毕业设计（论文）- II -ABSTRACTWhen the solution of the definite integral of some function values，because the original function is very complex and difficult to find the elementary function expression, the integral is difficult to accurately calculate, only managed to find the approximate value, and the case is small that allows to direct interface with the Newton - Leibniz formula to calculate the definite integral. Numerical integration is an effective method to solve such problems. The numerical integration is an important branch of numerical analysis; therefore, exploring the approximate calculation of the numerical integration method has obvious practical significance. This article departure from the numerical integration problem, described in detail some important numerical integration methods.This paper has introduced detail the Newton - Coates quadrature formula, and in order to improve the calculation accuracy of numerical integration formulas, More precise formulas have Romberg quadrature formulas and the Gauss - Legendre quadrature formula. In addition to the study of these numerical integration algorithm theory, the article also involve what these numerical integration algorithm be programmed by matlab software on the computer, and an example is calculated with a variety of quadrature formulas, finally analysis and comparison to various quadrature formulas calculation error.【Key words】 Numerical integration Newton-Cotes quadrature formula High-precision quadrature formula Matlab software哈尔滨工业大学本科毕业设计（论文）- III -目 录前 言 .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哈尔滨工业大学本科毕业设计（论文）- 1 -前 言对于定积分 ，在求某函数的定积分时，在一定条件下，虽然有()bafxd牛顿-莱布里茨公式 可以计算定积分的值，但在很()IfFba多情况下 的原函数不易求出或非常复杂。被积函数 的原函数很难用()f ()fx初等函数表达出来，例如 等；有的函数 的原函数 存2sin(),xfe ()Fx在，但其表达式太复杂，计算量太大，有的甚至无法有解析表达式。因此能够借助牛顿-莱布尼兹公式计算定积分的情形是不多的。另外，许多实际问题中的被积函数往往是列表函数或其他形式的非连续函数，对这类函数的定积分，也不能用不定积分方法求解，只能设法求其近似值。因此，探讨近似计算的数值积分方法是有明显的实际意义的，即有必要研究定积分的数值计算方法，以解决定积分的近似计算。而数值积分就是解决此类问题的一种有效的方法，它的特点是利用被积函数 在一些节点上的信息求出定积分的近()fx似值。微积分的发明是人类科学史上一项伟大的成就，在科学技术中，积分是经常遇到的一个重要计算环节。数值积分是数学上重要的课题之一，是数值分析中重要的内容之一，也是应用数学研究的重点。随着计算机的出现，近几十年来，对于数值积分问题的研究已经成为一个很活跃的研究领域。现在，数值积分在计算机图形学，积分方程，工程计算，金融数学等应用科学领域都有着相当重要的应用，所以研究数值积分问题有着很重要的意义。国内外众多学者在数值积分应用领域也提出了许多新方法。在很多实际应用中，只能知道积分函数在某些特定点的取值，比如天气测量中的气温、湿度、气压等，医学测量中的血压、浓度等等。通过这个课题的研究，我们将会更好地掌握运用数值积分算法求特殊积分函数的定积分的一些基本方法、理论基础；并且通过 matlab 软件编程的实现，应用于实际生活中。哈尔滨工业大学本科毕业设计（论文）- 2 -第一章 牛顿-科特斯求积公式第一节 数值求积公式的构造大多数实际问题的积分是需要用数值积分方法求出近似结果的。数值积分原则上可以用于计算各种被积函数的定积分，无论被积函数是解析解形式还是数表形式，其基本原理都是用多项式函数近似代替被积函数，用多项式的积分结果近似代替对被积函数的积分。由于所选多项式形式的不同，可以有许多种数值积分方法。而利用插值多项式来构造数值求积公式是最常用的一种方法。对于积分 ，用一个容易积分的函数 去代替被积函数 ，这样()bafxd ()x()fx的 自然以多项式 为最佳，因为多项式能很好的逼近任何连续函数，()()nL而且容易求出其原函数。一、 求积公式的推导在积分区间 上取有限个点 ，作 的 次插,ab01naxxb ()fxn值多项式 ，其中， 为 次插值基函数。0()()nkLxflx(),)kl用 近似代替被积函数 ，n f则得 0()()()(nbb bnkaa akfxdLxdfxldx(1.)若记011()()()()()bbkknkkaak kAlx xxx (.2)则得数值求积公式0()()nbkafdAf (1.3)哈尔滨工业大学本科毕业设计（论文）- 3 -其中 称为求积系数， 称为求积节点。则称该求积公式为插值型求积公式。kAkx知道了插值型求积公式以及其构造方法。为了便于计算与应用，常将积分区间的等分点作为求积节点，这样构造出来的插值型求积公式就称为牛顿-科特斯（Newton-Cotes）求积公式。在积分区间 上取 个等距节点 ，其中,ab1nkxah(0,1)n，做 次拉格朗日插值多项式 ，因为 ，所以bhn()nL)(fxLRx()()()bbnnaaafxdLxRxd(1)10 ()()!bbnk naakfl fxd记 1()()(bbnkkaakkxAlxd (1.4)(1)1()!bnn naRffd (.5)截去第二项得 0()nbkafxdAfx显然 与 无关，只与节点 有关。令 ，则当kA()f ,1kn xath时， ，于是,xb,tn111()()()2()nnahtt (1.6)而 0111()kkkkknxxxx !()!n从而得 0(1)1(1)()!kkhAttktktnd 记 () 0()()()()!)nknkCttktkt (1.7)哈尔滨工业大学本科毕业设计（论文）- 4 -则 ()(nkkAbaC故求积公式 可写成 (1.3)()0()nbkafxdafx (1.8)这就是牛顿-科特斯求积公式，其中 称为科特斯系数。()nkC部分科特斯系数取值如下表 1.1科特斯系数具有以下特点 1（1） （2） ()01niC()()nnii（3）当 8 时，出现负数，稳定性得不到保证。而且当 较大时，由于nRunge 现象，收敛性也无法保证。故一般不采用高阶的牛顿 -科特斯求积公式。（4）当 7 时，牛顿-科特斯公式是稳定的。n表 1.1 部分科特斯系数表知道了什么是牛顿-科特斯求积公式，下面我们来看它的误差估计，首先来看看牛顿-科特斯求积公式的截断误差。我们知道牛顿-科特斯求积公式是哈尔滨工业大学本科毕业设计（论文）- 5 -一个插值型数值求积公式，当用插值多项式 代替 进行积分时，其截()nLx()f断误差 即积分真值 和近似值之差，推导如下RfI，由插值多项()()()()bbbbnnaaaafIxdfxxdfxd式的误差估计可知，用 次拉格朗日多项式 逼近函数时产生的误差为nnL(1)1() ()!nbnnaffxLxd (1.9)其中 。对上式两边从 到 作定积分，便可得10(),(,)niinab出它的截断误差(1)1()()!bnn naRffxd (1.0)二、几个低次牛顿-科特斯求积公式从上面的讨论可知，用多项式近似代替被积函数进行数值积分时，虽然最高次数可是 8，但是 8 次多项式的计算式非常繁杂的。常用的是下面介绍的几种低次多项式。1、 矩形求积公式定义 1.1 在牛顿-科特斯求积公式中，如果取 ，用零次多项式（即0n常数）代替被积函数，即用矩形面积代替曲边梯形的面积，则有(0)0 0()()()(bbaafxdLxdacfxbafx (1.)称式 为矩形求积公式1.根据牛顿-科特斯求积公式的误差理论式 ，矩形求积公式的误差估(1.0)计为 (01)001()()!bafRxdfba哈尔滨工业大学本科毕业设计（论文）- 6 -2、梯形求积公式定义 1.21 在牛顿-科特斯求积公式中，如果取 ，用一次多项式代1n替被积函数，即用梯形面积代替曲边梯形的面积，则有 (1)(1)10()()bbaafxdLxacfxf其中， , 查表可得 代入上式得出0ff()()121()()2bbaaxxdfb (1.2)称式 为梯形求积公式(1.2)由于用一次多项式 近似代替被积函数 ，所以它的精度是 1。也1()Lx()fx就是说，只有当被积函数是一次多项式时，梯形求积公式才是准确的。根据牛顿-科特斯求积公式的误差理论式 ，梯形求积公式的误差估(1.0)计为 (2) 31 ()(! 2baf baRxdxf是被积函数 二阶导数在 点的取值，()f()f,3、辛浦生求积公式定义 1.32 在牛顿-科特斯求积公式中，如果取 ，用二次多项式代 2n替被积函数，即曲边用抛物线代替，则有 (2)()()2012()() bbaafxdLxacfxfcfx其中， ， ， ,查表可得 ， ，代入01()(2)6()13上式得出21()()()(63bbaa abfxdLxafff (.)称式 为辛浦生求积公式，也称抛物线求积公式。它的几何意义是：1.3用过 3 个点 ， ， 的抛物线和 ， 构(,)f(,()2bf(,)fxab成的曲边梯形面积，近似地代替了被积函数 形成的曲边和 ，x哈尔滨工业大学本科毕业设计（论文）- 7 -构成的曲线梯形面积。xb下面对辛浦生求积公式的误差进行估计。由于辛浦生求积公式是用二次多项式逼近被积函数推得的，原则上它的代数精度为 2.但因多项式次数是偶数，根据定理 1.1 可知，它的代数精度为 3过 ， 和 3 个点，构造一个 的三次(,)af(,()2baf(,)bf ()fxLagrange 插值多项式 ,且使 。根据 Lagrange 插值余3Lx3()2ab项定理得(4)(4) 23() )()(!fffx xx,ab对上式两边从 到 进行积分，即可得到ab(4) 2231()()()(!baabRffLxdf d (1.4)根据定积分中值定理可知，在 上总有一点满足下述关系：,(4) 2)()(bafxxb(4) 2)()(abf xd通过变量代换 , ，很容易求得tdt 52()()(120ba ax把这个结果代入式 ，便得出辛浦生求积公式的误差估计式1.455() (4)2()!20280bbaRfff,ab(1.5)4、科特斯求积公式由于 和 时具有相同的迭代精度，但是 时计算量小，故3n22n的 Newton-Cotes 积分公式用的很少。定义 1.4 在牛顿-科特斯求积公式中，如果取 时，牛顿科特斯公4式为 01234()7()32()()()7()9baafxdfxffxffx (1.6)哈尔滨工业大学本科毕业设计（论文）- 8 -称式 为科特斯求积公式。同理根据式 可求得其误差估计式(1.6) (1.0)6()42(954baRff(,ab(1.7)三、求积公式的代数精度如果被积函数 为任意一个次数不高于 次的多项式时，数值求积公()fxn式一般形式的截断误差 ；而当它是 次多项式时， ，则0Rf(1)0Rf说明数值求积公式具有 次代数精度。一个数值求积公式的代数精度越高，n表示用它求数值积分时所需逼近被积函数的多项式次数越高。定理 1.13 牛顿-科特斯求积公式的代数精度等于 ，当 为偶数时，n牛顿-科特斯求积公式的代数精度等于 1n证明 如果被积函数 是一个不大于 次的多项式，则 ，()fx (1)0nfx即 ；而当时任意一个 次多项式时， ，故 0Rf (1)0nfxR所以，按照代数精度的定义可知，一般情况下，牛顿- 科特斯求积公式的代数精度等于 。当 为 次多项式时， 牛顿-科特斯求n()fxn(1)nf(,)ab积公式的代数精度至少等于 。若设是一个 次多项式，这时 为(1)nf一常数，而 (1)01()()()!nb bn na afRffxLdxxd因此，只要证明在 是偶数时， ，定理就可得证。为此，设01()()0b nax Rf，令 ，1,2iih ,xathndxht于是 101 0()()()2()b naxxdtndt 由于 为偶数，不妨设 ， 为正整数，则 ，于是nk,k20 0(1)2()(1)(1)(1)(2ttthtttttkt 再引进变换 ，则 ， ， ，代入上式右侧，得出ukukutd,k哈尔滨工业大学本科毕业设计（论文）- 9 -0(1)2()ntttnd 1()1(1)(kukuukdu 2222()()k d最后的积分中被积函数是奇函数，所以积分结果等于零，定理 1.1 得证。第二节 复化求积公式前面导出的误差估计式表明，用牛顿-科特斯公式计算积分近似值时，步长越小，截断误差越小。但缩小步长等于增加节点，亦即提高插值多项式的次数。龙格现象表明，这样做并不一定能提高精度。理论上已经证明，当时，牛顿-科特斯公式所求得的近似值不一定收敛于积分的准确值，而n且随着 的增大，牛顿-科特斯公式是不稳定的。因此，实际中不常用高阶牛顿-科特斯公式。为了提高计算精度，可考虑对被积函数用分段低次多项式插值，由此导出复化求积公式。一、复化梯形求积公式在实际应用中，若将积分区间分成若干个小区间，在各个小区间上采用低次的求积分式（梯形公式或辛浦生公式） ，然后再利用积分的区间可加性，把各区间上的积分加起来，便得到新的求积公式，这就是复化积分公式的基本思想。以梯形面积近似曲边梯形面积，即用梯形公式求小区间上积分的近似值。这样求得近似值显然比用梯形公式计算精度高。定积分存在定理表明，只要被积函数连续，当小区间的长度趋于零时，小梯形面积之和即就趋于曲边梯形面积的准确值，即定积分的准确值。定义 1.5 4 将积分区间 进行 等分，记为 ,abNbahN在每个小区间 上用梯形公式求和，得kxah1kx(0,1)1 100()()2kNNb kkaxkfdfdfxf哈尔滨工业大学本科毕业设计（论文）- 10 -若将所得的近似值记为 ，整理得NT1()()2()2Nb kNa khfxdfafbfxT (1.8)称式 为复化梯形公式。记为 (1.8) NT当 时， N101()()2Nkkkfxhfxh()()()bbbaaafdffxd即 收敛于NT()bafxd如果 ，则在小区间 上，梯形公式的截断误差为(2),C1,kx131()()()2kxkkkhhfdfxff 1(,)kkx因此 310()()()NbT kakRffxTf因为 区间 上连续，由介值定理知存在 ，使得f, ,kab10()()Nkkff从而有32()()()(121bTNahbaRffxdTfhf (1.9)这就是复化梯形公式的截断误差。下面讨论复化梯形公式的数值稳定性。设计算函数值 时产生的误差()kfx为 ，则用式 计算结果的误差为(0,1)k (1.8)00012max()maxNkk kNNhnhb因此，无论 为多大，复化梯形公式都是稳定的。哈尔滨工业大学本科毕业设计（论文）- 11 -二、复化辛浦生求积公式如果用分段二次插值函数近似被积函数，即在小区间上用 Simpson 公式计算积分近似值，就可导出复化 Simpson 公式。定义 1.65 将积分区间 分成 等分，分点为 ，,ab2Nmkxah在每个小区间 上。用(0,1)kN bh,(1,)kxNSimpson 公式求积分，则有 2222212()()4()()6kxkkkkxfdfxfxf 2212()3kkkhfff求和得 21()()kmbxafxdfd22121()4()()3kkkkhfxfxf整理后得到12211()()()4()mmbN kka kSfxdfafbfxfx (.0)式 就称为复化辛浦生求积公式。记为 (1.20) NS如果 ， 则由 Simpson 插值余项公式可得复化公式的截断(4),fxCab误差为 12211()()()()4()3mmbS kka khRffdffbfxfx5(4)1280mkf2,k因为 为连续，故存在 , 使得(4)fab哈尔滨工业大学本科毕业设计（论文）- 12 -(4)(4)1mkkff代入上式得5(4)4()12()80180mSkhbaRffhf(,)ab(1.2)式 表明，步长 越小，截断误差越小。与复化梯形公式的分析相类.2似，可以证明，当 时，用复化 Simpson 公式所求得的近似值收敛2N于积分值，而且算法具有数值稳定性。三、复化科特斯求积公式定义 1.7 将积分区间 等分为 个子区间 ，每个子区间的,abN4,kx中点 , ，子区间长度 ， 在每个子区间上用科特42kx(0,1)N bah斯公式求和，得()bafxd 4342117()2()()90Nkkkkhfafxfx414132()()7(NNkkkkfxffbNC(1.2)式 就称为复化科特斯求积公式，记为 ，式中， ，(1.) bahN4kxa(0,12,)k类似地可以推出复化科特斯公式的截断误差为() 6()42()945NbahRff(,)b(1.23)第三节 本章小结本章节开篇介绍了数值求积公式的构造，主要是用运用插值多项式。接着介绍了几个低次的牛顿-科特斯求积公式，即梯形公式、辛浦生公式、科特斯公式，以及牛顿-科特斯求积公式的改进复化求积公式，并对各个求积公式哈尔滨工业大学本科毕业设计（论文）- 13 -进行了相应的误差分析。第二章 高精度数值积分算法复化求积公式是提高精确度的一种行之有效的方法，但是在使用复化型求积公式之前必须先给出步长。步长太大精度难以保证，步长太小则又会导致计算量的增加，而事先给出一个合适的步长往往是困难的。在实际计算中常常采用变步长的计算方法，即在步长逐次减半的过程中，反复利用复化求积公式进行计算，并同时查看相继两次计算结果的误差是否达到要求，直到所求得的积分值满足要求为止。下面以梯形公式为例第一节 梯形法的递推在变步长的过程中探讨梯形法的计算规律如下：设将积分区间 分为 等分，则一共有 个等分点,abN1N， ，这里用 表示复化梯形法求得的积分值，,kxah01k T其下标 表示等分数。N由余项公式 可知，积分值为(1.9)21()(NbaITf 1()ab再将各子区间分半，使得区间成 等分。此时所得积分近似值记为 ，则N2NT再由余项公式 可知，积分值为(1.9)22()(NbaITf 2()ab假定 在 上变化不大，即有 ，于是得()fx, 12)ff哈尔滨工业大学本科毕业设计（论文）- 14 -，左式也可以写成为24NIT222211()()34NNNITT(2.1)这说明用 作为积分 的近似值时，其误差近似为 。计算过程中2I 2()3NT常用 是否满足作为控制计算精度的条件。如果满足，则取()NT作为 的近似值；如果不满足，则再将区间分半，直到满足要求为止。2I实际计算中的递推公式为 1()baTffb21(21)2NNjabafjN1(2,)k (2.)在给定控制参数 后，当满足 时，则以 作为积分 的近似值。2T2NTI通过类似的推导，还可以得到下面的结论对于辛浦生公式，假定 在 上变化不大，则有(4)fx,ab222211()()5NNNNISSS(2.3)对于科特斯公式，假定 在 上变化不大，则有(6)fx,ab2222311()()634NNNNICCC(2.4)第二节 龙贝格求积公式梯形法的算法简单，单精度低，收敛的速度缓慢。如何提高收敛速度以节省计算量，这是人们极为关心的课题。由此引出了龙贝格公式。由梯形的递推法可以看出，将积分区间等分时，用复化梯形公式计算的结果 作为积2NT哈尔滨工业大学本科毕业设计（论文）- 15 -分 的近似值，其误差近似值为 。可以设想，如果用这个误差作I 21()3NT为 的一种补偿，即将2NT作为积分的近似值，可望提高其精确度。22241()=31NT直接根据复化求积公式，不难验证2224()31NNNTS(2.5)这说明，将区间对分前后两次复化梯形公式的值，按式 作线性组合恰好(2.)等于复化辛浦生公式的值 ，它比 更接近于近似值。NS2NT同样，根据式 用 于 作线性组合会得到比 更精确的值，且(2.3)2 2NS通过直接验证可得22241()51NNNCSS(2.6)再由式 用 与 作线性组合，又可得到比 更精确的值，通常(.4)2 2C记为 ,即NR322241()631NNNC(2.7)式 就称为龙贝格求积公式 2。(2.7)上述用若干个积分近似值推算出更为精确的积分近似值的方法，称为外推方法。我们将序列 , 和 分别称为梯形序列、辛浦生序列、NTSNCR科特斯序列和龙贝格序列。由龙贝格序列当然还可以继续进行外推，得到新的求积序列。但由于在新的求积序列中，其线性组合的系数分别为与 。因此，新的求积序列与前一个序列结果相差不41m0(4)m大。故通常外推到龙贝格序列为止。利用龙贝格序列求积的算法称为龙贝格算法。这种算法具有占用内存少、哈尔滨工业大学本科毕业设计（论文）- 16 -精确度高的优点。因此，成为实际中常用的求积方法。下面给出龙贝格求积算法的计算步骤：第 1 步：算出 和 的值，根据公式 计算 ；()fafb(2.)1T第 2 步：将 分半，算出 后，根据公式 计算 ,再根据,()af(.)2公式 计算 ；(.5)1S第 3 步：再将区间分半，算出 及 的值，并根据公()4bfa(3)4baf式 和 计算出 及 ，再由公式 计算出 ；(2.).4T22.61C第 4 步：将区间再次分半，计算 ， ， ，并由公式 计算 ；8T4S(2.7)1R第 5 步：将区间再次分半，类似上述过程计算 ， ， ， 。16T8S4重复上述过程可计算得到 ， ， ，一直算到龙贝格序列中前后两项之1R24差的绝对值不超过给定的误差限为止。定义 2.16 上述用若干个积分近似值算出更精确的积分近似值的方法，称之为外推法。上述计算步骤也可用表 2.1 表示表 2.1 龙贝格计算步骤表K 梯形序列 2KT辛浦生序列 1KS科特斯序列 2KC龙贝格序列 32KR0 11 2T12 42S1C3 8421R4 16T842哈尔滨工业大学本科毕业设计（论文）- 17 -5 32T16S8C4R 第三节 高斯求积公式前面介绍的 个节点的 Newton -Cotes 求积公式，其特征是节点是等1n距的。这种特点使得求积公式便于构造，复化求积公式易于形成。但同时也限制了公式的精度。 是偶数时，代数精度为 ， 是奇数时，代数精度为1n；我们知道 个节点的插值型求积公式的代数精确度不低于 。设想：能n1n n不能在区间 上适当选择个节点 使插值求积公式的代数精,ab012,nxx度高于 ？答案是肯定的，适当选择节点，可使公式的精度最高达到 ，这就是21n本节所要介绍的高斯型求积公式。例：101()()()fxdAfxf其中， 固定在 ， 可以适当选取，只有两个自由度，得到的是梯01,0,形公式，其代数精度只有 1。如果对求积节点 也进行适当的选取，将有01,x四个自由度，得到如下公式： ()()(3fxdff这个积分公式的代数精确度为 3，这就是高斯型求积公式，上面的求积节点 称为高斯点。1一、高斯型求积公式和高斯点对于含有 个参数， 的求积公式：2n,(0,1)KkAxn，适当选取这 个参数，可以使得数值积分公式的0()()bKkakfxdAfx 2哈尔滨工业大学本科毕业设计（论文）- 18 -代数精确度达到 ，我们称这一类求积公式为高斯型求积公式，称这类求21n积公式的积分节点为高斯点，系数 称为高斯系数。KA定义 2.27 如果 个求积节点的求积公式的代数精度为 ，则这1n21n个求积节点称为高斯点。1n因为高斯求积公式也是插值型求积公式，故有结论 ： 个节点的插值型求积公式的代数精度 ，满足1n d21nd二、高斯点的特征定理 2.17 设 是 个相异点，以这 个点为零点的01,nx 11n次多项式为 则 是高斯点的充1n1()()n nx 01,x要条件是对于任意不超过 次的多项式 ，成立q()0bnaqdx证明 1）必要性。设 为高斯点，则对任意不超过次 的01,nx 21n多项式 ，均有 ，则对于任意次数不超过 次多项式()fx0()()nbkafdAf， 是次数不超过 的多项式，且注意到 ，q1n211()0nkx，故 (0,)k 110()()0nbnknkaqxdqx2）充分性。设对于一切次数不超过 次的多项式 ，成立()qx，又设 是次数不超过 次的多项式，用 去1()0bnaqxd()fx21n1()n除 ，商 ，余 ，即 ，可知 和 均f()r()()qxr()xr是不超过 次的多项式，从而 ，1 ()bbbbnaaaafxdddx又因求积公式是插值多项式的构造导出的，由 的选取，其代数精0,nA确度可以达到 ，而 是次数不超过 次的多项式，因此成立n()rx哈尔滨工业大学本科毕业设计（论文）- 19 -，因 ，所以 ，0()()nbkarxdArx1()0nkx()()kkfxr，故而 也即,1k 0()()bkadAr 0()()nbkafdAx由于 是次数不超过 次的多项式，因此该积分公式的代数精确度至()fx21n少为 ，因而由定义 2.2 知节点 是高斯点。21n01,nx定义 2.3 对于关系 ，我们称之为正交性，即()0baqdx与任意 次多项式正交，而这样的多项式类称为正交多项式。1()nx三、高斯型求积公式的余项 (2)2101()()()!nb bnnk na aRfxdAfxfxd 其中， ， 。积分，可以计算得到,101()n nx234(2)3()!()n nbaRf (2.8)第四节 高斯-勒让德求积公式常用的高斯型求积公式有高斯-勒让德公式、高斯-切比雪夫公式、高斯-拉盖尔公式、高斯-埃尔米特公式，下面着重介绍高斯-勒让德公式一、Legendre 多项式， 21()(1)2!nnndPxx 0,2 (2.9)称式 为勒让德（Legendre）多项式，其具有前面提到的正交性性质，即.9哈尔滨工业大学本科毕业设计（论文）- 20 -对于任意次数不超过 的多项式 ，成立n()qx11()0nqxpd因此，多项式 的零点就是相应的高斯-勒让德求积公式的高斯点。1()px勒让德多项式的前几项如下：， ，0()px21()()dx222211()()(3)!dpxxx，33()5)2441()3503)8p5351()670),8pxx勒让德多项式的性质：性质 1 勒让德多项式的首项系数为 2()!na性质 2 当 为奇数时， 为奇函数；当 为偶数， 为偶函数n()npx()npx性质 3 对一切次数不高于 次的多项式 ，有1q1()0nqxpd二、高斯-勒让德求积公式1）当 时， ，其零点为 ，易得 ， 的高0n1()px0x02A1()fxd斯-勒让德求积公式是 ，其代数精确度为 12()fdf2）当 时， ，其零点为 ， ，此时1n2()31px03x1x， 的高斯-勒让德求积公式是：01,A1fd，其代数精确度为 31()()()3fxdff哈尔滨工业大学本科毕业设计（论文）- 21 -3）当 时， ，其零点为 和 0，2n331()5)2pxx35设高斯-勒让德求积公式为： 10123()()(0)()55fxdAffAf令 ，得到方程 f012 1令 ，得到方程 ()x 2令 ，得到方程 2f 023()53A 3联立方程 、 、 可解得 ， 1 2 3 025918A所以 ，其代数精确度为 51 3()5()8()9fxdfff4）当 时， ，由 ，得3n4241()50)pxx4()0px，故四个零点为 ，即 ，2150x303275， ，137523075x3275x相应的，可以解出， ，高斯-勒让德求积公式为031026A12036A其代数精确度为 71 23()()()()()fxdfxffxAf综上，高斯-勒让德求积公式的积分节点和积分系数表如 表 2.2表 2.2 积分节点和积分系数表哈尔滨工业大学本科毕业设计（论文）- 22 -三、高斯-勒让德求积公式的余项对 阶 求积公式（共有 个求积节点）nGausLegndr1n1(2)21()()!n nRfxd (2.10)其中， ，,101( )n nx为多项式 的零点。01,nx )p积分，可以计算得到： 234(2)()!nnnRf(2.1)四、一般区间上的高斯-勒让德求积公式前面讲解的高斯-勒让德求积公式是计算标准区间 上的定积分，对1,于一般的有界区间 上的定积分，可以通过变量代换转化为区间 的定,ab 1,积分。即在积分 中，令 ，当 时， ；当()fxd2baxtxat时，xb哈尔滨工业大学本科毕业设计（论文）- 23 -；且 。这个变换为线性变换，其反变换为1t2badxt 2batx于是成立 1()()2ba bafftdt然后，我们可以用 区间上的高斯-勒让德求积公式。,第五节 复化两点高斯-勒让德求积公式本节对两点 公式GausLegndr13()()(fxff (2.1)进行推广和复化，得到了新的数值积分公式11 1102233()()()66ni iihhhfxdfxfx (2.13)同时分析了它的积分误差和收敛阶定义 2.48 设 ， 是一种复化求积公式，若当 时成()baIfxdnI 0h立0lim,0nhIC(2.14)则称求积公式 是 阶收敛的。nIp例如，复化的 Trapezoid 公式和复化的 Simpson 公式分别具有二阶和四阶收敛性。定理 2.28 设 ， ，则复化两点高斯 -勒让德求(4),fxCabahn积公式为11102233()()()66nba i ikhhfxdfxfx (2.15)的余项表达式为哈尔滨工业大学本科毕业设计（论文）- 24 -4()5(,213bahRffab (2.16)该方法具有四阶收敛性。考虑积分 ，准确解为 0.74684204210xIed用复化两点高斯-勒让德公式可求得 =0.74680332 ，其绝对误差为I=0.00001960，与复化梯形公式，复化 Simpson 公式以及原两点 Gauss-RLegendre 公式相比，这里构造的复化两点高斯-勒让德求积公式方法在运算量和精度方面都有显著的改善。第六节 本章小结本章主要介绍了精确度比较高的两种数值求积公式，即龙贝格求积公式和高斯型求积公式，对其进行了相应的误差分析。其中高斯型求积公式主要介绍了高斯-勒让德求积公式，并对两点高斯-勒让德求积公式进行了复化。第三章 各种求积公式的 MATLAB 编程实现与应用MATLAB 是由 MathWorks 公式开发的一种主要用于数值计算及可视化图形处理的工程语言，是当今最优秀的科技应用软件之一。它将数值计算、矩阵运算、图形图像处理、信号处理和仿真等诸多强大的功能集成在较易使用的交互计算机环境中，为科学研究、工程应用提供了一种功能强、效率高的编程工具。下面我们将各种求积算法通过 MATLAB 软件编程实现，以下程序均用MATLAB7.0 编写，运行坏境：1、硬件环境 CPU（intel Core i3-2310M,2.1GHz） ，内存（2GB 昱联）,2、软件环境 windows7（32 位）操作系统。以下总共编写了六个算法程序，部分代码参考文献10-13,为了体现程序哈尔滨工业大学本科毕业设计（论文）- 25 -的正确性，以下程序都以 为例进行运算。原积分的精确值为10sinxId10sin.946837xId第一节 几个低次牛顿-科特斯求积公式的 MATLAB 实现先用 M 文件定义一个名为 f1.m 的函数：% i 是要调用第几个被积函数 g(i)，x 是自变量function f=f1(i,x) g(1)=sqrt(x);if x=0g(2)=1;elseg(2)=sin(x)/x;endg(3)=4/(1+x2);f=g(i);程序一：function C,g=NCotes(a,b,n,m)% a，b 分别为积分的上下限；% n 是子区间的个数；% m 是调用上面第几个被积函数；% 当 n=1 时计算梯形公式；当 n=2 时计算辛浦生公式，以此类推；i=n;h=(b-a)/i;z=0;for j=0:ix(j+1)=a+j*h;s=1;if j=0s=s