分析常微方程數(shù)值解法第八章常微分方程數(shù)值解法8-1第一頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-2第八章目錄§1歐拉（Euler）方法

1.1Euler法及其簡(jiǎn)單改進(jìn)1.2改進(jìn)的Euler法§2龍格庫(kù)塔（Runge-kutta）方法

2.1龍格-庫(kù)塔方法的基本思想2.2二階龍格-庫(kù)塔公式2.3高階R-K公式2.4變步長(zhǎng)R-K法§3線性多步法§4一階方程組與高階方程初值問題§5收斂性與穩(wěn)定性第二頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-3第八章序

許多科學(xué)技術(shù)問題，例如天文學(xué)中的星體運(yùn)動(dòng)，空間技術(shù)中的物體飛行，自動(dòng)控制中的系統(tǒng)分析，力學(xué)中的振動(dòng)，工程問題中的電路分析等，都可歸結(jié)為常微分方程的初值問題。

所謂初值問題，是函數(shù)及其必要的導(dǎo)數(shù)在積分的起始點(diǎn)為已知的一類問題，一般形式為：

我們先介紹簡(jiǎn)單的一階問題：第三頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-4第八章序由常微分方程的理論可知：上述問題的解唯一存在。

常微分方程求解求什么？應(yīng)求一滿足初值問題（8—1）的解函數(shù)y=y(x)，如對(duì)下列微分方程：《高等數(shù)學(xué)》中，微分方程求解，如對(duì)一階微分方程：y=f(x,y)是求解解函數(shù)y=y(x)，使?jié)M足上述方程。但能夠求出準(zhǔn)確的解析函數(shù)y(x)的微分方程是很少的，《高數(shù)》中研究微分方程的求解，是分門別類討論，對(duì)不同類型的微分方程，求解方法不一樣，因此，要求解微分方程，首先必須認(rèn)清類型。第四頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-5微分方程數(shù)值解而常微分方程初值問題的數(shù)值解法，是要尋求解函數(shù)y(x)在一系列點(diǎn)y(xi)（離散點(diǎn)）:上y(xi)的近似值yi（i=1,2,…,n），并且還可由這些（xi,yi)（i=1,2,…,n）構(gòu)造插值函數(shù)作為近似函數(shù)。上述離散點(diǎn)相鄰兩點(diǎn)間的距離hi=xi-1-xi稱為步長(zhǎng)，若hi都相等為一定數(shù)h,則稱為定步長(zhǎng)，否則為變步長(zhǎng)。由于在實(shí)際問題和科學(xué)研究中遇到的微分方程往往很復(fù)雜，絕大多數(shù)很難，甚至不可能求出解析函數(shù)y(x)，因此只能考慮求其數(shù)值解。本章重點(diǎn)討論如下一階微分方程：在此基礎(chǔ)上介紹一階微分方程組與高階微分方程的數(shù)值解法。第五頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-6§1歐拉（Euler）法

以Euler法及其改進(jìn)方法為例,說明常微分方程初值問題數(shù)值解法的一般概念，Euler法很簡(jiǎn)單，準(zhǔn)確度也不高，介紹此方法的目的，是由于對(duì)它的分析討論能夠比較清楚地顯示出方法的一些特點(diǎn)，而這些特點(diǎn)及基本方法反映了其它方法的特點(diǎn)。Euler法用于求解一階微分方程初值問題：第六頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-71.1Euler法及其簡(jiǎn)單改進(jìn)

Euler公式為：由x0出發(fā)x1,x2,…,xN，而利用此式可算出對(duì)應(yīng)的y1,y2,…,yN，式（8-2）稱為差分方程（序列{yn}滿足的方程）。下面是Euler公式的推導(dǎo)：一、從幾何意義出發(fā)：y=f(x,y)的解函數(shù)y=y(x)在xoy平面上是一族解曲線，而初值問題則是其中一條積分曲線，假定y=y(x)的曲線如圖8-1從給定的點(diǎn)P0(x0,y0)出發(fā)，以P0為切點(diǎn)，作切線，切線斜率為曲線y(x)的切線斜率

y

=f(x0,y0)，因此可得切線：（點(diǎn)斜式）P1P2y(x)P0

x2x1x0第七頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-8Euler公式的推導(dǎo)(續(xù)1）幾何意義：用折線近似解曲線y=y(x)，折線不會(huì)偏離太遠(yuǎn)，因?yàn)槊宽?xiàng)以f(x,y)（斜率）修正。切線與x=x1交于P1(x1,y1)，在[x0,x1]上以切線近似曲線，

第八頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-9Euler公式的推導(dǎo)（續(xù)2）二、利用Taylor級(jí)數(shù)：將y(x)在xn處展開：

第九頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-10Euler公式的推導(dǎo)（續(xù)3）公式（8-3）稱為后退Euler公式

所謂局部載斷誤差是指以yn代替y(xn)而得到y(tǒng)(xn+1)的近似值yn+1的誤差（只估計(jì)這一步的誤差）。三、利用數(shù)值微分公式：利用兩點(diǎn)公式并且Euler公式的局部截?cái)嗾`差為:后退Euler公式的局部截?cái)嗾`差為:第十頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-11Euler公式的推導(dǎo)（續(xù)4）第十一頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-12Euler公式的推導(dǎo)（續(xù)5）四、利用數(shù)值積分公式：在[xn,xn+1]上對(duì)y(x)=f(x,y(x))積分

對(duì)右端積分項(xiàng)采用不同的數(shù)值積分公式，便可得到各種不同的求解dE初值問題的計(jì)算公式。如，以矩形面積代替曲邊梯形面積1）以左矩形面積代替曲邊梯形面積如圖8-2，亦即以yf(x,y)xnxxn+1圖8-2

第十二頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-13yf(x,y)xnxxn+1圖8-3yf(x,y)xnxxn+1圖8-43）以梯形公式（面積）代替曲邊形如圖8-4則有

式（8-5）稱為求dE初值問題的梯形公式，梯形公式看作是以（xn,yn）(xn+1,yn+1)構(gòu)造的插值多項(xiàng)式代替被積函數(shù)得到的，而Euler公式則是以左端點(diǎn)函數(shù)值近似被積函數(shù)而得到，還可以用多個(gè)點(diǎn)做插值多項(xiàng)式近似被積函數(shù)構(gòu)造另一些精度更高的解微分方程的數(shù)值公式，梯形公式比Euler公式更準(zhǔn)確一些，誤差更小。Euler公式的推導(dǎo)（續(xù)6）2）以右矩形面積代替曲邊梯形：如圖8-3亦即以

第十三頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-14Euler公式注釋注1：Euler公式為顯式，右矩形，梯形公式為隱式；注2：Euler公式，梯形，右矩形公式為單步法，計(jì)算yn+1只用yn，而中點(diǎn)法公式為多步法（還可如上二所述，構(gòu)造多步法）即必須已知yn-1，yn才能計(jì)算yn+1，（求y0,y1不能用此公式。y0,y1稱為多步法的開始值，y0給定，而y1必須由其它公式算出，然后才能用中點(diǎn)法）；注3：前面已有Euler法的局部截?cái)嗾`差：后退Euler法的局部截?cái)嗾`差：

誤差階：如果局部截?cái)嗾`差則稱方法為P階的。第十四頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-15顯然，步長(zhǎng)h越小，階數(shù)P越高，局部截?cái)嗾`差越小，當(dāng)然計(jì)算精度越高；

注4：梯形法是幾階？梯形法精度比Euler法高，階數(shù)肯定比Euler法高，其實(shí)我們可以利用數(shù)值積分公式的誤差估計(jì)式，因?yàn)槲覀兪怯锰菪螖?shù)值積分公式計(jì)算因此由積分中梯形公式的誤差知此時(shí)的局部截?cái)嗾`差為：∴梯形法為2階方法！Euler法，后退Euler法為1階方法，而中點(diǎn)法為2階，Euler公式注釋（續(xù)）第十五頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-16關(guān)于Euler法的整體截?cái)嗾`差注釋注5：關(guān)于Euler法的整體截?cái)嗾`差：Euler方法的局部截?cái)嗾`差公式為：

實(shí)際計(jì)算時(shí)，yn是y(xn)的近似值，因此，計(jì)算過程中除每步所產(chǎn)生的局部截?cái)嗾`差外，還有因前面的計(jì)算不準(zhǔn)確而引起的誤差。在不考慮舍入誤差的情況下，稱y(xn+1)與yn+1之差為整體截?cái)嗾`差，記為：下面討論Euler方法的整體截?cái)嗾`差。

為簡(jiǎn)便起見，假定函數(shù)f(x,y)充分光滑，問題（8-1）解y(x)在[a,b]上二階連續(xù)可微，于是由式（8-6），局部截?cái)嗾`差有界，即存在M>0，使得對(duì)任意x[a,b]，都有|y′(x)|≤M，從而有：（緊接下屏）第十六頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-17式（8-7）表明，Euler方法的整體截?cái)嗾`差與h同階，當(dāng)h0時(shí)，en0。

關(guān)于Euler法的整體截?cái)嗾`差注釋（續(xù)）

反復(fù)遞推得：第十七頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-18結(jié)論

對(duì)于實(shí)際問題來說，由于L，M難以估計(jì)，因此（8-6）很難應(yīng)用，而且上述推導(dǎo)過程中一再放大了誤差上限，這樣的估計(jì)往往也很保守，遠(yuǎn)遠(yuǎn)大于實(shí)際的誤差，但是，從估計(jì)式（8-7）卻可以得到下面很有用處的結(jié)論。

1）當(dāng)h0時(shí)，en0即，亦即數(shù)值解yn，一致收斂于初值問題（8-1）的真解y(xn)，并且，Euler法的整體截?cái)嗾`差的階為O(h)與h同階，比局部截?cái)嗾`差低一階。2）舍入誤差局部截?cái)嗾`差對(duì)實(shí)際計(jì)算結(jié)果有影響，并且隨h減少而減少或增大。第十八頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-193）計(jì)算結(jié)果與解法的階數(shù)p，真解的導(dǎo)數(shù)y(p+1)有關(guān)，p越大，hp+1越小，|y(p+1)(ξ)|的上限越大，M也越大，因此為保證精度當(dāng)然應(yīng)選階數(shù)p較高的方法。但如果M很大，當(dāng)f(x,y)是分段連續(xù)的函數(shù)時(shí)，則應(yīng)采用低階的方法如用Euler法。結(jié)論（續(xù)）

4）計(jì)算結(jié)果還與開始值的精度有關(guān)，為使這種誤差的影響不致于超過局部截?cái)嗾`差，對(duì)多步法，應(yīng)采用跟多步法同階的方法計(jì)算開始值。第十九頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-201.2改進(jìn)的Euler法

梯形公式為二階方法，但卻是隱式格式，即若利用梯形公式求yn+1，就要求解方程（8-5）式，計(jì)算量較大，通常在實(shí)際計(jì)算時(shí)，將Euler法與梯形公式合起來使用，即先使用Euler公式,由（xn,yn）算出yn+1，記為yn+1(0)，稱為

預(yù)測(cè)值，然后用梯形公式去提高精度，將yn+1(0)

校正為較準(zhǔn)確的值：由于函數(shù)f(x,y)滿足Lipschitz條件，容易得出:第二十頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-21改進(jìn)的Euler法(續(xù)）第二十一頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-22預(yù)測(cè)──校正型公式

實(shí)際經(jīng)驗(yàn)表明，式（8-8）的迭代效果主要體現(xiàn)在第一次，由此構(gòu)成如下的預(yù)測(cè)──校正型公式:此式稱為改進(jìn)的Euler公式，為上機(jī)計(jì)算編程方便，常將式

（8-9）改寫為:下面分析改進(jìn)的Euler公式的局部截?cái)嗾`差：

第二十二頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-23改進(jìn)的Euler公式的局部截?cái)嗾`差分析假定yn=y(xn)，y(xn+1)的Taylor展式為:對(duì)于改進(jìn)的Euler公式，由于

這說明改進(jìn)的Euler法的局部截?cái)嗾`差為O(h3)，比Euler公式高一階，是二階方法。第二十三頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-24改進(jìn)的Euler公式舉例例1

這些結(jié)果在表8-1中，可見計(jì)算結(jié)果的精度，Euler法與后退Euler法差不多，與準(zhǔn)確值相比較Euler法偏小，而后退Euler法偏大；中點(diǎn)法與梯形法精度同為2階，但梯形法更好一些，這跟它們局部截?cái)嗾`差的符號(hào)，階數(shù)和系

數(shù)的大小是完全一致的。表見下屏：第二十四頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-25表格8-1表8-1

y=y,y(0)=1的數(shù)值解(h=0.1)

x精確解歐拉法后退歐拉中點(diǎn)法梯形法.1.904837.900000.909091.900000.904762.2.808731.810000.826446.820000.818594.3.740818.729000.751315.736000.740633.4.670320.656100.683013.627800.670096.5.060531.590490.620921.601440.606278.6.548812.531441.654474.552512.548537.7.496585.478298.5131458.490938.496295.8.449329.430467.466507.454324.449029.9.406570.387421.424098.400073.4062641.367879.348679.385543.374310.367573第二十五頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-26表格8-2

而表8-2是分別取了不同的h=0.1,h=0.01，h=0.001,h=0.0001，還是利用這些公式，經(jīng)過若干步的計(jì)算（h越小，計(jì)算量越大）算到y(tǒng)(1)的近似值，可見：隨著h的減小，y(1)的近似值的精度在提高，0.01比0.001差，即0.001比0.01時(shí)的y(1)準(zhǔn)確。Y′=-y,y(0)=1的解y(1)的近似值(y(1)=0.367879)h歐拉法后退歐拉法中點(diǎn)法梯形法0.1.348678.385543.374310.3675730.01.366033.369711.367944.3678770.001.367700.368052.367879.3678760.0001.367800.367800.367881.368020（緊接下屏）第二十六頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-27表8-2計(jì)算結(jié)果說明（續(xù)）

但h太小，到h=0.0001時(shí)卻又變得誤差大了，這與前面所說h越小，p階越高，應(yīng)該局部截?cái)嗾`差越小，因而計(jì)算精度更高矛盾了，為什么會(huì)產(chǎn)生這種情況呢？這是由于h太小而引起計(jì)算量大因而造成了舍入誤差和截?cái)嗾`差的積累，這種情況由于初值問題不同可能會(huì)影響更大，偏離更嚴(yán)重，如下面的例2。這種問題實(shí)際上是穩(wěn)定性問題，我們將會(huì)討論方法的穩(wěn)定性，由此得出對(duì)h有一定的要求的穩(wěn)定性制區(qū)域。第二十七頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-28例2求解初值問題y′=-20y，y(0)=1，計(jì)算y(1)的近似值。解：類似于例1，用歐拉法、后退歐拉法、中點(diǎn)法、梯形法求解，得到如下表8.3。

表8.3y=20,y(0)=1的解y(1)的近似值(y(1)=0.206116E8)h歐拉法后退歐拉法中點(diǎn)法梯形法.11.169351E-4.514229E+60.01.203704E-9.120746E-7.413244E+7.192743E-8.001.168300E-8.251090E-8.484136E+5.205979E-8.0001.202081E-8.210331E-8.475130+3.206103E-8

由表8.3可見，盡管中點(diǎn)法的階數(shù)與梯形法相同，比歐拉法和后退歐拉法的階數(shù)高，計(jì)算結(jié)果的精度卻很糟糕。此外，盡管歐拉法與后退歐拉法的階數(shù)相同，但歐拉法計(jì)算結(jié)果的精度，當(dāng)h=0.1時(shí)卻比后退歐拉法差。第二十八頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-29§2龍格-庫(kù)塔（Runge-kutta）方法2.1龍格-庫(kù)塔方法的基本思想：

因此只要對(duì)平均斜率k*提供一種算法，由（8-11）式便相應(yīng)地得到一種微分方程的數(shù)值計(jì)算公式。（緊接下屏）第二十九頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-30龍格-庫(kù)塔方法的基本思想(續(xù)）

改進(jìn)歐拉公式比歐拉公式精度高的原因，也就在于確定平均斜率時(shí)，多取了一個(gè)點(diǎn)的斜率值。因此它啟發(fā)我們，如果設(shè)法在[xi,xi+1]上多預(yù)報(bào)幾個(gè)點(diǎn)的斜率值，然后將它們加權(quán)平均作為k*的近似值，則有可能構(gòu)造出更高精度的計(jì)算公式，這是龍格-庫(kù)塔方法的基本思想。用這個(gè)觀點(diǎn)來研究歐拉公式與改進(jìn)歐拉公式，可以發(fā)現(xiàn)歐拉公式由于僅取xn一個(gè)點(diǎn)的斜率值f(xn,yn)作為平均斜率k*的近似值，因此精度很低。而改進(jìn)歐拉公式（8-10）卻是利用了xn與xn-1兩個(gè)點(diǎn)的斜率值k1=f(xn,yn)與k2=f(xn+1,yn+hk1)取算術(shù)平均作為平均斜率k*的近似值。其中k2是通過已知信息yn利用歐拉公式求得的。第三十頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-312.2二階龍格-庫(kù)塔公式

首先推廣改進(jìn)歐拉公式，考察區(qū)間[xn,xn+1]內(nèi)任一點(diǎn)：我們希望用xn和xn+1兩個(gè)點(diǎn)的斜率值k1和k2加權(quán)平均作為平均斜率k*的近似值：其中c1，c2為待定常數(shù)，同改進(jìn)歐拉公式一樣，這里仍取：?jiǎn)栴}在于怎樣預(yù)測(cè)xn+l處的斜率值k2。仿照改進(jìn)歐拉公式，先用歐拉公式提供y(xn+l)的預(yù)測(cè)值

然后再用預(yù)測(cè)值yn+l通過計(jì)算f產(chǎn)生斜率值k2=f(xn+l,yn+l)，這樣設(shè)計(jì)出的計(jì)算公式具有形式：第三十一頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-32二階龍格-庫(kù)塔公式（續(xù)1）公式（8-12）中含有三個(gè)待定參數(shù)c1,c2和l，我們希望適當(dāng)選取這些參數(shù)值，使得公式（8-12）具有二階精度，亦即使：現(xiàn)在仍假定yn=y(xn)，即yn是準(zhǔn)確的，將y(xn+1)與yn+1都在xi處作泰勒展開：第三十二頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-33二階龍格-庫(kù)塔公式（續(xù)2）代入（8-12）式，得：比較（8-13）與（8-14）兩式，要使公式具有二階精度，只有滿足下列條件：

這里一共有三個(gè)待定參數(shù)，但只需滿足兩個(gè)條件，因此有一個(gè)自由度，于是滿足條件（8-15）的參數(shù)不止一組，而是一族，相應(yīng)的公式（8-12）也有一族，這些公式統(tǒng)稱為二階龍格-庫(kù)塔公式，簡(jiǎn)稱二階R-K公式。第三十三頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-34

特別，當(dāng)l=1即xn+l=xn+1時(shí)，c1=c2=1/2，二階R-K公式就是改進(jìn)歐拉公式。

如果取l=1/2，則c1=0,c2=1，這時(shí)二階R-K公式稱為變形的歐拉公式，其形式見左邊：

從表面上看，變形的歐拉公式僅含一個(gè)斜率值k2，但k2是通過k1計(jì)算出來的，因此每完成一步，仍然需要兩次計(jì)算函數(shù)f的值，工作量和改進(jìn)歐拉公式相同。二階龍格-庫(kù)塔公式（續(xù)3）第三十四頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-35構(gòu)造二階R-K公式的主要步驟

綜上所述，構(gòu)造二階R-K公式主要由以下幾步產(chǎn)生：

在區(qū)間[xn,xn+1]上取二點(diǎn)，預(yù)報(bào)相應(yīng)點(diǎn)

的斜率值；2)對(duì)此兩斜率值加權(quán)平均作為平均斜率

值的近似值；3)將yn+1與y(xn+1)都在xi處作泰勒展開，

為使公式達(dá)到二階精度，比較相應(yīng)系

數(shù)，建立有關(guān)參數(shù)所應(yīng)滿足的方程組；4)解此方程組得一族二階R-K公式。第三十五頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-362.3高階R-K公式

為了進(jìn)一步提高精度，在[xn,xn+1]上除xn和xn+l外再增加一點(diǎn)xn+m=xn+mh(l

m1)，并用xn,xn+l,xn+m三處的斜率值k1,k2,k3加權(quán)平均作為k*的近似值，這時(shí)計(jì)算公式為：其中k1,k2仍用（8-12）式所取的形式。為了預(yù)測(cè)xn+m處的斜率值k3，要定出xn+m處所對(duì)應(yīng)的yn+m，可以看作在區(qū)間[xn,xn+m]上使用二階R-K公式，從而得到y(tǒng)(xn+m)的預(yù)測(cè)值：于是，再通過計(jì)算函數(shù)值f得到:第三十六頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-37高階R-K公式（續(xù)1）這樣設(shè)計(jì)出的計(jì)算公式具有形式:運(yùn)用泰勒展開方法，適當(dāng)選擇參數(shù)c1,c2,c3,l,m,及b1,b2使上述公式具有三階精度采用與(8-12)類似的處理方法，得到這些參數(shù)需要滿足的條件：(見下屏)第三十七頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-38高階R-K公式（續(xù)2）滿足5個(gè)條件的一族公式（8-18），統(tǒng)稱為三階R-K公式，其中常見的是庫(kù)塔公式：c1=1/6,c2=4/6,c3=1/6,l=1/2,m=1,b1=-1,b2=2第三十八頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-39

若需再將精度提高至四階，用上述類似處理方法，只是必須在[xi,xi+1]上用四個(gè)點(diǎn)處的斜率加權(quán)平均作為k*的近似值，構(gòu)成一族四階龍格-庫(kù)塔公式，由于推導(dǎo)復(fù)雜，這里從略，只將常用的公式介紹如下：四階公式中常用的還有下面的Gill公式：高階R-K公式（續(xù)3）一般稱它為標(biāo)準(zhǔn)四階R-K公式，其局部截?cái)嗾`差是O(h5)。第三十九頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-40Gill公式

在計(jì)算機(jī)上它具有節(jié)省存貯單元和控制舍入誤差增長(zhǎng)的優(yōu)點(diǎn)。第四十頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-41初值問題兩種方法舉例用改進(jìn)的Euler

法和四階標(biāo)準(zhǔn)

R-K公式解初值

問題：例3[解]用改進(jìn)的Euler公式，取h=0.1計(jì)算公式為：部分計(jì)算結(jié)果見表8-4：

表8-4xn改進(jìn)的歐拉法準(zhǔn)確解y(xn)誤差xn改進(jìn)的歐拉法準(zhǔn)確解y(xn)誤差01.100.61.4859561.4832400.0027160.21.1840961.1832160.0000880.81.6164761.6124520.0040240.41.3433601.3416410.0017191.01.7378691.7320510.005818第四十一頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-42對(duì)四階標(biāo)準(zhǔn)R-K法，取h=0.2，計(jì)算公式如下：例3（續(xù)1）第四十二頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-43例3（續(xù)2）表8-5

xn四階R-K法誤差xn四階R-K法誤差0110.61.4832810.0000410.21.1832290.0000130.81.6125130.0000610.41.3416670.000261.01.7321400.000089

計(jì)算結(jié)果見表8-5，與準(zhǔn)確解比較，至少有四位有效數(shù)字。而取步長(zhǎng)h=0.1時(shí)，用改進(jìn)的歐拉公式計(jì)算，最多只有三位有效數(shù)字（見表8-4），雖然改進(jìn)的歐拉公式每前進(jìn)一步只要計(jì)算兩次f值，而四階R-K法要計(jì)算4次f值，但由于后者步長(zhǎng)比前者大，所耗費(fèi)的計(jì)算總量還是差不多的。這個(gè)例子又一次顯示了選擇算法的重要意義。第四十三頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-44表8-6

計(jì)算f的次數(shù)與精度階數(shù)的關(guān)系

每步計(jì)算f的次數(shù)23456789精度的階數(shù)23445667表8-6的說明

一般，高精度的方法要求解有較好的光滑性，解的光滑性是由f(x,y)的可微性決定。如果解的光滑性差，則用四階R-K法求得的數(shù)值解反而不如用改進(jìn)的歐拉法好。因此，一定要針對(duì)具體問題的特點(diǎn)來選擇求解方法。例3（續(xù)3）第四十四頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-45表8-6的說明（續(xù)）

龍格-庫(kù)塔方法的推導(dǎo)是在用泰勒展開方法的基礎(chǔ)上進(jìn)行的，因而它要求所求的解具有較好的光滑性質(zhì)。假若解的光滑性差，那么使用四階R-K公式求得的數(shù)值解，其精度可能反而不如改進(jìn)歐拉公式。在實(shí)際計(jì)算時(shí)，應(yīng)當(dāng)針對(duì)問題的具體特點(diǎn)，選擇合適的算法。從理論上講，可以構(gòu)造任意高階的龍格-庫(kù)塔公式，但只要注意到精度的階數(shù)與計(jì)算函數(shù)值f(x,y)的次數(shù)之間的關(guān)系不是等量增加的（見表8-6），就可得出再提高公式階數(shù)已沒有多大意思了。因此更說明四階R-K公式是兼顧了精度及計(jì)算量的較理想的計(jì)算公式。

需要指出的是：第四十五頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-462.4變步長(zhǎng)R-K法

當(dāng)用數(shù)值方法解微分方程時(shí)，單從一步來看，步長(zhǎng)越小，局部截?cái)嗾`差就越小。但從整體來看，步長(zhǎng)越小，在求解范圍內(nèi)所要完成的步數(shù)就會(huì)越多。這一方面增加了計(jì)算量；另一方面也增大了舍入誤差的累積，因此有個(gè)如何合理選擇步長(zhǎng)的問題。在解決這個(gè)問題時(shí)，需要考慮兩個(gè)問題：

①怎樣衡量和檢驗(yàn)計(jì)算結(jié)果的精度？②如何依據(jù)獲得的精度信息及時(shí)處理步長(zhǎng)。

以四階標(biāo)準(zhǔn)R-K公式為例，從節(jié)點(diǎn)xn出發(fā)，先以步長(zhǎng)h求得xn+1處解的近似值，記作yn+1[h]，由于公式的局部截?cái)嗾`差為O(h5)，故有：第四十六頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-47變步長(zhǎng)R-K法（續(xù)1）

式中c與y(5)(x)在區(qū)間[xn,xn+1]內(nèi)的值有關(guān)，假設(shè)y(5)(x)在區(qū)間[xn,xn+1]內(nèi)變化不大，則可將c視作常數(shù)。將步長(zhǎng)折半，即以h/2為步長(zhǎng)，跨兩步到xn+1，再求得一個(gè)近似值每跨一步的局部截?cái)嗾`差為c(h/2)5，因此

第四十七頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-48具體做法

1.若<

且相差不多，則取，并以h為步長(zhǎng)計(jì)算xn+2對(duì)應(yīng)的yn+2。3.以h作基本步長(zhǎng)，當(dāng)>

時(shí)，反復(fù)將步長(zhǎng)折半，每折半一

次，由xn+1的步數(shù)增加一倍，直至2.如果<<

，又不需要節(jié)點(diǎn)xn+h處的數(shù)值解值，則可將步

長(zhǎng)加倍，只要保證所求節(jié)點(diǎn)處的y值符合精度要求即可。第四十八頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-49變步長(zhǎng)四階標(biāo)準(zhǔn)R-K法解初值問題舉例

例4解取基本步長(zhǎng)h=0.1，節(jié)點(diǎn)xk=kh，=106，按變步長(zhǎng)四階標(biāo)準(zhǔn)R-K法進(jìn)行計(jì)算，

其結(jié)果見表8-7，顯然，盡管每一步的步長(zhǎng)均只折半一次，但解的精度可達(dá)106，與準(zhǔn)確解：比較，說明以上分析是正確的。表8-7見下屏第四十九頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-50步長(zhǎng)折半次數(shù)xkyk準(zhǔn)確解y(xk)10.10.1049920.1049917110.20.2198680.2198680710.30.3443410.3443407710.40.4779530.4779534810.50.6201150.620114510.60.7701350.7701352810.70.9272700.9272702410.81.0907541.090753610.91.2598311.259830411.01.4337811.4337808表8-7第五十頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-51§3線性多步法

前兩節(jié)介紹的幾種差分格式有一共同的

特點(diǎn)，在計(jì)算yn+1時(shí)僅用前一步求得的yn，而

對(duì)前幾步的信息未予啟用，因而稱之為單步

法?，F(xiàn)討論多啟用一些已知信息的線性多步

法。線性多步法公式的一般形式為：

其中i、i均為與f,n無關(guān)的常數(shù)，r、r不同時(shí)為0。由于求yn+1時(shí)要到前r+1步的信息yn-r，…，yn及對(duì)應(yīng)的f值，因稱它為r+1步法。-1=0，式（8-21）為顯式，-10時(shí)，式（8-21）為隱式。第五十一頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-523.1阿當(dāng)姆斯Adams公式

當(dāng)式（8-21）中系數(shù)0=1,

1=…=r=0時(shí)，稱之為阿當(dāng)姆斯（Adams）公式，這類公式很容易借助數(shù)值積分導(dǎo)出。設(shè)y(x)是初值問題（8-1）的準(zhǔn)確解，則:兩邊從xn到xn+1積分，得:對(duì)此式右端應(yīng)用數(shù)值積分公式，則可導(dǎo)出不同的阿當(dāng)姆斯公式。例如，記fj=f(xj,y(xj))，取數(shù)據(jù)點(diǎn)(xn-1,fn-1)，(xn,fn)構(gòu)造f的線性插值多項(xiàng)式:其中：第五十二頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-53阿當(dāng)姆斯Adams公式（續(xù)）第五十三頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-54r+1步阿當(dāng)姆斯顯式公式

（又稱Adams-Bashforth公式）

類似地，取數(shù)據(jù)點(diǎn)（xn-r,fn-r）、(xn-r+1,fn-r+1)、…、(xn,fn)，構(gòu)造f(x,y(x))的r次插值多項(xiàng)式，可導(dǎo)出r+1步阿當(dāng)姆斯顯式公式（又稱Adams-Bashforth公式）及其局部截?cái)嗾`差:rArBr0Br1Br2Br3Br4r+1011

1/21231

5/1221223165

3/83245559379

251/7204720190127742616127425195/288其系數(shù)見如下表8-8，顯然具有r+1階精度。第五十四頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-55r+1步阿當(dāng)姆斯隱式公式

（又稱Adams-Moulton公式）

若改取數(shù)據(jù)點(diǎn)（xn-r+1,fn-r+1）、(xn-r+2,fnr+2)、…、(xn+1,fn+1)，構(gòu)造f(x,y(x))的r次插值多項(xiàng)式，則可以導(dǎo)出n+1步阿當(dāng)姆斯隱式公式（又稱Adams-Moulton公式）及其局部截?cái)嗾`差：第五十五頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-56

其系數(shù)見下表8-9，顯然也具有r+1階精度。

rArBr-1*Br0*Br1*Br2*Br3*r+1*011

1/21211

1/12212581

1/2432491951

19/7204720251646264106193/160可見，一階阿當(dāng)姆斯顯式公式（對(duì)應(yīng)r=0）即為歐拉公式；一階與二階阿當(dāng)姆斯隱式公式對(duì)應(yīng)r=0與r=1，分別為后退歐拉公式與梯形公式。r+1步阿當(dāng)姆斯隱式公式

（又稱Adams-Moulton公式）（續(xù)）

第五十六頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-57四階阿當(dāng)姆斯顯式公式

高階阿當(dāng)姆斯公式中最常用的是四階公式，對(duì)應(yīng)r=3。若選xn-3,xn-2，xn-1,xn作為插值節(jié)點(diǎn)，則四階阿當(dāng)姆斯顯式公式（又稱為外推公式）及其局部截?cái)嗾`差為：

此公式的優(yōu)點(diǎn)是精度高，其局部截?cái)嗾`差與四階R-K法同階，但每前進(jìn)一步只要計(jì)算一次函數(shù)值f(xn,yn)，其余的值f在前面的計(jì)算中已求出，計(jì)算量遠(yuǎn)小于四階R-K法。不足之處是不能從n=1開始使用，必須用其它方法（例如四階R-K法）求得出發(fā)值y1、y2、y3。此外，若要中途改變步長(zhǎng)也較麻煩。第五十七頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-58四階阿當(dāng)姆隱式公式內(nèi)插公式及其局部截?cái)嗾`差為（取xn-2,xn-1,xn,xn+1為插值節(jié)點(diǎn)）：

它也具有公式（8-26）的優(yōu)缺點(diǎn)，只是對(duì)yn+1必須提供一個(gè)預(yù)測(cè)值，然后通過迭代才能求得yn+1。將公式（8-26）與公式（8-27）聯(lián)合使用，可組成如下預(yù)測(cè)一校正系統(tǒng)：

四階阿當(dāng)姆斯顯式公式（續(xù)）第五十八頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-59阿當(dāng)姆斯公式舉例例5分別用四階阿當(dāng)姆斯顯式公式（8-26）與阿當(dāng)姆斯預(yù)測(cè)一校正公式（8-28）解例3中的初值問題。解

取步長(zhǎng)h=0.1，由于公式（8-9）、（8-28）均為四步法，因此必須用其它方法求出發(fā)值。現(xiàn)采用四階標(biāo)準(zhǔn)龍格-庫(kù)塔法求y1,y2,y3‘然后再分別用指定的方法進(jìn)行計(jì)算，計(jì)算結(jié)果見表8-10。表8-10xnyn四階標(biāo)準(zhǔn)R-K法顯式公式（8-25）預(yù)測(cè)-校正公式（8-28）準(zhǔn)確解0.11.095446

1.0954451150.21.183217

1.1832159560.31.264912

1.2649110460.4

1.3415521.3416411.3416407860.5

1.4140461.4142131.4142135620.6

1.4830191.4832391.4832396870.7

1.5489191.5491921.5481833380.8

1.6121161.6124501.6124515490.9

1.6729171.6733781.633200531.0

1.7315691.7320481.732050807第五十九頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-60阿當(dāng)姆斯預(yù)測(cè)—校正公式

為了進(jìn)一步提高上述預(yù)測(cè)一校正法的計(jì)算精度，可以利用預(yù)測(cè)公式與校正公式同階的特點(diǎn)，導(dǎo)出局部截?cái)嗾`差估計(jì)公式，對(duì)預(yù)測(cè)值與校正值分別用各自的局部截?cái)嗾`差估計(jì)值進(jìn)行修正，以此來提高解的精度。由式（8-26）、（8-27）可知，公式（8-28）中預(yù)測(cè)公式與校正公式的局部截?cái)嗾`差分別為：若記預(yù)測(cè)值為pn+1，校正值為cn+1，則：緊接下屏第六十頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-61于是得pn+1與cn+1的事后誤差估計(jì)：

阿當(dāng)姆斯預(yù)測(cè)—校正公式（續(xù)）分別代替pn+1與cn+1能提高解的精度。為簡(jiǎn)單計(jì)，在未求出cn+1以前，用pncn代替pn+1cn+1。第六十一頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-62阿當(dāng)姆斯帶誤差修正的預(yù)測(cè)—校正公式據(jù)此，得到下述帶誤差修正的預(yù)測(cè)一校正公式：式中初始值p3與c3可取為0。

第六十二頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-633.2哈明(Hamming)公式

一般的線性多步公式（8-21）還可用待系數(shù)法并借助泰勒展開導(dǎo)出。例如考慮確定如下形式的四階公式：其中yn表示fn。只需求系數(shù)j、j，使其局部截?cái)嗾`差T=O(h5)即可。局部截?cái)嗾`差為：將其在xn點(diǎn)進(jìn)行泰勒展開，得：

（緊接下屏）第六十三頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-64哈明(Hamming)公式（續(xù)1）令hk(k=0,1,…4)的系數(shù)為0,

得關(guān)于j、j的方程組:方程組有六個(gè)未知數(shù)，而只有五個(gè)方程，因此有無窮多解，取1為自由參量，則有:第六十四頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-65哈明(Hamming)公式（續(xù)2）當(dāng)1=0時(shí)，得公式：

這是一個(gè)四階隱式公式。完全類似，可導(dǎo)出如下形式的四階顯式公式：其中最常用的是：第六十五頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-66Hamming公式（預(yù)測(cè)一校正公式）

用顯式公式（8-34）與隱式公式（8-33）聯(lián)合可組成預(yù)測(cè)一校正公式

此式稱為哈明（Hamming）公式，大量數(shù)值實(shí)驗(yàn)表明，此公式的數(shù)值穩(wěn)定性在同類公式中幾乎是最好的。為了進(jìn)一步改善哈明公式計(jì)算結(jié)果的精度而又不致增加過多的計(jì)算量，可仿照帶誤差修正的阿當(dāng)姆斯預(yù)測(cè)一校正公式，建立一個(gè)帶誤差修正的哈明公式：（緊接下屏）第六十六頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-67

前面介紹的四套預(yù)測(cè)一校正公式（8-28）、（8-29）、（8-32）、（8-33）都是實(shí)際計(jì)算中常用的方法，它們的共同特點(diǎn)是計(jì)算量省而精度高，但都必須用別的單步法為其提供出發(fā)值，由于解題過程是遞推進(jìn)行的，所以出發(fā)值的精度對(duì)以后的計(jì)算影響很大，故一般采用四階龍格一庫(kù)塔法求其出發(fā)值。Hamming公式（預(yù)測(cè)一校正公式）（續(xù)）

第六十七頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-68§4一階方程組與高階方程初值問題

4.1一階方程組

下面以兩個(gè)方程的情形為例，介紹一階方程組的數(shù)值解法。設(shè)初值問題：若采用向量記號(hào)，記：則初值問題（8-34）可表示為：

它與初值問題（8-1）形式上完全一致，前面介紹的解初值問題（8-1）的各種數(shù)值方法均可用于解一階方程組初值問題（8-35），只要將公式中的y與f換成向量y與f即可。第六十八頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-69一階方程組（續(xù)1）例如：解初值問題（8-35）的歐拉公式為：

即：寫成分量形式即為：又如解初值問題（8-35）的四階標(biāo)準(zhǔn)龍格-庫(kù)塔公式為：第六十九頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-70一階方程組（續(xù)2）寫成分量形式為：

第七十頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-714.2高階方程初值問題

高階方程初值問題常轉(zhuǎn)化為一階方程組求解。例如：

引進(jìn)新變量z=y，則此問題等價(jià)于

對(duì)其使用四階標(biāo)龍格-庫(kù)塔法，計(jì)算公式為

第七十一頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-72高階方程初值問題舉例例5用四階龍格-庫(kù)塔方法在[0，1]上取步長(zhǎng)h=0.2，求解二階方程初值問題：解先將二階方程化為一階方程組。令z=y，則得方程組：使用四階龍格-庫(kù)塔公式，其相應(yīng)的形式為：其中

第七十二頁(yè)，共八十一頁(yè)，編輯于2023年，星期日第八章常微分方程數(shù)值解法8-73例5（續(xù)）表8-11

y-3y+2y=0解的近擬值xiyizik1l1k2l2k3l3k4l401.00001.00001.00001.00001.10001.10001.11001.11001.22201.22200.21.22141.22141.22141.22141.34351.34351.35581.35581.49261.49260.41.19181.19181.19181.19181.64101.64101.65591.65591.82301.82300.61.