Romberg積分及初值問題

1、.11 5.3 Romberg算法算法 綜合前幾節(jié)的內(nèi)容,我們知道梯形公式,Simpson公式,Cotes公式的代數(shù)精度分別為1次,3次和5次復(fù)合梯形、復(fù)合Simpson、復(fù)合Cotes公式的收斂階分別為2階、4階和6階無論從代數(shù)精度還是收斂速度,復(fù)合梯形公式都是較差的有沒有辦法改善梯形公式呢?.22一、復(fù)合梯形公式的遞推化等份分割為的積分區(qū)間將定積分nbadxxfIba,)(njjhaxk, 1 ,0,nabh各節(jié)點為 )()(2)(211njjnbfxfafnabT復(fù)合梯形(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三種公式之間的關(guān)系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二、外推加速公式由復(fù)合梯形公式的余項公式)(3122nnnTTTInnTTI31342可得nnjjn

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

5、S)(1kTnS2當然)1(2kT)1(151)(151611kTkT令nC自己證明-(6)nC-(7).1010)1(2kT12kCnC-(8)即)(2kTnC2當然同樣由復(fù)合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算法將上述結(jié)果綜合后.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計算積分前側(cè)矩形公式 z1 = 0.99212545660563 z11 = 0.99212545660563后側(cè)矩形公式 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如何構(gòu)造Romberg算法.153 龍貝格(Romberg)積分方法 我們已經(jīng)

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

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

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

11、來是利用所謂外推法構(gòu)造出的一種計算積分的方法。為了避免從外推引入而帶來理論上的麻煩,我們將直接從構(gòu)造一個T數(shù)表開始。 首先將a,b依次作20,21,22,等分,記,0,1,2,2iibahi.24 按復(fù)合梯形公式(520)算得的值相應(yīng)地記為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 我們可按照類似的方法繼續(xù)進行改進,也即由S2m與Sm的線性組合得到改進值,依次記為T(k)2(k=0,1,2,),即 2(1)( )(

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

13、(0)1,即(0)(0)(0)101441TTT若 T(0)1-T(0)0, (為預(yù)給的精度) 則停止計算;否則繼續(xù)往下計算;.28 (3)依次分別算出T(2)0,T(1)1,T(0)2,這一行地往下推算,每一行算完,就得驗證T(0)m(m=1,2,)是否滿足預(yù)給的精度,即若(0)(0)1mmTT則停止計算;否則繼續(xù)進行下一行。為了便于在計算機上實現(xiàn),可運用下列公式編制程序:.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簡單的數(shù)值方法與基本概念簡單的數(shù)值方法與基本概念 1. 簡單簡單歐拉法（歐拉法（Euler） 2后退的歐拉法后退的歐拉法 3梯形法梯形法 4改進改進EulerEuler法法 2、初值問題的數(shù)值解法、初值問題的數(shù)值解法單步法單步法.371. 簡單的歐拉簡單的歐拉(Euler)方法方法考慮模型：00( , ) (1.1)() (1.2)yf x yaxby xy 在精度要求不高時通過歐拉方法的討論歐拉方法歐拉方法.382. 歐拉方法的導(dǎo)出歐拉方法的導(dǎo)出把區(qū)間a,b分為n個小區(qū)間步長為要計算出解函數(shù) y(x) 在一系列節(jié)點()iiyy x)/()( iiixaihhhban,一般取即等距節(jié)點

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右端積分用右端積分用左矩形數(shù)值左矩形數(shù)值求積公式求積公式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或用向前差商近似導(dǎo)數(shù)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依次計算下去，依次計算下去，部分計算結(jié)果部分計算結(jié)果見下表見下表. xy21 與準確解與準確解 相比，可看出歐拉公式的計算結(jié)相比，可看出歐拉公式的計算結(jié)果精度很差果精度很差. xn 歐拉公式數(shù)值解歐拉公式數(shù)值解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、. .假假設(shè)設(shè)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 定義定義在假設(shè) yi = y(xi)，即第 i 步計算是精確的前提下，考慮的截斷誤差 Ri = y(xi+1) yi+1 稱為局部截斷誤差 /* local truncation error */。.465. 歐拉公式的改進歐拉公式的改進： 隱式歐拉法隱式歐拉法 /* implicit Euler method */向后差商近似導(dǎo)數(shù)向后差商近似導(dǎo)數(shù)hxyxyxy)()()(011 x0 x1)(,()(1101xyxfhyxy )1,., 0(),(111 niyxfhyyiiii由于未知數(shù)由于未知數(shù) yi+1

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

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

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

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