湖南大学科学计算导论数值积分

1、第第7 7章章 数值积分数值积分/* Numerical Integration */引言引言在科学研究与工程技术应用中，经常需要计算定积分在科学研究与工程技术应用中，经常需要计算定积分( )baf x dx而在实际应用中常出现如下几种情况：而在实际应用中常出现如下几种情况： 被积函数的原函数不能用初等函数表示，如被积函数的原函数不能用初等函数表示，如 ， ；10sin xdxx210 xedx 的原函数的表达式过于复杂；的原函数的表达式过于复杂；( )f x 是以表格的形式给出的。是以表格的形式给出的。( )f x 对于上述情况，都需要建立定积分的近似计算方法，即所谓对于上述情况，都需要建立

2、定积分的近似计算方法，即所谓的积分的数值计算方法。的积分的数值计算方法。7.1 机械求积公式机械求积公式 /* mechanical quadrature fomula*/ 基本思想基本思想ab( )f x 几何观点几何观点直线段直线段(1) ),0,1yfab近似代替近似代替( )( , )yf xxa b(a+b)/2当当 时时12 ( /2/2)yf ab( )() ( /2/2)baf x dxba f ab中矩形公式中矩形公式左矩形公式？右矩形公式？梯形公式？左矩形公式？右矩形公式？梯形公式？若用过点若用过点 的直线段的直线段( ,( ),( ,( )A a f aB b f b(

3、)( )( )(), , f bf ayf axaxa bba近似代替曲线段近似代替曲线段 则可得到梯形积分公式则可得到梯形积分公式( )( , )yf xxa b( )( ( )( )2babaf x dxf af b类似过点类似过点 的抛物线段的抛物线段( ,( ),(,(),( ,( )22ababA a f aCfB b f b( )( ( )4 ()( )62babaabf x dxf aff b称为称为Simpson公式公式2 , 0,2),1(2)(020tabhttftfafy( )( /)2(2)baf x dxfbaba( )( (2)baf x dxfbaaf b( )(

4、 ( )46()( )2baabf x dxf abbaff以上几个公式可以写成一般形式以上几个公式可以写成一般形式0( )()nbkkakf x dxA f xkxkA求积节点求积节点求积系数求积系数(与被积函数有关与被积函数有关?) 此求积公式具有通用性，其特点是将积分求值问题转化为求此求积公式具有通用性，其特点是将积分求值问题转化为求函数值的问题，避免了求原函数的困难。而这类数值积分方法函数值的问题，避免了求原函数的困难。而这类数值积分方法通常称为机械求积。通常称为机械求积。 待定系数法待定系数法 根据根据Weierstrass多项式逼近定理多项式逼近定理p209th6.6p209th6

5、.6可知，对于闭区间可知，对于闭区间上的连续函数，都可以用多项式去一致地逼近它。换句话说，上的连续函数，都可以用多项式去一致地逼近它。换句话说，任一连续函数都可以用多项式作为它的最简单的近似函数，一般任一连续函数都可以用多项式作为它的最简单的近似函数，一般说来，多项式的次数取得越高，用它们来近似连续函数的程度也说来，多项式的次数取得越高，用它们来近似连续函数的程度也就越高。这自然使我们想到利用多项式的次数去规定求积公式的就越高。这自然使我们想到利用多项式的次数去规定求积公式的精确性程度（所谓代数精度）。精确性程度（所谓代数精度）。 定义定义若某个求积公式对若某个求积公式对任意任意 k n 阶阶

6、的多项式成立，且对的多项式成立，且对某某个个 n+1 阶多项式公式不精确成立，则称此求积公式的阶多项式公式不精确成立，则称此求积公式的代数精度代数精度为为n或称此公式具有或称此公式具有n次代数精度次代数精度 。例：例：对于对于a, b上上1次插值，有次插值，有)()()(1bfafxLabaxbabx )()()(2221bfafdxxfAAabbaab 考察其代数精度。考察其代数精度。f(x)abf(a)f(b)梯形公式梯形公式/* trapezoidal rule*/解：解：逐次检查公式是否精确成立逐次检查公式是否精确成立代入代入 P0 = 1： baabdx111 2 ab=代入代入 P

7、1 = x ：=代入代入 P2 = x2 ： 222abbadxx 2baab 3233abbadxx 222baab 代数精度代数精度 = 1一般地，要使求积公式具有一般地，要使求积公式具有 次代数精度，只要令它对于次代数精度，只要令它对于m( )(0,1,)if xx im都能精确成立，即要求都能精确成立，即要求110,0,1,1iinikkkbaA ximi当节点当节点 给定且互异时，有可能够通过上式确定给定且互异时，有可能够通过上式确定 而这样的以代数精度为标准构造求积公式的方法称为待定系数而这样的以代数精度为标准构造求积公式的方法称为待定系数法。见法。见p237例例7.1(0,1,

8、)kxknkA 插值型求积公式插值型求积公式 可以从不同的角度出发通过各种途径来构数值求积公式。可以从不同的角度出发通过各种途径来构数值求积公式。但常用的一个方法是，利用插值多项式来构造数值求积公式。但常用的一个方法是，利用插值多项式来构造数值求积公式。具体做法如下：具体做法如下： 在积分区间在积分区间 a，b 上取一组点上取一组点作f (x)的的n 次插值多项式次插值多项式：为为n 次插值基函数。用次插值基函数。用 Ln(x) 近似代替被积函数近似代替被积函数f (x)，则得，则得其中其中若记若记(7.1)(7.2)得数值求积公式得数值求积公式(7.3)形如形如（7.3）的求积公式根据之前的

9、内容可以称为的求积公式根据之前的内容可以称为机械求积公式。机械求积公式。 若求积公式若求积公式（7.3）中的求积系数中的求积系数 Ak 是由是由（7.2）确定确定的，则称该求积公式为的，则称该求积公式为插值型求积公式插值型求积公式。本章主要讨论插值型求积公式。本章主要讨论插值型求积公式。 求积公式的余项求积公式的余项 与由某求积公式给出的近似之差，称为该求积公式的余项，与由某求积公式给出的近似之差，称为该求积公式的余项，记作记作R(f ) 。例如，求积公式例如，求积公式（7.3）的余项为的余项为积分的真值积分的真值（7.4） 如果求积公式如果求积公式（7.3）是插值型的，则由上知是插值型的，则

10、由上知于是，由插值余项公式得于是，由插值余项公式得其中其中含有含有n +1个节点个节点 xk (k=0,1,n )的求积公式的求积公式（7.3）是插是插值型（精确成立）的充要条件是该公式值型（精确成立）的充要条件是该公式的代数精度至少的代数精度至少为为n 。 定理定理分析：分析：由于对次数不超过由于对次数不超过n次的任意多项式次的任意多项式 来说，来说， 所以所以n+1个结点的插值型求积公式的代数精度大于等于个结点的插值型求积公式的代数精度大于等于n；反之，反之，容易证明，代数精度大于等于容易证明，代数精度大于等于n的的n+1个结点的求积公式一定是插个结点的求积公式一定是插值型求积公式。事实上

11、，因值型求积公式。事实上，因LagrangeLagrange插值基函数插值基函数 ，而，而因求积公式的代数精度大于等于因求积公式的代数精度大于等于n，因此该求积公式对，因此该求积公式对 必精必精确成立，即确成立，即( )f x( )0nRf 1knlxP klx 1,1, .nbkj kjkajlx dxA lxA kn求积公式的收敛性与稳定性1.收敛性定义：2稳定性定义（定义7.3）3.结论：若Ak0,则求积公式稳定。（证明略）nkbakkhndxxfxfA00)()(limn0kkk| )(A|,|)f(x|, 0, 0kkkfxff就有只要思思路路利用利用插值多项式插值多项式 则积分易算

12、。则积分易算。)()(xfxPn 在在a, b上取上取 a x0 x1 xn b，做，做 f 的的 n 次插值次插值多项式多项式 ，即得到，即得到 nkkknxlxfxL0)()()( babaknkkdxxlxfdxxf)()()(0Ak bakjxxxxkdxAjkj)()(由由 决定，决定，与与 无关。无关。节点节点 f (x)插值型求积公式插值型求积公式/*interpolatory quadrature*/0(1)0( )( )() ( )( )( )()()(1)!nnbkkakbbnnaannbxkakRff x dxA f xf xLx dxRx dxfxxdxn 误差误差 7

13、.2 Newton-Cotes公式公式 当节点当节点等距分布等距分布时：时：ninabhhiaxi,., 1, 0, dxxxxxAnxxijjiji 0)()( njiinnjidtjtininabdthhjihjt00)()!( !) 1()()()(令令htax Cotes系数系数)(niC注：注：Cotes 系数仅取决于系数仅取决于 n 和和 i，可查表得到。与可查表得到。与 f (x) 及区及区间间a, b均无关。均无关。 ( )()niiAba C( )0( )()()nbniiaif x dxbaCf x（7.5）21,21)1(1)1(0 CCn = 1:)()(2)(bfaf

14、abdxxfba Trapezoidal RuledxbxaxffRbax)(!2)( /* 令令 x = a+th, h = b a, 用中用中值定理值定理 */1, , )(1213abhbafh 代数精度代数精度 = 1n = 2:61,32,61)2(2)2(1)2(0 CCC)()(4)(6)(2bffafabdxxfbaba Simpsons Rule代数精度代数精度 = 32,),(,)(901)4(5abhbafhfR n = 3: Simpsons 3/8-Rule, 代数精度代数精度 = 3,)(803)5(5 fhfR n = 4: Cotes Rule, 代数精度代数精

15、度 = 5,)(9458)6(7 fhfR n 为为偶数阶偶数阶的的Newton-Cotes 公式至少有公式至少有 n+1 次代数精度。次代数精度。n123456Cn)(0Cn)(1Cn)(2Cn)(3Cn)(4Cn)(5Cn)(62161819072881984041216483451696253599017962528098145161442510534618315214425280928819935984041为了便于应用，部分柯特斯系数列见下表为了便于应用，部分柯特斯系数列见下表 利用这张柯特斯系数表由利用这张柯特斯系数表由（7.5）可以直接写出当可以直接写出当n =1,2, 6 时的

16、牛顿时的牛顿- -柯特斯公式柯特斯公式。例如，当例如，当 n =1时有时有两两点公式点公式(7.6)当当n =2时有时有三点公式三点公式(7.7)当当n =4 时有时有五点公式五点公式(7.8)其中其中., 1 , 0nk求积公式求积公式（7.6）就是梯形公式梯形公式。当当n =3时有时有4点公式点公式)()32(3)32(3)(8)(bfabfabfafabdxxfba求积公式求积公式（7.7）称为称为辛普生辛普生（Simpson）公式公式。其几何其几何意义就是通过意义就是通过A, B, C 三点的抛物线三点的抛物线 y =L2(x) 围成的曲边围成的曲边梯形面积近似地代替原曲边梯形面积梯形

17、面积近似地代替原曲边梯形面积（见图图）。 因此，因此，求积公式求积公式（7.7）又名又名抛物线公式。抛物线公式。求积公式求积公式（7.8）称为称为柯特斯公式柯特斯公式。 梯形公式、辛普生公式和柯特斯公式，是三个最基梯形公式、辛普生公式和柯特斯公式，是三个最基本、最常用的等距节点下的求积公式。本、最常用的等距节点下的求积公式。 下述定理给出了下述定理给出了这些求积公式的余项为：这些求积公式的余项为：yx0a2)(ba bACB)(xfy 误差估计定理定理7.3：n 为为偶数阶偶数阶的的Newton-Cotes 公式至少有公式至少有 n+1 次代数精度。次代数精度。证明：只证对证明：只证对f(x)

18、=xn+1,R(f)=0.nninbanbanndtithdxxdxxnf00211)1(0)()()()!1()(R(f)因为：022241)2)(1t(2n23420dttt时，例 若若 f (x) 在在 a,b 上连续，则梯形公式上连续，则梯形公式（7.6）的余项为的余项为 定理定理3 33 2 1)(12)(2)(3)(2)()()(2)()(2)()()(21)(abfababfdxbaaxaxfdxbxaxfdxbxaxxffRbababa证明：由插值余项定理知由插值余项定理知等式两边积分得等式两边积分得由于由于f(x) C2a,b ，且，且(x-a)(x-b)在在a,b上非正上非

19、正 ( (不变不变号号),),故根据积分中值定理知故根据积分中值定理知, ,至少存在一点至少存在一点 a,b, ,使使),()()(21)()(1babxaxfxLxf 1( )( )()()d2bTaRffxa xbx311( )( )()()()( )212bTaRffxa xb dxbaf 定理证明定理证明若若 f（4）(x) 在在 a,b 上连续，则辛普生公式上连续，则辛普生公式（7.7）的余项为：的余项为：定理定理)(2880)()()2)(!4)()()2)()(f!41)()()2)(4!)(fP(x)-f(x)P(x),f(x)33.7)4(52)4(2(4)2(4)fabdx

20、bxbaxaxfdxbxbaxaxfRbxbaxaxbaba则的逼近多项式为：设。知其代数精度为至少证明：由定理抛物线求积公式的代数精度为抛物线求积公式的代数精度为3,为此构造三次多项式为此构造三次多项式P3(x)，满足满足 P3(a)=f(a), 则则等式两边从等式两边从a到到b积分得积分得由于由于P3(x)是三次多项式，故抛物线求积公式对它准确成立，是三次多项式，故抛物线求积公式对它准确成立，即即)2()2(),()(),2()2(333bafbaPbfbPbafbaP),()()2)()(! 41)()(2)4(3babxbaxaxfxPxf(4)231( )d( )d( )()() (

21、)d4!2bbbaaaabf xxP xxfxa xxbx这样这样由于由于f(x) C4a,b ，且，且 在在a,b上非上非正正( (不变号不变号),),故根据积分中值定理知故根据积分中值定理知, ,至少存在一点至少存在一点 a,b, ,使使)()2(4)(6)()2(4)(6d3333bfbafafabbPbaPaPabbax(x)P(4)21( )( )()() ()d4!2bsaabRffxa xxbx)()2)(2bxbaxax5(4)( )( )d( )4 ()( )62()( ) , 2880bsabaabRff xxf aff bbafa b 解 (1)将f(x)在a处展开，得两

22、边在a, b上积分，得由于x-a在a, b上不变号，故有a,b, 使从而得 21( )d() ( )( )()2baf xxba f afba31( )d() ()( )()224baabf xxba ffba( )( )( )(),( , )f xf afxaa x( )d( )d( )()d() ( )( )()bbbaaabaf xxf axfxaxba f afxa dx( )d() ( )( )()dbbaaf xxba f afxax21( )d() ( )( )() , , 2baf xxba f afbaa b练习：练习：推导下列矩形求积公式推导下列矩形求积公式(1)(2)将f

23、(x)在(a+b)/2处展开，得两边积分，得 由于 在a, b上不变号，故有a,b, 使 21( )()()()( )()( , )22222ababababf xffxfxa b21( )d() ()()()d( )() d22222bbbaaaababababf x xba ffxxfxx21( )d() ()( )() d222bbaaababf xxba ffxx31() ()( )() , , 224abba ffbaa b2)2(bax 7.3 复合求积公式复合求积公式 /* Composite Quadrature */ 从从Newton-Cotes求积公式的余项可知，被积函数所

24、用的插求积公式的余项可知，被积函数所用的插值多项式的次数越高，相应的求积公式的精度也越高，而高次值多项式的次数越高，相应的求积公式的精度也越高，而高次插值多项式稳定性很差，从而导致高次插值求积公式很不稳定。插值多项式稳定性很差，从而导致高次插值求积公式很不稳定。因此在实际应用中，为了既能提高结果的精度，又使算法简便因此在实际应用中，为了既能提高结果的精度，又使算法简便且易在电子计算机上实现，往往采用复合求积的方法。所谓复且易在电子计算机上实现，往往采用复合求积的方法。所谓复合求积，就是先将积分区间分成几个小区间，并在每个小区间合求积，就是先将积分区间分成几个小区间，并在每个小区间上用低阶牛顿上

25、用低阶牛顿柯特斯公式计算积分的近似值，然后对这些近柯特斯公式计算积分的近似值，然后对这些近似值求和，从而得到所求积分的近似值。由此得到的一些具有似值求和，从而得到所求积分的近似值。由此得到的一些具有更大实用价值的数值求积公式，统称为更大实用价值的数值求积公式，统称为复合求积公式。复合求积公式。例如例如，先将区间先将区间 a，b n 等分，记分点为等分，记分点为其中其中称为称为步长步长，然后在每个小区间然后在每个小区间 xi-1，xi 上应用梯形公式即上应用梯形公式即余项为：余项为：)(max*)()(min*)(12)( 10 103fnffnfhfRknkknk就可导出就可导出复复合合梯形公

26、式梯形公式若将所得积分近似值记成若将所得积分近似值记成Tn ，并注意到并注意到x0=a，xn=b，则则上式即为上式即为)(12)()( 23fnabfR余项为余项为仿上，可得仿上，可得复复合合辛普生公式辛普生公式也可以像教材上介绍的将区间也可以像教材上介绍的将区间a,b分为分为2n等分，得到等分，得到辛普生公式辛普生公式1010212)(2)(4)()(3)(nknkkkbaxfxfbfafhdxxf例例7.4.利用梯形公式和利用梯形公式和辛普生公式求辛普生公式求dxdxex1011xsinx2和和和复复合合柯特斯公式柯特斯公式例例7.4。解：。解：4835. 1)0()5 . 0()5 .

27、0(4) 1 () 1(35 . 0112fffffdxex若若 f (x) 在积分区间在积分区间 a，b 上分别具有二阶、四阶和六阶上分别具有二阶、四阶和六阶连续导数连续导数, ,则复则复合合求积公式求积公式的余项分别为的余项分别为)(12)( 2fhabTdxxfban)()2(180)()4(4fhabSdxxfban)()2(945)(2)()6(6fhabCdxxfban定理定理证明证明 只对复只对复合合梯形公式证明余项公式。梯形公式证明余项公式。由于由于 在在 上连续，故由定理知，对每个小区间上积上连续，故由定理知，对每个小区间上积分分 使用梯形公式时，所得近似值的误差为使用梯形公

28、式时，所得近似值的误差为故故 即即 )(xfba,kkxxdxxf1)(),)(12113kkkkxxfhbannfffhTdxxf)()()(121)(213bannfffnhabTdxxf)()()(112)(212)( 12)()(1223 2fnabfhab由上式再根据已知条件容易判定复由上式再根据已知条件容易判定复合合梯形公式是收敛的。梯形公式是收敛的。而事实上只要而事实上只要 就可得到收敛性。就可得到收敛性。( ) , f xC a b因为因为由介值定理知，在由介值定理知，在 中必有点中必有点 ，使，使故余项公式成立。故余项公式成立。 ba,)()()(1)(21nfffnf)()

29、()(1)(21,minnbaxfffnxf)( max,xfbax)(12)( 2fhabTdxxfban复化复化Simpson公式证明余项公式公式证明余项公式,见见p345公式公式7.26 7.4 变步长求积公式变步长求积公式 /* Composite Quadrature */ 复合求积公式的截断误差随着复合求积公式的截断误差随着n的增大而减小，但是对于一的增大而减小，但是对于一个给定的积分，选定了某种积分方法之后，如何选择适当的个给定的积分，选定了某种积分方法之后，如何选择适当的n，使得计算结果达到预先选定的精度？由误差的要求可以定出使得计算结果达到预先选定的精度？由误差的要求可以定出n ，但事先计算时很难估计：应边算边估计进而但事先计算时很难估计：应边算边估计进而加密。加密。 接下来以接下来以复合梯形求积公式复合梯形求积公式为例，探讨如何避免事先误差为例，探讨如何避免事先误差估计中的高阶导数的计算，以及事后误差估计中步长的选择问估计中的高阶导数的计算，以及事后误差估计中步长的选择问题。题。研究对象：研究对象： 上的定积分上的定积分 。 , a b( )baf x dx取取2n等分计算等分计算2nT2222n( ,)() ( , )12( )2nnbaR ffa bhT 2n( ,)() ( , )12nnbaR ffa