數(shù)值分析第九章

1、數(shù)值分析,第九章 常微分方程的數(shù)值解,一、 Euler方法,三、 單步法的收斂性和穩(wěn)定性,二、 Runge-Kutta方法,四、 線性多步法,很多科學(xué)技術(shù)和工程問題常用常微分方程的形式建立數(shù)學(xué)模型.,但是對(duì)于絕大多數(shù)的微分方程問題，很難或者根本不可能得到它的解析解.,本章重點(diǎn)考察一階方程的初值問題,的數(shù)值解法，就是尋求解y(x)在一系列離散點(diǎn),處的近似值 的方法.,相鄰兩個(gè)節(jié)點(diǎn)間的距離

2、 稱為步長.,一、 Euler方法,1 歐拉公式,由初值條件 表示積分曲線從,出發(fā)，并在 處的切線斜率為,因此可以設(shè)想積分曲線在 x=x0 附近可以用切,線近似的代替曲線.,切線方程為,當(dāng)x=x1時(shí)，代入有,這樣得到y(tǒng)(x1)的近似值y1的方法.,重復(fù)上述方法，當(dāng) x=x2 時(shí),依次可以計(jì)算出x3, x4, …處的近似值y3, y4, …,由此得到Euler公式：,由于

3、用折線近似代替方程的解析解，所以Euler方法也稱為Euler折線法.,例 用Euler法計(jì)算初值問題的解在x=0.3時(shí)的近似值，取步長h=0.1 .,解：,Euler公式的截?cái)嗾`差,局部截?cái)嗾`差：一步Euler公式產(chǎn)生的誤差;,總體截?cái)嗾`差：Euler公式的累積總誤差;,歐拉法的局部截?cái)嗾`差：,所以歐拉法具有 1 階精度.,Lipschitiz條件：若存在正數(shù)L,使得對(duì)一切,x, y1, y2有,則稱f(x, y)滿足Lipschit

4、iz條件.,歐拉法的總體截?cái)嗾`差：,那么,設(shè),為局部截?cái)嗾`差，所以,特別當(dāng)n=m-1時(shí)，有,總體誤差與h是同階的.,上式還說明，當(dāng) 時(shí)，有 即,也就是說，ym收斂到方程的準(zhǔn)確解,后退Euler公式,（隱式歐拉法）,（隱式歐拉公式 ）,利用向后差商近似導(dǎo)數(shù),由于未知數(shù)yn+1同時(shí)出現(xiàn)在等式的兩邊，不能直接得到，故稱為隱式歐拉公式，而前者稱為顯式歐拉公式.,一般先用顯式計(jì)算一個(gè)初值，再迭代求解.,隱式

5、歐拉法的局部截?cái)嗾`差：,即隱式歐拉公式具有1階精度.,2 梯形公式和改進(jìn)Euler方法,梯形公式,設(shè)y=y(x)是 的解, 故,由此得到,,用yn來近似y(xn),用yn+1來近似y(xn+1), 得,梯形公式,梯形公式是隱式的，可以用迭代法求解.,具有2階精度.,梯形公式的局部截?cái)嗾`差,中點(diǎn)歐拉公式,假設(shè) , 則可以導(dǎo)出即中點(diǎn)公式

6、具有 2 階精度.,,需要2個(gè)初值 y0和 y1來啟動(dòng)遞推過程，這樣的算法稱為雙步法 ，而前面的三種算法都是單步法.,簡單,精度低,穩(wěn)定性最好,精度低, 計(jì)算量大,精度提高,計(jì)算量大,精度提高, 顯式,多一個(gè)初值, 可能影響精度,有沒有一種方法，既有這些方法的優(yōu)點(diǎn)，而沒有它們的缺點(diǎn)？,改進(jìn)歐拉法,(1)先用顯式歐拉公式作預(yù)測，算出,,(2)再將 代入梯形公式的右邊作校正,得到,注：此法亦稱為預(yù)測-校正法.可以證明該算法具有

7、2階精度，同時(shí)可以看到它是個(gè)單步遞推格式，比隱式公式的迭代求解過程簡單.,例 用梯形公式求解初值問題(步長h=0.2),解：,梯形公式為,于是,整理得,由 y(1) = y0 = 2 依次可得 y1, y2, y3, y4, y5.,例 用改進(jìn)歐拉法求解初值問題,要求步長h=0.2，并計(jì)算y(1.2)和y(1.4),解：,改進(jìn)歐拉法公式為,即,由 y(1) = y0 =1 計(jì)算得,二、Runge-Kutta方法,建立高精度的單步遞推

8、格式.,單步遞推法的基本思想是從(xn , yn )點(diǎn)出發(fā)，以某一斜率沿直線達(dá)到(xn+1 , yn+1 )點(diǎn). 歐拉法及其各種變形所能達(dá)到的最高精度為2階.,考察改進(jìn)的歐拉法，可以將其改寫為：,,斜率一定取K1 K2 的平均值嗎？,步長一定是一個(gè)h 嗎？,首先希望能確定系數(shù) ?1、?2、p，使得到的算法格式有2階精度，即在 的前提假設(shè)下，使得,將改進(jìn)歐拉法推廣為：,(1) 將K2在(xn ,

9、yn )點(diǎn)作 Taylor 展開,1 二階Runge-Kutta方法,(2) 將 K2 代入第1式，得到,(3) 將yn+1與y( xn+1 )在xn點(diǎn)的泰勒展開作比較,要求 , 則必須有：,所以存在無窮多個(gè)解！,所有滿足上式的統(tǒng)稱為2階Runge-Kutta格式.,若 則,改進(jìn)的歐拉方法,若

10、 則,中點(diǎn)公式,2 四階Runge-Kutta方法,其中?i (i = 1, …, m)，?i (i =2, …, m) 和?ij (i = 2, …, m; j =1, …, i?1 ) 均為待定系數(shù)，確定這些系數(shù)的步驟與前面相似.,由于方程的個(gè)數(shù)少于未知量的個(gè)數(shù)，所以方程有無窮多個(gè)解，可以根據(jù)情況得到幾種常用的解，即得到相應(yīng)的四階公式.,最常用為四階經(jīng)典龍格-庫塔法,也稱為標(biāo)準(zhǔn)四階龍

11、格-庫塔公式,Gill公式,(2) 龍格-庫塔法的導(dǎo)出基于泰勒展開，故精度主要受解函數(shù)的光滑性影響. 對(duì)于光滑性不太好的解, 最好采用低階算法而將步長h 取小.,注：,(1) 龍格-庫塔法的主要運(yùn)算在于計(jì)算Ki 的值,即計(jì)算f的值. Butcher于1965年給出了計(jì)算量與可達(dá)到的最高精度階數(shù)的關(guān)系：,例 用標(biāo)準(zhǔn)四階Runge-Kutta法求初值問題,在x=0.1處的近似值，取步長為h=0.1 .,解：,所以,那么,例 用標(biāo)準(zhǔn)四階Run

12、ge-Kutta法求初值問題,在x=0.4處的近似值，取步長為h=0.2 .,解：,所以,而,所以,1 單步法的收斂性,三、 單步法的收斂性和穩(wěn)定性,單步法是在計(jì)算yn+1時(shí)只用到前一步的信息yn .,顯式單步法的共同特征是它們都是將yn加上某種形式的增量，得出yn+1,計(jì)算公式如下：,增量函數(shù),Euler方法的增量函數(shù),改進(jìn)Euler方法的增量函數(shù),則稱,為顯式單步法在xn+1處的局部截?cái)嗾`差 .,例：考察歐拉顯式格式的收斂性:,解：

13、該問題的精確解為,歐拉公式為,,對(duì)任意固定的 x=xi =ih，有,,?,?,Tn+1按h展開的第一項(xiàng)，又稱為主項(xiàng).,若局部截?cái)嗾`差的展開式寫成,則稱 為局部截?cái)嗾`差的主項(xiàng),單步法的收斂定理,設(shè)單步法 具有p階精度,其增量函數(shù) 關(guān)于y滿足Lipschitz條件,即存在常數(shù)L，使對(duì)任何的 及

14、任意的x有,又設(shè)初值y0是準(zhǔn)確的，即 則總體,截?cái)嗾`差是p階的，也就是,特別的當(dāng) 時(shí), 不論n為何值, 總有,即方法收斂.,在f(x,y)對(duì)y滿足Lipschitz條件下, Euler法，改進(jìn)Euler法和Runge-Kutta法的增量函數(shù),都對(duì)y滿足Lipschitz條件，所以上述結(jié)論對(duì)這些方法都成立.,例 設(shè),是求解微分方程的單步法，試求其局部截?cái)嗾`差的主項(xiàng)，并說出它具有幾階精度.,

15、解:,考慮 在xn處的Taylor展式,所以,該方法的局部截?cái)嗾`差的主項(xiàng)是,具有一階精度.,解:,考慮 在xn處的Taylor展式,所以,該方法的局部截?cái)嗾`差的主項(xiàng)是,具有二階精度.,2 單步法的穩(wěn)定性,收斂性是在假定每一步計(jì)算都準(zhǔn)確的前提下, 討論步長 時(shí)，方法的總體截?cái)嗾`差是否趨于零的問題.,穩(wěn)定性是討論舍入誤差的積

16、累能否對(duì)計(jì)算結(jié)果有嚴(yán)重的影響.,例：考察初值問題 在區(qū)間,1.0000?2.0000 4.0000?8.0000 1.6000?101 ?3.2000?101,1.00002.5000?10?1 6.2500?10?21.5625?10?23.9063?10?39.7656?10?4,1.00002.50006.25001.5626?1013.

17、9063?1019.7656?101,1.00004.9787?10?22.4788?10?31.2341?10?46.1442?10?63.0590?10?7,[0, 0.5]上的解，分別用歐拉顯、隱式格式和改進(jìn)的歐拉格式計(jì)算數(shù)值解.,一般分析時(shí)為簡單起見，只考慮試驗(yàn)方程,λ 為復(fù)數(shù)且Re(λ)<0,設(shè) 在節(jié)點(diǎn)值yn處有擾動(dòng) 令,那么,于是,反復(fù)應(yīng)用可得,為使,

18、則,可得如下定義,下面討論已知的幾種方法的絕對(duì)穩(wěn)定區(qū)間和絕對(duì)穩(wěn)定區(qū)域.,顯式歐拉法:,,在復(fù)平面上的絕對(duì)穩(wěn)定區(qū)域是,即以-1為中心，1為半徑的圓域,所以,相應(yīng)的絕對(duì)穩(wěn)定區(qū)間是,隱式歐拉法(后退歐拉法)：,在復(fù)平面上的絕對(duì)穩(wěn)定區(qū)域是,是以1為中心，1為半徑的圓的外域,所以,相應(yīng)的絕對(duì)穩(wěn)定區(qū)間是,即,如果只考慮λ <0的實(shí)數(shù)，則相應(yīng)的絕對(duì)穩(wěn)定區(qū)間對(duì)于任意的h都成立，所以是無條件穩(wěn)定的,梯形公式：,在復(fù)平面上的絕對(duì)穩(wěn)定區(qū)域是,也是無條件

19、穩(wěn)定的,相應(yīng)的絕對(duì)穩(wěn)定區(qū)間是,龍格-庫塔法：,而顯式1~ 4 階方法的絕對(duì)穩(wěn)定區(qū)域?yàn)?例 設(shè),是求解微分方程的單步法，分析它的穩(wěn)定性.,解:,所以絕對(duì)穩(wěn)定區(qū)域是,即 為復(fù)平面的左半平面.,在實(shí)數(shù)域上是無條件穩(wěn)定的.,解：,將 代入得,即,例 討論求解初值問題 的求解,公式：,的穩(wěn)定性. (λ>0為實(shí)

20、數(shù)),所以絕對(duì)穩(wěn)定區(qū)域是,所以,因此是條件穩(wěn)定的.,四、 線性多步法,在逐步推進(jìn)的求解過程中，計(jì)算yn+1之前已經(jīng)求出了一系列的近似值y0, y1, …, yn ,如果充分利用前面信息來預(yù)測yn+1，則可期望會(huì)獲得較高的精度，這就是線性多步法的基本思想.,1 線性多步法的一般公式,最常用的線性多步法公式為,其中 為常數(shù)，yn-k為y(xn-k)的近似值,fn-k = f(xn-k ,yn-k),特別的當(dāng)

21、 時(shí), 上式為顯式, 否則是隱式.,若 則稱該方法具有p階精度.,若 則稱局部截?cái)嗾`差的主項(xiàng)為 為誤差常數(shù).,例 設(shè) yn+1 = yn-1+2hf(xn, yn) 為求解常微分初值問題的線性二步法，試求該二步公式的局部截?cái)嗾`差主項(xiàng)，和精度.,解:,

22、由局部截?cái)嗾`差的定義可知,考慮 在xn處的Taylor展式,代入可得,所以局部截?cái)嗾`差的主項(xiàng)為,具有二階精度.,解:,局部截?cái)嗾`差為,考慮 在xn處的Taylor展式,于是,為使方法具有二階精度則,解得,因此該方法為,局部截?cái)嗾`差的主項(xiàng)為,解:,局部截?cái)嗾`差為,考慮

23、 在xn處的Taylor展式,所以,解得,局部截?cái)嗾`差的主項(xiàng)為,2 Adams外推公式,考慮用r+1個(gè)點(diǎn)(xn-k, f(xn-k, yn-k))構(gòu)造一個(gè)r次多項(xiàng)式來近似,的被積函數(shù)f(x, y(x)), 這里用yn-k作為y(xn-k)的近似值，令fn-k=f(xn-k, yn-k), 用Newton向后插值公式來構(gòu)造r次多項(xiàng)式, 即,這里,而,將yn+1作為y(xn+1)的近似值，因此得到,Adams外推公式,顯式格式,當(dāng)r=0

24、時(shí), 為Euler公式, 最常用的是r=3情況,Adams外推系數(shù),3 Adams內(nèi)插公式,用xn+1, xn ,…, xn-r+1為插值節(jié)點(diǎn), 構(gòu)造f(x,y(x))的r次多項(xiàng)式，得到內(nèi)插公式,隱式格式,最常用的是r=3情況,Adams外推系數(shù),4 預(yù)報(bào)-校正公式,不論單步法或多步法，隱式公式比顯式公式穩(wěn)定性好，但在實(shí)際使用隱式公式時(shí)，都會(huì)遇到兩個(gè)問題：一個(gè)是隱式公式如何能方便地進(jìn)行計(jì)算；另一個(gè)是實(shí)際計(jì)算步長取多大.,如隱式梯形公式，

25、每往前推進(jìn)一步，不必進(jìn)行多次迭代，而是采用一階顯式Euler公式預(yù)測，二階隱式梯形公式校正一次，構(gòu)成顯式改進(jìn)Euler公式，能達(dá)到與梯形公式同階的精度，即二階精度.,預(yù)測-校正技術(shù)即保證了計(jì)算精度，又使隱式計(jì)算顯式化，克服了隱式公式要反復(fù)迭代的困難.,類似改進(jìn)Euler方法，將Adams外推公式和內(nèi)插公式結(jié)合起來，得到實(shí)用的預(yù)報(bào)-校正公式,預(yù)報(bào),校正,該方法的精度是四階的,解：,為運(yùn)用Adams方法進(jìn)行計(jì)算，先用四階,Runge-Kut