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

在假設yn=y(xn)，即第i步計算是精確的前提下，考慮的截斷誤差Rn=y(xn+1)

yn+1稱為局部截斷誤差.定義歐拉法的局部截斷誤差：所以歐拉法具有1階精度.

若某算法的局部截斷誤差為O(hp+1)，則稱該算法有p

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

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

,yn+1

)點.歐拉法及其各種變形所能達到的最高精度為2階.考察改進的歐拉法，可以將其改寫為：斜率一定取K1K2的平均值嗎？步長一定是一個h

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

1、

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

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

將K2在(xn

,yn

)點作Taylor展開1二階Runge-Kutta方法第二十三頁，共69頁。(2)

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

i(i=1,…,m)，

i(i=2,…,m)

和

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

)

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

取小.注：(1)龍格-庫塔法的主要運算在于計算Ki

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

若某算法對于任意固定的x=xi=x0+ih,當h0

(同時i

)時有yi

y(xi

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

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

設y(x)是微分方程初值問題的準確解，定義則稱為顯式單步法在xn+1處的局部截斷誤差.第三十五頁，共69頁。例：考察歐拉顯式格式的收斂性:解：該問題的精確解為

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

第三十六頁，共69頁。

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

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

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

歐拉隱式歐拉顯式

節(jié)點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]上的解，分別用歐拉顯、隱式格式和改進的歐拉格式計算數(shù)值解.第四十三頁，共69頁。

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

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

若中，則稱單步法是絕對穩(wěn)定的.在復平面上，λh滿足

的區(qū)域稱為方法的絕對穩(wěn)定區(qū)域,它與實軸的交稱為絕對穩(wěn)定區(qū)間定義下面討論已知的幾種方法的絕對穩(wěn)定區(qū)間和絕對穩(wěn)定區(qū)域.第四十五頁，共69頁。顯式歐拉法:0-1-2ReImg在復平面上的絕對穩(wěn)定區(qū)域是即以-1為中心，1為半徑的圓域所以相應的絕對穩(wěn)定區(qū)間是第四十六頁，共69頁。隱式歐拉法(后退歐拉法)：210ReImg在復平面上的絕對穩(wěn)定區(qū)域是是以1為中心，1為半徑的圓的外域所以相應的絕對穩(wěn)定區(qū)間是即如果只考慮λ

<0的實數(shù)，則相應的絕對穩(wěn)定區(qū)間對于任意的h都成立，所以是無條件穩(wěn)定的第四十七頁，共69頁。梯形公式：在復平面上的絕對穩(wěn)定區(qū)域是是復平面的左半平面即也是無條件穩(wěn)定的相應的絕對穩(wěn)定區(qū)間是第四十八頁，共69頁。龍格-庫塔法：而顯式1~4階方法的絕對穩(wěn)定區(qū)域為k=1k=2k=3k=4-1-2-3---123ReImg第四十九頁，共69頁。例設是求解微分方程的單步法，分析它的穩(wěn)定性.解:所以絕對穩(wěn)定區(qū)域是即為復平面的左半平面.在實數(shù)域上是無條件穩(wěn)定的.第五十頁，共69頁。解：將代入得即例討論求解初值問題的求解公式：的穩(wěn)定性.(λ>0為實數(shù))所以絕對穩(wěn)定區(qū)域是所以因此是條件穩(wěn)定的.第五十一頁，共69頁。四、線性多步法在逐步推進的求解過程中，計算yn+1之前已經(jīng)求出了一系列的近似值y0,y1,…,yn,如果充分利用前面信息來預測yn+1，則可期望會獲得較高的精度，這就是線性多步法的基本思想.1線性多步法的一般公式最常用的線性多步法公式為其中為常數(shù)，yn-k為y(xn-k)的近似值fn-k=f(xn-k,yn-k)第五十二頁，共69頁。特別的當時,上式為顯式,否則是隱式.

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