常微分方程數(shù)值解

常微分方程數(shù)值解第一頁(yè)，共一百五十一頁(yè)，2022年，8月28日1第二頁(yè)，共一百五十一頁(yè)，2022年，8月28日29.1

引言考慮一階常微分方程的初值問題（1.1）（1.2）如果存在實(shí)數(shù)，使得（1.3）則稱關(guān)于滿足利普希茨(Lipschitz）條件，稱為的利普希茨常數(shù)（簡(jiǎn)稱Lips.常數(shù)）.第三頁(yè)，共一百五十一頁(yè)，2022年，8月28日3

定理1設(shè)在區(qū)域上連續(xù)，關(guān)于滿足利普希茨條件，則對(duì)任意，常微分方程(1.1),(1.2)式當(dāng)時(shí)存在唯一的連續(xù)可微解.關(guān)于解對(duì)擾動(dòng)的敏感性，有以下結(jié)論.

定理2設(shè)在區(qū)域（如定理1所定義）上連續(xù)，且關(guān)于滿足利普希茨條件，設(shè)初值問題的解為，則第四頁(yè)，共一百五十一頁(yè)，2022年，8月28日4兩者的區(qū)別:1.問題第五頁(yè)，共一百五十一頁(yè)，2022年，8月28日5第六頁(yè)，共一百五十一頁(yè)，2022年，8月28日6第七頁(yè)，共一百五十一頁(yè)，2022年，8月28日7左矩形右矩形梯形公式單步法:對(duì)初值問題,計(jì)算yn+1時(shí)只用到前一點(diǎn)的值yn,即yn+1=f(yn)k步法:計(jì)算yn+1時(shí)需要用到前k點(diǎn)的值yn,,yn-1,…,yn-k+1,即yn+1=f(yn,yn-1,…,yn-k+1)對(duì)方程離散化，建立求數(shù)值解的遞推公式.描述這類算法，只要給出用已知信息計(jì)算的遞推公式.第八頁(yè)，共一百五十一頁(yè)，2022年，8月28日89.2簡(jiǎn)單的數(shù)值方法

求解一階微分方程初值問題:第九頁(yè)，共一百五十一頁(yè)，2022年，8月28日9幾何意義:它是用一條自點(diǎn)(x0,y0)出發(fā)的折線段去逼近積分曲線y=y(x)如下圖9-1P280第十頁(yè)，共一百五十一頁(yè)，2022年，8月28日10

例1求解初值問題

（2.2）

解歐拉公式的具體形式為

取步長(zhǎng)，計(jì)算結(jié)果見表9-1.

初值問題（2.2）的解為，按這個(gè)解析式子算出的準(zhǔn)確值同近似值一起列在表9-1中，兩者相比較可以看出歐拉方法的精度很差.第十一頁(yè)，共一百五十一頁(yè)，2022年，8月28日11還可以通過幾何直觀來考察歐拉方法的精度.假設(shè)，即頂點(diǎn)落在積分曲線上，那么，按歐拉方法做出的折線便是過點(diǎn)的切線（圖9-2）.第十二頁(yè)，共一百五十一頁(yè)，2022年，8月28日12圖9-2從圖形上看，這樣定出的頂點(diǎn)顯著地偏離了原來的積分曲線，可見歐拉方法是相當(dāng)粗糙的.誤差分析：為了分析計(jì)算公式的精度,通?？捎锰├照归_將在處展開，則有第十三頁(yè)，共一百五十一頁(yè)，2022年，8月28日13在的前提下，稱為此方法的局部截?cái)嗾`差.于是可得歐拉法（2.1）的誤差（2.3）（2.1）估算=精確第十四頁(yè)，共一百五十一頁(yè)，2022年，8月28日14（2.5）稱為后退的歐拉法(隱式歐拉公式).歐拉公式是關(guān)于的一個(gè)直接的計(jì)算公式，這類公式稱作是顯式的；后退歐拉公式的右端含有未知的，它是關(guān)于的一個(gè)函數(shù)方程，這類公式稱作是隱式的.

第十五頁(yè)，共一百五十一頁(yè)，2022年，8月28日15隱式方程通常用迭代法求解，而迭代過程的實(shí)質(zhì)是逐步顯示化.設(shè)用歐拉公式給出迭代初值，用它代入（2.5）式的右端，使之轉(zhuǎn)化為顯式，直接計(jì)算得然后再用代入（2.5）式，又有第十六頁(yè)，共一百五十一頁(yè)，2022年，8月28日16如此反復(fù)進(jìn)行，得（2.6）由于對(duì)滿足利普希茨條件(1.3).由(2.6)減(2.5)得由此可知，只要迭代法(2.6)就收斂到解.第十七頁(yè)，共一百五十一頁(yè)，2022年，8月28日17

梯形方法（2.7）稱為梯形方法.梯形方法是隱式單步法，可用迭代法求解.第十八頁(yè)，共一百五十一頁(yè)，2022年，8月28日18為了分析迭代過程的收斂性，將（2.7）與（2.8）式相減，得（2.8）同后退的歐拉方法一樣，仍用歐拉方法提供迭代初值，則梯形法的迭代公式為（2.7）第十九頁(yè)，共一百五十一頁(yè)，2022年，8月28日19如果選取充分小，使得則當(dāng)時(shí)有，這說明迭代過程(2.8)是收斂的.于是有式中為關(guān)于的利普希茨常數(shù).第二十頁(yè)，共一百五十一頁(yè)，2022年，8月28日20

改進(jìn)歐拉公式梯形方法雖然提高了精度，但其算法復(fù)雜.在應(yīng)用迭代公式（2.8）進(jìn)行實(shí)際計(jì)算時(shí)，每迭代一次，都要重新計(jì)算函數(shù)的值.為了控制計(jì)算量，通常只迭代一兩次就轉(zhuǎn)入下一步的計(jì)算，這就簡(jiǎn)化了算法.具體地，先用歐拉公式求得一個(gè)初步的近似值,而迭代又要反復(fù)進(jìn)行若干次，計(jì)算量很大，而且往往難以預(yù)測(cè).稱之為預(yù)測(cè)值，第二十一頁(yè)，共一百五十一頁(yè)，2022年，8月28日21這樣建立的預(yù)測(cè)-校正系統(tǒng)通常稱為改進(jìn)的歐拉公式：預(yù)測(cè)值的精度可能很差，再用梯形公式（2.7）將它校正一次，即按（2.8）式迭代一次得，這個(gè)結(jié)果稱校正值.預(yù)測(cè)校正（2.9）也可以表為下列平均化形式（2.7）（2.8）第二十二頁(yè)，共一百五十一頁(yè)，2022年，8月28日22

例2用改進(jìn)的歐拉方法求解初值問題（2.2）.

解這里改進(jìn)的歐拉公式為（2.2）第二十三頁(yè)，共一百五十一頁(yè)，2022年，8月28日23仍取，計(jì)算結(jié)果見表9-2.同例1中歐拉法的計(jì)算結(jié)果比較，改進(jìn)歐拉法明顯改善了精度.第二十四頁(yè)，共一百五十一頁(yè)，2022年，8月28日24

單步法的局部截?cái)嗾`差與階初值問題(1.1)，(1.2)的單步法可用一般形式表示為（2.10）其中多元函數(shù)與有關(guān)，當(dāng)含有時(shí)，方法是隱式的，若不含則為顯式方法，（2.11）稱為增量函數(shù)，所以顯式單步法可表示為例如對(duì)歐拉法（2.1）有它的局部截?cái)嗾`差已由(2.3)給出.（1.1）（1.2）（2.1）（2.3）第二十五頁(yè)，共一百五十一頁(yè)，2022年，8月28日25對(duì)一般顯式單步法則可如下定義.

定義1設(shè)是初值問題(1.1)，(1.2)的準(zhǔn)確解，稱（2.12）為顯式單步法（2.11）的局部截?cái)嗾`差.之所以稱為局部的，是假設(shè)在前各步?jīng)]有誤差.當(dāng)時(shí)，計(jì)算一步，則有（1.1）（1.2）（2.11）第二十六頁(yè)，共一百五十一頁(yè)，2022年，8月28日26在前一步精確的情況下用公式（2.11）計(jì)算產(chǎn)生的公式誤差.根據(jù)定義，歐拉法的局部截?cái)嗾`差即為（2.3）的結(jié)果.這里稱為局部截?cái)嗾`差主項(xiàng).局部截?cái)嗾`差可理解為用方法（2.11）計(jì)算一步的誤差，即顯然（2.11）（2.3）第二十七頁(yè)，共一百五十一頁(yè)，2022年，8月28日27

定義2設(shè)是初值問題(1.1)，(1.2)的準(zhǔn)確解，若存在最大整數(shù)使顯式單步法(2.11)的局部截?cái)嗾`差滿足

（2.13）則稱方法（2.11）具有階精度.若將（2.13）展開式寫成則稱為局部截?cái)嗾`差主項(xiàng).以上定義對(duì)隱式單步法（2.10）也是適用的.（1.1）（1.2）（2.11）（2.10）第二十八頁(yè)，共一百五十一頁(yè)，2022年，8月28日28對(duì)后退歐拉法（2.5）其局部截?cái)嗾`差為這里，是1階方法，局部截?cái)嗾`差主項(xiàng)為.（2.5）第二十九頁(yè)，共一百五十一頁(yè)，2022年，8月28日29對(duì)梯形法（2.7）有所以梯形方法是二階的，其局部誤差主項(xiàng)為（2.7）第三十頁(yè)，共一百五十一頁(yè)，2022年，8月28日309.3

龍格-庫(kù)塔方法第三十一頁(yè)，共一百五十一頁(yè)，2022年，8月28日31一、Taylor展開法取等式右邊前p+1項(xiàng)第三十二頁(yè)，共一百五十一頁(yè)，2022年，8月28日32例取h=0.1,用三階Taylor展開法求解第三十三頁(yè)，共一百五十一頁(yè)，2022年，8月28日33從計(jì)算高階導(dǎo)數(shù)的公式知道,方法的截?cái)嗾`差提高一階,需要增加的計(jì)算量很大.下面我們用區(qū)間上若干點(diǎn)的導(dǎo)數(shù)f,而不是高階導(dǎo)數(shù),將它們作線性組合得到平均斜率,將其與解的Taylor展開相比較,使前面若干項(xiàng)吻合，從而得到提高階的方法第三十四頁(yè)，共一百五十一頁(yè)，2022年，8月28日342龍格-庫(kù)塔法(Runge-Kutta法)第三十五頁(yè)，共一百五十一頁(yè)，2022年，8月28日35（2）龍格-庫(kù)塔法的一般形式第三十六頁(yè)，共一百五十一頁(yè)，2022年，8月28日36第三十七頁(yè)，共一百五十一頁(yè)，2022年，8月28日37第三十八頁(yè)，共一百五十一頁(yè)，2022年，8月28日38將以上結(jié)果代入局部截?cái)嗾`差公式則有要使公式（3.6）具有階，必須使（3.6）第三十九頁(yè)，共一百五十一頁(yè)，2022年，8月28日39即非線性方程組（3.9）的解是不唯一的.令，則得這樣得到的公式稱為二階R-K方法，如取，則這就是改進(jìn)歐拉法（3.1）.（3.9）第四十頁(yè)，共一百五十一頁(yè)，2022年，8月28日40若取，則得計(jì)算公式.

稱為中點(diǎn)公式，相當(dāng)于數(shù)值積分的中矩形公式.

（3.10）也可表示為（3.10）第四十一頁(yè)，共一百五十一頁(yè)，2022年，8月28日41

三階與四階顯式R-K方法要得到三階顯式R-K方法，必須.（3.11）其中及均為待定參數(shù).此時(shí)（3.4）,（3.5）的公式表示為公式（3.11）的局部截?cái)嗾`差為（3.4）（3.5）第四十二頁(yè)，共一百五十一頁(yè)，2022年，8月28日42只要將按二元函數(shù)泰勒展開，使，可得待定參數(shù)滿足方程（3.12）第四十三頁(yè)，共一百五十一頁(yè)，2022年，8月28日43這是8個(gè)未知數(shù)6個(gè)方程的方程組，解也不是唯一的.所以這是一簇公式.滿足條件（3.12）的公式（3.11）統(tǒng)稱為三階R-K公式.一個(gè)常見的公式為此公式稱為庫(kù)塔三階方法.第四十四頁(yè)，共一百五十一頁(yè)，2022年，8月28日44繼續(xù)上述過程，經(jīng)過較復(fù)雜的數(shù)學(xué)演算，可以導(dǎo)出各種四階龍格-庫(kù)塔公式，下列經(jīng)典公式是其中常用的一個(gè)：可以證明其截?cái)嗾`差為.四階龍格-庫(kù)塔方法的每一步需要計(jì)算四次函數(shù)值，（3.13）第四十五頁(yè)，共一百五十一頁(yè)，2022年，8月28日45

謝謝！第四十六頁(yè)，共一百五十一頁(yè)，2022年，8月28日46第四十七頁(yè)，共一百五十一頁(yè)，2022年，8月28日47

變步長(zhǎng)的龍格-庫(kù)塔方法單從每一步看，步長(zhǎng)越小，截?cái)嗾`差就越小，但隨著步長(zhǎng)的縮小，在一定求解范圍內(nèi)所要完成的步數(shù)就增加了.步數(shù)的增加不但引起計(jì)算量的增大，而且可能導(dǎo)致舍入誤差的嚴(yán)重積累.因此同積分的數(shù)值計(jì)算一樣，微分方程的數(shù)值解法也有個(gè)選擇步長(zhǎng)的問題.在選擇步長(zhǎng)時(shí)，需要考慮兩個(gè)問題：1°怎樣衡量和檢驗(yàn)計(jì)算結(jié)果的精度？第四十八頁(yè)，共一百五十一頁(yè)，2022年，8月28日482°如何依據(jù)所獲得的精度處理步長(zhǎng)？考察經(jīng)典的四階龍格-庫(kù)塔公式（3.13）從節(jié)點(diǎn)出發(fā)，先以為步長(zhǎng)求出一個(gè)近似值，第四十九頁(yè)，共一百五十一頁(yè)，2022年，8月28日49（3.14）然后將步長(zhǎng)折半，即取為步長(zhǎng)從跨兩步到，再求得一個(gè)近似值，每跨一步的截?cái)嗾`差是，因此有（3.15）比較（3.14）式和（3.15）式我們看到，步長(zhǎng)折半后，由于公式的局部截?cái)嗾`差為，故有誤差大約減少到，第五十頁(yè)，共一百五十一頁(yè)，2022年，8月28日50由此易得下列事后估計(jì)式這樣，可以通過檢查步長(zhǎng)，折半前后兩次計(jì)算結(jié)果的偏差即有來判定所選的步長(zhǎng)是否合適.具體地說，將區(qū)分以下兩種情況處理：第五十一頁(yè)，共一百五十一頁(yè)，2022年，8月28日511.對(duì)于給定的精度，如果，反復(fù)將步長(zhǎng)折半進(jìn)行計(jì)算，直至為止.這時(shí)取最終得到的作為結(jié)果；2.如果，反復(fù)將步長(zhǎng)加倍，直到為止，這種通過加倍或折半處理步長(zhǎng)的方法稱為變步長(zhǎng)方法.這時(shí)再將步長(zhǎng)折半一次，就得到所要的結(jié)果.表面上看，為了選擇步長(zhǎng)，每一步的計(jì)算量增加了，但總體考慮往往是合算的.第五十二頁(yè)，共一百五十一頁(yè)，2022年，8月28日529.4

單步法的收斂性與穩(wěn)定性

收斂性與相容性

數(shù)值解法的基本思想是通過某種離散化手段將微分方程轉(zhuǎn)化為差分方程，如單步法(2.11)，即

（4.1）它在處的解為，而初值問題（1.1）,（1.2）在處的精確解為，記稱為整體截?cái)嗾`差.（1.1）（1.2）第五十三頁(yè)，共一百五十一頁(yè)，2022年，8月28日53收斂性就是討論當(dāng)固定且時(shí)的問題.

定義3若一種數(shù)值方法對(duì)于固定的,當(dāng)時(shí)有，其中是(1.1),(1.2)的準(zhǔn)確解，則稱該方法是收斂的.

顯然數(shù)值方法收斂是指.對(duì)單步法（4.1）有下述收斂性定理：（1.1）（1.2）（4.1）第五十四頁(yè)，共一百五十一頁(yè)，2022年，8月28日54

定理3假設(shè)單步法（4.1）具有階精度，且增量函數(shù)關(guān)于滿足利普希茨條件（4.2）又設(shè)初值是準(zhǔn)確的，即，則其整體截?cái)嗾`差

（4.3）

證明設(shè)以表示取用公式（4.1）求得的結(jié)果，即（4.4）則為局部截?cái)嗾`差，（4.1）（4.1）第五十五頁(yè)，共一百五十一頁(yè)，2022年，8月28日55由于所給方法具有階精度，按定義2，存在定數(shù)，使又由式（4.4）與（4.1），得利用假設(shè)條件（4.2），有從而有（4.2）（4.1）（4.4）第五十六頁(yè)，共一百五十一頁(yè)，2022年，8月28日56即對(duì)整體截?cái)嗾`差成立下列遞推關(guān)系式（4.5）反復(fù)遞推，可得（4.6）再注意到當(dāng)時(shí)最終得下列估計(jì)式（4.7）第五十七頁(yè)，共一百五十一頁(yè)，2022年，8月28日57由此可以斷定，如果初值是準(zhǔn)確的，即，則（4.3）式成立.依據(jù)這一定理，判斷單步法（4.1）的收斂性，歸結(jié)為驗(yàn)證增量函數(shù)能否滿足利普希茨條件（4.2）.對(duì)于歐拉方法，由于其增量函數(shù)就是，故當(dāng)關(guān)于滿足利普希茨條件時(shí)它是收斂的.再考察改進(jìn)的歐拉方法，其增量函數(shù)給出，這時(shí)有（4.3）（4.2）（4.1）第五十八頁(yè)，共一百五十一頁(yè)，2022年，8月28日58假設(shè)關(guān)于滿足利普希茨條件，記利普希茨常數(shù)為,設(shè)為定數(shù)），上式表明關(guān)于的利普希茨常數(shù)則由上式推得因此改進(jìn)的歐拉方法也是收斂的.第五十九頁(yè)，共一百五十一頁(yè)，2022年，8月28日59類似地，也可驗(yàn)證其他龍格-庫(kù)塔方法的收斂性.定理3表明時(shí)單步法收斂，并且當(dāng)是初值問題（1.1），（1.2）的解，（4.1）具有階精度時(shí)，有展開式所以的充要條件是，（4.1）（1.1）（1.2）第六十頁(yè)，共一百五十一頁(yè)，2022年，8月28日60而，于是可給出如下定義：

定義4若單步法（4.1）的增量函數(shù)滿足

則稱單步法（4.1）與初值問題（1.1），（1.2）相容.相容性是指數(shù)值方法逼近微分方程(1.1)，即微分方程(1.1)離散化得到的數(shù)值方法當(dāng)時(shí)可得到第六十一頁(yè)，共一百五十一頁(yè)，2022年，8月28日61

定理4階方法(4.1)與初值問題(1.1)，(1.2)相容的充分必要條件是由定理3可知單步法（4.1）收斂的充分必要條件是(4.1)是相容的.

以上討論表明階方法（4.1）當(dāng)時(shí)與(1.1),(1.2)相容，反之相容方法至少是1階的.

第六十二頁(yè)，共一百五十一頁(yè)，2022年，8月28日62

9.4.2絕對(duì)穩(wěn)定性與絕對(duì)穩(wěn)定域

定義5若一種數(shù)值方法在節(jié)點(diǎn)值上大小為的擾動(dòng)，于以后各節(jié)點(diǎn)值上產(chǎn)生的偏差均不超過，則稱該方法是穩(wěn)定的.以歐拉法為例考察計(jì)算穩(wěn)定性.

例4考察初值問題

其準(zhǔn)確解是一個(gè)按指數(shù)曲線衰減得很快的函數(shù)，如圖9-3所示.第六十三頁(yè)，共一百五十一頁(yè)，2022年，8月28日63若取，則歐拉公式的具體形式為計(jì)算結(jié)果列于表9-4的第2列.可以看到，歐拉方法的解（圖9-3中用×號(hào)標(biāo)出）在準(zhǔn)確值的上下波動(dòng)，計(jì)算過程明顯地不穩(wěn)定.圖9-3用歐拉法解方程得第六十四頁(yè)，共一百五十一頁(yè)，2022年，8月28日64再考察后退的歐拉方法，取時(shí)計(jì)算公式為計(jì)算結(jié)果列于表9-4的第3列（圖9-3中標(biāo)以·號(hào)），這時(shí)計(jì)算過程是穩(wěn)定的.但若取則計(jì)算過程穩(wěn)定.第六十五頁(yè)，共一百五十一頁(yè)，2022年，8月28日65這表明穩(wěn)定性不但與方法有關(guān)，也與步長(zhǎng)的大小有關(guān)，當(dāng)然也與方程中的有關(guān).

為了只考察數(shù)值方法本身，通常只檢驗(yàn)將數(shù)值方法用于解模型方程的穩(wěn)定性，（4.8）其中為復(fù)數(shù).例如在的鄰域，可展開為模型方程為對(duì)一般方程可以通過局部線性化化為這種形式.第六十六頁(yè)，共一百五十一頁(yè)，2022年，8月28日66略去高階項(xiàng)，再做變換即可得到的形式.對(duì)于個(gè)方程的方程組，也可線性化為，這里為的雅可比矩陣.若有個(gè)特征值，則還可能是復(fù)數(shù)，為保證微分方程本身的穩(wěn)定性，還應(yīng)假定.先研究歐拉方法的穩(wěn)定性.模型方程的歐拉公式為所以，為了使模型方程結(jié)果能推廣到方程組，方程中應(yīng)為復(fù)數(shù).第六十七頁(yè)，共一百五十一頁(yè)，2022年，8月28日67（4.9）設(shè)在節(jié)點(diǎn)值上有一擾動(dòng)值，它的傳播使節(jié)點(diǎn)值產(chǎn)生大小為的擾動(dòng)值，假設(shè)用按歐拉公式得出的計(jì)算過程不再有新的誤差，則擾動(dòng)值滿足可見擾動(dòng)值滿足原來的差分方程（4.9）.如果差分方程的解是不增長(zhǎng)的，即有則它就是穩(wěn)定的.第六十八頁(yè)，共一百五十一頁(yè)，2022年，8月28日68即圖9-4顯然，為要保證差分方程（4.9）的解不增長(zhǎng)，只要選取充分小，（4.10）在的復(fù)平面上，這是以為圓心，1為半徑的單位圓內(nèi)部（圖9-4）.這個(gè)圓域稱為歐拉法的絕對(duì)穩(wěn)定域，一般情形可由下面定義.

定義6單步法（4.1）用于解模型方程（4.8），若得到的解，滿足，則稱方法（4.1）是絕對(duì)穩(wěn)定的.（4.9）使（4.1）（4.8）第六十九頁(yè)，共一百五十一頁(yè)，2022年，8月28日69在的平面上，使的變量圍成的區(qū)域，稱為絕對(duì)穩(wěn)定域，對(duì)歐拉法，給出，絕對(duì)穩(wěn)定區(qū)間為.它與實(shí)軸的交稱為絕對(duì)穩(wěn)定區(qū)間.其絕對(duì)穩(wěn)定域由在例5中，即為絕對(duì)穩(wěn)定區(qū)間.當(dāng)取時(shí)例4中取故它是不穩(wěn)定的，它是穩(wěn)定的.第七十頁(yè)，共一百五十一頁(yè)，2022年，8月28日70故絕對(duì)穩(wěn)定域由得到.絕對(duì)穩(wěn)定區(qū)間為，即.類似可得三階及四階的R-K方法的分別為用二階R-K方法解模型方程可得到第七十一頁(yè)，共一百五十一頁(yè)，2022年，8月28日71由可得到相應(yīng)的絕對(duì)穩(wěn)定域.當(dāng)為實(shí)數(shù)時(shí)則得絕對(duì)穩(wěn)定區(qū)間.分別為三階顯式R-K方法：即四階顯式R-K方法：即圖9-5給出了R-K方法到的絕對(duì)穩(wěn)定域.從以上討論可知顯式的R-K方法的絕對(duì)穩(wěn)定域均為有限域，都對(duì)步長(zhǎng)有限制.如果不在所給的絕對(duì)穩(wěn)定區(qū)間內(nèi)，方法就不穩(wěn)定.圖9-5第七十二頁(yè)，共一百五十一頁(yè)，2022年，8月28日72例5分別取，及用經(jīng)典的四階R-K方法（3.13）計(jì)算.

解本例分別為及.前者在絕對(duì)穩(wěn)定區(qū)間內(nèi)，后者則不在.經(jīng)典的四階R-K方法為（3.13）第七十三頁(yè)，共一百五十一頁(yè)，2022年，8月28日73以上結(jié)果看到，如果步長(zhǎng)不滿足絕對(duì)穩(wěn)定條件，誤差增長(zhǎng)很快.用四階R-K方法計(jì)算其誤差見下表：這里第七十四頁(yè)，共一百五十一頁(yè)，2022年，8月28日74對(duì)隱式單步法，可以同樣討論方法的絕對(duì)穩(wěn)定性，例如對(duì)后退歐拉法，用它解模型方程可得故由可得絕對(duì)穩(wěn)定域?yàn)?它是以為圓心，1為半徑的單位圓外部，故絕對(duì)穩(wěn)定區(qū)間為.第七十五頁(yè)，共一百五十一頁(yè)，2022年，8月28日75當(dāng)時(shí)，則，即對(duì)任何步長(zhǎng)均為穩(wěn)定的.故對(duì)有，故絕對(duì)穩(wěn)定域?yàn)榈淖蟀肫矫?，?duì)隱式梯形法，解模型方程得第七十六頁(yè)，共一百五十一頁(yè)，2022年，8月28日76絕對(duì)穩(wěn)定區(qū)間為，即時(shí)梯形法均是穩(wěn)定的.

定義7如果數(shù)值方法的絕對(duì)穩(wěn)定域包含了那么稱此方法是A-穩(wěn)定的.隱式歐拉法與梯形方法的絕對(duì)穩(wěn)定域均為在具體計(jì)算中步長(zhǎng)的選擇只需考慮計(jì)算精度及迭代收斂性要求而不必考慮穩(wěn)定性，具有這種特點(diǎn)的方法特別重要.由定義知A-穩(wěn)定方法對(duì)步長(zhǎng)沒有限制.第七十七頁(yè)，共一百五十一頁(yè)，2022年，8月28日779.5

線性多步法在逐步推進(jìn)的求解過程中，計(jì)算之前事實(shí)上已經(jīng)求出了一系列的近似值，如果充分利用前面多步的信息來預(yù)測(cè)，則可以期望會(huì)獲得較高的精度.這就是構(gòu)造所謂線性多步法的基本思想.本節(jié)主要介紹基于泰勒展開的構(gòu)造方法.構(gòu)造多步法的主要途徑是基于數(shù)值積分方法和基于泰勒展開方法，前者可直接由方程兩端積分后利用插值求積公式得到.第七十八頁(yè)，共一百五十一頁(yè)，2022年，8月28日78

線性多步法的一般公式一般的線性多步法公式可表示為（5.1）其中為的近似，計(jì)算時(shí)需先給出前面?zhèn)€近似值,再由（5.1）逐次求出.如果計(jì)算時(shí)，除了使用的值，還用到的值，則稱此方法為線性多步法.為常數(shù).及不全為零，則稱為線性步法.若第七十九頁(yè)，共一百五十一頁(yè)，2022年，8月28日79這時(shí)可直接由（5.1）算出；如果，稱（5.1）為顯式步法，如果，則（5.1）稱為隱式步法，求解時(shí)與梯形法（2.7）相同，要用迭代法方可算出.（5.1）（2.7）第八十頁(yè)，共一百五十一頁(yè)，2022年，8月28日80

定義8設(shè)是初值問題(1.1)，(1.2)的準(zhǔn)確解，線性多步法（5.1）在上的局部截?cái)嗾`差為

（5.2）（5.1）中系數(shù)及可根據(jù)方法的局部截?cái)嗾`差及階確定，其定義為：（5.1）（1.1）（1.2）若，則稱方法（5.1）是階的，則稱方法（5.1）與方程（1.1）是相容的.第八十一頁(yè)，共一百五十一頁(yè)，2022年，8月28日81由定義8，對(duì)在處做泰勒展開，由于代入（5.2）得（5.3）（5.2）第八十二頁(yè)，共一百五十一頁(yè)，2022年，8月28日82其中（5.4）若在公式（5.1）中選擇系數(shù)及，使它滿足由定義可知此時(shí)所構(gòu)造的多步法是階的.（5.1）第八十三頁(yè)，共一百五十一頁(yè)，2022年，8月28日83（5.5）稱右端第一項(xiàng)為局部截?cái)嗾`差主項(xiàng)，稱為誤差常數(shù).根據(jù)相容性定義，,即，故方法（5.1）與微分方程（1.1）相容的充分必要條件是（5.6）成立.且（5.6）由（5.4）得（5.1）第八十四頁(yè)，共一百五十一頁(yè)，2022年，8月28日84當(dāng)時(shí)，若，則由（5.6）可求得此時(shí)公式（5.1）為即為歐拉法.從（5.4）可求得，故方法為1階精度，這和第2節(jié)給出的定義及結(jié)果是一致的.且局部截?cái)嗾`差為（5.6）第八十五頁(yè)，共一百五十一頁(yè)，2022年，8月28日85對(duì),若，方法為隱式公式.為了確定系數(shù)，可由解得于是得到公式即為梯形法.由（5.4）可求得,故，所以梯形法是二階方法，其局部截?cái)嗾`差主項(xiàng)是.這與第9.2節(jié)中的討論也是一致的.第八十六頁(yè)，共一百五十一頁(yè)，2022年，8月28日86對(duì)的多步法公式都可利用（5.4）確定系數(shù)，并由（5.5）給出局部截?cái)嗾`差.（5.5）第八十七頁(yè)，共一百五十一頁(yè)，2022年，8月28日87

阿當(dāng)姆斯顯式與隱式公式考慮形如（5.7）的步法，稱為阿當(dāng)姆斯（Adams)方法.為顯式方法，為隱式方法，通常稱為阿當(dāng)姆斯顯式與隱式公式，也稱Adams-Bashforth公式與Adams-Monlton公式.這類公式可直接由對(duì)方程兩端從到積分求得.第八十八頁(yè)，共一百五十一頁(yè)，2022年，8月28日88也可以利用(5.4)由推出，對(duì)比與可知此時(shí)系數(shù).顯然成立，下面只需確定系數(shù)，可令，求得.若，則令來求得.第八十九頁(yè)，共一百五十一頁(yè)，2022年，8月28日89以為例，由，根據(jù)若，則由前三個(gè)方程解得得到的阿當(dāng)姆斯顯式公式是（5.8）第九十頁(yè)，共一百五十一頁(yè)，2022年，8月28日90由（5.4）求得，所以（5.8）是三階方法，局部截?cái)嗾`差是若，則可解得于是得的阿當(dāng)姆斯隱式公式為（5.9）它是四階方法，局部截?cái)嗾`差是（5.8）第九十一頁(yè)，共一百五十一頁(yè)，2022年，8月28日91（5.10）類似的方法可求得阿當(dāng)姆斯顯式方法和隱式方法的公式，表9-5及表9-6分別列出了時(shí)的阿當(dāng)姆斯顯式公式與阿當(dāng)姆斯隱式公式，其中為步數(shù)，為方法的階，為誤差常數(shù).第九十二頁(yè)，共一百五十一頁(yè)，2022年，8月28日92第九十三頁(yè)，共一百五十一頁(yè)，2022年，8月28日93第九十四頁(yè)，共一百五十一頁(yè)，2022年，8月28日94

例6用四階阿當(dāng)姆斯顯式和隱式方法解初值問題

取步長(zhǎng).

解本題.從四階阿當(dāng)姆斯隱式公式得到從四階阿當(dāng)姆斯顯式公式得到第九十五頁(yè)，共一百五十一頁(yè)，2022年，8月28日95由此可直接解出而不用迭代，得到計(jì)算結(jié)果見表9-8，其中顯式方法中的初值及隱式方法中的均用準(zhǔn)確解得到，對(duì)一般方程，可用四階R-K方法計(jì)算初始近似.(表略，見教材.)第九十六頁(yè)，共一百五十一頁(yè)，2022年，8月28日96從以上例子看到同階的阿當(dāng)姆斯方法，隱式方法要比顯式方法誤差小.這可以從兩種方法的局部截?cái)嗾`差主項(xiàng)的系數(shù)大小得到解釋.這里分別為及.第九十七頁(yè)，共一百五十一頁(yè)，2022年，8月28日97

米爾尼方法與辛普森方法考慮另一個(gè)的顯式公式其中為待定常數(shù).由（5.4）可知，再令得到可根據(jù)使公式的階盡可能高這一條件來確定其數(shù)值.第九十八頁(yè)，共一百五十一頁(yè)，2022年，8月28日98解此方程組得于是得到四步顯式公式（5.11）稱為米爾尼(Milne）方法.由于，故方法為4階，其局部截?cái)嗾`差為（5.12）米爾尼方法也可以通過方程兩端積分第九十九頁(yè)，共一百五十一頁(yè)，2022年，8月28日99得到.右端積分通過辛普森求積公式就有（5.13）稱為辛普森方法.（5.14）它是隱式二步四階方法，其局部截?cái)嗾`差為若將方程從到積分，可得第一百頁(yè)，共一百五十一頁(yè)，2022年，8月28日100

漢明方法辛普森公式是二步方法中階數(shù)最高的，但它的穩(wěn)定性差，為了改善穩(wěn)定性，考察另一類三步法公式

其中系數(shù)及為常數(shù).若希望導(dǎo)出的公式是四階的，則系數(shù)中至少有一個(gè)自由參數(shù).若取，則可得到辛普森公式.若取，仍利用泰勒展開，由（5.4），令第一百零一頁(yè)，共一百五十一頁(yè)，2022年，8月28日101解此方程組得于是有（5.15）則可得到第一百零二頁(yè)，共一百五十一頁(yè)，2022年，8月28日102稱為漢明(Hamming)方法.由于，故方法是四階的，且局部截?cái)嗾`差（5.16）第一百零三頁(yè)，共一百五十一頁(yè)，2022年，8月28日103

預(yù)測(cè)-校正方法對(duì)于隱式的線性多步法，計(jì)算時(shí)要進(jìn)行迭代，計(jì)算量較大.為了避免進(jìn)行迭代，通常采用顯式公式給出的一個(gè)初始近似，記為，稱為預(yù)測(cè)(predictor).接著計(jì)算的值(evaluation)，再用隱式公式計(jì)算，稱為校正（corrector）.

在（2.13）中用歐拉法做預(yù)測(cè)，再用梯形法校正，得到改進(jìn)歐拉法，它就是一個(gè)二階預(yù)測(cè)-校正方法.（2.13）第一百零四頁(yè)，共一百五十一頁(yè)，2022年，8月28日104例如用四階的阿當(dāng)姆斯顯式方法做預(yù)測(cè)，再用四階阿當(dāng)姆斯隱式公式做校正，得到以下格式：一般情況下，預(yù)測(cè)公式與校正公式都取同階的顯式方法與隱式方法相匹配.預(yù)測(cè)P：求值E：校正C：求值E：此公式稱為阿當(dāng)姆斯四階預(yù)測(cè)-校正格式（PECE).第一百零五頁(yè)，共一百五十一頁(yè)，2022年，8月28日105依據(jù)四階阿當(dāng)姆斯公式的截?cái)嗾`差，對(duì)于PECE的預(yù)測(cè)步P有對(duì)校正步C有兩式相減得于是有下列事后誤差估計(jì)第一百零六頁(yè)，共一百五十一頁(yè)，2022年，8月28日106容易看出比更好.（5.17）但在的表達(dá)式中是未知的，因此計(jì)算時(shí)用上一步代替，從而構(gòu)造一種修正預(yù)測(cè)-校正格式(PMECME):第一百零七頁(yè)，共一百五十一頁(yè)，2022年，8月28日107P：M：E：C：M：E：在PMECME格式中已將（5.17）的及分別改為及.第一百零八頁(yè)，共一百五十一頁(yè)，2022年，8月28日108利用米爾尼公式（5.11）和漢明公式（5.15）相匹配，并利用截?cái)嗾`差（5.12），（5.16）改進(jìn)計(jì)算結(jié)果，可類似地建立四階修正米爾尼-漢明預(yù)測(cè)-校正格式PMECME）：P：M：E：C：M：E：（5.15）（5.11）（5.16）（5.12）第一百零九頁(yè)，共一百五十一頁(yè)，2022年，8月28日109

構(gòu)造多步法公式的注記和例構(gòu)造多步法公式有基于數(shù)值積分和泰勒展開兩種途徑，數(shù)值積分方法只能用于可轉(zhuǎn)化為等價(jià)的積分方程的微分方程，它是有局限性的.而用泰勒展開則可構(gòu)造任意多步法公式，其做法是根據(jù)多步法公式的形式，直接在處做泰勒展開.確定多步法（5.1）的系數(shù)及時(shí)不必套用系數(shù)公式（5.4），因?yàn)槎嗖椒ü讲灰欢ㄈ纾?.1）的形式，而且套用公式容易記錯(cuò).（5.1）第一百一十頁(yè)，共一百五十一頁(yè)，2022年，8月28日110

例7解初值問題

用顯式二步法其中試確定參數(shù)使方法階數(shù)盡可能高,并求局部截?cái)嗾`差.

解根據(jù)局部截?cái)嗾`差定義，用泰勒展開確定參數(shù)滿足的方程.第一百一十一頁(yè)，共一百五十一頁(yè)，2022年，8月28日111由于第一百一十二頁(yè)，共一百五十一頁(yè)，2022年，8月28日112為求參數(shù)使方法階數(shù)盡量高，可令即得方程組第一百一十三頁(yè)，共一百五十一頁(yè)，2022年，8月28日113解得，此時(shí)公式為三階，即為所求局部截?cái)嗾`差.而所得二步法為而且第一百一十四頁(yè)，共一百五十一頁(yè)，2022年，8月28日114

例8證明存在的一個(gè)值，使線性多步法

是四階的.

證明只要證明局部截?cái)嗾`差，則方法仍用泰勒展開，由于為四階.第一百一十五頁(yè)，共一百五十一頁(yè)，2022年，8月28日115第一百一十六頁(yè)，共一百五十一頁(yè)，2022年，8月28日116當(dāng)時(shí)，，故方法是四階的.第一百一十七頁(yè)，共一百五十一頁(yè)，2022年，8月28日1179.6線性多步法的收斂性與穩(wěn)定性

相容性及收斂性

線性多步法(5.1)式的相容性：在定義8中給出的局部截?cái)嗾`差(5.2)中，若稱步法(5.1)與微分方程(1.1)相容，它等價(jià)于對(duì)多步法(5.1)可引入多項(xiàng)式（6.1）（6.2）第一百一十八頁(yè)，共一百五十一頁(yè)，2022年，8月28日118和（6.3）分別稱為線性多步法(5.1)的第一特征多項(xiàng)式和第二特征多項(xiàng)式.可以看出，如果(5.1)式給定，則和也完全確定.反之也成立.第一百一十九頁(yè)，共一百五十一頁(yè)，2022年，8月28日119

定理5線性多步法(5.1)式與微分方程(1.1)相容的充分必要條件是關(guān)于多步法(5.1)的收斂性：由于用多步法求數(shù)值解需要個(gè)初值，而微分方程(1.1)只給出一個(gè)初值，因此還要給出個(gè)初值才能用多步法(5.1)進(jìn)行求解，即（6.4）（6.5）其中由微分方程初值給定，可由相應(yīng)單步法給出.第一百二十頁(yè)，共一百五十一頁(yè)，2022年，8月28日120設(shè)由(6.5)式在處得到的數(shù)值解為，這里為固定點(diǎn)，，于是有下面定義.

定義9設(shè)初值問題(1.1),(1.2)有精確解.如果初始條件滿足條件的線性步法(6.5)在處的解有則稱線性步法是收斂的.

定理6設(shè)線性多步法(6.5)是收斂的，則它是相容的.此定理的逆定理是不成立的.第一百二十一頁(yè)，共一百五十一頁(yè)，2022年，8月28日121

例9用線性二步法解初值問題

解此初值問題精確解，而由(6.6)式知故有，故方法(6.6)是相容的，但方法(6.6)的解并不收斂，在方法(6.6)中取初值

（6.6）（6.7）此時(shí)方法(6.6)為二階差分方程第一百二十二頁(yè)，共一百五十一頁(yè)，2022年，8月28日122其特征方程為解得其根為及.于是可求得(6.8)的解為故方法不收斂.多步法是否收斂與的根有關(guān).（6.8）第一百二十三頁(yè)，共一百五十一頁(yè)，2022年，8月28日123

定義10如果多步法（5.1）式的第一特征多項(xiàng)式的根都在單位圓內(nèi)或圓上，且在單位圓上的根為單根，則稱線性多步法（5.1）滿足根條件.

定理7線性多步法是相容的，則線性多步法（6.5）收斂的充分必要條件是線性多步法（5.1）滿足根條件.第一百二十四頁(yè)，共一百五十一頁(yè)，2022年，8月28日124

穩(wěn)定性與絕對(duì)穩(wěn)定性穩(wěn)定性主要研究初始條件擾動(dòng)與差分方程右端項(xiàng)擾動(dòng)對(duì)數(shù)值解的影響，假設(shè)多步法（6.5）有擾動(dòng)，則經(jīng)過擾動(dòng)后的解為，它滿足方程（6.9）第一百二十五頁(yè)，共一百五十一頁(yè)，2022年，8月28日125

定義11對(duì)初值問題(1.1)，(1.2)，由方法（6.5）得到的差分方程解，由于有擾動(dòng)，使得方程（6.9）的解為，若存在常數(shù)及，使對(duì)所有，當(dāng)時(shí)，有則稱多步法（5.1）是穩(wěn)定的或稱為零穩(wěn)定的.研究零穩(wěn)定性就是研究時(shí)差分方程（6.5）解的穩(wěn)定性.它表明當(dāng)初始擾動(dòng)或右端擾動(dòng)不大時(shí)，解的誤差也不大.第一百二十六頁(yè)，共一百五十一頁(yè)，2022年，8月28日126對(duì)多步法（5.1），當(dāng)時(shí)對(duì)應(yīng)差分方程的特征方程為.故有

定理8線性多步法（5.1）是穩(wěn)定的充分必要條件是它滿足根條件.關(guān)于絕對(duì)穩(wěn)定性只要將多步法（5.1）用于解模型方程（4.8），得到線性差分方程（6.10）第一百二十七頁(yè)，共一百五十一頁(yè)，2022年，8月28日127利用線性多步法的第一、第二特征多項(xiàng)式，令此式稱為線性多步法的穩(wěn)定多項(xiàng)式，它是關(guān)于的次多項(xiàng)式.如果它的所有零點(diǎn)滿足，則(6.10)式的解當(dāng)時(shí)，有.（6.11）第一百二十八頁(yè)，共一百五十一頁(yè)，2022年，8月28日128

定義12對(duì)于給定的，如果穩(wěn)定多項(xiàng)式(6.11)的零點(diǎn)滿足，則稱線性多步法(5.1)關(guān)于此值是絕對(duì)穩(wěn)定的.若在的復(fù)平面的某個(gè)區(qū)域中所有值線性多步法(5.1)都是絕對(duì)穩(wěn)定的，而在區(qū)域外，方法是不穩(wěn)定的，則稱為多步法(5.1)的絕對(duì)穩(wěn)定域.與實(shí)軸的交集稱為線性多步法(5.1)的絕對(duì)穩(wěn)定區(qū)間.為實(shí)數(shù)時(shí)，可以只討論絕對(duì)穩(wěn)定區(qū)間.由于線性多步法的絕對(duì)穩(wěn)定域較為復(fù)雜，通常采用根軌跡法.第一百二十九頁(yè)，共一百五十一頁(yè)，2022年，8月28日129阿當(dāng)姆斯顯式方法和隱式方法的絕對(duì)穩(wěn)定域圖形分別為圖9-6和圖9-7，絕對(duì)穩(wěn)定區(qū)間見表9-9.圖9-6圖9-7第一百三十頁(yè)，共一百五十一頁(yè)，2022年，8月28日130

例10討論辛普森方法的穩(wěn)定性

解辛普森方法的第第二特征多項(xiàng)式為的根分別為-1及1，它滿足根條件，故方法是零穩(wěn)定的.但它的穩(wěn)定多項(xiàng)式為求絕對(duì)穩(wěn)定域的邊界軌跡.若，則可令，在平面域的邊界軌跡為第一百三十一頁(yè)，共一百五十一頁(yè)，2022年，8月28日131可看出在虛軸上，且對(duì)全部從而可知為虛軸上從到的線段，故辛普森公式的絕對(duì)穩(wěn)定域?yàn)榭占?即步長(zhǎng)，此方法都不是絕對(duì)穩(wěn)定的，故它不能用于求解.第一百三十二頁(yè)，共一百五十一頁(yè)，2022年，8月28日1329.7

一階方程組與剛性方程

一階方程組

前面研究了單個(gè)方程的數(shù)值解法，只要把和理解為向量，那么，所提供的各種計(jì)算公式即可應(yīng)用到一階方程組的情形.考察一階方程組的初值問題，初始條件給為第一百三十三頁(yè)，共一百五十一頁(yè)，2022年，8月28日133則上述方程組的初值問題可表示為（7.1）求解這一初值問題的四階龍格-庫(kù)塔公式為式中若采用向量的記號(hào)，記第一百三十四頁(yè)，共一百五十一頁(yè)，2022年，8月28日134考察兩個(gè)方程的特殊情形：第一百三十五頁(yè)，共一百五十一頁(yè)，2022年，8月28日135這時(shí)四階龍格-庫(kù)塔公式具有形式其中（7.2）（7.3）第一百三十六頁(yè)，共一百五十一頁(yè)，2022年，8月28日136這是一步法，利用節(jié)點(diǎn)上的值，由（7.3）式順序計(jì)算，然后代入（7.2）式即可求得節(jié)點(diǎn)上的.（7.3）第一百三十七頁(yè)，共一百五十一頁(yè)，2022年，8月28日137

化高階方程為一階方程組高階微分方程（或方程組）的初值問題，原則上總可以歸結(jié)為一階方程組來求解.例如，考察下列階微分方程（7.4）初始條件為（7.5）只要引進(jìn)新的變量第一百三十八頁(yè)，共一百五