精選數(shù)值分析第九章講義

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

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

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

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

階精度.定義現(xiàn)在是7頁\一共有68頁\編輯于星期日Lipschitiz條件：若存在正數(shù)L,使得對一切x,y1,y2有則稱f(x,y)滿足Lipschitiz條件.歐拉法的總體截斷誤差：那么設(shè)為局部截斷誤差，所以現(xiàn)在是8頁\一共有68頁\編輯于星期日現(xiàn)在是9頁\一共有68頁\編輯于星期日特別當(dāng)n=m-1時，有總體誤差與h是同階的.上式還說明，當(dāng)時，有即也就是說，ym收斂到方程的準(zhǔn)確解現(xiàn)在是10頁\一共有68頁\編輯于星期日后退Euler公式（隱式歐拉法）（隱式歐拉公式）利用向后差商近似導(dǎo)數(shù)現(xiàn)在是11頁\一共有68頁\編輯于星期日由于未知數(shù)yn+1同時出現(xiàn)在等式的兩邊，不能直接得到，故稱為隱式歐拉公式，而前者稱為顯式歐拉公式.一般先用顯式計算一個初值，再迭代求解.隱式歐拉法的局部截斷誤差：即隱式歐拉公式具有1階精度.現(xiàn)在是12頁\一共有68頁\編輯于星期日2梯形公式和改進Euler方法梯形公式設(shè)y=y(x)是的解,故由此得到用梯形公式近似現(xiàn)在是13頁\一共有68頁\編輯于星期日用yn來近似y(xn),用yn+1來近似y(xn+1),得梯形公式梯形公式是隱式的，可以用迭代法求解.現(xiàn)在是14頁\一共有68頁\編輯于星期日具有2階精度.梯形公式的局部截斷誤差現(xiàn)在是15頁\一共有68頁\編輯于星期日中點歐拉公式中心差商近似導(dǎo)數(shù)假設(shè),則可以導(dǎo)出即中點公式具有2階精度.需要2個初值y0和y1來啟動遞推過程，這樣的算法稱為雙步法，而前面的三種算法都是單步法.現(xiàn)在是16頁\一共有68頁\編輯于星期日方法顯式歐拉隱式歐拉梯形公式中點公式簡單精度低穩(wěn)定性最好精度低,計算量大精度提高計算量大精度提高,顯式多一個初值,可能影響精度有沒有一種方法，既有這些方法的優(yōu)點，而沒有它們的缺點？現(xiàn)在是17頁\一共有68頁\編輯于星期日改進歐拉法(1)先用顯式歐拉公式作預(yù)測，算出(2)再將代入梯形公式的右邊作校正,得到注：此法亦稱為預(yù)測-校正法.可以證明該算法具有2階精度，同時可以看到它是個單步遞推格式，比隱式公式的迭代求解過程簡單.現(xiàn)在是18頁\一共有68頁\編輯于星期日例用梯形公式求解初值問題(步長h=0.2)解：梯形公式為于是整理得由y(1)=y0=2依次可得y1,y2,y3,y4,y5.現(xiàn)在是19頁\一共有68頁\編輯于星期日例用改進歐拉法求解初值問題要求步長h=0.2，并計算y(1.2)和y(1.4)解：改進歐拉法公式為即現(xiàn)在是20頁\一共有68頁\編輯于星期日由y(1)=y0=1計算得現(xiàn)在是21頁\一共有68頁\編輯于星期日二、Runge-Kutta方法建立高精度的單步遞推格式.單步遞推法的基本思想是從(xn,yn)點出發(fā)，以某一斜率沿直線達到(xn+1

,yn+1

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

嗎？現(xiàn)在是22頁\一共有68頁\編輯于星期日首先希望能確定系數(shù)1、2、p，使得到的算法格式有2階精度，即在的前提假設(shè)下，使得

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

將K2在(xn

,yn

)點作Taylor展開1二階Runge-Kutta方法現(xiàn)在是23頁\一共有68頁\編輯于星期日(2)

將K2代入第1式，得到現(xiàn)在是24頁\一共有68頁\編輯于星期日(3)將yn+1與y(xn+1)在xn點的泰勒展開作比較要求,則必須有：這里有個未知數(shù)，個方程。32所以存在無窮多個解！現(xiàn)在是25頁\一共有68頁\編輯于星期日所有滿足上式的統(tǒng)稱為2階Runge-Kutta格式.若則改進的歐拉方法若則中點公式現(xiàn)在是26頁\一共有68頁\編輯于星期日2四階Runge-Kutta方法其中i(i=1,…,m)，i(i=2,…,m)

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

)

均為待定系數(shù)，確定這些系數(shù)的步驟與前面相似.現(xiàn)在是27頁\一共有68頁\編輯于星期日由于方程的個數(shù)少于未知量的個數(shù)，所以方程有無窮多個解，可以根據(jù)情況得到幾種常用的解，即得到相應(yīng)的四階公式.最常用為四階經(jīng)典龍格-庫塔法也稱為標(biāo)準(zhǔn)四階龍格-庫塔公式現(xiàn)在是28頁\一共有68頁\編輯于星期日Gill公式現(xiàn)在是29頁\一共有68頁\編輯于星期日753可達到的最高精度642每步須算Ki的個數(shù)(2)龍格-庫塔法的導(dǎo)出基于泰勒展開，故精度主要受解函數(shù)的光滑性影響.對于光滑性不太好的解,最好采用低階算法而將步長h

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

的值,即計算f的值.Butcher于1965年給出了計算量與可達到的最高精度階數(shù)的關(guān)系：現(xiàn)在是30頁\一共有68頁\編輯于星期日例用標(biāo)準(zhǔn)四階Runge-Kutta法求初值問題在x=0.1處的近似值，取步長為h=0.1.解：所以現(xiàn)在是31頁\一共有68頁\編輯于星期日那么例用標(biāo)準(zhǔn)四階Runge-Kutta法求初值問題在x=0.4處的近似值，取步長為h=0.2.現(xiàn)在是32頁\一共有68頁\編輯于星期日解：所以而所以現(xiàn)在是33頁\一共有68頁\編輯于星期日

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

(同時i)時有yi

y(xi

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

.顯式單步法的共同特征是它們都是將yn加上某種形式的增量，得出yn+1,計算公式如下：增量函數(shù)現(xiàn)在是34頁\一共有68頁\編輯于星期日Euler方法的增量函數(shù)改進Euler方法的增量函數(shù)

設(shè)y(x)是微分方程初值問題的準(zhǔn)確解，定義則稱為顯式單步法在xn+1處的局部截斷誤差.現(xiàn)在是35頁\一共有68頁\編輯于星期日例：考察歐拉顯式格式的收斂性:解：該問題的精確解為

歐拉公式為對任意固定的x=xi=ih，有現(xiàn)在是36頁\一共有68頁\編輯于星期日

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

總有即方法收斂.現(xiàn)在是38頁\一共有68頁\編輯于星期日在f(x,y)對y滿足Lipschitz條件下,Euler法，改進Euler法和Runge-Kutta法的增量函數(shù)

都對y滿足Lipschitz條件，所以上述結(jié)論對這些方法都成立.例設(shè)是求解微分方程的單步法，試求其局部截斷誤差的主項，并說出它具有幾階精度.解:現(xiàn)在是39頁\一共有68頁\編輯于星期日考慮在xn處的Taylor展式所以該方法的局部截斷誤差的主項是具有一階精度.現(xiàn)在是40頁\一共有68頁\編輯于星期日例設(shè)試求出它具有幾階精度.解:考慮在xn處的Taylor展式現(xiàn)在是41頁\一共有68頁\編輯于星期日所以該方法的局部截斷誤差的主項是具有二階精度.2單步法的穩(wěn)定性收斂性是在假定每一步計算都準(zhǔn)確的前提下,討論步長時，方法的總體截斷誤差是否趨于零的問題.穩(wěn)定性是討論舍入誤差的積累能否對計算結(jié)果有嚴(yán)重的影響.現(xiàn)在是42頁\一共有68頁\編輯于星期日例：考察初值問題在區(qū)間0.00.10.20.30.40.5精確解改進歐拉法

歐拉隱式歐拉顯式

節(jié)點xi

1.00002.00004.00008.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ù)值解.現(xiàn)在是43頁\一共有68頁\編輯于星期日

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

為復(fù)數(shù)且Re(λ)<0設(shè)在節(jié)點值yn處有擾動令那么于是現(xiàn)在是44頁\一共有68頁\編輯于星期日反復(fù)應(yīng)用可得為使則可得如下定義

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

的區(qū)域稱為方法的絕對穩(wěn)定區(qū)域,它與實軸的交稱為絕對穩(wěn)定區(qū)間定義下面討論已知的幾種方法的絕對穩(wěn)定區(qū)間和絕對穩(wěn)定區(qū)域.現(xiàn)在是45頁\一共有68頁\編輯于星期日顯式歐拉法:0-1-2ReImg在復(fù)平面上的絕對穩(wěn)定區(qū)域是即以-1為中心，1為半徑的圓域所以相應(yīng)的絕對穩(wěn)定區(qū)間是現(xiàn)在是46頁\一共有68頁\編輯于星期日隱式歐拉法(后退歐拉法)：210ReImg在復(fù)平面上的絕對穩(wěn)定區(qū)域是是以1為中心，1為半徑的圓的外域所以相應(yīng)的絕對穩(wěn)定區(qū)間是即如果只考慮λ

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

設(shè)y(x)是微分方程初值問題的準(zhǔn)確解，定義稱為線性多步法在xn+1上的局部截斷誤差.若則稱該方法具有p階精度.若則稱局部截斷誤差的主項為為誤差常數(shù).現(xiàn)在是53頁\一共有68頁\編輯于星期日例設(shè)yn+1=yn-1+2hf(xn,yn)為求解常微分初值問題的線性二步法，試求該二步公式的局部截斷誤差主項，和精度.解:由局部截斷誤差的定義可知考慮在xn處的Taylor展式現(xiàn)在是54頁\一共有68頁\編輯于星期日代入可得所以局部截斷誤差的主項為具有二階精度.例試建立求解為微分方程初值問題具有如下形式的線性二步法，并使該方法具有二階精度，同時求其局部截斷誤差的主項.解:局部截斷誤差為現(xiàn)在是55頁\一共有68頁\編輯于星期日考慮在xn處的Taylor展式于是為使方法具有二階精度則現(xiàn)在是56頁\一共有68頁\編輯于星期日解得因此該方法為局部截斷誤差的主項為例試建立求解為微分方程初值問題具有如下形式的線性二步法，并使該方法具有三階精度，同時求其局部截斷誤差的主項.解:局部截斷誤差為現(xiàn)在是57頁\一共有68頁\編輯于星期日考慮在xn處的Taylor展式現(xiàn)在是58頁\一共有68頁\編輯于星期日所以解得局部截斷誤差的主項為現(xiàn)在是59頁\一共有68頁\編輯于星期日2Adams外推公式考慮用r+1個點(xn-k,f(xn-k,yn-k))構(gòu)造一個r次多項式來近似的被積函數(shù)f(x,y(x)),這里用yn-k作為y(xn-k)的近似值，令fn-k=f(xn-k,yn-