《數(shù)值分析》第九章常微分方程的數(shù)值解同步課程

數(shù)值分析第九章第一頁(yè)，共69頁(yè)。第九章常微分方程的數(shù)值解一、Euler方法三、單步法的收斂性和穩(wěn)定性二、Runge-Kutta方法四、線(xiàn)性多步法第二頁(yè)，共69頁(yè)。很多科學(xué)技術(shù)和工程問(wèn)題常用常微分方程的形式建立數(shù)學(xué)模型.但是對(duì)于絕大多數(shù)的微分方程問(wèn)題，很難或者根本不可能得到它的解析解.本章重點(diǎn)考察一階方程的初值問(wèn)題的數(shù)值解法，就是尋求解y(x)在一系列離散點(diǎn)處的近似值的方法.相鄰兩個(gè)節(jié)點(diǎn)間的距離稱(chēng)為步長(zhǎng).第三頁(yè)，共69頁(yè)。一、Euler方法1歐拉公式由初值條件表示積分曲線(xiàn)從出發(fā)，并在處的切線(xiàn)斜率為因此可以設(shè)想積分曲線(xiàn)在x=x0附近可以用切線(xiàn)近似的代替曲線(xiàn).切線(xiàn)方程為當(dāng)x=x1時(shí)，代入有這樣得到y(tǒng)(x1)的近似值y1的方法.第四頁(yè)，共69頁(yè)。重復(fù)上述方法，當(dāng)x=x2時(shí)依次可以計(jì)算出x3,x4,…處的近似值y3,y4,…由此得到Euler公式：由于用折線(xiàn)近似代替方程的解析解，所以Euler方法也稱(chēng)為Euler折線(xiàn)法.例用Euler法計(jì)算初值問(wèn)題的解在x=0.3時(shí)的近似值，取步長(zhǎng)h=0.1.第五頁(yè)，共69頁(yè)。解：Euler公式的截?cái)嗾`差局部截?cái)嗾`差：一步Euler公式產(chǎn)生的誤差;總體截?cái)嗾`差：Euler公式的累積總誤差;第六頁(yè)，共69頁(yè)。

在假設(shè)yn=y(xn)，即第i步計(jì)算是精確的前提下，考慮的截?cái)嗾`差Rn=y(xn+1)

yn+1稱(chēng)為局部截?cái)嗾`差.定義歐拉法的局部截?cái)嗾`差：所以歐拉法具有1階精度.

若某算法的局部截?cái)嗾`差為O(hp+1)，則稱(chēng)該算法有p

階精度.定義第七頁(yè)，共69頁(yè)。Lipschitiz條件：若存在正數(shù)L,使得對(duì)一切x,y1,y2有則稱(chēng)f(x,y)滿(mǎn)足Lipschitiz條件.歐拉法的總體截?cái)嗾`差：那么設(shè)為局部截?cái)嗾`差，所以第八頁(yè)，共69頁(yè)。第九頁(yè)，共69頁(yè)。特別當(dāng)n=m-1時(shí)，有總體誤差與h是同階的.上式還說(shuō)明，當(dāng)時(shí)，有即也就是說(shuō)，ym收斂到方程的準(zhǔn)確解第十頁(yè)，共69頁(yè)。后退Euler公式（隱式歐拉法）（隱式歐拉公式）利用向后差商近似導(dǎo)數(shù)第十一頁(yè)，共69頁(yè)。由于未知數(shù)yn+1同時(shí)出現(xiàn)在等式的兩邊，不能直接得到，故稱(chēng)為隱式歐拉公式，而前者稱(chēng)為顯式歐拉公式.一般先用顯式計(jì)算一個(gè)初值，再迭代求解.隱式歐拉法的局部截?cái)嗾`差：即隱式歐拉公式具有1階精度.第十二頁(yè)，共69頁(yè)。2梯形公式和改進(jìn)Euler方法梯形公式設(shè)y=y(x)是的解,故由此得到用梯形公式近似第十三頁(yè)，共69頁(yè)。用yn來(lái)近似y(xn),用yn+1來(lái)近似y(xn+1),得梯形公式梯形公式是隱式的，可以用迭代法求解.第十四頁(yè)，共69頁(yè)。具有2階精度.梯形公式的局部截?cái)嗾`差第十五頁(yè)，共69頁(yè)。中點(diǎn)歐拉公式中心差商近似導(dǎo)數(shù)假設(shè),則可以導(dǎo)出即中點(diǎn)公式具有2階精度.需要2個(gè)初值y0和y1來(lái)啟動(dòng)遞推過(guò)程，這樣的算法稱(chēng)為雙步法，而前面的三種算法都是單步法.第十六頁(yè)，共69頁(yè)。方法

顯式歐拉隱式歐拉梯形公式中點(diǎn)公式簡(jiǎn)單精度低穩(wěn)定性最好精度低,計(jì)算量大精度提高計(jì)算量大精度提高,顯式多一個(gè)初值,可能影響精度有沒(méi)有一種方法，既有這些方法的優(yōu)點(diǎn)，而沒(méi)有它們的缺點(diǎn)？第十七頁(yè)，共69頁(yè)。改進(jìn)歐拉法(1)先用顯式歐拉公式作預(yù)測(cè)，算出(2)再將代入梯形公式的右邊作校正,得到注：此法亦稱(chēng)為預(yù)測(cè)-校正法.可以證明該算法具有2階精度，同時(shí)可以看到它是個(gè)單步遞推格式，比隱式公式的迭代求解過(guò)程簡(jiǎn)單.第十八頁(yè)，共69頁(yè)。例用梯形公式求解初值問(wèn)題(步長(zhǎng)h=0.2)解：梯形公式為于是整理得由y(1)=y0=2依次可得y1,y2,y3,y4,y5.第十九頁(yè)，共69頁(yè)。例用改進(jìn)歐拉法求解初值問(wèn)題要求步長(zhǎng)h=0.2，并計(jì)算y(1.2)和y(1.4)解：改進(jìn)歐拉法公式為即第二十頁(yè)，共69頁(yè)。由y(1)=y0=1計(jì)算得第二十一頁(yè)，共69頁(yè)。二、Runge-Kutta方法建立高精度的單步遞推格式.單步遞推法的基本思想是從(xn,yn)點(diǎn)出發(fā)，以某一斜率沿直線(xiàn)達(dá)到(xn+1

,yn+1

)點(diǎn).歐拉法及其各種變形所能達(dá)到的最高精度為2階.考察改進(jìn)的歐拉法，可以將其改寫(xiě)為：斜率一定取K1K2的平均值嗎？步長(zhǎng)一定是一個(gè)h

嗎？第二十二頁(yè)，共69頁(yè)。首先希望能確定系數(shù)

1、

2、p，使得到的算法格式有2階精度，即在的前提假設(shè)下，使得

將改進(jìn)歐拉法推廣為：(1)

將K2在(xn

,yn

)點(diǎn)作Taylor展開(kāi)1二階Runge-Kutta方法第二十三頁(yè)，共69頁(yè)。(2)

將K2代入第1式，得到第二十四頁(yè)，共69頁(yè)。(3)將yn+1與y(xn+1)在xn點(diǎn)的泰勒展開(kāi)作比較要求,則必須有：這里有個(gè)未知數(shù)，個(gè)方程。32所以存在無(wú)窮多個(gè)解！第二十五頁(yè)，共69頁(yè)。所有滿(mǎn)足上式的統(tǒng)稱(chēng)為2階Runge-Kutta格式.若則改進(jìn)的歐拉方法若則中點(diǎn)公式第二十六頁(yè)，共69頁(yè)。2四階Runge-Kutta方法其中

i(i=1,…,m)，

i(i=2,…,m)

和

ij(i=2,…,m;j=1,…,i1

)

均為待定系數(shù)，確定這些系數(shù)的步驟與前面相似.第二十七頁(yè)，共69頁(yè)。由于方程的個(gè)數(shù)少于未知量的個(gè)數(shù)，所以方程有無(wú)窮多個(gè)解，可以根據(jù)情況得到幾種常用的解，即得到相應(yīng)的四階公式.最常用為四階經(jīng)典龍格-庫(kù)塔法也稱(chēng)為標(biāo)準(zhǔn)四階龍格-庫(kù)塔公式第二十八頁(yè)，共69頁(yè)。Gill公式第二十九頁(yè)，共69頁(yè)。753可達(dá)到的最高精度642每步須算Ki的個(gè)數(shù)(2)龍格-庫(kù)塔法的導(dǎo)出基于泰勒展開(kāi)，故精度主要受解函數(shù)的光滑性影響.對(duì)于光滑性不太好的解,最好采用低階算法而將步長(zhǎng)h

取小.注：(1)龍格-庫(kù)塔法的主要運(yùn)算在于計(jì)算Ki

的值,即計(jì)算f的值.Butcher于1965年給出了計(jì)算量與可達(dá)到的最高精度階數(shù)的關(guān)系：第三十頁(yè)，共69頁(yè)。例用標(biāo)準(zhǔn)四階Runge-Kutta法求初值問(wèn)題在x=0.1處的近似值，取步長(zhǎng)為h=0.1.解：所以第三十一頁(yè)，共69頁(yè)。那么例用標(biāo)準(zhǔn)四階Runge-Kutta法求初值問(wèn)題在x=0.4處的近似值，取步長(zhǎng)為h=0.2.第三十二頁(yè)，共69頁(yè)。解：所以而所以第三十三頁(yè)，共69頁(yè)。

若某算法對(duì)于任意固定的x=xi=x0+ih,當(dāng)h0

(同時(shí)i

)時(shí)有yi

y(xi

)，則稱(chēng)該算法是收斂的.定義1單步法的收斂性三、單步法的收斂性和穩(wěn)定性單步法是在計(jì)算yn+1時(shí)只用到前一步的信息yn

.顯式單步法的共同特征是它們都是將yn加上某種形式的增量，得出yn+1,計(jì)算公式如下：增量函數(shù)第三十四頁(yè)，共69頁(yè)。Euler方法的增量函數(shù)改進(jìn)Euler方法的增量函數(shù)

設(shè)y(x)是微分方程初值問(wèn)題的準(zhǔn)確解，定義則稱(chēng)為顯式單步法在xn+1處的局部截?cái)嗾`差.第三十五頁(yè)，共69頁(yè)。例：考察歐拉顯式格式的收斂性:解：該問(wèn)題的精確解為

歐拉公式為對(duì)任意固定的x=xi=ih，有

第三十六頁(yè)，共69頁(yè)。

設(shè)y(x)是微分方程初值問(wèn)題的準(zhǔn)確解，定義若存在最大整數(shù)p，使顯式單步法的局部截?cái)嗾`差滿(mǎn)足則稱(chēng)該方法具有p階精度，或稱(chēng)為p方法.Tn+1按h展開(kāi)的第一項(xiàng)，又稱(chēng)為主項(xiàng).若局部截?cái)嗾`差的展開(kāi)式寫(xiě)成則稱(chēng)為局部截?cái)嗾`差的主項(xiàng)第三十七頁(yè)，共69頁(yè)。單步法的收斂定理設(shè)單步法具有p階精度其增量函數(shù)關(guān)于y滿(mǎn)足Lipschitz條件即存在常數(shù)L，使對(duì)任何的及任意的x有又設(shè)初值y0是準(zhǔn)確的，即則總體截?cái)嗾`差是p階的，也就是特別的當(dāng)時(shí),不論n為何值,

總有即方法收斂.第三十八頁(yè)，共69頁(yè)。在f(x,y)對(duì)y滿(mǎn)足Lipschitz條件下,Euler法，改進(jìn)Euler法和Runge-Kutta法的增量函數(shù)

都對(duì)y滿(mǎn)足Lipschitz條件，所以上述結(jié)論對(duì)這些方法都成立.例設(shè)是求解微分方程的單步法，試求其局部截?cái)嗾`差的主項(xiàng)，并說(shuō)出它具有幾階精度.解:第三十九頁(yè)，共69頁(yè)?？紤]在xn處的Taylor展式所以該方法的局部截?cái)嗾`差的主項(xiàng)是具有一階精度.第四十頁(yè)，共69頁(yè)。例設(shè)試求出它具有幾階精度.解:考慮在xn處的Taylor展式第四十一頁(yè)，共69頁(yè)。所以該方法的局部截?cái)嗾`差的主項(xiàng)是具有二階精度.2單步法的穩(wěn)定性收斂性是在假定每一步計(jì)算都準(zhǔn)確的前提下,討論步長(zhǎng)時(shí)，方法的總體截?cái)嗾`差是否趨于零的問(wèn)題.穩(wěn)定性是討論舍入誤差的積累能否對(duì)計(jì)算結(jié)果有嚴(yán)重的影響.第四十二頁(yè)，共69頁(yè)。例：考察初值問(wèn)題在區(qū)間0.00.10.20.30.40.5精確解改進(jìn)歐拉法

歐拉隱式歐拉顯式

節(jié)點(diǎn)xi

1.0000

2.00004.0000

8.00001.6000101

3.2000101

1.00002.5000101

6.25001021.56251023.90631039.76561041.00002.50006.25001.56261013.90631019.76561011.00004.97871022.47881031.23411046.14421063.0590107[0,0.5]上的解，分別用歐拉顯、隱式格式和改進(jìn)的歐拉格式計(jì)算數(shù)值解.第四十三頁(yè)，共69頁(yè)。

若一種數(shù)值解法僅在節(jié)點(diǎn)值yn上有大小為δ的擾動(dòng)，于以后各節(jié)點(diǎn)值ym(m>n)上，僅由δ所引起的擾動(dòng)都不超過(guò)δ時(shí)，稱(chēng)該方法是穩(wěn)定的.定義一般分析時(shí)為簡(jiǎn)單起見(jiàn)，只考慮試驗(yàn)方程λ

為復(fù)數(shù)且Re(λ)<0設(shè)在節(jié)點(diǎn)值yn處有擾動(dòng)令那么于是第四十四頁(yè)，共69頁(yè)。反復(fù)應(yīng)用可得為使則可得如下定義

若中，則稱(chēng)單步法是絕對(duì)穩(wěn)定的.在復(fù)平面上，λh滿(mǎn)足

的區(qū)域稱(chēng)為方法的絕對(duì)穩(wěn)定區(qū)域,它與實(shí)軸的交稱(chēng)為絕對(duì)穩(wěn)定區(qū)間定義下面討論已知的幾種方法的絕對(duì)穩(wěn)定區(qū)間和絕對(duì)穩(wěn)定區(qū)域.第四十五頁(yè)，共69頁(yè)。顯式歐拉法:0-1-2ReImg在復(fù)平面上的絕對(duì)穩(wěn)定區(qū)域是即以-1為中心，1為半徑的圓域所以相應(yīng)的絕對(duì)穩(wěn)定區(qū)間是第四十六頁(yè)，共69頁(yè)。隱式歐拉法(后退歐拉法)：210ReImg在復(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ú)條件穩(wěn)定的第四十七頁(yè)，共69頁(yè)。梯形公式：在復(fù)平面上的絕對(duì)穩(wěn)定區(qū)域是是復(fù)平面的左半平面即也是無(wú)條件穩(wěn)定的相應(yīng)的絕對(duì)穩(wěn)定區(qū)間是第四十八頁(yè)，共69頁(yè)。龍格-庫(kù)塔法：而顯式1~4階方法的絕對(duì)穩(wěn)定區(qū)域?yàn)閗=1k=2k=3k=4-1-2-3---123ReImg第四十九頁(yè)，共69頁(yè)。例設(shè)是求解微分方程的單步法，分析它的穩(wěn)定性.解:所以絕對(duì)穩(wěn)定區(qū)域是即為復(fù)平面的左半平面.在實(shí)數(shù)域上是無(wú)條件穩(wěn)定的.第五十頁(yè)，共69頁(yè)。解：將代入得即例討論求解初值問(wèn)題的求解公式：的穩(wěn)定性.(λ>0為實(shí)數(shù))所以絕對(duì)穩(wěn)定區(qū)域是所以因此是條件穩(wěn)定的.第五十一頁(yè)，共69頁(yè)。四、線(xiàn)性多步法在逐步推進(jìn)的求解過(guò)程中，計(jì)算yn+1之前已經(jīng)求出了一系列的近似值y0,y1,…,yn,如果充分利用前面信息來(lái)預(yù)測(cè)yn+1，則可期望會(huì)獲得較高的精度，這就是線(xiàn)性多步法的基本思想.1線(xiàn)性多步法的一般公式最常用的線(xiàn)性多步法公式為其中為常數(shù)，yn-k為y(xn-k)的近似值fn-k=f(xn-k,yn-k)第五十二頁(yè)，共69頁(yè)。特別的當(dāng)時(shí),上式為顯式,否則是隱式.

設(shè)y(x)是微分方程初值問(wèn)題的準(zhǔn)確解，定義稱(chēng)為線(xiàn)性多步法在xn+1上的局部截?cái)嗾`差.若則稱(chēng)該方法具有p階精度.若則稱(chēng)局部截?cái)嗾`差的主項(xiàng)為為誤差常數(shù).第五十三頁(yè)，共69頁(yè)。例設(shè)yn+1=yn-1+2hf(xn,yn)為求解常微分初值問(wèn)題的線(xiàn)性二步法，試求該二步公式的局部截?cái)嗾`差主項(xiàng)，和精度.解:由局部截?cái)嗾`差的定義可知考慮在xn處的Taylor展式第五十四頁(yè)，共69頁(yè)。代入可得所以局部截?cái)嗾`差的主項(xiàng)為具有二階精度.例試建立求解為微分方程初值問(wèn)題具有如下形式的線(xiàn)性二步法，并使該方法具有二階精度，同時(shí)求其局部截?cái)嗾`差的主項(xiàng).解:局部截?cái)嗾`差為第五十五頁(yè)，共69頁(yè)。考慮在xn處的Taylor展式于是為使方法具有二階精度則第五十六頁(yè)，共69頁(yè)。解得因此該方法為局部截?cái)嗾`差的主項(xiàng)為例試建立求解為微分方程初值問(wèn)題具有如下形式的線(xiàn)性二步法，并使該方法具有三階精度，同時(shí)求其局部截?cái)嗾`差的主項(xiàng).解:局部截?cái)嗾`差為第五十七頁(yè)，共69頁(yè)?？紤]在xn處的Taylor展式第五十八頁(yè)，共69頁(yè)。所以解得局部截?cái)嗾`差的主項(xiàng)為第五十九頁(yè)，共69頁(yè)。2Adams外推公式考慮用r+1個(gè)點(diǎn)(xn-k,f(xn-k,yn-k))構(gòu)造一個(gè)r次多項(xiàng)式來(lái)近似的被積函數(shù)f(x,y(x)),這里用yn-k作為y(xn-k)的近似值，令fn-k=f(xn-k,yn-