数值计算方法第07章数值微分与数值积分

1、1第七章第七章 数值微分与数值积分数值微分与数值积分1 数值微分数值微分2 NewtonCotes求积公式求积公式3 复化求积公式复化求积公式4 Romberg求积公式求积公式5 Gauss型求积公式型求积公式23 利用离散点上函数的信息求函数导数近似值利用离散点上函数的信息求函数导数近似值的方法的方法, 称为称为数值微分数值微分. 差商型求导公式差商型求导公式 插值型求导公式插值型求导公式1 数值微分数值微分4由导数定义由导数定义hxfhxfxfh)()(lim)(0 当当h很小时很小时, 可用可用差商差商近似导数近似导数.5 差商型求导公式差商型求导公式 (3)中心差商公式中心差商公式0,

2、)()()( hhxfhxfxf,)()()(hhxfxfxf .2)()()(hhxfhxfxf (1) 向前差商公式向前差商公式 (2) 向后差商公式向后差商公式6 几何意义几何意义hx xhx ABChxfhxfkBC)()( hhxfxfkAB)()( hhxfhxfkAC2)()( B点切线斜率点切线斜率)(xf 从几何直观看从几何直观看: 中心差商效果最好中心差商效果最好7 截断误差截断误差)(2)( )()()( 1hOhhxfhxfhxfxf )(2)( )()()( 2hOhhxfhhxfxfxf )(12)()(2)()()( 223)3(3)3(hOhhxfhxfhhxf

3、hxfxf 其中其中1,0321 由由Taylor公式可得公式可得8 二阶导数的中心差商公式二阶导数的中心差商公式2)()(2)()( hhxfxfhxfxf 截断误差截断误差)(12)()(2)()( )4(22 fhhhxfxfhxfxf 9101112数值积分数值积分一、数值积分的必要性一、数值积分的必要性本章主要讨论如下形式的一元函数积分本章主要讨论如下形式的一元函数积分badxxffI)()( )( )( )( )baI ff x dxF bF a在微积分里，按在微积分里，按Newton-Leibniz公式公式求定积分求定积分 f x F x13 F x f x实际问题实际问题142

4、 sinf xx这个问题就是要求由函数这个问题就是要求由函数0 x 48x 15dxxdxxfL48024802)(cos1)(1 16类似的，下列函数也不存在由初等函数表示的原函数类似的，下列函数也不存在由初等函数表示的原函数:2,1,ln1,sin,cos,sin322xexxxxxx173222xx)322ln(21693216332412222xxxxxx18 f x f xx1423454.5688.5原来通过原函数来原来通过原函数来计算积分有它的局计算积分有它的局限性。那限性。那怎么办呢？怎么办呢？呵呵呵呵这就需要积这就需要积分的数值方法来帮分的数值方法来帮忙啦。忙啦。19二、数值

5、积分的基本思想二、数值积分的基本思想1、定积分的几何意义、定积分的几何意义badxxffI)()(abxyo f x202、数值积分的理论依据、数值积分的理论依据 f x依据依据积分中值定理积分中值定理, 对于连续函数对于连续函数 ,在在 内存在一点内存在一点 ,使得使得, a b)()()()(fabdxxffIba称称 为为 在区间在区间 上的平均高度上的平均高度. f, a b ?f f x213、求积公式的构造、求积公式的构造 若简单选取区间端点或中点的函数值作为平均高度，则若简单选取区间端点或中点的函数值作为平均高度，则可得一点求积公式如下：可得一点求积公式如下：左矩形公式：左矩形公

6、式：中矩形公式：中矩形公式：右矩形公式：右矩形公式： Iff aba 2abIffba Iff bba22xyOab f x f a左矩形公式：左矩形公式： Iff aba23xyOab f x2abf2ab中矩形公式：中矩形公式： 2abIffba24xyOab f x f b右矩形公式：右矩形公式： Iff bba25 若取若取 两点，并令两点，并令 ，则可得梯形公式，则可得梯形公式（两点求积公式）（两点求积公式）, a b 2f af bf 2f af bI fba26xyOab f x f a f b27则可得则可得Simpson公式公式(三点求积公式三点求积公式), ,2aba b

7、c 46f af cf bf 若取三点，若取三点， 并令并令 46f af cf bIfba28 一般地一般地 ，取区间，取区间 内内 个点个点, a b1n ,0,1,2,.,ixin ,0,1,.,if xin f处的高度处的高度通过通过加权平均加权平均的方法近似地得出平均高度的方法近似地得出平均高度这类求积方法称为这类求积方法称为机械求积机械求积:)()()(0ibaniixfabdxxf29 或写成或写成: :数值积分公式数值积分公式求积系数求积系数 求积节点求积节点 )()(0kbankkxfAdxxf30记记0( )()nnkkkIfA f x0( )( )( )( )(),nbn

8、kkakR fI fIff x dxA f x称称为数值为数值求积公式求积公式称为求积公称为求积公式余项式余项(误误差差).（1）（2）31构造或确定一个求积公式，要解决的问题包括构造或确定一个求积公式，要解决的问题包括:(i) 确定求积系数确定求积系数 和求积节点和求积节点 kAkx；(iii) 求积公式的误差估计和收敛性分析求积公式的误差估计和收敛性分析.(ii) 确定衡量求积公式好坏的标准；确定衡量求积公式好坏的标准；3233 依据积分中值定理依据积分中值定理 就是说，底为就是说，底为b a 而高为而高为 f ( ) 的矩形面积恰恰等的矩形面积恰恰等于所求曲边梯形于所求曲边梯形 f (x

9、)的面积的面积.)()(abfdxxfba 取取a, b内若干个节点内若干个节点xk 处的高度处的高度 f (xk ), 通过加通过加权平均的方法生成平均高度权平均的方法生成平均高度 f ( ), 这类求积公式称这类求积公式称机械求积公式机械求积公式式中式中 xk 称为称为求积节点求积节点, Ak 称为称为求积系数求积系数, 亦称伴随亦称伴随节点的权节点的权. nkkkbaxfAdxxf0)()(数值积分基本思想数值积分基本思想342 Newton-Cotes 公式公式基本思想基本思想: 利用利用插值多项式插值多项式).()(xfxLn banbadxxLdxxffI.)()()(其中其中Ln

10、(x)是是n阶阶Lagrange插值多项式，用插值多项式，用Ln (x)的的积分近似积分近似 f (x)的积分，即的积分，即插值型求积公式插值型求积公式35 bankjjjkjkdxxxxxA,0)()(由由 决定决定,与与 无关无关.节点节点 f (x) 在在a, b上取上取 a x0 x1 0, 使得使得|,|max0knk 则称该求积公式是则称该求积公式是稳定稳定的的. 求积公式的稳定性求积公式的稳定性72 若求积公式是稳定的若求积公式是稳定的, 则则 f (x)的观察值的较小的的观察值的较小的误差引起的求积结果的误差也是较小的误差引起的求积结果的误差也是较小的. 求积公式求积公式没有把

11、没有把 f (x)的误差的误差“放大放大”很多很多.73), 1 , 0()(nkfxfkk nkkkknnfxfAfIfI0)()()(证明证明因此复化梯形公式是数值稳定的因此复化梯形公式是数值稳定的.当当).(2)1(2abhhnh 定理定理 复化梯形公式是复化梯形公式是数值稳定数值稳定的的. nkkkkfxfA0)( nkkA0 74x0 x2xf (x)x4hhxn 2hxnmnnabh2, .hx3x1xn 1 复化复化Simpson公式公式分片二次多项式近似分片二次多项式近似75 将积分区间将积分区间a, b划分为划分为n=2m等分等分, 步长步长 h=( b a )/n, 分点分

12、点 xk= a+kh ( k=0, 1, , n). 在每个在每个小区间小区间 x2k 2 , x 2k ( k=1, , m)上用上用Simpson公式：公式： kkxxkkkkkxfxfxfxxdxxf222)()(4)(6)(21222222 )()(4)(321222kkkxfxfxfh 复化复化Simpson公式公式k=1, , m76 111122)(4)(2)()(3mkmkkkxfxfbfafh= Sn( f ) mkxxbadxxfdxxffIkk1222)()()( )()(4)(3212221kkkmkxfxfxfh 复化复化Simpson公式公式77)()()(fSfI

13、fRnn 当当 f (x)在在a, b上具有四阶连续导数时上具有四阶连续导数时, ),(),()(1)4()4(bamffmkk ),(222kkkxx 故得故得),(),(180)()(90)()4(4)4(5bafhabfmhfRn 复化复化Simpson公式的截断误差公式的截断误差 mkkfh1)4(5),(90 78 由复化由复化Simpson公式的截断误差知公式的截断误差知, 误差阶为误差阶为 h4, 收敛性是显然的收敛性是显然的, 事实上事实上,只要只要 f (x) Ca, b则则可得到可得到收敛性收敛性, 即即.)()(limdxxffSbann 由于求积系数均为正由于求积系数均

14、为正, 与复化梯形公式一样的与复化梯形公式一样的证法可得复化证法可得复化 Simpson公式是公式是数值稳定数值稳定的的.79例：例：计算计算dxx 10214 解：解： )1()(2)0(161718fxffTkk8kxk 其中其中= 3.138988494 )1()(2)(4)0(241oddeven8fxfxffSkk8kxk 其中其中= 3.141592502运算量运算量基本相同基本相同 显然用复化显然用复化Simpson公式计算精度较高公式计算精度较高, 这与它们这与它们的误差阶的结论是相符的的误差阶的结论是相符的.80例例 对于函数对于函数,sin)(xxxf 给出给出n=8的函数

15、表的函数表, 试用试用复化梯形公式及复化复化梯形公式及复化Simpson公式计算积分公式计算积分 10.sindxxxI解解ix)(ixf0.00.1250.250.3750.50.6250.750.8751.01.00.99739780.98961580.97672670.95885100.93615560.90885160.84147090.8771925应用复化梯形公式求得应用复化梯形公式求得T8=0.9456909应用复化应用复化Simpson公式求得公式求得S8=0.9460832准确值准确值 I=0.9460831两者运算量基本相同两者运算量基本相同8126.26.利用下面数据表，

16、利用下面数据表， 1. 用复化梯形公式计算积分用复化梯形公式计算积分的近似值；的近似值；2. 用复化用复化Simpson公式计算积分公式计算积分的近似值。的近似值。( (要求计算结果保留到小数点后六位要求计算结果保留到小数点后六位). 2.61.8( )If x dx2.61.8( )If x dxx1.82.02.22.42.6F(x) 3.120144.42569 6.04241 8.03014 10.4667582trapz: 复化梯形公式求积分复化梯形公式求积分.用法用法: trapz(X, Y), 其中其中X, Y为相同维数的向量为相同维数的向量.例例: X=0.125:0.125:

17、1.0;Y=sin(X)./X;X=0,X;Y=1,Y;trapz(X,Y)ans = 0.94569086358270Matlab函数函数83例例 若用复化求积公式计算积分若用复化求积公式计算积分dxeIx 10的近似值的近似值, 若要求计算结果有若要求计算结果有4位有效数字位有效数字, n应取多大应取多大?解解, 1110 dxeIex.105 . 04 1 , 0, 1| )(|)( xexfxk复化梯形公式的误差复化梯形公式的误差)( )(12|2 fabhRT .83.40 n若用复化梯形公式求积分若用复化梯形公式求积分, n取取41能达到精度要求能达到精度要求.2121n 4102

18、1 84故应取故应取n=4. 该例表明该例表明, 为达到相同的精度为达到相同的精度, 用复化用复化Simpson公式所需的计算量比复化梯形公式要少公式所需的计算量比复化梯形公式要少, 这也说明这也说明了复化了复化Simpson公式的精度高公式的精度高.复化复化Simpson公式的误差公式的误差)()(180|)4(4 fabhRS .25. 3 n41801n 41021 857.3.3 逐次分半算法逐次分半算法(变步长方法变步长方法) 86 梯形法的梯形法的递推公式递推公式 abhbfafabT),()(211、)2(22)()2(2)(4)2)()2(21(2)2(212)(22,2,2,

19、212bafabTbfbafafabbfbafabbafafabTabhbbabaaba分半为、将 ,)d( baxxfI对87 复化复化梯形公式的逐次分半算法梯形公式的逐次分半算法将区间将区间a, b分成分成n=2m等分等分, 记记,2mmabh ), 2 , 1 , 0( m883122312422)(2)()(2)(21)()(21(44 , 3 , 2 , 1 , 0 ,4,4)(3kkkkhafbfafhbfkhafafabTkkabakhaxabh再分半、加密一次区间)3()(21 )4)( 3()4(421)4)( 3()4(4)4)( 2(2)(2)(4222222241212

20、hafhafhTTabafabafabTabafabafababafbfafabTT，亦即即, 2 , 1 , 0 )(2)()(21212mkhafbfafhTmmkmm梯形公式的逐次分半算法（续梯形公式的逐次分半算法（续1）89梯形公式的逐次分半算法梯形公式的逐次分半算法1,kkxx)(21121kkkxxx1,kkxx)()(2)(411kkkxfxfxfhnabh90梯形公式的逐次分半算法（续梯形公式的逐次分半算法（续2）21)-(7 ) 1(2(21112122mmmkmmhkafhTT )( 2 )()(41 -n0211012knkkknxfhxfxfhT91 mhTTmm212

21、2所有新增加节点的函数值之和所有新增加节点的函数值之和其中其中.2mmabh 复化复化梯形公式的逐次分半算法梯形公式的逐次分半算法92以以n=8, m=3为例为例. 记记 fk= f (xk)x0 x2x4x6x3x1x5x7x8 )(2167654321808fffffffffabT )(21664280fffffab 75318ffffab 24T 所有新增加节点的函数值之和所有新增加节点的函数值之和. 3h93 复化梯形公式余项的后验估计复化梯形公式余项的后验估计);,(),(12112bafhabTIn );,(),(2122222bafhabTIn f ( 1 ), f ( 2 )

22、分别是分别是 f (x) 在在a, b上的上的n个点与个点与 2n 个点处的算术平均值个点处的算术平均值 (每个小区间上取一个点每个小区间上取一个点). 当当n较大时较大时, 有有.)( 1)()(21dxxfabffba 94因此因此, 若事先给定误差限若事先给定误差限 , 则当则当.3|2 nnTT时时, 就可停止计算就可停止计算, 并认为并认为 T2n是满足精度要求的近是满足精度要求的近似值似值.;412 nnTITI);()(21 ff ).(3122nnnTTTI 95 复化复化Simpson公式的逐次分半算法公式的逐次分半算法将区间将区间a, b分成分成 n=2m 等分等分, 记记

23、,2mmabh evenodd2)(2)(4)()(3kmkmmkhafkhafbfafhSm, 2 , 1 m称称 为为Simpson序列序列.2mS96;1612 nnSISI);,(),(18011)4(4bafhabSIn );,(),(218022)4(42bafhabSIn );()(2)4(1)4( ff ).(15122nnnSSSI 因此因此, 若事先给定误差限若事先给定误差限 , 则当则当.15|2 nnSS时可停止计算时可停止计算, 取取 S2n为满足精度要求的近似值为满足精度要求的近似值. 复化复化Simpson公式余项的后验估计公式余项的后验估计9722、已知函数f(

24、x)在若干点处的值： 试计算积分 的梯形值 以及Simpson值 。X-1-0.500.51F(x)02.12532.125011( )f x dx124,T T T12,S S表表达达式式知知：由由复复化化梯梯形形公公式式的的余余项项，但但收收敛敛慢慢，精精度度低低。复复化化梯梯形形公公式式算算法法简简单单).,(),( )2(12);,(),( 12222bafhabTIbafhabTInn )( )( ff 假假定定412 nnTITI 7.4 7.4 龙贝格算法龙贝格算法).(3122nnnTTTI 1443134)(31:2222nnnnnnnTTTTTTTI于于是是的的误误差差很很

25、小小很很小小，可可保保证证若若nnnTTT22事后估计法利用计算结果来估计误差的方法利用计算结果来估计误差的方法龙贝格算法龙贝格算法 当当n=1 =1 时时, ,我们计算上式右端我们计算上式右端 )(21)(21)(31)(21)2()(212)(3414412bfafabbfbafafabTTT1)(61)2(64)(61)(Sbfbafafab 这恰好是辛普森公式的结果，即有这恰好是辛普森公式的结果，即有121141144TTS 比梯形公式有比梯形公式有更好的精确度更好的精确度龙贝格算法龙贝格算法类似地可验证：类似地可验证：nnnTTSTTS1411441411442242 即即nnnTT

26、S31342 辛普森积分值辛普森积分值_nS龙贝格算法龙贝格算法式式的的余余项项注注意意复复化化辛辛普普森森求求积积公公)()(21804)4(4hOfhabSIRnn ) (4180)4(422 fhabSIRnn )()()4()4( ff 假假定定144)(15122222 nnnnnSSSSSI)(,215121612nnnSISISSSInn 龙贝格算法龙贝格算法可以验证可以验证nnnSSC15115162 复复化化柯柯特特斯斯积积分分值值_nC事实上事实上 C1=(42S2-S1)/(42-1) =16S2/15-S1/15=(7y0+32y1+12y2+32y3+7y4)/90

27、恰为柯斯特公式。恰为柯斯特公式。同理，同理， C2=(42S4-S2)/(42-1)，. )()4(945)(2)6(6 fhabCIRnn 余项，余项，注意柯斯特积分公式的注意柯斯特积分公式的) ()8(945)(2)6(622 fhabCIRnn 龙贝格算法龙贝格算法 即即, R, R1 1=(4=(43 3C C2 2-C-C1 1)/(4)/(43 3-1), -1), R R2 2 =(4 =(43 3C C4 4-C-C2 2)/(4)/(43 3-1), . -1), . Rn= (4 Rn= (43 3 C2n -Cn )/ (4C2n -Cn )/ (43 3-1);-1);

28、 上式即为龙贝格公式，得龙贝格值序列上式即为龙贝格公式，得龙贝格值序列)()()6()6( ff 假假定定.63163642nnCCI ,6412 nnCICI.144323 nnnCCR龙龙贝贝格格积积分分值值_nR龙贝格算法龙贝格算法105nnnnnTTTTTI3134)(31222 4 Romberg求积公式求积公式启示启示: 是否用是否用 复化梯形公式余项的后验估计表明复化梯形公式余项的后验估计表明nnTT31342 逼近逼近 I ( f ) 比用比用 T2n要好要好. 事实上有事实上有.313422nnnSTT 即梯形值序列的巧妙线性组合得到即梯形值序列的巧妙线性组合得到Simpso

29、n序列序列!106以以n=4为例加以说明为例加以说明. 记记 fk= f (xk), )(2167654321808fffffffffabT )(28642804fffffabT 8abkaxk )(422847654321808fffffffffabT .3)4(488TTS )( 2)( 48642753180fffffffffab 484TT 107).(15122nnnSSSI .1516)(151222nnnnnSSSSSI 逼近逼近 I ( f ) 比用比用 S2n要好要好.回答回答: 是的是的, 记记,15)16(22nnnSSC 15)16(2nnSS 则它恰为复化则它恰为复化

30、Cotes公式公式; 且有如下误差估计式且有如下误差估计式).(| )(|62hOfICn 复化复化Simpson公式余项的后验估计表明公式余项的后验估计表明 问题问题: 是否用是否用108.631636422nnnRCC 类似地可以得到类似地可以得到2nR其中其中被称为被称为Romberg序列序列.).(| )(|82hOfIRn 截断误差截断误差:龙贝格龙贝格积分积分 /* Romberg Integration */例：例：计算计算dxx 10142 已知对于已知对于 = 10 6 须将区间对分须将区间对分 9 次，得到次，得到 T512 = 3.14159202由由 来计算来计算 I

31、效果是否好些？效果是否好些？nnnnTTTTI313414422 考察考察412 nnTITI483134TT = 3.141592502 = S4一般有：一般有：nnnSTT 1442nnnCSS 144222nnnRCC 144323Romberg 序列序列 Romberg 算法：算法： ? ? ? T1 =)0(0T T8 =)3(0T T4 =)2(0T T2 =)1(0T S1 =)0(1T R1 =)0(3T S2 =)1(1T C1 =)0(2T C2 =)1(2T S4 =)2(1T1101T2T4T8T16T32T2S4S8S16S32S4C8C16C32C8R16R32R3

32、422nnnTTS 151622nnnSSC 636422nnnCCR 停机准则停机准则:梯形值序列梯形值序列Simpson序列序列Cotes序列序列Romberg序列序列 Romberg求积公式求积公式 用若干个积分近似值推算出更为精确的积分近用若干个积分近似值推算出更为精确的积分近似值的方法，称为似值的方法，称为外推方法外推方法。 列和龙贝格序列。列和龙贝格序列。辛普森序列、柯特斯序辛普森序列、柯特斯序分别称为梯形序列、分别称为梯形序列、序列序列NNNNRCST,龙龙贝贝格格序序列列为为止止。故故通通常常到到及及续续外外推推，但但因因由由龙龙贝贝格格序序列列还还可可以以继继),4(0141

33、1144 mmmm龙贝格算法龙贝格算法), 2 , 1 , 0, 2 , 1(144)(1)1(1)( kmTTTmkmkmmkm，记记2)(1kSTk是是辛辛普普森森值值序序列列2)(2kCTk是是柯柯特特斯斯值值序序列列2)(0kTTk就就是是梯梯形形值值序序列列则则2)(3kRTk是是龙龙贝贝格格序序列列计算过程见表：计算过程见表：龙贝格算法龙贝格算法 T1 T2 S1 T4 S2 C1 T8 S4 C2 R1 T16 S8 C4 R2 T32 S16 C8 R4.1442 TTSnnn144323 CCRnnn144222 SSCnnn 122 nnRR上面是上面是RombergRom

34、berg的计算表的计算表若若 则计算停止则计算停止龙贝格算法龙贝格算法61021 要求误差不超过要求误差不超过9207355. 0)1()0(21)0(0 ffT解：解：9397932. 0)5 . 0(2121)0(0)1(0 fTT9461457. 034)0(0)1(0)0(1 TTT9445136. 0)75. 0(25. 0(4121)1(0)2(0 ffTT 用用RombergRomberg方法计算积分方法计算积分 10sindxxxI近似值近似值例例9460870. 034)1(0)2(0)1(1 TTT9460830. 01516)0(1)1(1)0(2 TTT9456910.

35、 0)875. 0()625. 0()375. 0(125. 0(8121)2(0)3(0 ffffTT9460835. 034)2(0)3(0)2(1 TTT9460832. 01516)1(1)2(1)1(2 TTT9460832. 06364)0(2)1(2)0(3 TTT67)0(2)0(31021102 TT9460832. 0sin10 dxxx117例例 计算计算.sin10dxxxI 2)1()0(1ffT =0.9207355)5 . 0(212112fTT =0.9397933 )43()41(412124ffTT=0.9445135=0.9456909 )87()85()

36、83()81(812148ffffTT解解先求梯形值序列先求梯形值序列第二种写法118nnTnSnCnR24180.92073550.93979330.94451350.94569090.94614590.95608690.94608330.94608300.94608310.94608313422nnnTTS 151622nnnSSC .636422nnnCCR 利用只有两三位有效数字利用只有两三位有效数字的的T1, ,T8 经过三次外推得经过三次外推得到到7位有效数字位有效数字. 可见加速的可见加速的效果十分显著效果十分显著.用用Romberg算法计算如下算法计算如下fx_:=Sinx/x

37、;fx_:=Sinx/x;a=0;b=1;f0=1;f1=0.8414709;a=0;b=1;f0=1;f1=0.8414709;hk_:=(b-a)/2k;T0,0=(b-a)/2hk_:=(b-a)/2k;T0,0=(b-a)/2* *(fa+fb);(fa+fb);T0_,k_:=1/2T0_,k_:=1/2* *T0,k-1+hkT0,k-1+hk* *Sumfa+hkSumfa+hk* *(2i-(2i-1),i,1,2(k-1);1),i,1,2(k-1);NTableT0,k,k,0,3,6;NTableT0,k,k,0,3,6;MatrixForm%MatrixForm%Tm_

38、,k_:=(4mTm_,k_:=(4m* *Tm-1,k+1-Tm-1,k)/(4m-1);Tm-1,k+1-Tm-1,k)/(4m-1);NTableTm,k,m,1,3,k,0,2,7;NTableTm,k,m,1,3,k,0,2,7;MatrixForm%;NTableT1,k,k,0,2,7;MatrixForm%;NTableT1,k,k,0,2,7;MatrixForm%;NTableT2,k,k,0,1,7;MatrixForm%;NTableT2,k,k,0,1,7;MatrixForm%;NT3,0,7MatrixForm%;NT3,0,7算算 法法 badxxfI)(设设k

39、kabhba22, 等等分分，将将区区间间)(hTTn记记为为近近似似值值为为用用复复化化梯梯形形公公式式求求得得的的记记.)(634221 hahahaIhT则则)()( 12)(:22hOfhabhTI注意到注意到 .212121)2(636424212 hahahaIhT将余项泰将余项泰勒展开勒展开)2(,)0()(lim20hTTIThTnh 而而 外推法的一般讨论外推法的一般讨论.3)()2(4)(62411 hbhbIhThThT)()(41hOIhT 即即,),2(),(211nnSShThT就就是是辛辛普普森森序序列列如如此此构构造造序序列列)(151)2(1516)(,16)

40、(11262411hThThThhIhT 若若令令由由 外推法外推法82612)(hhIhT ,)(22nnCChT就就是是柯柯特特斯斯序序列列这这样样构构造造的的序序列列)()(62hOIhT )2(22)1(21)(), 2 , 1(2mmmhhIhTmm 列列形形式式：次次加加速速后后，余余项项便便取取下下经经过过阶阶，高高一一次次，误误差差的的量量级级便便提提如如此此继继续续下下去去，每每加加速速此处理方法称为理查森此处理方法称为理查森(Richardson)(Richardson)外推加速方法外推加速方法。外推法外推法123 理论依据理论依据: 复化梯形公式的余项展开复化梯形公式的余

41、项展开.记记),(hTTn 定理定理 设设,)(baCxf 则则 kkhhhIhT24221)( 其中系数其中系数 k ( k=0, 1, )是与是与 h 无关的常数无关的常数. T (h) 逼近逼近 I 的速度是的速度是 O ( h2 )阶阶.124,3)()2(4)(262411 kkhhhIhThThT 当区间当区间a, b 2n等分时等分时, 则有则有),2(2hTTn 在定理中以在定理中以 h/2 代替代替 h 得得,2164224221 kkhhhIhT 上式乘以上式乘以4减去减去 T(h) 再除以再除以3, 记之为记之为 T1(h), 得得T1 (h) 逼近逼近 I 的速度是的速

42、度是 O ( h4 )阶阶, 效果比效果比 T (h)好好, 它它不是别的不是别的, 就是就是Simpson序列序列.125 类似地类似地 kkhhhIhThThT2826111215)()2(16)( ,162262411 kkhhhIhT 上式乘以上式乘以16减去减去 T1(h) 再除以再除以15, 记之为记之为 T2(h), 得得T2 (h) 逼近逼近 I 的速度是的速度是 O ( h6 )阶阶, 效果比效果比 T1 (h)好好, 它不是别的它不是别的, 就是就是Cotes公式序列公式序列.126 对对Cotes公式序列进行同样处理得到公式序列进行同样处理得到Romberg公式序列公式序

43、列.Richardson外推加速方法外推加速方法 也称为也称为Romberg求积算法求积算法 收敛性说明收敛性说明: 如果如果 f (x) 充分光滑充分光滑, 那么梯形公那么梯形公式序列式序列, Simpson公式序列公式序列, Cotes公式序列公式序列, Romberg公式序列均收敛到所求的积分值公式序列均收敛到所求的积分值. 对于对于 f (x)不充分光滑的函数也可用不充分光滑的函数也可用Romberg算算法计算法计算, 只是收敛慢一些只是收敛慢一些. 也可以直接使用复化也可以直接使用复化Simpson公式计算公式计算.127例例 用用Romberg算法计算积分算法计算积分.1023dx

44、xI 解解23)(xxf 在在0, 1上仅是一次连续可微上仅是一次连续可微用用Romberg算法计算结果见下表算法计算结果见下表nnTnSnCnR241816320.50.4267770.4070180.4018120.4004630.4001180.4023690.4004320.4000770.4000140.4000020.4003020.4000540.4000090.4000020.4000500.4000090.400002128 23、用Romberg方法计算积分 的近似值，要求误差不超过12041dxx510129. )()(1 nkkkbaxfAdxxf5 Gauss型求积公

45、式型求积公式 基本思想基本思想 设计求积公式设计求积公式:在节点数在节点数 n 固定时固定时, 适当地选取求积节点适当地选取求积节点 xk 与求与求积系数积系数 Ak , 使求积公式具有使求积公式具有最高的代数精确度最高的代数精确度.130例例 确定确定x1, x2, A1, A2, 使求积公式使求积公式)()()(221111xfAxfAdxxf 具有最高次的代数精确度具有最高次的代数精确度.x2x11 选取选取 (A1 , A2 , x1 , x2)使该求积公式对使该求积公式对 f (x) = 1, x, x2, x3 时等号成立时等号成立. 1131)()()(221111xfAxfAd

46、xxf 313111 0 32 0 21 12121322311113322221111222211112111xxAAxAxAdxxxfxAxAdxxxfxAxAxdxxfAAdxf)31()31()(11ffdxxfI 对对 f = 1, x, x2, x3 积分精确成立积分精确成立 四个方程四个未知数四个方程四个未知数132梯形公式与梯形公式与Gauss求积公式的比较求积公式的比较 对对1, x求积公式精确求积公式精确成立成立(1(1次代数精确度次代数精确度) ) 对对1, x, x2, x3求积公式精求积公式精确成立确成立(3 次代数精确度次代数精确度)133)()()()( :333

47、221111xfAxfAxfAdxxfn x3x11x2 1 选取选取(A1, A2 , A3 , x1, x2 , x3) 使该求积公式对使该求积公式对 f (x) = x0, x1, x2, x3, x4, x5 时等号成立时等号成立.区间区间 1, 1上的上的Gauss求积公式求积公式134 3/503/55/98/95/9 0 52 0 32 0 21 132132153352112511554334211241144333321123113323322112211223321121131121xxxAAAxAxAxAdxxxfxAxAxAdxxxfxAxAxAdxxxfxAxAxAd

48、xxxfxAxAxAxdxxfAAAdxf)53(95)0(98)53(95)(11fffdxxfI 对对 f = x0, x1, x2, x3, x4, x5 求积公式等号成立求积公式等号成立135)()()()(221111nnxfAxfAxfAdxxf 选取选取(A1, A2, , An , x1, x2 , , xn)使该求积公式对使该求积公式对f (x) = x0, x1, x2, , x2n 1时等号成立时等号成立. .区间区间 1, 1上的上的Gauss求积公式求积公式136 0 121 32 0 21 11212211212111212222221122211222222211

49、2211222112111121 nnnnnnnnnnnnnnnnnnnxAxAxAdxxxfxAxAxAndxxxfxAxAxAdxxxfxAxAxAxdxxfAAAdxf 2n个方程个方程2n个未知数个未知数, 非线性方程非线性方程, 解是否存在唯解是否存在唯一？若存在一？若存在, 如何求解？如何求解？ 对对 f = x0, x1, x2, , x2n 1求积公式等号成立求积公式等号成立137. )()()(1 nkkkbaxfAdxxfx 研究最一般情形的带权积分研究最一般情形的带权积分:具有具有2n1次代数精确度次代数精确度, 则称这组节点则称这组节点 xk为为 Gauss 点点, 上

50、述公式称为带权函数上述公式称为带权函数 (x)的的Gauss型型求求积公式积公式.定义定义 如果一组节点如果一组节点x1, x2, , xn a, b能使求积公式能使求积公式 nkkkbaxfAdxxfx1)()()( Gauss型型求积公式求积公式138. 12, 1 , 0,)(1 nmxAdxxxmknkkbam 由上述由上述 2n个方程确定全部个方程确定全部 的的2n个待定参数个待定参数 xk , Ak ( k=1, , n), 使求积公式至少具有使求积公式至少具有 2n 1次代数精次代数精确度确度. 但上述方程组是非线性方程组但上述方程组是非线性方程组, 求解十分困求解十分困难难.

51、一般利用一般利用 正交多项式正交多项式来求出来求出Gauss点与求积系点与求积系数数. Gauss型型求积公式求积公式 对对 f = x0, x1, x2, , x2n 1求积公式等号成立求积公式等号成立139 21、试确定常数A，B，C 和 ，使得数值积分公式 有尽可能高的代数精度，试问所得的数值积分公式代数精度是多少？它是否为Gauss型的？22( )()(0)( )f x dxAfBfCf 几种常用的几种常用的Gauss型型求积公式：求积公式： Gauss-Legendre 求积公式：求积公式：1)( x 定义在定义在 1, 1上，上，21( )(1)2!kkkkkdPxxkdx满足：满

52、足：2210(,)klkklP PklxPP 10, 1 ，递推公式：，递推公式：11(1)(21)kkkkPkxPkP其中，求积节点为其中，求积节点为 Pn+1 的根（求积系数通过解线性方程组的根（求积系数通过解线性方程组得到）。得到）。Legendre多项式：多项式：110( )()nkkkf x dxA f xGauss-Legendre 公式：公式：140区间区间a,b上的上的Gauss-Legendre 公式：公式：110( )222222bankkkbababaf x dxftdtbababaAft其中，其中， 为为 n+1次次Legendre多项式多项式Pn+1 的根。的根。kt141 Gauss-Chebyshev 求积公式：求积公式：211)(xx 定义在定义在 1, 1上，上，) arccos( cos)(xkxTk Tn+1 的根为：的根为： 2212cosnkxkk = 0, , n以此为节点构造公式以此为节点构造公式12100( )()()11nnkkkkkf xdxA f xf xnx称为称为 Gauss-Chebyshev 公式公式。注意到积分端点注意到积分端点 1 可能是积分可能是积分的的奇点奇点，用普通，用普通Newton-Cotes公式在端点会出问题。而公式在端点会出问题。而Gauss公式可能避免此问题的发生。