定积分的数值计算方法[含论文、综述、开题-可编辑]

设计（20 届）定积分的数值计算方法所在学院专业班级信息与计算科学学生姓名学号指导教师职称完成日期年月摘要:数值计算是许多科学与工程计算的核心.定积分的数值计算方法有很多，其中一些常用的计算方法有牛顿-科茨求积公式，梯形求积公式，辛普森求积公式，复合求积公式，龙贝格积分法，高斯求积公式，切比雪夫求积法等．本篇论文主要介绍定积分数值计算的多种方法，并对其中几种做了比较评述，最后给出了梯形求积公式，龙贝格积分法在Matlab环境中的编程实现.关键词：牛顿-科茨求积公式；复合求积公式；高斯求积公式Some numerical methods of definite integralAbstract: Numerical calculation is the core of many science and engineering calculation. There are many numerical calculation methods, including some commonly used numerical methods are Newton – Cotes Quadrature formula, Trapezoidal Quadrature formula, Simpson formula,Composite Quadrature formula, Romberg Quadrature method, Gaussian Quadrature formula, chebyshev Quadrature formula, and so on. This theies mainly introduces Some numerical methods of definite integral and compare several of these methods, finally gives the Trapezoidal Quadrature formula, Romberg Quadrature method in the Matlab environment for programming realize.Key words:Newton – Cotes Quadrature formula; Composite Quadrature formula; Gaussian Quadrature formula目录1 绪论 (1)1.1 问题的背景 (1)2 牛顿-科茨求积公式 (3)2.1 公式的一般形式 (3)2．2 梯形公式 (4)2．3 辛普森公式 (4)3 复化求积公式 (6)3．1 复化梯形求积公式 (6)3.2 复化辛普森求积公式 (6)4 龙贝格求积公式 (8)4．1 递推梯形法则 (8)4．2 龙贝格算法 (8)5 高斯求积公式 (9)5．1 高斯求积公式 (9)5．2 高斯—勒让德（Gauss-Legendre）求积公式 (10)5．3 高斯—埃尔米特求积公式（Gauss-Hermite） (11)5．4 高斯—切比雪夫（Gauss-Chebyshev）求积公式 (12)5.5 递推型高斯求积 (13)6 几种数值积分方法的比较评述和Matlab 实例 (14)致谢 (18)参考文献 (19)1 绪论1.1 问题的背景在科学与工程计算中，经常要计算定积分()()()baI f f x dx a b =-∞≤≤≤∞⎰． （1．1）这个积分的计算似乎很简单，只要求出f 的原函数F 就可以得出积分（1．1）的值，即()()()I f F b F a =-．（1．2） 如果原函数F 非常简单又便于使用，那么式（1．2）就提供了计算起来最快的积分法．但是，积分过程往往将导出新的超越函数，例如，简单积分1dx x ⎰可引出对数函数，它已不是代数函数了；而积分2x edx -⎰，将引出一个无法用有限个代数运算、对数运算或指数运算组合表示的函数．有些积分虽然容易求解，并且原函数仍然是一个初等函数，但可能过于复杂，以致于人们采用（1．2）来计算之前还得三思而行[1]．例如411dx C x =++⎰ ． （1．3） 采用式（1．3）这种“精确”表达式时，所需运算次数是个根本问题．由式（1．3）看出，需计算对数和反正切，因此只能计算到一定的近似程度．因此可以看出，这类表面上是“精确”的方法，实际上也是近似的．因此，我们常常需要探讨一些近似计算定积分的数值方法[2]．通过人们的研究和发现，得出了很多数值计算的方法，比如利用牛顿-科茨求积公式，复合求积公式，龙贝格积分法，高斯求积公式，切比雪夫求积法等来解决定积分的数值计算问题．构造数值积分公式最通常的方法是用积分区间上的n 次插值多项式代替被积函数，由此导出的求积公式称为插值型求积公式．特别在节点分布等距的情形称为牛顿-柯茨公式，例如梯形公式与抛物线公式就是最基本的近似公式．但它们的精度较差．龙贝格算法是在区间逐次分半过程中，对梯形公式的近似值进行加权平均获得准确程度较高的积分近似值的一种方法，它具有公式简练、计算结果准确、使用方便、稳定性好等优点，因此在等距情形宜采用龙贝格求积公式．当用不等距节点进行计算时，常用高斯型求积公式计算，它在节点数目相同情况下，准确程度较高，稳定性好，而且还可以计算无穷积分[3]．各种定积分的数值计算方法的出现和发展，加快和简化了求解定积分的效率和步骤．各种数值积分的方法——牛顿-科茨求积公式，复合求积公式，龙贝格积分法，高斯求积公式等，每种方法都有各自的优缺点，针对不同的积分函数采用不同的方法，所以在实际计算时，要做适当的采取．相信随着理论分析和研究的日益深入，求定积分的数值计算方法将更加简单和完善，为我们的计算带来前所未有的方便，在数学领域也将会更上一层楼．2 牛顿-科茨求积公式2.1 公式的一般形式将积分（1．1）中的积分区间[],a b 分成n 等分，其节点k x 为1,()k x a kh h b a n=+=- (0,1,,)k n =．对于给定的函数f ，在节点k x (0,1,,)k n =上的值()k f x 为已知．那么f 在n+1个节点01,,,n x x x 上的n 次代数插值多项式为00()()n nj n k k j k j j k x x p x f x x x ==≠⎡⎤-⎢⎥=⎢⎥-⎢⎥⎣⎦∑∏． 如果记x a th =+，则上式可以写为00()()n nn k k j j k t j p x f x k j ==≠⎡⎤-⎢⎥=⎢⎥-⎢⎥⎣⎦∑∏．（2．1） 在积分（1．1）中的被积函数f 用其n+1个节点的代数插值多项式()n p x 来代替，可 得 ()()()()bbn n aaI f f x dx I f p x dx =≈=⎰⎰．多项式的积分是容易求出的，因此把上式写为()()()nn n k k I f I f A f x =≈=∑， （2．2）其中 ()00(),n n n k k j j kb a t j A dt b ac n k j=≠--==--∏⎰ (2．3) ()00(1)()!()!n kn n n kj j kct j dt k n k n -=≠-=--∏⎰．(2．4) 公式（2．2）称为牛顿-科茨求积公式[4]或称为等距节点求积公式，k A 称为求积公式系数，()n kc 称为科茨求积系数．牛顿-科茨求积公式的误差估计()n E f ()()n I f I f =-，由下面定理给出 定理2．1 (1) 如果n 为偶数，(2)n f +在[],a b 上连续，则有[]3(2)()(),,n n n n E f c hf a b ηη++=∈， （2．5）其中 201(1)(2)()(2)!n n c t t t t n dt n =---+⎰．(2) 如果n 为奇数，(1)n f+在[],a b 上连续，则有[]2(1)()(),,n n n n E f c h f a b ηη++=∈， （2．6）其中 01(1)(2)()(1)!n n c t t t t n dt n =---+⎰．定义2．1 如果求积公式()()nbk k ak f x dx A f x =≈∑⎰对所有次数不高于n 的代数多项式等式精确成立，但存在n+1次的代数多项式使等式不成立，则称上式求积公式具有n 次代数精度．由定理2．1可知，牛顿-科茨求积公式（2．2）的代数精度至少是n 次，而当n 是偶数时，（2．2）的代数精度可达n+1次．2．2 梯形公式在牛顿-科茨公式（2．2）中，取n=1时(1)(1)011,2c c ==所以有 []1()()()()2b aI f I f f a f b -≈=+．（2．7） 公式（2．7）称为梯形公式[5]，如果用连接(),()a f a 和(),()b f b 的直线来逼近f ，并对这线性函数进行积分可得到1()I f ．再用1()I f 来逼近()I f ． 定理 2．2 若[]2,f Ca b ∈，则梯形公式（2．7）的误差为[]3111()()()()''(),,12E f I f I f b a f a b ηη=-=--∈． 2．3 辛普森公式在牛顿-科茨公式（2．2）中,取n=2，则有220011(1)(2),46c t t dt =--=⎰ 221014(2),26c t t dt =--=⎰222011(1),46c t t dt =-=⎰由此得到2()()()4()().32h a b I f I f f a f f b +⎡⎤≈=++⎢⎥⎣⎦（2．8）其中1()2h b a =-．式（2．8）称为辛普森公式[6]． 定理2．3 若[]4,f C a b ∈，则辛普森公式（2．8）的误差为 []5(4)221()()()(),,90E f I f I f h f a b ηη=-=-∈．3 复化求积公式上面已经给出了计算积分()()baI f f x dx =⎰的3个基本的求积公式：梯形公式，辛普森公式，牛顿-科茨公式，并给出了它们误差的表达式．由这些表达式可知其截断误差依赖于求积区间的长度．若积分区间的长度是小量的话，则这些求积公式的截断误差是该长度的高阶小量．但若积分区间的长度比较大，直接使用这些公式，则精度难以保证．为了提高计算积分的精度，可把积分区间分为若干个小区间，()I f 等于这些小区间上的积分和，然后对每个小区间上的积分应用上述求积公式，并把每个小区间上的结果累加，所得到的求积公式称为复化求积公式[7]．将积分区间[],a b 作n 等分，并记,,0,1,,k b ah x a kh k n n-==+=，于是11()()k kn x x k I f f x dx +-==∑⎰．3．1 复化梯形求积公式如果需要求出一个已知函数()f x 在一个很大区间[],a b 上的积分，那么我们可以把区间分成n 个长度为x h ∆=的小区间，对每一个小区间用梯形法则，然后再把这些小区间上的积分值相加．于是就得到了计算定积分的复化梯形公式[8]：1101210()()(222)22n bi i n n ai h hf x dx f f f f f f f -+-=≈+=+++++∑⎰． (3.1)整体积分误差等于n 个小区间上的积分误差之和：整体误差= []312''()''()''()12n h f f f ξξξ-+++，其中i ξ是第i 个小区间上的某一点．如果''()f x 在区间[],a b 上连续，那么由连续函数的性质可知，在区间[],a b 上存在点ξ使得''()i f ξ的平均值等于()f ξ．于是由于nh b a =-,有整体误差= 322''()''()()1212nh b a f h f O h ξξ--=-=， 局部误差是3()O h ，整体误差是2()O h ．3.2 复化辛普森求积公式对于积分()baf x dx ⎰，将[],a b 等分，每个小区间长度b ah n-=，节点记为 (0,1,2,,)k x a kh k n =+=，第k 个小区间记为[]1,(1,2,,)k k x x k n -=．记[]1,k k x x -的中点为1121()2k k k xx x --=+，则复化辛普森公式[9]为 1112()()()4()()6n bk k ak k h f x dx S h f x f x f x --=⎡⎤≈=++⎢⎥⎣⎦∑⎰．现在要介绍用龙贝格（Romberg ）命名的一个算法，龙贝格首先给出了这种算法的递推形式，假设需要积分()baI f x dx =⎰ （4.1）的近似值．在讨论过程中函数()f x 和区间[],a b 将保持不变．4．1 递推梯形法则设()T n 表示在长度是()/h b a n =-的n 个子区间上积分I 的梯形法则．根据()''()nbai f x dx h f a ih =≈+∑⎰，我们有 00()()''()''()nn n i i b a b a T h f a ih f a i n n ==--=+=+∑∑， （4.2） 这里求和符号中的两撇表示和式中第一项和最后一项减半．4．2 龙贝格算法在龙贝格算法[10]中使用上述公式．设(,0)R n 表示具有2n个子区间的梯形估计，我们有[]1211(0,0)()()()21(,0)(1,0)((21))2n n n i R b a f a f b R n R n hf a i h -=⎧=-+⎪⎪⎨⎪=-++-⎪⎩∑ ， （4.3） 对于一个适度的M 值，计算(0,0),(1,0),(2,0),,(,0)R R R R M ，并且其中没有重复的函数值的计算．在龙贝格算法的其余部分中，还要计算附加值(,)R n m ．所有这些都可以被理解为积分I 的估计．计算出(,0)R M 后，不再需要被积函数f 值的计算．根据公式 []1(,)(,1)(,1)(1,1)41mR n m R n m R n m R n m =-+-----， （4.4） 对于1n ≥和1m ≥构造R 阵列的各列．定理 4.1（龙贝格算法收敛性定理）[10]若[],f C a b ∈，则龙贝格阵列中每一列都收敛于f 的积分．因此，对每个m ，lim (,)()baR n m f x dx =⎰．前面研究的求积公式都是事先确定了n 个节点，然后按使求积公式阶数达到最大的原则选取最佳权．由于自由参数为n 个，所以阶数一般为n-1，但如果节点的位置也自由选择，则自由参数的个数将变为2n ，因此求积公式的阶数可达到2n-1．高斯求积公式就是通过选择最佳的节点和权，使求积公式的阶数最大化．一般地，对每个n ，n 点高斯公式都是唯一的，而且阶数为2n-1．因而，对一定的节点个数，高斯求积公式的精度是最高的．但它的求得比牛顿—柯特斯公式要困难得多．虽然它的节点和权也可由待定系数法确定，但得到的方程是非线性的．5．1 高斯求积公式为说明高斯求积公式，推导区间[]1,1-上的两点公式1112221()()()()()I f f x dx w f x w f x G f -=≈+=⎰，其中的节点1x 、2x 及权1w 、2w 按使求积公式阶数最大化的原则选取．令公式对前四个单项式精确成立，得力矩方程组112111122112221122113331122112,0,2,30.w w dx w x w x xdx w x w x x dx w x w x x dx ----⎧+==⎪⎪+==⎪⎪⎨⎪+==⎪⎪+==⎪⎩⎰⎰⎰⎰这个非线性方程组的一个解为12121,1,x x w w =-===另一个解可通过改变1x ，2x 的符号而得到．这样，两点高斯求积公式为2()(G f f f =-+，阶数为3．另外，高斯求积公式的节点也可以由正交多项式得到．若p 是n 次多项式，且满足()0,0,,1,bk ap x x dx k n ==-⎰则p 与[],a b 区间上所有次数小于n 的多项式正交，容易证明：1． p 的n 个零点都是实的、单的，且位于开区间(,)a b ．2． 区间[],a b 上以p 的零点为节点的n 点插值型求积公式的阶数为2n-1，是唯一的n 点高斯公式[11]．定义5.1[12] 如果1n +个节点的求积公式()()()nbk k ak x f x dx A f x ρ=≈∑⎰（5.1）的代数精度达到21n +，则称式（5.1）为高斯型求积公式，此时称节点k x 为高斯点，系数k A 称为高斯系数．定理5.1[12] 以01,,,n x x x 为高斯点的插值型求积公式具有21n +次代数精确度的充要条件是以这些节点为零点的多项式101()()()()n n x x x x x x x ω+=---，与任意次数不超过n 的多项式()p x 带权()x ρ均在区间[],a b 上正交，即1()()()0bn ax p x x dx ρω+=⎰． （5.2）定理5.2 高斯公式()()nbi i ai f x dx A f x =≈∑⎰（5.3）的求积系数k A 全为正，且 2()(),0,1,,bbk k k aaA l x dx l x dx k n ===⎰⎰． （5.4）定理5.3 对于高斯公式（5.1），其余项为 (22)211()()()()(22)!b n n a R f f x x dx n ξρω++=+⎰ ， （5.5） 其中[]101,,()()()().n n a b x x x x x x x ξω+∈=---5．2 高斯—勒让德（Gauss-Legendre ）求积公式对于任意求积区间[],a b ，通过变换22a b b ax t +-=+，可化为区间[]1,1-，这时11()()222bab a a b b af x dx f t dt --+-=+⎰⎰．因此，不失一般性，可取1,1a b =-=，考查区间[]1,1-上的高斯公式11()()ni i i f x dx A f x -==∑⎰． （5.6）[]423(22)32![](),11,(23)[(22)!]n n n R f f n n ηη++=-≤≤++ (5.7)我们知道，勒让德多项式1211111()(1)2(1)!n n n n n d L x x n dx +++++⎡⎤=-⎣⎦+， （5.8） '12(1)()()i n i n i A n L x L x +=+． (5.9)是区间[]1,1-上的正交多项式，因此，1()n L x +的n+1个零点就是高斯公式（5.6）的n+1个节点．特别地，称1()n L x +的零点为高斯点，形如（5.6）的高斯公式称为高斯—勒让德公式[13]．以上这些公式中的节点和求积系数可查表得到．表5.1 高斯-勒让德求积公式的高斯点与高斯系数5．3 高斯—埃尔米特求积公式（Gauss-Hermite ）Gauss-Hermite 求积公式[14]2()()nx i k ef x dx A f x +∞--∞≈∑⎰， （5.10）其余项为(22)1([]()2(22)!n n n n R f f n η+++=+． ,η-∞<<+∞ （5.11） 其中高斯点为在区间（-∞，+∞）上带权函数2()x x e ω-=正交的埃尔米特多项式221111()(1)()n n x x n n d H x e e dx ++-++=- （5.12）的零点，高斯系数为1i n i m i A += （5.13）表5.2 高斯－埃尔米特求积公式的高斯点与高斯系数5．4 高斯—切比雪夫（Gauss-Chebyshev ）求积公式区间为[]1,1-，权函数()x ρ=Gauss 型求积公式，其节点k x 是Chebyshev多项式1()n T x +的零点，即21cos (0,1,,)2(1)k k x k n n π⎡⎤+==⎢⎥+⎣⎦，而(0,1,,)1k A k n n π==+，于是得到1021cos 12(1)nk k f n n ππ-=⎡⎤+≈⎢⎥++⎣⎦∑⎰， （5.14）称为Gauss-Chebyshev 求积公式[15]，公式的余项为 (22)2(1)2()(),(1,1)2(22)!n n n R f f n πηη++=∈-+ ， （5.15） 这种求积公式可用于计算奇异积分．5.5 递推型高斯求积高斯求积公式不具有递推性：当节点个数一定时，如果自由选择所有的节点和权以达到最高的阶数，则节点个数不同的公式一般没有公共节点，这意味着与一组节点对应的积分值，在用另外一组节点计算积分值时不能被利用．Kronrod 求积公式避免了这种工作量的增加，这类公式是对称的，n 点高斯公式n G 与2n+1个点Kronrod 公式21n K +对应．21n K +节点的约束条件为：以n G 的节点作为21n K +的节点，按求积公式达到最高阶数的要求确定21n K +中剩下的n+1个节点及2n+1个权（其中包括n G 的节点的权）．这样，求积公式的阶数可达到3n+1，而真正2n+1个点高斯公式应该是4n+1阶的，所以精度和效率是一对矛盾．使用两个节点个数不同的求积公式的主要原因是可以用它们的差估计积分近似值的误差．使用Gauss-Kronrod 公式对时，若以21n K +的值作为积分的近似值，则一半基于理论，一半基于经验，可以得到关于误差的保守估计： 1.521(200)n n G K +-．Gauss-Kronrod 公式不仅有效地提供了较高的精度，还给出了可靠误差估计，所以它被认为是最有效的求积公式之一，并且构成了主要软件库中求积程序的基础，特别地，公式715(,)G K 已被普遍使用．6 几种数值积分方法的比较评述和Matlab 实例由于实际问题计算时要考虑到精度要求和计算量，因而不同求积公式各有特点．梯形求积公式和辛普森求积公式虽然计算简单、使用方便，但是精度较差，但对于光滑性较差的被积函数有时比高精度方法更为有效．尤其梯形公式对被积函数是周期函数的效果更为突出．7n>时，Newton-Cotes公式是不稳定的，而复化梯形公式和复化辛普森公式在保留了低阶公式的特点，又能获得较高的精度，因此在实际计算中应用的最为广泛．龙贝格方法的算法简单，易于编程实现．当节点加密提高积分近似程度时，前面计算的结果可以为后面所用，对减少计算两很有好处，并且有比较简单的误差估计，能得到若干积分序列，如果在做收敛性控制时，同时检查各行各列，对于不同性态的函数可以用其中最快的收敛序列来逼近积分．高斯型求积公式的精度较高，特别对计算无穷积分方面是其他方法所不能比拟的．但因为其节点不规则，当节点增加时，前面计算的函数值不能为后面所用，上机计算时还要预先存入不同N 值的节点值和系数表，比较麻烦．对于以上几种方法的实现，本文应用Matlab工具来实现．Matlab是数学、自然科学和工程学的标准指导语言.它有一个由许多实现特殊应用任务的功能集所组成的工具箱.对于本论文的目的而言，我们只需要使用一般的Matlab功能集.[16]例6.1用梯形求积公式求22xπ-⎰．function T=rctrap(fun,a,b,m)%fun 函数，a 积分上限b积分下限m 递归次数n=1;h=b-a; T=zeros(1,m+1); x=a;T(1)=h*(feval(fun,a)+feval(fun,b))/2;for i=1:mh=h/2; n=2*n; s=0;for k=1:n/2x=a+h*(2*k-1); s=s+feval(fun,x);endT(i+1)=T(i)/2+h*s;endT=T(1:m);function y=fun(x)%y = sin(5.*x^2-21);y=exp((-x^.2./2)./(sqrt(2*pi)));>> T=rctrap(@fun,0,pi/2,14), syms tT =Columns 1 through 71.4168 1.3578 1.3313 1.3195 1.3142 1.3119 1.3109 Columns 8 through 141.3105 1.3103 1.3102 1.3102 1.3101 1.3101 1.3101例6.2用龙贝格数值积分求221.5x-⎰．function [RT,R,wugu,h]=romberg(fun,a,b, wucha,m)%fun被积函数a，b积分上下限wucha两次相邻迭代绝对差值m 龙贝格积分表最大行数%RT 龙贝格积分表R 数值积分结果wucha 误差估计h 最小步长n=1;h=b-a; wugu=1; x=a;k=0; RT=zeros(4,4);RT(1,1)=h*(feval(fun,a)+feval (fun,b))/2;while((wugu>wucha)&(k<m)|(k<4))k=k+1; h=h/2; s=0;for j=1:nx=a+h*(2*j-1); s=s+feval(fun,x);endRT(k+1,1)= RT(k,1)/2+h*s; n=2*n;for i=1:kRT(k+1,i+1)=((4^i)*RT(k+1,i)-RT(k,i))/(4^i-1);endwugu=abs(RT(k+1,k)-RT(k+1,k+1));endR=RT(k+1,k+1);function y=fun(x)%y = sin(5.*x^2-21);y=exp((-x^.2./2)./(sqrt(2*pi)));>> syms x>> [RT,R,wugu,h]=romberg(@fun,0,1.5,1.e-8,13)RT =Columns 1 through 81.3541 0 0 0 0 0 0 01.2983 1.2797 0 0 0 0 0 01.2732 1.2648 1.2638 0 0 0 0 01.2620 1.2583 1.2578 1.2577 0 0 0 01.2570 1.2554 1.2552 1.2552 1.2551 0 0 01.2549 1.2541 1.2541 1.2540 1.2540 1.2540 0 01.2539 1.2536 1.2535 1.2535 1.2535 1.2535 1.2535 01.2535 1.2533 1.2533 1.2533 1.2533 1.2533 1.2533 1.25331.2533 1.2532 1.2532 1.2532 1.2532 1.2532 1.2532 1.2532 Column 91.2532R =1.2532wugu =1.4243e-009h =0.0059参考文献[1] 孙志忠，吴宏伟，袁慰平，闻震初．计算方法与实习（第4版）[M]．南京：东南大学出版社， 2009，2：128~129．[2]Micheal T．Heath．张威，贺华，冷爱萍译．科学计算导论（第2版）[M]．北京：清华大学出版社，2005，10：396~297．[3]李桂成．计算方法[M]．北京：电子工业出版社，2005，10：186．[4] 现代应用数学手册编委会．现代应用数学手册——计算与数值分析卷[M]．北京：清华大学出版社，2005，1：163~168．[5] 林成森．数值计算方法（上）[M]．北京：科学出版社，2004，5：220~221．[6]冯康．数值计算方法[M]．北京：国防工业出版社，1978，12：45~47．[7]孙志忠，袁慰平，闻震初．数值分析（第2版）[M]．南京：东南大学出版社，2002，1：191~194．[8] （美）柯蒂斯F．杰拉尔德帕特里克O．惠特莱．应用数值分析（第7版）[M]．北京：机械工业出版社，2006，8： 222~225．[9]夏爱生，胡宝安，孙利民，夏凌辉．复化Simpson数值求积公式的外推算法[J]．军事交通学院学报．2006，第8卷（第1期）： 66~68．[10]（美）David Kincaid, Ward Cheney ．王国荣，俞耀明，徐兆亮译．数值分析（原书第三版）[M]．北京：机械工业出版社，2005，9：400~403．[11]M．T．Heath． Scientific Computing:An Introductory Survey, Sscond Edition[M]．清华大学出版社．英文影印版． 2001，10：351~355．[12]封建湖，车刚明，聂玉．数值分析原理[M]．北京：科学出版社，2001，9：111~114．[13]杨大地，涂光裕．数值分析[M]．重庆：重庆大学出版社,，2006，9：139~142．[14] 曾金平．数值计算方法[M]．长沙：湖南大学出版社，2004，8：253~256．[15]Jeffery J．Leader． Numerical Analysis and Scientific Computation[M]．英文影印本．北京：清华大学出版社，2005，8：342~349．[16]（美）拉克唐瓦尔德（Recktenwald,G）．伍为国，万群，张辉等译．数值方法和MATLAB实现与应用[M]．北京：机械工业出版社，2004，9：422~430．文献综述定积分的数值计算方法一、 前言部分在科学与工程计算中，经常要计算定积分()()().ba I f f x dx ab =-∞≤≤≤∞⎰ （1．1） 这个积分的计算似乎很简单，只要求出f 的原函数F 就可以得出积分（1．1）的值，即()()().I f F b F a =- （1．2）如果原函数F 非常简单又便于使用，那么式（1．2）就提供了计算起来最快的积分法．但是，积分过程往往将导出新的超越函数，例如，简单积分1dx x ⎰可引出对数函数，它已不是代数函数了；而积分2x e dx -⎰，将引出一个无法用有限个代数运算、对数运算或指数运算组合表示的函数．有些积分虽然容易求解，并且原函数仍然是一个初等函数，但可能过于复杂，以致于人们采用（1．2）来计算之前还得三思而行[1]．例如411dx C x =++⎰， （1．3） 采用式（1．3）这种“精确”表达式时，所需运算次数是个根本问题．由式（1．3）看出，需计算对数和反正切，因此只能计算到一定的近似程度．因此可以看出，这类表面上是“精确”的方法，实际上也是近似的．因此，我们常常需要探讨一些近似计算定积分的数值方法[2]．通过人们的研究和发现，得出了很多数值计算的方法，比如利用牛顿-科茨求积公式，复合求积公式，龙贝格积分法，高斯求积公式，切比雪夫求积法等来解决定积分的数值计算问题．构造数值积分公式最通常的方法是用积分区间上的n 次插值多项式代替被积函数，由此导出的求积公式称为插值型求积公式．特别在节点分布等距的情形称为牛顿-柯茨公式，例如梯形公式与抛物线公式就是最基本的近似公式．但它们的精度较差．龙贝格算法是在区间逐次分半过程中，对梯形公式的近似值进行加权平均获得准确程度较高的积分近似值的一种方法，它具有公式简练、计算结果准确、使用方便、稳定性好等优点，因此在等距情形宜采用龙贝格求积公式．当用不等距节点进行计算时，常用高斯型求积公式计算，它在节点数目相同情况下，准确程度较高，稳定性好，而且还可以计算无穷积分[3]．二、 主题部分2．1 牛顿-科茨求积公式[4]2．1．1 公式的一般形式[4]将积分（1．1）中的积分区间[],a b 分成n 等分，其节点k x 为1,()k x a kh h b a n =+=- (0,1,,)k n =．对于给定的函数f ，在节点k x (0,1,,)k n =上的值()k f x 为已知．那么f 在n+1个节点01,,,n x x x 上的n 次代数插值多项式为00()().n nj n k k j k j j k x x p x f x x x ==≠⎡⎤-⎢⎥=⎢⎥-⎢⎥⎣⎦∑∏ 如果记x a th =+，则上式可以写为00()().n nn k k j j k t j p x f x k j ==≠⎡⎤-⎢⎥=⎢⎥-⎢⎥⎣⎦∑∏ （2．1） 在积分（1．1）中的被积函数f 用其n+1个节点的代数插值多项式()n p x 来代替，可得 ()()()()bbn n a a I f f x dx I f p x dx =≈=⎰⎰． 多项式的积分是容易求出的，因此把上式写为()()()nn n k k I f I f A f x =≈=∑， （2．2）其中 ()00(),n n n k k j j kb a t j A dt b ac n k j=≠--==--∏⎰ (2．3)()00(1)().!()!n kn n n kj j k c t j dt k n k n -=≠-=--∏⎰ (2．4) 公式（2．2）称为牛顿-科茨求积公式或称为等距节点求积公式，k A 称为求积公式系数，()n k c 称为科茨求积系数．牛顿-科茨求积公式的误差估计()n E f ()()n I f I f =-，由下面定理给出定理2．1 (1) 如果n 为偶数，(2)n f +在[],a b 上连续，则有[]3(2)()(),,n n n n E f c h f a b ηη++=∈， （2．5）其中 201(1)(2)()(2)!n n c t t t t n dt n =---+⎰．(2) 如果n 为奇数，(1)n f +在[],a b 上连续，则有[]2(1)()(),,n n n n E f c h f a b ηη++=∈， （2．6）其中 01(1)(2)()(1)!n n c t t t t n dt n =---+⎰．定义2．1 如果求积公式0()()nbk k a k f x dx A f x =≈∑⎰对所有次数不高于n 的代数多项式等式精确成立，但存在n+1次的代数多项式使等式不成立，则称上式求积公式具有n 次代数精度．由定理2．1可知，牛顿-科茨求积公式（2．2）的代数精度至少是n 次，而当n 是偶数时，（2．2）的代数精度可达n+1次．2．1．2 梯形公式[5]在牛顿-科茨公式（2．2）中，取n=1时(1)(1)011,2c c ==所以有 []1()()()().2b a I f I f f a f b -≈=+ （2．7） 公式（2．7）称为梯形公式，如果用连接(),()a f a 和(),()b f b 的直线来逼近f ，并对这线性函数。