第五章数值微分与数值积分

1、1第五章第五章 数值微分与数值积分数值微分与数值积分2数值微分数值微分hhxfhxfhhxfxfhxfhxfxfhhh2)()(lim)()(lim)()(lim)( 0001. 函数函数f(x)以离散点列给出时，而要求我们给出导数值，以离散点列给出时，而要求我们给出导数值，2. 函数函数f(x)过于复杂过于复杂这两种情况都要求我们用数值的方法求函数的导数值这两种情况都要求我们用数值的方法求函数的导数值微积分中，关于导数的定义如下：微积分中，关于导数的定义如下：自然，而又简单的方法就是，取极限的近似值，即差商自然，而又简单的方法就是，取极限的近似值，即差商3hxfhxfxf)()()( 000

2、由由Taylor展开展开hxxfhxhfxfhxf002000),( ! 2)( )()(因此，有误差因此，有误差)()( ! 2)()()( )(000hOfhhxfhxfxfxR向前差商向前差商4hhxfxfxf)()()( 000由由Taylor展开展开hxxfhxhfxfhxf002000),( ! 2)( )()(因此，有误差因此，有误差)()( ! 2)()()( )(000hOfhhhxfxfxfxR向后差商向后差商5hhxfhxfxf2)()()( 000由由Taylor展开展开23000010102300002020()()()()( ),2!3!()()()()(),2!3

3、!hhf xhf xhfxfxfxxhhhf xhf xhfxfxfxhx因此，有误差因此，有误差)()( 6)( )( 12 2)()()( )(22212000hOfhffhhhxfhxfxfxR中心差商中心差商6 由误差表达式，由误差表达式，h越小，误差越小，但同时舍入误差增越小，误差越小，但同时舍入误差增大，所以，有个最佳步长大，所以，有个最佳步长我们可以用事后误差估计的方法来确定我们可以用事后误差估计的方法来确定设设D(h),D(h/2)分别为步长为分别为步长为h,h/2的差商公式。则的差商公式。则)2()(hDhD时的步长时的步长h/2就是合适的步长就是合适的步长( )( )( )

4、( )( /2)( /2)fxD hO hfxD hO h( )( )( )2( )( /2)( /2)fxD hO hfxD hO h( )( )2( )2 ( /2)fxD hfxD h( )( /2)( )( /2)fxD hD hD h7f(x)=exp(x)h hf(1.15)f(1.15)R(xR(x) )h hf(1.15)f(1.15) R(xR(x) )0.100.103.16303.1630-0.0048-0.00480.050.053.15903.1590-0.0008-0.00080.090.093.16223.1622-0.0040-0.00400.040.043.1

5、5883.1588-0.0006-0.00060.080.083.16133.1613-0.0031-0.00310.030.033.15833.1583-0.0001-0.00010.070.073.16073.1607-0.0025-0.00250.020.023.15753.1575-0.0007-0.00070.060.063.16003.1600-0.0018-0.00180.010.013.15503.1550-0.0032-0.0032例：8 插值是建立逼近函数的手段，用以研究原函数的性质。插值是建立逼近函数的手段，用以研究原函数的性质。因此，可以用插值函数的导数近似为原函数的导

6、数因此，可以用插值函数的导数近似为原函数的导数)()()()(xLxfknk误差误差)()()()!1()()()1(xLxfxnfxRnnnn)()!1()()()1()(xnfdxdxRnnkkkn插值型数值微分插值型数值微分9给定点列给定点列10( ,( )iiixf x且且10 xxh，求，求10(),()fxfx解：解：01101()()( )()()xxxxL xf xf xhh-例例1101()() ( )f xf xLxh余项余项10010()()() ()( )2f xf xhfxLxfh10111()()() ()( )2f xf xhfxLxfhTaylor展开分析，可以

7、知道，它们都是展开分析，可以知道，它们都是( )O h称为称为两点公式两点公式10给定点列给定点列20)(,(iiixfx且且hxxxx0112，求，求)( ),( ),( 012xfxfxf解：解：)(2)()()()(2)()(2210122002212xfhxxxxxfhxxxxxfhxxxxxL)(2)()()()(2)()(2210122002212xfhxxxxxfhxxxxxfhxxxxxL例例21120200121() ()3 ()4 ()()( )23hfxLxf xf xf xfh2121021() ()()()( )26hfxLxf xf xfh22220121() ()

8、()4 ()3 ()( )23hfxLxf xf xf xfhTaylor展开分析，可以知道，它们都是展开分析，可以知道，它们都是)(2hO称为称为三点公式三点公式2(4)0200121221() ()()2 ()()( )()6hfxLxf xf xf xhffh 2(4)11201221() ()()2 ()()( )12hfxLxf xf xf xfh2(4)2220121221() ()()2 ()()( )()6hfxLxf xf xf xhffh12例例3已知函数已知函数y=ex的下列数值的下列数值x2.52.62.72.82.9y=ex12.182513.463714.87971

9、6.444618.1741试用两点数值微分公式和三点数值微分公式计算试用两点数值微分公式和三点数值微分公式计算x=2.7处函处函数的一、二阶导数值。数的一、二阶导数值。解：取解：取h=0.2，x0=2.5,x1=2.7,x2=2.9，计算结果如下：，计算结果如下：114.8797 12.1825(2.7) (2.7)13.4860.2fL（两点）218.1741 12.1825(2.7) (2.7)14.9792 0.2fL（三点）121(2.7) (2.7)12.18252 14.8797 18.174114.930(0.2)fL 13 取取h=0.1，x0=2.6,x1=2.7,x2=2.

10、8，计算结果如下：，计算结果如下：114.8797 13.4637(2.7) (2.7)14.1600.1fL（两点）216.4446 13.4637(2.7) (2.7)14.9052 0.1fL（三点）121(2.7) (2.7)13.46372 14.8797 16.444614.890(0.1)fL f(2.7)与与f(2.7)的真值都是的真值都是14.87973，上面计算表明，上面计算表明，三点公式三点公式 比两点公式准确，步长比两点公式准确，步长h越小结果也越准确。越小结果也越准确。14)()()(aFbFdxxfba关于积分，有关于积分，有Newton-Leibniz公式公式但是

11、，在很多情况下，还是要数值积分：但是，在很多情况下，还是要数值积分：1、函数、函数f(x)的积分存在，但的积分存在，但F(x)不能用初等函数表示，不能用初等函数表示，例如：例如：2、被积函数、被积函数f(x)表达式未知，表达式未知，f(x)是用表格形式给出的。是用表格形式给出的。3、用微积分中的换元积分、分部积分等方法能积出、用微积分中的换元积分、分部积分等方法能积出f(x)的原函数，但过程复杂（比如复杂的有理函数积分）。的原函数，但过程复杂（比如复杂的有理函数积分）。11200sin,sinxdxx dxx数值积分数值积分15积分中值定理告诉我们积分中值定理告诉我们 若若f(x)在区间在区间

12、a,b上连续，则在上连续，则在a,b上至少存在一点上至少存在一点 ， 满足满足( )() ( )baf x dxba f几何意义：几何意义：b( )yfxyxOa16梯形公式：梯形公式：( )( )2baTf af b矩形公式：矩形公式：()2abRba f几何意义：几何意义：（5.1）（5.2）yxOab( )yf x17更一般地，定义数值积分是离散点上的函数值的线性组合更一般地，定义数值积分是离散点上的函数值的线性组合0()()()nbiinaifx dxa fxIf称为求积系数，与称为求积系数，与f(x)无关，与积分区间和求积节点无关，与积分区间和求积节点xi有关有关这类数值积分方法通常

13、称为这类数值积分方法通常称为机械求积机械求积。两个问题：两个问题：1、系数、系数ai如何选取，即选取原则如何选取，即选取原则2、若节点可以自由选取，取什么点好？、若节点可以自由选取，取什么点好？（5.3）18求积公式的求积公式的余项余项（截断误差）：（截断误差）： 0( )( )nbniiaiRff x dxa f x19代数精度代数精度)()(0iniinxfafI为数值积分，为数值积分，badxxffI)()(为积分，则称数值为积分，则称数值积分有积分有m阶代数精度阶代数精度是指：是指：11()(),0,;()()iimmnnIxI xim IxI x定义定义对任意次数不高于对任意次数不高

14、于m次的多项式次的多项式f(x)，数值积分没有误差，数值积分没有误差 可以验证，梯形公式（可以验证，梯形公式（5.1）和矩形公式（）和矩形公式（5.2）均具）均具有一次代数精度。有一次代数精度。20例例4、确定求积公式、确定求积公式的代数精度。的代数精度。解：引入记号，记求积公式解：引入记号，记求积公式(5.3）的左边为）的左边为I(f)，当，当f(x)分别取分别取1,x,x2,x3,x4时，计算如下：时，计算如下：1111( )33f x dxff11111222133133314414441(1)12,(1)1 12.11( )0,( )0.332112(),().333311()0,()

15、0.332112(),().5933nnnnnIdxII xxdxIxI xx dxIxI xx dxIxI xx dxIx 21 可见，当可见，当f(x)分别为分别为1,x,x2,x3时，求积公式（时，求积公式（5.3）的）的左右两边精确相等，而当左右两边精确相等，而当f(x)= x4 时，左右两边不等，所时，左右两边不等，所以该求积公式具有以该求积公式具有3次代数精度。次代数精度。显然，一个求积公式的代数精度越高，用它来进行积分的显然，一个求积公式的代数精度越高，用它来进行积分的近似计算越具有好的实际计算意义。近似计算越具有好的实际计算意义。22用插值函数的积分，作为数值积分用插值函数的积

16、分，作为数值积分)()()()()()(00inibaiibaniibannxfdxxldxxfxldxxLfI ia代数精度代数精度:由由Lagrange插值的余项表达式插值的余项表达式(1)1( )( )( )(1)!nnnfR xxn于是，余项于是，余项(1)1()()()( )( )( ,)(1)!nbbnnnxaafIfIfRx dxx dxa bn可以看出，至少可以看出，至少n阶代数精度阶代数精度nkxxfxfkn,)(, 0)()1(插值型插值型（5.4）23 反之，如果求积公式反之，如果求积公式(5.3)至少具有至少具有n次代数精次代数精度，则它必是插值型的。度，则它必是插值型

17、的。 事实上，如果求积公式事实上，如果求积公式(5.3)至少具有至少具有n次代数精度，则次代数精度，则当当f(x)分别取分别取l0(x),l1(x),ln(x)（注意到它们都是（注意到它们都是n次多项式）次多项式）时，求积公式均精确成立。比如取时，求积公式均精确成立。比如取f(x)=li(x)代入得代入得因此，求积公式是插值型的。因此，求积公式是插值型的。0011( )()()()(0,1, )biiin inail x dxa l xal xa l xain定理定理1、形如、形如(5.3)的求积公式至少具有的求积公式至少具有n次代数精度的充分必次代数精度的充分必要条件是，它是插值型的。要条件

18、是，它是插值型的。24例例5、验证求积公式、验证求积公式是插值型求积公式。是插值型求积公式。解：从求积公式中可以看出，求积节点为解：从求积公式中可以看出，求积节点为求积系数求积系数a0=1,a1=1.要验证求积公式是插值型的，就是要验要验证求积公式是插值型的，就是要验证证计算知计算知1111( )33f x dxff0111,33xx 11( )0,1iial x dxi11001113( )11133xlx dxdxa 2511111113( )11133xl x dxdxa 所以所给求积公式是插值型求积公式。所以所给求积公式是插值型求积公式。26例例6、给定求积节点、给定求积节点x0=0,

19、x1=1,试推导出积分试推导出积分的插值型求积公式，并写出其截断误差。的插值型求积公式，并写出其截断误差。解：设求积公式为解：设求积公式为要使其成为插值型，则要使其成为插值型，则所以求积公式为所以求积公式为11( )f x dx1011( )(0)(1)f x dxa fa f1100111111111( )20 10( )01 0 xalx dxdxxal x dxdx11( )2 (0)f x dxf27该插值型求积公式的该插值型求积公式的截断误差截断误差为为 111111( )1()(0)(1)2xRfR x dxfxxdx其中其中( 1,1)x 28Newton-Cotes 积分积分

20、将积分区间将积分区间a,bn等分，分点记为等分，分点记为 a=x0 x1xn=b, 则有则有xi=a+ih, (i=0,1,2,n) 在求积节点等距情形下，构造的插值型求积公式称为在求积节点等距情形下，构造的插值型求积公式称为牛牛顿柯特斯（顿柯特斯（NewtonCotes）公式）公式。bahn29dtntititttnininhdthinintititttdxxlaninninthaxbaii00)() 1)(1() 1()!( !) 1( ) 1()!( !)() 1)(1() 1()()()(niiCaba)( niC(b-a)仅与步长仅与步长n有关，可以预先求出有关，可以预先求出30于是

21、，插值型求积公式于是，插值型求积公式(5.4)可写成可写成( )0( )()( )( )nbniinaif x dxbaCf xIx（5.5）式（式（5.5）称为）称为Newton-Cotes求积公式求积公式， 叫做叫做柯特斯系数柯特斯系数，其求积其求积截断误差截断误差为为( )niC (1)()1(1)!2(1)0( )() (1)()( , )(1)!nxbfnnnannnxxRfx dxhft ttn dta bn（5.6）31( )01nniiC（权质性）.数(n)(n)i i柯柯特特斯斯系系C 的C 的性性：证数为对积将n nb b( (n n) )i ii ia ai i= =0

22、0n nb b( (n n) )i ii ia ai i= =0 0n n( (n n) )i ii ii i= =0 0明明：由由于于N Ne ew wt to on nC Co ot te es s公公式式代代精精度度至至少少n n次次，故故f f( (x x) )= = 1 1, ,求求公公式式f f( (x x) )d dx x( (b b - - a a) )C Cf f( (x x ) )精精确确成成立立，f f( (x x) )= = 1 1代代入入得得1 1d dx x = =( (b b - - a a) )C Cf f( (x x ) )即即C Cf f( (x x ) )

23、= = 1 132N1时2121) 1(10)1(110)1(0dttCdttC1( )( )( )2baIff af b梯形公式33N2时时61) 1(4164)2(2161)2)(1(4120)2(220)2(120)2(0dtttCdtttCdtttC2( )( )4 ()( )62babaIff aff bSimpson公式34N4时，计算可得柯特斯系数为时，计算可得柯特斯系数为4( )7 ( )32 ( ) 12 ( )32 ( )7 ( )90baIff af df cf ef b柯特斯公式柯特斯公式(4)(4)(4)(4)(4)01234732123279090909090CCC

24、CC其中其中3,424baabdaceaba351、梯形公式、梯形公式31( )( )()( )()()()()( )2!2!12bbaaffb aE fx a x bdxx a x bdxf 此处用了积分中值定理此处用了积分中值定理误差误差362、Simpson公式公式2332(4)2(4)5(4)( )( )( )( )()()( )( ) ()()4!2( )() ()()( )4!22880babaEfI fS fI fI PI PS ffabxaxxb dxfabbaxaxxb dxf 注意到，注意到，Simpson公式有公式有3次代数精度，因此为了对误差有次代数精度，因此为了对误差

25、有更精确地估计，我们用更精确地估计，我们用3次多项式估计误差次多项式估计误差)2( )2( ),2()2(),()(),()(3333bafbaPbafbaPbfbPafaP为037 猜想：猜想：N-C积分，对偶数有积分，对偶数有n1次代数次代数精度，而奇数为精度，而奇数为n次代数精度。次代数精度。38定理定理3、当等分数、当等分数n为偶数时，为偶数时，NewtonCotes公式公式(5.5)至少至少具有具有n+1次代数精度。次代数精度。证明：由定理证明：由定理1可知，可知， (5.5)至少具有至少具有n次代数精度。下面只次代数精度。下面只需证明，需证明，当当n为偶数时，式为偶数时，式(5.5

26、)对对f(x)=xn+1精确成立精确成立。也就。也就是截断误差是截断误差Rnf=0. 由于由于f(x)=xn+1,所以所以f(n+1)(x)=(n+1)!，由式，由式(5.6)得得令积分变量变换令积分变量变换 ，因为，因为n为偶数，故为偶数，故 为整数，为整数， 20(1)()nnnRfht ttn dt2ntz2n39 222222222222()(1)22(1) (1)(1)()22(1 )(2 )20nnnnnnnnnRfhzznnzz zzzdznhz zzzdz证毕。证毕。40复化求积公式复化求积公式 复化梯形公式复化梯形公式 将区间将区间a,b划分为划分为n等分，分点等分，分点在每

27、个子区间在每个子区间xk,xk+1(k=0,1,n-1)上采用梯形公式，则上采用梯形公式，则,0,1,.,kb axakh hknn110110( )( )()()( )2kknbxaxknkknkIf x dxf x dxhf xf xRf41记记11011()()2( )2()( )2nnkkknkkhTf xf xhf af xf b（5.7）称为称为复化梯形公式复化梯形公式，其余项为，其余项为3110( )() ,(,)12nnnkkkkkhRfITfxx42 项为当长时敛n-1n-1kkkkkk0 k n-10 k n-10 k n-10 k n-1k=0k=0n-1n-1k kk=

28、0k=0n nn na a2 2b b2 21 1min f (min f ()f ()f ()max f ()max f () )n n1 1f (f ()=f ()=f (由由于于f(x) C a,b,且f(x) C a,b,且所所以以 (a,b)使(a,b)使于于是是复复化化梯梯形形公公式式余余由由上上式式可可知知，步步h h 0即0即n n +，R (R (复复) )n nb-ab-aR (f)= -R (f)= -化化梯梯形形公公式式是是收收于于f(f(f)f) 0 0 x)dxx)dx, ,h f (h f (1212故故) )的的。（5.8）43 复化复化Simpson公式公式

29、将区间将区间a,b划分为划分为n等分，在每个子区间等分，在每个子区间xk,xk+1 上采用上采用Simpson公式，若记公式，若记xk+1/2=xk+h/2，则，则11011/210( )( )()4 ()()( )6kknbxaxknkkknkIf x dxf x dxhf xf xf xRf44记记11/210111/201()4 ()()6( )4()2()( )6nnkkkknnkkkkhSf xf xf xhf af xf xf b(5.9)称为称为复化复化Simpson公式公式，其余项为，其余项为41(4)10( )(),(,)180 2nnnkkkkkhhRfISfxx 4544

30、(4)( ) ,( )02(18nbahRfff xC a b 于是当时，类似复化梯形公式余项有与复化梯形公式的复化Simps讨论类似，可以on公式证明，是收敛的。(5.10)46例例7、已知某河宽为、已知某河宽为20m，测得水深，测得水深f(x)如下表（单位：如下表（单位：m）：）：xk02468101214161820f(xk)1.01.51.83.02.82.53.02.82.01.81.4分别用复化梯形公式及复化分别用复化梯形公式及复化Simpson公式计算河水的截面公式计算河水的截面积积200( )f x dx47解解：用复化梯形公式。取用复化梯形公式。取n=10，步长，步长h=2,

31、于是于是92010100912( )()()2( )2()( )221.02(1.5 1.83.02.82.523.02.82.0 1.8) 1.444.8()kkkkkhf x dx Tf xf xhf af xf bm48（2）用复化）用复化Simpson公式。公式。n=5,h=4,这时这时12051/2100441/2012( )()4 ()()64( )4()2()( )641.04(1.53.02.52.8 1.8)62(1.82.83.02.0) 1.445.3()nkkkkkkkkhf x dxSf xf xf xf af xf xf bm49计积时断 误绝 对过问 应 将 区

32、间1 1- - x x0 0- - 4 4例例 8 8 、 用用 复复 化化 梯梯 形形 公公 式式算算分分e ed d x x 的的 近近 似似 值值，1 1要要 求求 截截差差 的的值值 不不 超超 1 1 0 0，2 2 0 0 , ,1 1 分分 成成 多多 少少 等等 份份 ？01222424,( ),max |( ) |1,1| |()()1212112111| |10102122140.841)xxxefxefxbaRfhfhfhRfhnh 解 ： f(x)=求 出 导 数于 是用 复 化 梯 形 公 式 时 ， 根 据 其 余 项 估 计 式要 使， 只 需因 此h0.0244

33、9从 而n=（ 取50龙贝格（龙贝格（RombergRomberg）求积公式）求积公式v 复化梯形公式逐次分半算法复化梯形公式逐次分半算法11211121 , 2 (0,1,2,)(0,1,2,2 ).2( )2()( )2mmmmmkmmmmka bnmbahxakhkhTf af akhf b将积分区间分成等份，步长记为，则分点根据复化梯形公式，可以计算有当，于时是m-1m-1mm1 12 2n = 2T的n = 2T的值值51现将区间将时个区间个间递计mmm mm mm mm m1 1mm-1mm-12 22 2 2 2m m2 2k 1k 12 21 1m m= =a,b分a,b分成成

34、2 等2 等份份，也也就就是是2等2等分分的的每每小小逐逐分分半半，下下面面根根据据复复化化梯梯形形公公式式，算算T Th hT=f(a)+2f(a+kh )+f(b)T=f(a)+2f(a+kh )+f(b)下下面面建建立立T T2 2与与T 之T 之的的推推公公式式。52于是m mm mm m- -1 1m m- -1 1m m- -1 1m m- -1 1m mm m- -1 12 2 1 1m mm m2 2k k= =1 12 21 1m mm mk k= =1 12 2m mk k= =1 12 21 1m mm m- -1 1k k= =1 12 2m mm mk k= =1 1

35、m mm m2 22 2k k= =h hT Tf f( (a a) )+ + 2 2f f( (a a + + k kh h ) )+ + f f( (b b) )2 2h hf f( (a a) )+ + 2 2f f( (a a + + 2 2k kh h ) )2 2+ +2 2f f( (a a + +( (2 2k k - -1 1) )h h ) )+ + f f( (b b) )h hf f( (a a) )+ + 2 2f f( (a a + + k kh h) )+ + f f( (b b) )2 2+ +h hf f( (a a + +( (2 2k k - -1 1)

36、)h h ) )1 1T TT Th hf f( (a a + +( (2 2k k - -1 1) )h h ) )2 2(1,2,)m m m- -1 12 21 1(5.11)53式（式（5.11）称为复化梯形公式的）称为复化梯形公式的逐次分半计算公式逐次分半计算公式。注意注意：运用逐次分半算法时，一般预先给定误差限：运用逐次分半算法时，一般预先给定误差限，当，当|T2mT2m-1|时，停止计算，取时，停止计算，取T2m为所给积分的近似值。为所给积分的近似值。例例9、用复化梯形公式逐次分半算法计算、用复化梯形公式逐次分半算法计算的近似值，要求的近似值，要求|T2mT2m-1|10-512041dxx5421214241()0,1111(0)(1)(42)322120,15.1121113.12221330,1445.1111133.131176742444fxxTffxTTfxxTTff解 ： （ ）， 在 区 间上 用 梯 形 公 式（ ） 将分 成 两 等 份 ，是 新 分 点 ， 由 （） 式 得（ ） 将分 成 四 等 份 ， 这 时,是 新 分 点 ， 由（） 式 得这 样 不 断1248,TTTT-5512256512将 区 间 逐 次 分 半 ， 由 递 推 公 式 可 依 次 求 出