1、.11 5.3 Romberg算法算法 綜合前幾節的內容,我們知道梯形公式,Simpson公式,Cotes公式的代數精度分別為1次,3次和5次復合梯形、復合Simpson、復合Cotes公式的收斂階分別為2階、4階和6階無論從代數精度還是收斂速度,復合梯形公式都是較差的有沒有辦法改善梯形公式呢?.22一、復合梯形公式的遞推化等份分割為的積分區間將定積分nbadxxfIba,)(njjhaxk, 1 ,0,nabh各節點為 )()(2)(211njjnbfxfafnabT復合梯形(Trapz)公式為則不變等份，而分割為如果將,/)(2,nabhnba)()(2)(2)(41021112bfxfx

2、fafnabTnjjnjjn-(1)-(2).33)()(2)(2)(41021112bfxfxfafnabTnjjnjjnhjahxxjj)21(2121其中102111)(24)()(2)(4njjnjjxfnabbfxfafnab1021)(221njjnxfnabT10)21(221njnhjafnabT-(3)10)2)12(221njnnabjafnabT.44abhn 時，1則由(1)(2)(3)式,有)()(21bfafabT)21(22112hafabTT)0(0T)1(0T)1(0kTTn記12kn若,2 , 1kjhaxj12kabhhxxjj212112)21(kabj

3、a12kabjakabja2)12(.55因此(1)(2)(3)式可化為如下遞推公式)()(2bfafab)0(0T120001)2)12(2)1(21)(kjkkabjafabkTkT,2 , 1k(4)-上式稱為遞推的梯形公式.66三種公式之間的關系244 1nnnTTS222441nnnSSC)1()(4141)1(11kTkTkTmmmmm)0(0T)1(0T)0(1T)2(0T)1(1T)0(2TRomberg算法求解步驟（3，-3）（-1,2）（0,1）（1，-1）（3，-4）.77二、外推加速公式由復合梯形公式的余項公式)(3122nnnTTTInnTTI31342可得nnjjn

4、TxfnabTI31) )(221(341021由(3)式1021)(6)(431njjnxfnabT.8812kn設)1(31)(3400kTkT102111)(6)(4 )(2)()(231njjnjjxfnabxfbfafnabI )(4)(2)()(6102111njjnjjxfxfbfafnabnS復合Simpson公式nnTTI31342令引入),1(1kT)1(1kT)1(31)(3400kTkTnS12kS-(5)-(6).99)(15122nnnSSSI因此由復合Simpson公式的余項可得nnSSI15115162)1(1kT12kS即)1(151)(151611kTkTn

5、S)(1kTnS2當然)1(2kT)1(151)(151611kTkT令nC自己證明-(6)nC-(7).1010)1(2kT12kCnC-(8)即)(2kTnC2當然同樣由復合Cotes公式的余項)(63122nnnCCCInnCCI63163642)1(631)(636422kTkT得)1(3kT令)1(631)(636422kTkT-(9).1111)1(1kT)1(31)(3400kTkT)1(151)(151611kTkT)1(2kT)1(3kT)1(631)(636422kTkT)()(2bfafab)0(0T120001)2)12(2)1(21)(kjkkabjafabkTkT,

6、2 , 1k外推加速公式以上整個過程稱為Romberg算法將上述結果綜合后.1212)1()(4141)1(11kTkTkTmmmmm其中外推加速公式可簡化為-(9)0(0T)1(0T)0(1T)2(0T)1(1T)0(2T)3(0T)2(1T)1(2T)0(3T,2 , 1mm可以推廣到并且,2 , 1kRomberg算法求解步驟.1313romberg.m例:20sinxdx計算積分前側矩形公式 z1 = 0.99212545660563 z11 = 0.99212545660563后側矩形公式 z2 = 1.358 z22 = 1.358梯形公式 z3 = 0.9999794382396

7、1Simpson公式 z4 = 1.6138階simpson公式 z5 =1.000自選步長梯形公式 z6 = 0.99999921563419自選步長Simpson公式 z7 =1.471Romberg公式 z8 = 0.99999999999802Mote-Carlo算法 z9 = 0.99821071589516 -0.437 -0.437 0.358 0.358 -0.039 0.613 -0.000 -0.581 0.471 -0.198 -0.484Jifenbijiao.m積分法積分值絕對誤差.1414如何構造Romberg算法.153 龍貝格(Romberg)積分方法 我們已經

8、知道,當被積函數f(x)在區間a,b上連續時,要使得復合梯形公式或復合拋物線公式比較精確地代替定積分 可將分點(即基點)加密,也就是將區間a,b細分,然后利用復合梯形公式或復合拋物線公式求積。( )baf x dx.16 若用Tm表示把a,b作m等分并按復合梯形公式求積的結果,將每一小段再對分,令新的小段的長h=h/2,則T2m與Tm之間有如下關系: 2112(21) mmmkTThf akh(528) 其中 .17 另外,若用Sm表示把a,b分成m(偶數)個小段按復合拋物線公式計算的結果,那么只要把Sm中的m改為2m,h改為h就有20122221201222212024222(4224)3(

9、244442)3(222)3MmmmmmmmmhSffffffhffffffhfffff 從Tm的定義可得到關系式 22441mmmTTS(529) .18 我們再舉一個計算上半單位圓面積的例子(它的準確面積為/2)。現用內接正多邊形的逼近方法來計算。 如圖5.6,圖(a) 、 圖(b)是用同樣的內接正多邊形計算上半單位圓的面積。圖(a)是用梯形方法計算其面積,圖(b)是用三角形方法計算其面積。 .19 圖 5.6 .20 設正多邊形邊數為n=2k,則由圖(b)利用三角形公式算得面積為 ( )201235111133524(1)350112sincossin2sin242112()()23!

10、25! 222()()23!25!222()()23!25!2kkkkkkkkkkkkkkkknnTnnnT同理 .21 如果組合一下,就會得到更精確的結果,即 (1)( )( )00154441225!(2 )kkkkTTT同理 (2)(1)(1)00151 4441225!(2)kkkkTTT.22 再以類似方法組合得 2(1)( )( )11164411()2(2 )kkkkTTTO 這樣繼續下去,其值越來越接近上半單位圓面積/2。這種方法可以用到計算定積分 ( )baf x dx.23 為了推廣公式(529)和上述計算上半單位圓面積的組合方法,我們引進龍貝格求積算法。 龍貝格求積算法本

11、來是利用所謂外推法構造出的一種計算積分的方法。為了避免從外推引入而帶來理論上的麻煩,我們將直接從構造一個T數表開始。 首先將a,b依次作20,21,22,等分,記,0,1,2,2iibahi.24 按復合梯形公式(520)算得的值相應地記為T(k)0(k=0,1,2,);把按式(529)算得的S2m依次記為T(k)1(k=0,1,2,崐),而這每一個S2m又理解為由T2m與Tm的線性組合得到的改進值,即(1)( )( )0014,0,1,2,41kkkTTTk 我們可按照類似的方法繼續進行改進,也即由S2m與Sm的線性組合得到改進值,依次記為T(k)2(k=0,1,2,),即 2(1)( )(

12、 )10224,0,1,2,41kkkTTTk.25 這樣就可構造出一個數表 (5-30).26 其中除第0列(即最左一列)的T(k)0是按復合梯形公式計算外,其余各列都按下述規則(對m) (1)( )( )111,2,4,0,1,41mkkkmmmmmTTTk(531) 遞推地計算出來。箭頭表示計算流程。其計算步驟為: (1)將區間a,b等分為20,用梯形公式計算T(0)0,即(0)0 ( )( )2baTf af b.27 (2)將區間a,b等分為21,用梯形公式算出T(1)0,即(1)0 ( )2 ()( )42baabTf aff b 再由T(0)0,T(1)0根據公式(531)算出T

13、(0)1,即(0)(0)(0)101441TTT若 T(0)1-T(0)0, (為預給的精度) 則停止計算;否則繼續往下計算;.28 (3)依次分別算出T(2)0,T(1)1,T(0)2,這一行地往下推算,每一行算完,就得驗證T(0)m(m=1,2,)是否滿足預給的精度,即若(0)(0)1mmTT則停止計算;否則繼續進行下一行。為了便于在計算機上實現,可運用下列公式編制程序:.291(0)02( )(1)001(1)( )( )11 ( )( )21(21)2221,2,4,0,1,41kkkkkiikkkiiiibaTf af bbabaTTf aiiTTTk.30 例 4 計算積分1204

14、1Idxx精確到10-4。 解(0)0(1)(0)00(1)(0)(0)00111(1) (0)(1)(42)322111(2)( )3.122243.13333341TffTTfTTT.31(2)(1)00(2)(1)(1)0012(1)(0)(0)11221113(3) ( )( )3.131176244443.1415694143.14211841TTffTTTTTT.32(3)(2)00(3)(2)(2)0012(2)(1)(1)11223(1)(0)(0)2233111357(4) ( )( )( )( )2888883.138988443141594414

15、43.14158641TTffffTTTTTTTTT.33(4)(3)00(4)(3)(3)0012(3)(2)(2)11223(2)(1)(1)223311135(5) ()()()21616161679111315()()()()()16161616163.14094243.1415934143.1415934143.14159341TTffffffffTTTTTTTTT.344(1)(0)(0)3344(0)34(0)1201102043.141593410.0000070.0000543.1415931443.14159261TTTTTdxxIdxarctgxx于是 由于 實際上 .3

16、5.36簡單的數值方法與基本概念簡單的數值方法與基本概念 1. 簡單簡單歐拉法（歐拉法（Euler） 2后退的歐拉法后退的歐拉法 3梯形法梯形法 4改進改進EulerEuler法法 2、初值問題的數值解法、初值問題的數值解法單步法單步法.371. 簡單的歐拉簡單的歐拉(Euler)方法方法考慮模型：00( , ) (1.1)() (1.2)yf x yaxby xy 在精度要求不高時通過歐拉方法的討論歐拉方法歐拉方法.382. 歐拉方法的導出歐拉方法的導出把區間a,b分為n個小區間步長為要計算出解函數 y(x) 在一系列節點()iiyy x)/()( iiixaihhhban,一般取即等距節點

17、處的近似值01 naxxxb1(-)iiihxxN N等分等分.3900( ,)(1.1)() (1.2)yfx yaxby xy 對微分方程(1.1)兩端從1nnxx到進行積分11( ,( )nnnnxxxxy dxfx y xdx11()()( , ( )nnxnnxy xy xf x y x dx.40右端積分用右端積分用左矩形數值左矩形數值求積公式求積公式2()( )()()()2bagg x dxba g aba( )( , ( )g xf x y x令1)1(, ()(, () nnnnxxnnf xy xnnyyf xy xh得.41x0 x11()()()(,)nnnnnny

18、xy xhy xyh f xy)1,., 0(),(1 niyxfhyyiiii亦稱為亦稱為歐拉折線法歐拉折線法 /* Eulers polygonal arc method*/ 每步計算只用到1ny或用向前差商近似導數1()()()nnny xy xyxh100021111(,)(,) (,)nnnnyyhf xyyyhf xyyyhf xy依上述公式逐次計算可得：ny故也稱Euler為單步法。公式右端只含有已知項所以又稱為顯格式的單步法。.42 例例1 用歐拉公式求解初值問題用歐拉公式求解初值問題)2 . 2(. 1)0(),10(2 yxyxyy 解解 取步長取步長h=0.1，歐拉公式的

19、具體形式為，歐拉公式的具體形式為)2(1nnnnnyxyhyy 其中其中xn=nh=0.1n (n=0,1,10), 已知已知y0 =1, 由此式可得由此式可得191818. 1)1 . 12 . 01 . 1 ( 1 . 01 . 1)2(1 . 11 . 01)2(1111200001 yxyhyyyxyhyy.43依次計算下去，依次計算下去，部分計算結果部分計算結果見下表見下表. xy21 與準確解與準確解 相比，可看出歐拉公式的計算結相比，可看出歐拉公式的計算結果精度很差果精度很差. xn 歐拉公式數值解歐拉公式數值解yn準確解準確解y(xn) 誤差誤差 0.2 0.4 0.6 0.8

20、 1.0 1.191818 1.358213 1.508966 1.649783 1.784770 1.183216 1.341641 1.483240 1.612452 1.732051 0.008602 0.016572 0.025726 0.037331 0.052719.44 歐拉公式具有明顯的幾何意義歐拉公式具有明顯的幾何意義, , 就是就是用折線近似用折線近似代替方程的解曲線代替方程的解曲線，因而常稱公式，因而常稱公式(2.1)為為歐拉折線法歐拉折線法. .( )yy xxynx1nxnp1np1np x 還可以通過幾何直觀來考察歐拉方法的精度還可以通過幾何直觀來考察歐拉方法的精度

21、. .假假設設yn=y(xn), ,即頂點即頂點Pn落在積分曲線落在積分曲線y=y(x)上，那么，上，那么，按歐拉方法做出的折線按歐拉方法做出的折線PnPn+1便是便是y=y(x)過點過點Pn的切線的切線. .從圖形上看從圖形上看, ,這這樣定出的頂點樣定出的頂點Pn+1顯著顯著地偏離了原來的積分曲地偏離了原來的積分曲線，可見歐拉方法是線，可見歐拉方法是相相當粗糙當粗糙的的. .45定義定義若某算法的局部截斷誤差為O(hp+1)，則稱該算法有p 階精度。Ri 的的主項主項/* leading term */ 歐拉法的局部截斷誤差：232()()hiyxO h歐拉法具有 1 階精度。),()()

22、()()()(32112iiiihiiiiiyxhfyhOxyxyhxyyxyR 定義定義在假設 yi = y(xi)，即第 i 步計算是精確的前提下，考慮的截斷誤差 Ri = y(xi+1) yi+1 稱為局部截斷誤差 /* local truncation error */。.465. 歐拉公式的改進歐拉公式的改進： 隱式歐拉法隱式歐拉法 /* implicit Euler method */向后差商近似導數向后差商近似導數hxyxyxy)()()(011 x0 x1)(,()(1101xyxfhyxy )1,., 0(),(111 niyxfhyyiiii由于未知數由于未知數 yi+1

23、同時出現在等式的兩邊，不能直接得到，故同時出現在等式的兩邊，不能直接得到，故稱為稱為隱式隱式 /* implicit */ 歐拉公式，而前者稱為歐拉公式，而前者稱為顯式顯式 /* explicit */ 歐拉公式。歐拉公式。一般先用顯式計算一個初值，再一般先用顯式計算一個初值，再迭代迭代求解。求解。 隱式隱式歐拉法的局部截斷誤差：歐拉法的局部截斷誤差：11)(iiiyxyR)()(322hOxyih 即隱式歐拉公式具有即隱式歐拉公式具有 1 階精度。階精度。.47 設用歐拉公式設用歐拉公式),() 0(1nnnnyxhfyy 給出迭代初值給出迭代初值 ，用它代入，用它代入(2.5)式的式的右端

24、，使之轉右端，使之轉化為顯式，直接計算得化為顯式，直接計算得) 0(1 ny),() 0(11) 1 (1 nnnnyxhfyy然后再用然后再用 代入代入(2.5)式，又有式，又有) 1 (1 ny).,() 1 (11) 2(1 nnnnyxhfyy如此反復進行，得如此反復進行，得) 6 . 2 ()., 1 , 0(),()(11) 1(1 kyxhfyyknnnkn.486.梯形公式梯形公式 /* trapezoid formula */ 顯、隱式兩種算法的顯、隱式兩種算法的平均平均)1,., 0(),(),(2111 niyxfyxfhyyiiiiii注：注：的確有局部截斷誤差的確有局

25、部截斷誤差 ， 即梯形公式具有即梯形公式具有2 階精度，比歐拉方法有了進步。階精度，比歐拉方法有了進步。但注意到該公式是但注意到該公式是隱式隱式公式，計算時不得不用到公式，計算時不得不用到迭代法，其迭代收斂性與歐拉公式相似。迭代法，其迭代收斂性與歐拉公式相似。)()(311hOyxyRiii 中點歐拉公式中點歐拉公式 /* midpoint formula */中心差商近似導數中心差商近似導數hxyxyxy2)()()(021 x0 x2x1)(,(2)()(1102xyxfhxyxy 1,., 1),(211 niyxfhyyiiii假設假設 ，則可以導出，則可以導出即中點公式具有即中點公式

26、具有 2 階精度。階精度。)(),(11iiiixyyxyy )()(311hOyxyRiii 需要需要2個初值個初值 y0和和 y1來啟動遞推來啟動遞推過程，這樣的算法稱為過程，這樣的算法稱為雙步法雙步法 /* double-step method */，而前面的三種算法都是，而前面的三種算法都是單步法單步法 /* single-step method */。.49方方 法法 顯式歐拉顯式歐拉隱式歐拉隱式歐拉梯形公式梯形公式中點公式中點公式簡單簡單精度低精度低穩定性最好穩定性最好精度低精度低, 計算量大計算量大精度提高精度提高計算量大計算量大精度提高精度提高, 顯式顯式多一個初值多一個初值, 可能影響精度可能影響精度.5