數(shù)值計算方法第07章數(shù)值微分與數(shù)值積分

1、1第七章第七章 數(shù)值微分與數(shù)值積分?jǐn)?shù)值微分與數(shù)值積分1 數(shù)值微分?jǐn)?shù)值微分2 NewtonCotes求積公式求積公式3 復(fù)化求積公式復(fù)化求積公式4 Romberg求積公式求積公式5 Gauss型求積公式型求積公式23 利用離散點上函數(shù)的信息求函數(shù)導(dǎo)數(shù)近似值利用離散點上函數(shù)的信息求函數(shù)導(dǎo)數(shù)近似值的方法的方法, 稱為稱為數(shù)值微分?jǐn)?shù)值微分. 差商型求導(dǎo)公式差商型求導(dǎo)公式 插值型求導(dǎo)公式插值型求導(dǎo)公式1 數(shù)值微分?jǐn)?shù)值微分4由導(dǎo)數(shù)定義由導(dǎo)數(shù)定義hxfhxfxfh)()(lim)(0 當(dāng)當(dāng)h很小時很小時, 可用可用差商差商近似導(dǎo)數(shù)近似導(dǎo)數(shù).5 差商型求導(dǎo)公式差商型求導(dǎo)公式 (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 二階導(dǎo)數(shù)的中心差商公式二階導(dǎo)數(shù)的中心差商公式2)()(2)()( hhxfxfhxfxf 截斷誤差截斷誤差)(12)()(2)()( )4(22 fhhhxfxfhxfxf 9101112數(shù)值積分?jǐn)?shù)值積分一、數(shù)值積分的必要性一、數(shù)值積分的必要性本章主要討論如下形式的一元函數(shù)積分本章主要討論如下形式的一元函數(shù)積分badxxffI)()( )( )( )( )baI ff x dxF bF a在微積分里，按在微積分里，按Newton-Leibniz公式公式求定積分求定積分 f x F x13 F x f x實際問題實際問題142

4、 sinf xx這個問題就是要求由函數(shù)這個問題就是要求由函數(shù)0 x 48x 15dxxdxxfL48024802)(cos1)(1 16類似的，下列函數(shù)也不存在由初等函數(shù)表示的原函數(shù)類似的，下列函數(shù)也不存在由初等函數(shù)表示的原函數(shù):2,1,ln1,sin,cos,sin322xexxxxxx173222xx)322ln(21693216332412222xxxxxx18 f x f xx1423454.5688.5原來通過原函數(shù)來原來通過原函數(shù)來計算積分有它的局計算積分有它的局限性。那限性。那怎么辦呢？怎么辦呢？呵呵呵呵這就需要積這就需要積分的數(shù)值方法來幫分的數(shù)值方法來幫忙啦。忙啦。19二、數(shù)值

5、積分的基本思想二、數(shù)值積分的基本思想1、定積分的幾何意義、定積分的幾何意義badxxffI)()(abxyo f x202、數(shù)值積分的理論依據(jù)、數(shù)值積分的理論依據(jù) f x依據(jù)依據(jù)積分中值定理積分中值定理, 對于連續(xù)函數(shù)對于連續(xù)函數(shù) ,在在 內(nèi)存在一點內(nèi)存在一點 ,使得使得, a b)()()()(fabdxxffIba稱稱 為為 在區(qū)間在區(qū)間 上的平均高度上的平均高度. f, a b ?f f x213、求積公式的構(gòu)造、求積公式的構(gòu)造 若簡單選取區(qū)間端點或中點的函數(shù)值作為平均高度，則若簡單選取區(qū)間端點或中點的函數(shù)值作為平均高度，則可得一點求積公式如下：可得一點求積公式如下：左矩形公式：左矩形公

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 一般地一般地 ，取區(qū)間，取區(qū)間 內(nèi)內(nèi) 個點個點, a b1n ,0,1,2,.,ixin ,0,1,.,if xin f處的高度處的高度通過通過加權(quán)平均加權(quán)平均的方法近似地得出平均高度的方法近似地得出平均高度這類求積方法稱為這類求積方法稱為機(jī)械求積機(jī)械求積:)()()(0ibaniixfabdxxf29 或?qū)懗苫驅(qū)懗? :數(shù)值積分公式數(shù)值積分公式求積系數(shù)求積系數(shù) 求積節(jié)點求積節(jié)點 )()(0kbankkxfAdxxf30記記0( )()nnkkkIfA f x0( )( )( )( )(),nbn

8、kkakR fI fIff x dxA f x稱稱為數(shù)值為數(shù)值求積公式求積公式稱為求積公稱為求積公式余項式余項(誤誤差差).（1）（2）31構(gòu)造或確定一個求積公式，要解決的問題包括構(gòu)造或確定一個求積公式，要解決的問題包括:(i) 確定求積系數(shù)確定求積系數(shù) 和求積節(jié)點和求積節(jié)點 kAkx；(iii) 求積公式的誤差估計和收斂性分析求積公式的誤差估計和收斂性分析.(ii) 確定衡量求積公式好壞的標(biāo)準(zhǔn)；確定衡量求積公式好壞的標(biāo)準(zhǔn)；3233 依據(jù)積分中值定理依據(jù)積分中值定理 就是說，底為就是說，底為b a 而高為而高為 f ( ) 的矩形面積恰恰等的矩形面積恰恰等于所求曲邊梯形于所求曲邊梯形 f (x

9、)的面積的面積.)()(abfdxxfba 取取a, b內(nèi)若干個節(jié)點內(nèi)若干個節(jié)點xk 處的高度處的高度 f (xk ), 通過加通過加權(quán)平均的方法生成平均高度權(quán)平均的方法生成平均高度 f ( ), 這類求積公式稱這類求積公式稱機(jī)械求積公式機(jī)械求積公式式中式中 xk 稱為稱為求積節(jié)點求積節(jié)點, Ak 稱為稱為求積系數(shù)求積系數(shù), 亦稱伴隨亦稱伴隨節(jié)點的權(quán)節(jié)點的權(quán). nkkkbaxfAdxxf0)()(數(shù)值積分基本思想數(shù)值積分基本思想342 Newton-Cotes 公式公式基本思想基本思想: 利用利用插值多項式插值多項式).()(xfxLn banbadxxLdxxffI.)()()(其中其中Ln

10、(x)是是n階階Lagrange插值多項式，用插值多項式，用Ln (x)的的積分近似積分近似 f (x)的積分，即的積分，即插值型求積公式插值型求積公式35 bankjjjkjkdxxxxxA,0)()(由由 決定決定,與與 無關(guān)無關(guān).節(jié)點節(jié)點 f (x) 在在a, b上取上取 a x0 x1 0, 使得使得|,|max0knk 則稱該求積公式是則稱該求積公式是穩(wěn)定穩(wěn)定的的. 求積公式的穩(wěn)定性求積公式的穩(wěn)定性72 若求積公式是穩(wěn)定的若求積公式是穩(wěn)定的, 則則 f (x)的觀察值的較小的的觀察值的較小的誤差引起的求積結(jié)果的誤差也是較小的誤差引起的求積結(jié)果的誤差也是較小的. 求積公式求積公式?jīng)]有把

11、沒有把 f (x)的誤差的誤差“放大放大”很多很多.73), 1 , 0()(nkfxfkk nkkkknnfxfAfIfI0)()()(證明證明因此復(fù)化梯形公式是數(shù)值穩(wěn)定的因此復(fù)化梯形公式是數(shù)值穩(wěn)定的.當(dāng)當(dāng)).(2)1(2abhhnh 定理定理 復(fù)化梯形公式是復(fù)化梯形公式是數(shù)值穩(wěn)定數(shù)值穩(wěn)定的的. nkkkkfxfA0)( nkkA0 74x0 x2xf (x)x4hhxn 2hxnmnnabh2, .hx3x1xn 1 復(fù)化復(fù)化Simpson公式公式分片二次多項式近似分片二次多項式近似75 將積分區(qū)間將積分區(qū)間a, b劃分為劃分為n=2m等分等分, 步長步長 h=( b a )/n, 分點分

12、點 xk= a+kh ( k=0, 1, , n). 在每個在每個小區(qū)間小區(qū)間 x2k 2 , x 2k ( k=1, , m)上用上用Simpson公式：公式： kkxxkkkkkxfxfxfxxdxxf222)()(4)(6)(21222222 )()(4)(321222kkkxfxfxfh 復(fù)化復(fù)化Simpson公式公式k=1, , m76 111122)(4)(2)()(3mkmkkkxfxfbfafh= Sn( f ) mkxxbadxxfdxxffIkk1222)()()( )()(4)(3212221kkkmkxfxfxfh 復(fù)化復(fù)化Simpson公式公式77)()()(fSfI

13、fRnn 當(dāng)當(dāng) f (x)在在a, b上具有四階連續(xù)導(dǎo)數(shù)時上具有四階連續(xù)導(dǎo)數(shù)時, ),(),()(1)4()4(bamffmkk ),(222kkkxx 故得故得),(),(180)()(90)()4(4)4(5bafhabfmhfRn 復(fù)化復(fù)化Simpson公式的截斷誤差公式的截斷誤差 mkkfh1)4(5),(90 78 由復(fù)化由復(fù)化Simpson公式的截斷誤差知公式的截斷誤差知, 誤差階為誤差階為 h4, 收斂性是顯然的收斂性是顯然的, 事實上事實上,只要只要 f (x) Ca, b則則可得到可得到收斂性收斂性, 即即.)()(limdxxffSbann 由于求積系數(shù)均為正由于求積系數(shù)均

14、為正, 與復(fù)化梯形公式一樣的與復(fù)化梯形公式一樣的證法可得復(fù)化證法可得復(fù)化 Simpson公式是公式是數(shù)值穩(wěn)定數(shù)值穩(wěn)定的的.79例：例：計算計算dxx 10214 解：解： )1()(2)0(161718fxffTkk8kxk 其中其中= 3.138988494 )1()(2)(4)0(241oddeven8fxfxffSkk8kxk 其中其中= 3.141592502運算量運算量基本相同基本相同 顯然用復(fù)化顯然用復(fù)化Simpson公式計算精度較高公式計算精度較高, 這與它們這與它們的誤差階的結(jié)論是相符的的誤差階的結(jié)論是相符的.80例例 對于函數(shù)對于函數(shù),sin)(xxxf 給出給出n=8的函數(shù)

15、表的函數(shù)表, 試用試用復(fù)化梯形公式及復(fù)化復(fù)化梯形公式及復(fù)化Simpson公式計算積分公式計算積分 10.sindxxxI解解ix)(ixf0.00.1250.250.3750.50.6250.750.8751.01.00.99739780.98961580.97672670.95885100.93615560.90885160.84147090.8771925應(yīng)用復(fù)化梯形公式求得應(yīng)用復(fù)化梯形公式求得T8=0.9456909應(yīng)用復(fù)化應(yīng)用復(fù)化Simpson公式求得公式求得S8=0.9460832準(zhǔn)確值準(zhǔn)確值 I=0.9460831兩者運算量基本相同兩者運算量基本相同8126.26.利用下面數(shù)據(jù)表，

16、利用下面數(shù)據(jù)表， 1. 用復(fù)化梯形公式計算積分用復(fù)化梯形公式計算積分的近似值；的近似值；2. 用復(fù)化用復(fù)化Simpson公式計算積分公式計算積分的近似值。的近似值。( (要求計算結(jié)果保留到小數(shù)點后六位要求計算結(jié)果保留到小數(shù)點后六位). 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: 復(fù)化梯形公式求積分復(fù)化梯形公式求積分.用法用法: trapz(X, Y), 其中其中X, Y為相同維數(shù)的向量為相同維數(shù)的向量.例例: X=0.125:0.125:

17、1.0;Y=sin(X)./X;X=0,X;Y=1,Y;trapz(X,Y)ans = 0.94569086358270Matlab函數(shù)函數(shù)83例例 若用復(fù)化求積公式計算積分若用復(fù)化求積公式計算積分dxeIx 10的近似值的近似值, 若要求計算結(jié)果有若要求計算結(jié)果有4位有效數(shù)字位有效數(shù)字, n應(yīng)取多大應(yīng)取多大?解解, 1110 dxeIex.105 . 04 1 , 0, 1| )(|)( xexfxk復(fù)化梯形公式的誤差復(fù)化梯形公式的誤差)( )(12|2 fabhRT .83.40 n若用復(fù)化梯形公式求積分若用復(fù)化梯形公式求積分, n取取41能達(dá)到精度要求能達(dá)到精度要求.2121n 4102

18、1 84故應(yīng)取故應(yīng)取n=4. 該例表明該例表明, 為達(dá)到相同的精度為達(dá)到相同的精度, 用復(fù)化用復(fù)化Simpson公式所需的計算量比復(fù)化梯形公式要少公式所需的計算量比復(fù)化梯形公式要少, 這也說明這也說明了復(fù)化了復(fù)化Simpson公式的精度高公式的精度高.復(fù)化復(fù)化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 復(fù)化復(fù)化梯形公式的逐次分半算法梯形公式的逐次分半算法將區(qū)間將區(qū)間a, b分成分成n=2m等分等分, 記記,2mmabh ), 2 , 1 , 0( m883122312422)(2)()(2)(21)()(21(44 , 3 , 2 , 1 , 0 ,4,4)(3kkkkhafbfafhbfkhafafabTkkabakhaxabh再分半、加密一次區(qū)間)3()(21 )4)( 3()4(421)4)( 3()4(4)4)( 2(2)(2)(4222222241212

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

21、2所有新增加節(jié)點的函數(shù)值之和所有新增加節(jié)點的函數(shù)值之和其中其中.2mmabh 復(fù)化復(fù)化梯形公式的逐次分半算法梯形公式的逐次分半算法92以以n=8, m=3為例為例. 記記 fk= f (xk)x0 x2x4x6x3x1x5x7x8 )(2167654321808fffffffffabT )(21664280fffffab 75318ffffab 24T 所有新增加節(jié)點的函數(shù)值之和所有新增加節(jié)點的函數(shù)值之和. 3h93 復(fù)化梯形公式余項的后驗估計復(fù)化梯形公式余項的后驗估計);,(),(12112bafhabTIn );,(),(2122222bafhabTIn f ( 1 ), f ( 2 )

22、分別是分別是 f (x) 在在a, b上的上的n個點與個點與 2n 個點處的算術(shù)平均值個點處的算術(shù)平均值 (每個小區(qū)間上取一個點每個小區(qū)間上取一個點). 當(dāng)當(dāng)n較大時較大時, 有有.)( 1)()(21dxxfabffba 94因此因此, 若事先給定誤差限若事先給定誤差限 , 則當(dāng)則當(dāng).3|2 nnTT時時, 就可停止計算就可停止計算, 并認(rèn)為并認(rèn)為 T2n是滿足精度要求的近是滿足精度要求的近似值似值.;412 nnTITI);()(21 ff ).(3122nnnTTTI 95 復(fù)化復(fù)化Simpson公式的逐次分半算法公式的逐次分半算法將區(qū)間將區(qū)間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 因此因此, 若事先給定誤差限若事先給定誤差限 , 則當(dāng)則當(dāng).15|2 nnSS時可停止計算時可停止計算, 取取 S2n為滿足精度要求的近似值為滿足精度要求的近似值. 復(fù)化復(fù)化Simpson公式余項的后驗估計公式余項的后驗估計9722、已知函數(shù)f(

24、x)在若干點處的值： 試計算積分 的梯形值 以及Simpson值 。X-1-0.500.51F(x)02.12532.125011( )f x dx124,T T T12,S S表表達(dá)達(dá)式式知知：由由復(fù)復(fù)化化梯梯形形公公式式的的余余項項，但但收收斂斂慢慢，精精度度低低。復(fù)復(fù)化化梯梯形形公公式式算算法法簡簡單單).,(),( )2(12);,(),( 12222bafhabTIbafhabTInn )( )( ff 假假定定412 nnTITI 7.4 7.4 龍貝格算法龍貝格算法).(3122nnnTTTI 1443134)(31:2222nnnnnnnTTTTTTTI于于是是的的誤誤差差很很

25、小小很很小小，可可保保證證若若nnnTTT22事后估計法利用計算結(jié)果來估計誤差的方法利用計算結(jié)果來估計誤差的方法龍貝格算法龍貝格算法 當(dāng)當(dāng)n=1 =1 時時, ,我們計算上式右端我們計算上式右端 )(21)(21)(31)(21)2()(212)(3414412bfafabbfbafafabTTT1)(61)2(64)(61)(Sbfbafafab 這恰好是辛普森公式的結(jié)果，即有這恰好是辛普森公式的結(jié)果，即有121141144TTS 比梯形公式有比梯形公式有更好的精確度更好的精確度龍貝格算法龍貝格算法類似地可驗證：類似地可驗證：nnnTTSTTS1411441411442242 即即nnnTT

26、S31342 辛普森積分值辛普森積分值_nS龍貝格算法龍貝格算法式式的的余余項項注注意意復(fù)復(fù)化化辛辛普普森森求求積積公公)()(21804)4(4hOfhabSIRnn ) (4180)4(422 fhabSIRnn )()()4()4( ff 假假定定144)(15122222 nnnnnSSSSSI)(,215121612nnnSISISSSInn 龍貝格算法龍貝格算法可以驗證可以驗證nnnSSC15115162 復(fù)復(fù)化化柯柯特特斯斯積積分分值值_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求積公式求積公式啟示啟示: 是否用是否用 復(fù)化梯形公式余項的后驗估計表明復(fù)化梯形公式余項的后驗估計表明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 則它恰為復(fù)化則它恰為復(fù)化

30、Cotes公式公式; 且有如下誤差估計式且有如下誤差估計式).(| )(|62hOfICn 復(fù)化復(fù)化Simpson公式余項的后驗估計表明公式余項的后驗估計表明 問題問題: 是否用是否用108.631636422nnnRCC 類似地可以得到類似地可以得到2nR其中其中被稱為被稱為Romberg序列序列.).(| )(|82hOfIRn 截斷誤差截斷誤差:龍貝格龍貝格積分積分 /* Romberg Integration */例：例：計算計算dxx 10142 已知對于已知對于 = 10 6 須將區(qū)間對分須將區(qū)間對分 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 停機(jī)準(zhǔn)則停機(jī)準(zhǔn)則:梯形值序列梯形值序列Simpson序列序列Cotes序列序列Romberg序列序列 Romberg求積公式求積公式 用若干個積分近似值推算出更為精確的積分近用若干個積分近似值推算出更為精確的積分近似值的方法，稱為似值的方法，稱為外推方法外推方法。 列和龍貝格序列。列和龍貝格序列。辛普森序列、柯特斯序辛普森序列、柯特斯序分別稱為梯形序列、分別稱為梯形序列、序列序列NNNNRCST,龍龍貝貝格格序序列列為為止止。故故通通常常到到及及續(xù)續(xù)外外推推，但但因因由由龍龍貝貝格格序序列列還還可可以以繼繼),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 利用只有兩三位有效數(shù)字利用只有兩三位有效數(shù)字的的T1, ,T8 經(jīng)過三次外推得經(jīng)過三次外推得到到7位有效數(shù)字位有效數(shù)字. 可見加速的可見加速的效果十分顯著效果十分顯著.用用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)(設(shè)設(shè)k

39、kabhba22, 等等分分，將將區(qū)區(qū)間間)(hTTn記記為為近近似似值值為為用用復(fù)復(fù)化化梯梯形形公公式式求求得得的的記記.)(634221 hahahaIhT則則)()( 12)(:22hOfhabhTI注意到注意到 .212121)2(636424212 hahahaIhT將余項泰將余項泰勒展開勒展開)2(,)0()(lim20hTTIThTnh 而而 外推法的一般討論外推法的一般討論.3)()2(4)(62411 hbhbIhThThT)()(41hOIhT 即即,),2(),(211nnSShThT就就是是辛辛普普森森序序列列如如此此構(gòu)構(gòu)造造序序列列)(151)2(1516)(,16)

40、(11262411hThThThhIhT 若若令令由由 外推法外推法82612)(hhIhT ,)(22nnCChT就就是是柯柯特特斯斯序序列列這這樣樣構(gòu)構(gòu)造造的的序序列列)()(62hOIhT )2(22)1(21)(), 2 , 1(2mmmhhIhTmm 列列形形式式：次次加加速速后后，余余項項便便取取下下經(jīng)經(jīng)過過階階，高高一一次次，誤誤差差的的量量級級便便提提如如此此繼繼續(xù)續(xù)下下去去，每每加加速速此處理方法稱為理查森此處理方法稱為理查森(Richardson)(Richardson)外推加速方法外推加速方法。外推法外推法123 理論依據(jù)理論依據(jù): 復(fù)化梯形公式的余項展開復(fù)化梯形公式的余

41、項展開.記記),(hTTn 定理定理 設(shè)設(shè),)(baCxf 則則 kkhhhIhT24221)( 其中系數(shù)其中系數(shù) k ( k=0, 1, )是與是與 h 無關(guān)的常數(shù)無關(guān)的常數(shù). T (h) 逼近逼近 I 的速度是的速度是 O ( h2 )階階.124,3)()2(4)(262411 kkhhhIhThThT 當(dāng)區(qū)間當(dāng)區(qū)間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公式序列進(jìn)行同樣處理得到公式序列進(jìn)行同樣處理得到Romberg公式序列公式序

43、列.Richardson外推加速方法外推加速方法 也稱為也稱為Romberg求積算法求積算法 收斂性說明收斂性說明: 如果如果 f (x) 充分光滑充分光滑, 那么梯形公那么梯形公式序列式序列, Simpson公式序列公式序列, Cotes公式序列公式序列, Romberg公式序列均收斂到所求的積分值公式序列均收斂到所求的積分值. 對于對于 f (x)不充分光滑的函數(shù)也可用不充分光滑的函數(shù)也可用Romberg算算法計算法計算, 只是收斂慢一些只是收斂慢一些. 也可以直接使用復(fù)化也可以直接使用復(fù)化Simpson公式計算公式計算.127例例 用用Romberg算法計算積分算法計算積分.1023dx

44、xI 解解23)(xxf 在在0, 1上僅是一次連續(xù)可微上僅是一次連續(xù)可微用用Romberg算法計算結(jié)果見下表算法計算結(jié)果見下表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、式型求積公式 基本思想基本思想 設(shè)計求積公式設(shè)計求積公式:在節(jié)點數(shù)在節(jié)點數(shù) n 固定時固定時, 適當(dāng)?shù)剡x取求積節(jié)點適當(dāng)?shù)剡x取求積節(jié)點 xk 與求與求積系數(shù)積系數(shù) Ak , 使求積公式具有使求積公式具有最高的代數(shù)精確度最高的代數(shù)精確度.130例例 確定確定x1, x2, A1, A2, 使求積公式使求積公式)()()(221111xfAxfAdxxf 具有最高次的代數(shù)精確度具有最高次的代數(shù)精確度.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 積分精確成立積分精確成立 四個方程四個未知數(shù)四個方程四個未知數(shù)132梯形公式與梯形公式與Gauss求積公式的比較求積公式的比較 對對1, x求積公式精確求積公式精確成立成立(1(1次代數(shù)精確度次代數(shù)精確度) ) 對對1, x, x2, x3求積公式精求積公式精確成立確成立(3 次代數(shù)精確度次代數(shù)精確度)133)()()()( :333

47、221111xfAxfAxfAdxxfn x3x11x2 1 選取選取(A1, A2 , A3 , x1, x2 , x3) 使該求積公式對使該求積公式對 f (x) = x0, x1, x2, x3, x4, x5 時等號成立時等號成立.區(qū)間區(qū)間 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時等號成立時等號成立. .區(qū)間區(qū)間 1, 1上的上的Gauss求積公式求積公式136 0 121 32 0 21 11212211212111212222221122211222222211

49、2211222112111121 nnnnnnnnnnnnnnnnnnnxAxAxAdxxxfxAxAxAndxxxfxAxAxAdxxxfxAxAxAxdxxfAAAdxf 2n個方程個方程2n個未知數(shù)個未知數(shù), 非線性方程非線性方程, 解是否存在唯解是否存在唯一？若存在一？若存在, 如何求解？如何求解？ 對對 f = x0, x1, x2, , x2n 1求積公式等號成立求積公式等號成立137. )()()(1 nkkkbaxfAdxxfx 研究最一般情形的帶權(quán)積分研究最一般情形的帶權(quán)積分:具有具有2n1次代數(shù)精確度次代數(shù)精確度, 則稱這組節(jié)點則稱這組節(jié)點 xk為為 Gauss 點點, 上

50、述公式稱為帶權(quán)函數(shù)上述公式稱為帶權(quán)函數(shù) (x)的的Gauss型型求求積公式積公式.定義定義 如果一組節(jié)點如果一組節(jié)點x1, x2, , xn a, b能使求積公式能使求積公式 nkkkbaxfAdxxfx1)()()( Gauss型型求積公式求積公式138. 12, 1 , 0,)(1 nmxAdxxxmknkkbam 由上述由上述 2n個方程確定全部個方程確定全部 的的2n個待定參數(shù)個待定參數(shù) xk , Ak ( k=1, , n), 使求積公式至少具有使求積公式至少具有 2n 1次代數(shù)精次代數(shù)精確度確度. 但上述方程組是非線性方程組但上述方程組是非線性方程組, 求解十分困求解十分困難難.

51、一般利用一般利用 正交多項式正交多項式來求出來求出Gauss點與求積系點與求積系數(shù)數(shù). Gauss型型求積公式求積公式 對對 f = x0, x1, x2, , x2n 1求積公式等號成立求積公式等號成立139 21、試確定常數(shù)A，B，C 和 ，使得數(shù)值積分公式 有盡可能高的代數(shù)精度，試問所得的數(shù)值積分公式代數(shù)精度是多少？它是否為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其中，求積節(jié)點為其中，求積節(jié)點為 Pn+1 的根（求積系數(shù)通過解線性方程組的根（求積系數(shù)通過解線性方程組得到）。得到）。Legendre多項式：多項式：110( )()nkkkf x dxA f xGauss-Legendre 公式：公式：140區(qū)間區(qū)間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以此為節(jié)點構(gòu)造公式以此為節(jié)點構(gòu)造公式12100( )()()11nnkkkkkf xdxA f xf xnx稱為稱為 Gauss-Chebyshev 公式公式。注意到積分端點注意到積分端點 1 可能是積分可能是積分的的奇點奇點，用普通，用普通Newton-Cotes公式在端點會出問題。而公式在端點會出問題。而Gauss公式可能避免此問題的發(fā)生。