龍貝格求積法

1、4.4 外推原理與外推原理與Romberg求積方法求積方法4.4.1 外推原理外推原理 在科學(xué)與工程計(jì)算中，很多算法與步長(zhǎng)h有關(guān)，特別是數(shù)值積分、數(shù)值微分和微分方程數(shù)值解的問(wèn)題。對(duì)于這些算法，我們可以通過(guò)外推技巧提高計(jì)算精度。例1 計(jì)算的近似值。sin x由函數(shù) 的Taylor展開(kāi)式有3524sin3!5!nnnn,( )6sin66hF h若記 ，則有24( ),6120F hhh2411( ),26 412016hFhh由此構(gòu)造新的表達(dá)式414 ( )( )12( )3120 4hFF hF hh 可見(jiàn)計(jì)算的近似值的算法F(h)的截?cái)嗾`差是O (h2),而算法F1(h)的截?cái)嗾`差是O (h

2、4)。外推一次，精度就提高了。這就是外推法的基本思想。若重復(fù)以上過(guò)程，不斷外推，即不斷折半步長(zhǎng)h, 得到計(jì)算的算法序列Fk(h)。隨著k的增加，算法的截?cái)嗾`差階越來(lái)越高，計(jì)算精度越來(lái)越好。 可將上述外推思想推廣到一般情況。設(shè)F(h)是計(jì)算F(0)的一種近似算式，帶截?cái)嗾`差的表達(dá)式為( )(0)(),pspF hFa hO hsppa( )(0)( )(),psphhFFaO hqq1( )( )( )(0)(),1psphq FF hqF hFO hq().sO h其中， 與h無(wú)關(guān)。如果我們用h和h/q (q1)兩種步長(zhǎng)分別計(jì)算F(h)和F(h /q)，則有削去截?cái)嗾`差的主項(xiàng)，得新的算法我們稱

3、這個(gè)計(jì)算過(guò)程為Richardson外推法。 F1(h)逼近F(0)的截?cái)嗾`差是 只要知道F(h)的更加完整的關(guān)于h冪的展開(kāi)式，而無(wú)需知道展開(kāi)式中各個(gè)系數(shù)的具體數(shù)值，就能重復(fù)使用Richardson外推法，直到截?cái)嗾`差達(dá)到容許誤差。用歸納法可以證明下面更一般的定理。定理4.5 假設(shè)F(h)逼近F(0)的余項(xiàng)為312123( )(0),pppF hFa ha ha h其中，123,(0,1,2,)kpppa k是與h無(wú)關(guān)的非零常數(shù)。則由01()( )( )( ),( ),0,1,2,1kkpkkkphq FF hpF hF h Fhkq定義的序列Fk(h)有123( )( )( )123( )(0

4、),nnnpppnnnknnnF hFahahah( )(1,2,3,)nn kakh其中與無(wú)關(guān)，q1. Richardson外推法應(yīng)用非常廣泛且有效，下面介紹應(yīng)用于數(shù)值積分的情形。 變步長(zhǎng)梯形求積法變步長(zhǎng)梯形求積法算法簡(jiǎn)單，但精度較差，收斂速度較慢，算法簡(jiǎn)單，但精度較差，收斂速度較慢，但可以利用梯形法算法簡(jiǎn)單的優(yōu)點(diǎn)，形成一個(gè)新算法，這就是龍但可以利用梯形法算法簡(jiǎn)單的優(yōu)點(diǎn)，形成一個(gè)新算法，這就是龍貝格求積公式。貝格求積公式。龍貝格公式龍貝格公式又稱又稱逐次分半加速法逐次分半加速法。 Romberg算法算法是在積分區(qū)間逐次分半的過(guò)程中，對(duì)用復(fù)化梯是在積分區(qū)間逐次分半的過(guò)程中，對(duì)用復(fù)化梯形法形法產(chǎn)

5、生的近似值產(chǎn)生的近似值進(jìn)行進(jìn)行加權(quán)平均加權(quán)平均，以獲得準(zhǔn)確度較高的近似值的，以獲得準(zhǔn)確度較高的近似值的一種方法，具有公式簡(jiǎn)練，使用方便，結(jié)果較可靠的優(yōu)點(diǎn)。一種方法，具有公式簡(jiǎn)練，使用方便，結(jié)果較可靠的優(yōu)點(diǎn)。二、二、Romberg算法算法 根據(jù)梯形法的誤差公式，可知積分值根據(jù)梯形法的誤差公式，可知積分值 的截?cái)嗾`差大致與的截?cái)嗾`差大致與 成正比，因此當(dāng)步長(zhǎng)二分后，截?cái)嗾`差將減至原有誤差的成正比，因此當(dāng)步長(zhǎng)二分后，截?cái)嗾`差將減至原有誤差的1/4，即有即有 將上式移項(xiàng)整理，可得將上式移項(xiàng)整理，可得nT2h（4.4.3） 由此可見(jiàn)，只要二分前后的兩個(gè)積分值由此可見(jiàn)，只要二分前后的兩個(gè)積分值 與與 相當(dāng)

6、接近，就可以保相當(dāng)接近，就可以保證證 計(jì)算結(jié)果的誤差很小。這樣直接用計(jì)算結(jié)果來(lái)估計(jì)誤差的方法計(jì)算結(jié)果的誤差很小。這樣直接用計(jì)算結(jié)果來(lái)估計(jì)誤差的方法通常稱作通常稱作誤差的事后估計(jì)法誤差的事后估計(jì)法。 nT2nT2nT按式按式(4.4.3),積分近似值與積分近似值與 的誤差大致等于的誤差大致等于 因此如果因此如果用這個(gè)誤差值作為用這個(gè)誤差值作為 的的一種補(bǔ)償一種補(bǔ)償，可以期望所得到的，可以期望所得到的 213nnTT2nT2nT（4.4.4）可能是更好的結(jié)果。可能是更好的結(jié)果。nnnTTT31342(4.4.5) 有可能比有可能比 更好地接近于積分更好地接近于積分 的真值的真值 InT2 dxxf

7、ba)(3122nnnnTTTT即即這就是說(shuō)，用梯形法二分前后兩個(gè)積分值這就是說(shuō)，用梯形法二分前后兩個(gè)積分值 和和 作線性組合作線性組合nTnT2與與這表明在收斂緩慢的梯形數(shù)列這表明在收斂緩慢的梯形數(shù)列 的基礎(chǔ)上，若將的基礎(chǔ)上，若將 kT2nT2nT按按(4.4.6)作線性組合就可產(chǎn)生收斂速度較快的作線性組合就可產(chǎn)生收斂速度較快的Simpson序列：序列：kS21S：、2S、關(guān)于關(guān)于(4.4.3.6)的證明：由的證明：由(4.4.1)可知：可知：這種新近似值這種新近似值 實(shí)質(zhì)上又是什么呢？可以驗(yàn)證：實(shí)質(zhì)上又是什么呢？可以驗(yàn)證：nTnnST 即：即：144313422nnnnnTTTTS (4.

8、4.6)4.S31211223412nknnnnabkafnabTTT 故（故（4.4.6）式成立）式成立. 1111122322nnkkkhfafxf bhfakh 111102246nnknkkkhfafxfxf bS(由上節(jié)梯形公式由上節(jié)梯形公式)（由上節(jié)（由上節(jié)Simpson公式）公式）nnnSSSI22151144151151622222nnnnnSSSSC同理，由上節(jié)近似式同理，由上節(jié)近似式類似推導(dǎo)可得：類似推導(dǎo)可得：（4.4.7）即將即將Simpson序列序列 按按(4.4.7) 作線性組合就可產(chǎn)生收斂作線性組合就可產(chǎn)生收斂kS2速度更快的新序列速度更快的新序列1242:,kCC

9、 C C -Cotes序列序列即在即在Cotes序列序列 的基礎(chǔ)上，產(chǎn)生了一個(gè)稱為的基礎(chǔ)上，產(chǎn)生了一個(gè)稱為RombergkC2 :2kR1R、2R、4R.序列序列的新序列的新序列（4.4.8）14463163643232nnnnnCCCCR 由近似式 ， 類似推導(dǎo)可得：)(63122nnnCCCI 我們?cè)谧儾介L(zhǎng)的過(guò)程中運(yùn)用公式我們?cè)谧儾介L(zhǎng)的過(guò)程中運(yùn)用公式(4.4.4-4.4.8)(也稱它們?yōu)橐卜Q它們?yōu)榧铀俟郊铀俟? ，就能將粗糙的梯形值，就能將粗糙的梯形值Tn逐步加工成精度較高的辛逐步加工成精度較高的辛普森值普森值Sn 、柯特斯值、柯特斯值Cn和龍貝格值和龍貝格值Rn .可以證明：可以證明

10、：當(dāng)當(dāng)f（x）滿足一定條件時(shí)，）滿足一定條件時(shí)，Romberg序列序列 kR2比比Cotes序列序列 能更快地收斂到積分能更快地收斂到積分 的真值的真值I。kC2 dxxfba越精確的近似值越精確的近似值 也就是將收斂速度緩慢也就是將收斂速度緩慢綜上可知：在積分區(qū)間逐次分半過(guò)程中利用公式綜上可知：在積分區(qū)間逐次分半過(guò)程中利用公式可以將粗糙的近似值可以將粗糙的近似值 逐步地逐步地“加工加工”成越來(lái)成越來(lái)的梯形序列的梯形序列 逐步地逐步地“加工加工”成收斂速度越來(lái)越快的新成收斂速度越來(lái)越快的新 kT2新序列新序列 4,3,2,nnnRCS,222kkkRCS(4.4.3-4.4.8) 這種加速的方

11、法稱為這種加速的方法稱為 Romberg算法。其算法。其“加工加工”過(guò)程如下圖，過(guò)程如下圖，其中圓圈中號(hào)碼表示計(jì)算順序。其中圓圈中號(hào)碼表示計(jì)算順序。 1T2T4T8T16T1S2S4S8S1C2C4C1R2R22241161641,3315156363 1,2,4,8,.nnnnnnnnnSTT CSSRCCnT1,T2、T4、T6 由梯由梯形公式、復(fù)化梯形公式形公式、復(fù)化梯形公式的遞推化公式求得的遞推化公式求得kk2kT212 kS22 kC32 kR0 20=1 T11 21=2 T2 S12 22=4 T4 S2 C13 23=8 T8 S4 C2 R14 24=16 T16 S8 C4

12、 R25 25=32 T32 S16 C8 R4 區(qū)間等分?jǐn)?shù)區(qū)間等分?jǐn)?shù) 梯形序列梯形序列 辛普森序列辛普森序列 柯特斯序列柯特斯序列 龍貝格序列龍貝格序列 龍貝格求積算法也可用下表來(lái)表示：龍貝格求積算法也可用下表來(lái)表示： 例例1 1: 利用利用Romberg算法計(jì)算算法計(jì)算.,8421TTTT214)(, 1, 0 xxfba3)24(21) 1 ()0(211ffT1 . 351621321)(21212112fTT解解: :由題意由題意12041dxx計(jì)算到計(jì)算到R1.121413.13333333STT3142441111( )( )3.1(3.7647062.56)3.13117724

13、24TTff242413.14156933STT357184888811( )( )( )( )3.13898928TTffff1211613.1421181515CSS484413.14159333STT121413.13333333STT242413.14156933STT2421613SS1216413.1415862926363RCC 例例2 用用Romberg算法計(jì)算算法計(jì)算 得到的梯形值，計(jì)得到的梯形值，計(jì)算結(jié)果見(jiàn)算結(jié)果見(jiàn) 表表4-5(k代表二分次數(shù)代表二分次數(shù))。計(jì)算值的誤差不超過(guò)。計(jì)算值的誤差不超過(guò)0.5 10-6.表表4-5我們看到，這里利用二分我們

14、看到，這里利用二分3次的數(shù)據(jù)次的數(shù)據(jù)(它們的精度都很差，只有二三位它們的精度都很差，只有二三位是有效數(shù)字是有效數(shù)字)，通過(guò)三次加速求得，通過(guò)三次加速求得 =0.9460831，這個(gè)結(jié)果的每一位，這個(gè)結(jié)果的每一位數(shù)字都是有效數(shù)字，可見(jiàn)加速的效果是十分顯著的。三次外推后達(dá)到數(shù)字都是有效數(shù)字，可見(jiàn)加速的效果是十分顯著的。三次外推后達(dá)到6位有效數(shù)字。位有效數(shù)字。10sin.xIdxx注意注意 教材中介紹的教材中介紹的Richardson外推法外推法，為便于上機(jī)計(jì)算，引用記，為便于上機(jī)計(jì)算，引用記號(hào)號(hào) 來(lái)表示各近似值，其中來(lái)表示各近似值，其中k仍代表積分區(qū)間的二分次數(shù)，而仍代表積分區(qū)間的二分次數(shù)，而下標(biāo)

15、下標(biāo)m則指出了近似值則指出了近似值 所在序列的性質(zhì)。如所在序列的性質(zhì)。如m=1在梯形序列中，在梯形序列中，m=2在在Simpson序列中，序列中，m=3在在Cotes序列中，序列中， 引入上面記號(hào)引入上面記號(hào)后，后，Romberg算法所用到的各個(gè)計(jì)算公式可統(tǒng)一化為：算法所用到的各個(gè)計(jì)算公式可統(tǒng)一化為： kmT kmT三、三、Romberg 算法計(jì)算公式的簡(jiǎn)化算法計(jì)算公式的簡(jiǎn)化(0)1( )( )2baTf af b211111121() ,1,2,.222iiiiijbajTTf abai( )(1)(1)14,1,2,.,1,2,.,41mkkkmmmmTTTmki由此可逐行構(gòu)造由此可逐行構(gòu)

16、造 出一個(gè)三角形數(shù)表出一個(gè)三角形數(shù)表-稱為稱為T(mén)數(shù)表數(shù)表 ki0123Romberg算法中止準(zhǔn)則，一般取同列或同行相鄰兩數(shù)值的誤差算法中止準(zhǔn)則，一般取同列或同行相鄰兩數(shù)值的誤差絕對(duì)絕對(duì)實(shí)際計(jì)算中常常只計(jì)算到第實(shí)際計(jì)算中常常只計(jì)算到第4列列,只使用只使用3次理查遜外推法。次理查遜外推法。 注：注：值小于事先給定的精度要求。取最后一次的數(shù)值積分值作為積值小于事先給定的精度要求。取最后一次的數(shù)值積分值作為積分的分的近似值。近似值。 外推次數(shù)外推次數(shù)分分半半次次數(shù)數(shù)復(fù)化梯復(fù)化梯形序列形序列Simpson序列序列Cotes序列序列Romberg序列序列 T1 = T8 = T4 = T2 = ? ? ?

17、 S1 = R1 = S2 = C1 = C2 = S4 =1kT234 kkkTTT01112131TTTT021222TTT0313TT04T1、在上面、在上面“加工加工”過(guò)程中的系數(shù)過(guò)程中的系數(shù) 和和 ，當(dāng)，當(dāng) m 4 時(shí)，時(shí)， 而另一個(gè)而另一個(gè) 系數(shù)則接近于系數(shù)則接近于1，也就是，也就是144mm141m2551141m新公式與原公式差別不大，新公式與原公式差別不大， 但工作量卻大增。但工作量卻大增。 kmkmTT1因此，在實(shí)際計(jì)算中常規(guī)定因此，在實(shí)際計(jì)算中常規(guī)定 m 3，即計(jì)算到出現(xiàn)，即計(jì)算到出現(xiàn)Romberg序列為止。序列為止。 2、 可用二維數(shù)組來(lái)存放并參加運(yùn)算，也可用一維數(shù)組。

18、可用二維數(shù)組來(lái)存放并參加運(yùn)算，也可用一維數(shù)組。 kmT四、幾點(diǎn)說(shuō)明：四、幾點(diǎn)說(shuō)明：3、 對(duì)于積分限為無(wú)窮的積分對(duì)于積分限為無(wú)窮的積分 ，可利用變量代換化成有限區(qū)，可利用變量代換化成有限區(qū) 間的積分然后再進(jìn)行計(jì)算。例如：間的積分然后再進(jìn)行計(jì)算。例如：4、若被積函數(shù)有奇異點(diǎn)（間斷點(diǎn)）存在于積分區(qū)間內(nèi)，則、若被積函數(shù)有奇異點(diǎn)（間斷點(diǎn)）存在于積分區(qū)間內(nèi)，則 可將積分可將積分 區(qū)間分成小部分，使間斷點(diǎn)在子區(qū)間的端點(diǎn)處。區(qū)間分成小部分，使間斷點(diǎn)在子區(qū)間的端點(diǎn)處。 也可用變量代換法處理。也可用變量代換法處理。 1321xxdxtx1令令則則代入得：代入得：1030132213211111tttdttttdxxxdx21dxdtt 例例3 用龍貝格方法計(jì)算橢圓用龍貝格方法計(jì)算橢圓 x2/4 + y2 l 的周長(zhǎng)，使結(jié)果的周長(zhǎng)，使結(jié)果具有五位有效數(shù)字具有五位有效數(shù)字 分析分析 為便于計(jì)算，先將橢圓方程采用參數(shù)形式表示為便于計(jì)算，先將橢圓方程采用參數(shù)形式表示, ,再根再根據(jù)弧長(zhǎng)公式將橢圓周長(zhǎng)用積分形式表示由于計(jì)算結(jié)果要求具據(jù)弧長(zhǎng)公式將橢圓周長(zhǎng)用積分形式表示由于計(jì)算結(jié)果要求具有五位有效數(shù)字，因此需要估計(jì)所求積分值有幾位整數(shù)，從而有五位有效數(shù)字，因此需要估計(jì)所求積分值有幾位整數(shù)，從而確定所求積分值的絕對(duì)誤差限最后再應(yīng)