《計(jì)算方法》課件第7章

第7章常微分方程的數(shù)值解法7.1引言7.2歐拉方法及改進(jìn)的歐拉方法7.3龍格－庫(kù)塔方法7.4線性多步法 7.1引言

函數(shù)是客觀事物的內(nèi)部聯(lián)系在數(shù)量方面的反映，利用函數(shù)關(guān)系式可以對(duì)客觀事物的規(guī)律性進(jìn)行研究，因此尋求函數(shù)關(guān)系式在實(shí)踐中具有重要意義。然而工程實(shí)踐和科學(xué)研究中，許多問題僅能夠提供包含有要找的函數(shù)及其導(dǎo)數(shù)的關(guān)系式，也就是微分方程。通過微分方程的求解可以得到函數(shù)關(guān)系式，這在高等數(shù)學(xué)中已經(jīng)給出了一些典型微分方程求解的基本方法。但實(shí)際中遇到的微分方程往往比較復(fù)雜，在很多情況下，不能給出解的解析表達(dá)式。有時(shí)候即使能用解析表達(dá)式表示，又因?yàn)橛?jì)算量太大而不實(shí)用。更重要的是，許多實(shí)際問題，并不需要方程解的表達(dá)式，而僅僅需要得到一些點(diǎn)上的近似值即可。因此，有必要學(xué)習(xí)和研究微分方程的數(shù)值解法。在微分方程中，如果未知函數(shù)是一元函數(shù)，則微分方程稱做常微分方程，未知函數(shù)是多元函數(shù)的稱為偏微分方程。而微分方程中，出現(xiàn)的未知函數(shù)最高階導(dǎo)數(shù)的階數(shù)稱為微分方程的階。如果給定了微分方程，并給出了初始條件，則被稱為微分方程初值問題;如果給出的是邊界條件，則稱之為邊值問題。本章研究微分方程的數(shù)值解法，主要討論微分方程中最簡(jiǎn)單的一類問題，即一階常微分方程初值問題。一階常微分方程初值問題描述如下:(7.1)假設(shè)式(7.1)在區(qū)間［a，b］上存在唯一且足夠光滑的解y(x)。式(7.1)的數(shù)值求解就是尋求解y(x)在一系列離散點(diǎn)a≤x0<x1<…<xn≤b(也稱為節(jié)點(diǎn))處的近似函數(shù)值yi=f(xi)(i=0，1，2，…，n)。相鄰兩個(gè)節(jié)點(diǎn)間的距離——節(jié)點(diǎn)間距，也稱為步長(zhǎng)，用hi表示，即hi=xi+1－xi(i=0，1，2，…，n－1)。通常采用等距節(jié)點(diǎn)，即hi=h(常數(shù))，此時(shí)可以記xi=a+ih(i=0，1，2，…，n)。7.2歐拉方法及改進(jìn)的歐拉方法

7.2.1歐拉方法

1.歐拉公式

歐拉方法是解一階常微分方程初值問題的基本方法，可以利用數(shù)值微分的概念來推導(dǎo)歐拉公式。

關(guān)于一階微分方程初值問題，其目標(biāo)是計(jì)算函數(shù)在x1，x2，…，xn點(diǎn)上的函數(shù)值y1，y2，…，yn。為此，利用數(shù)值微分中的向前差商近似導(dǎo)數(shù)的方法可得由此可得根據(jù)式(7.1)可以得到y(tǒng)′(x0)=f(x0，y0)及y(x0)=y(a)=y0，代入式(7.2)得并將其記為(7.2)(7.3)式中,x0、y0、h及f(x，y)的表達(dá)式都是已知的，所以可以根據(jù)式(7.3)計(jì)算出節(jié)點(diǎn)x1處的函數(shù)近似值y1。用同樣的方法，可以得到節(jié)點(diǎn)x2處的函數(shù)近似值y2=y1+hf(x1，y1)，…，直到節(jié)點(diǎn)xn處的函數(shù)近似值yn=yn－1+hf(xn－1，yn－1)，則計(jì)算過程的通用表達(dá)式為式(7.4)即為求解一階常微分方程初值問題的歐拉公式。

例7－1取步長(zhǎng)h=0.02，用歐拉公式求一階微分方程初值問題：解根據(jù)歐拉公式及已知條件可得由于h=0.02，區(qū)間為［0，0.1］，所以節(jié)點(diǎn)為x0=0，x1=0.02，x2=0.04，x3=0.06，x4=0.08，x5=0.1，計(jì)算過程略，計(jì)算結(jié)果如表7.1所示。表7.1例7－1的計(jì)算結(jié)果

2.歐拉公式的幾何意義

設(shè)一階常微分方程的解函數(shù)y(x)為如圖7.1中所示的曲線，根據(jù)一階常微分方程初值問題的已知條件，曲線y(x)上的初始點(diǎn)(x0，y0)已知，故過該初始點(diǎn)以f(x0，y0)為斜率,可以得到一條直線，且該直線的方程為

y=y0+f(x0，y0)(x－x0)

該直線與直線x=x1的交點(diǎn)設(shè)為P1，則P1的縱坐標(biāo)為記為y1，即此與利用歐拉公式得到的值相同。然后過點(diǎn)(x1，y1)以f(x1，y1)為斜率，可以得到一條直線，且該直線方程為同樣可以得到該直線與直線x=x2的交點(diǎn)P2的縱坐標(biāo)，記為y2，即此與利用歐拉公式得到的值也相同。依此類推，得到圖7.1歐拉公式的幾何意義

3.歐拉方法的局部截?cái)嗾`差

定義7.1假設(shè)yi=y(xi)，即第i步計(jì)算是精確的前提下，考慮的截?cái)嗾`差Ri+1=y(xi+1)－yi+1，稱為局部截?cái)嗾`差，可用圖7.2表示，圖中|APi+1|即為歐拉方法在xi+1點(diǎn)的局部截?cái)嗾`差。圖7.2歐拉方法的局部截?cái)嗾`差示意圖

定義7.2如果數(shù)值方法得到的局部截?cái)嗾`差為O(hp+1)，則稱該方法具有p階精度或稱為p階的。

根據(jù)局部截?cái)嗾`差的定義，可以計(jì)算歐拉公式的局部截?cái)嗾`差。將y(xi+1)在xi處泰勒展開可得(7.5)由歐拉公式可得(7.6)根據(jù)局部截?cái)嗾`差的定義有則由式(7.5)和式(7.6)可得(7.7)式(7.7)即為歐拉方法的局部截?cái)嗾`差，故其具有一階精度。7.2.2改進(jìn)的歐拉方法

從歐拉方法的局部截?cái)嗾`差和歐拉方法的幾何意義都可以看出，歐拉方法的誤差較大，為此研究其改進(jìn)方法，以獲得精度較好的計(jì)算方法。

一階微分方程初值問題的微分方程為對(duì)其在小區(qū)間［xi，xi+1］(i=0，1，…，n－1)上進(jìn)行積分，可得即(7.8)根據(jù)問題的初始條件，如果能夠計(jì)算出式(7.8)右端的積分值，則可以獲得微分方程的數(shù)值解,故求解一階微分方程初值問題變成了求解函數(shù)積分的問題。

1.應(yīng)用矩形公式

對(duì)式(7.8)右端的積分應(yīng)用左矩形公式可得將其代入式(7.8)可得將y(xi)記做yi，y(xi+1)記做yi+1，則有可以看出式(7.9)與歐拉公式完全相同，這說明用數(shù)值積分的方法也可以推導(dǎo)得到歐拉公式。

2.應(yīng)用梯形公式

對(duì)式(7.8)右端的積分應(yīng)用梯形公式可得將其代入式(7.8)可得將y(xi)記做yi，y(xi+1)記做yi+1，則有(7.10)式(7.10)是應(yīng)用了數(shù)值積分中的梯形公式推導(dǎo)得到的，因此也被稱為梯形公式。根據(jù)局部截?cái)嗾`差的定義可以得到梯形公式的局部截?cái)嗾`差為(7.11)同樣，將y(xi+1)在xi處泰勒展開可得(7.12)其中，將式(7.12)代入式(7.11)有(7.13)將y′(xi+1)在xi處泰勒展開可得(7.14)其中,xi<τi<xi+1。將式(7.14)代入式(7.13)得(7.15)

3.改進(jìn)的歐拉公式

歐拉公式計(jì)算簡(jiǎn)單，但誤差較大。梯形公式雖可以提高精度，但為隱式公式，缺點(diǎn)也比較明顯(計(jì)算量增大了)。為此，將歐拉公式和梯形公式聯(lián)合使用，克服梯形公式的缺點(diǎn)，由此得到了改進(jìn)的歐拉公式。改進(jìn)的歐拉公式就是先用歐拉公式得出一個(gè)y(xi+1)的近似值，用yi+1來表示，稱為預(yù)估值。然后對(duì)預(yù)估值yi+1使用梯形公式進(jìn)行調(diào)整，得到較為精確的近似值yi+1，稱之為校正值。整個(gè)計(jì)算公式為(7.16)式(7.16)中第一式為預(yù)估算式，第二式為校正算式，整個(gè)公式稱為改進(jìn)的歐拉公式。為了便于編寫程序，常將式(7.16)改寫為(7.17)

例7－2取步長(zhǎng)h=0.1，用歐拉方法和改進(jìn)的歐拉方法求如下一階微分方程初值問題：

解由于h=0.1，區(qū)間為［0，1.5］，x0=0，所以需計(jì)算函數(shù)值的節(jié)點(diǎn)為xk=0+0.1k(k=1，2，…，15)。

(1)用歐拉方法求解，其計(jì)算公式為

(2)用改進(jìn)的歐拉方法求解，其計(jì)算公式為將已知數(shù)據(jù)代入上面的計(jì)算公式，可以得到微分方程的數(shù)值解，計(jì)算過程略，計(jì)算結(jié)果見表7.2。 由于該微分方程的解析解為 ，為了比較算法計(jì)算結(jié)果的誤差，將真實(shí)值y(xi)(i=0，1，2，…，15)也列入表7.2中。從表中數(shù)據(jù)可以看出，改進(jìn)的歐拉方法的計(jì)算結(jié)果更接近真實(shí)值。表7.2例7－2的計(jì)算結(jié)果

7.3龍格－庫(kù)塔方法

7.3.1龍格－庫(kù)塔方法的基本思想

對(duì)微分方程y′=f(x，y)，在小區(qū)間［xi，xi+1］上應(yīng)用拉格朗日中值定理，可得(7.18)令ξ=xi+θh，θ∈(0，1)，則式(7.18)變?yōu)?7.19)根據(jù)式(7.19)，如果能夠計(jì)算出f(xi+θh，y(xi+θh))，則可以根據(jù)一階微分方程初值問題的條件遞推得到問題的數(shù)值解，而f(xi+θh，y(xi+θh))在幾何意義上表示曲線上點(diǎn)的切線斜率，因此令K*=f(xi+θh，y(xi+θh))(稱其為區(qū)間［xi，xi+1］上的平均斜率)，則式(7.19)變?yōu)?7.20)因此只要對(duì)K*提供一種算法，就可以求得微分方程的數(shù)值解。下面根據(jù)該觀點(diǎn)對(duì)歐拉方法及改進(jìn)的歐拉方法進(jìn)行分析。(7.21)將y(xi)記做yi，y(xi+1)記做yi+1，則有(7.22)該式即歐拉公式，或者說歐拉公式以點(diǎn)(xi，yi)上的斜率來計(jì)算K*。當(dāng)取θ=0和θ=1時(shí)，可以得到兩個(gè)點(diǎn)(xi，yi)和(xi+1，y

i+1)上的斜率，并分別令其為K1和K2，則有(7.23)其中,點(diǎn)(xi+1，yi+1)處的斜率K2是由已知信息預(yù)測(cè)得到的。 如果取點(diǎn)(xi，yi)和點(diǎn)(xi+1，yi+1)上的斜率的算術(shù)平均值作為區(qū)間［xi，xi+1］上的平均斜率，即(7.24)則由式(7.20)、式(7.23)和式(7.24)，并將y(xi)和y(xi+1)分別記作yi和yi+1，得(7.25)

這顯然是改進(jìn)的歐拉公式，或者說改進(jìn)的歐拉公式以兩個(gè)點(diǎn)(xi，yi)和(xi+1，yi+1)上的斜率的算術(shù)平均值來計(jì)算K*。

已經(jīng)知道改進(jìn)的歐拉公式具有二階精度，而歐拉公式只有一階精度，結(jié)合上面的分析可知，用兩個(gè)點(diǎn)斜率的算術(shù)平均作為平均斜率比只取一個(gè)點(diǎn)的斜率作為平均斜率的計(jì)算精度要高。由此，如果設(shè)法在區(qū)間［xi，xi+1］內(nèi)多預(yù)測(cè)幾個(gè)點(diǎn)的斜率值，然后將它們加權(quán)平均作為平均斜率K*，則有可能構(gòu)造出具有更高精度的計(jì)算公式，這就是龍格－庫(kù)塔方法的基本思想。7.3.2二階龍格－庫(kù)塔方法

在區(qū)間［xi，xi+1］上，改進(jìn)的歐拉方法使用了xi、xi+1兩個(gè)節(jié)點(diǎn)上的斜率值?，F(xiàn)將其推廣到一般情況，取區(qū)間［xi，xi+1］上的點(diǎn)xi、xi+p，其中xi+p=xi+ph(0<p≤1)。設(shè)xi

和xi+p兩個(gè)點(diǎn)上的斜率分別為K1和K2，取其加權(quán)平均作為平均斜率K*的近似值，即(7.26)其中λ1和λ2為權(quán)值。顯然K1=f(xi，yi)，而K2=f(xi+p，yi+p)，對(duì)于其中的yi+p，應(yīng)用歐拉公式計(jì)算得(7.27)由式(7.20)、式(7.26)和式(7.27)可得(7.28)式(7.28)中，λ1、λ2、p為待定的未知數(shù)，如果能夠確定λ1、λ2、p的值，就可以計(jì)算微分方程的數(shù)值解。為了確定未知數(shù)的值，可以利用式(7.28)具有二階精度的條件，為此，首先計(jì)算式(7.28)的局部截?cái)嗾`差。根據(jù)泰勒公式有根據(jù)式(7.28)有(7.30)對(duì)進(jìn)行泰勒展開可得將其代入式(7.30)得(7.31)由式(7.29)和式(7.31)可得式(7.28)的局部截?cái)嗾`差為(7.32)式(7.28)具有二階精度，則其局部截?cái)嗾`差Ri+1=O(h3)，所以根據(jù)式(7.32)有(7.33)方程組(7.33)為三個(gè)未知數(shù)，兩個(gè)方程，因此方程組有無窮多個(gè)解，每一個(gè)解都可以得到一個(gè)具有二階精度的龍格－庫(kù)塔公式。特別地，當(dāng)λ1=λ2=1/2，p=1時(shí),式(7.28)可寫為(7.34)式(7.34)顯然是改進(jìn)的歐拉公式。此外，當(dāng)p=1/2時(shí)，可以得到λ1=0，λ2=1/2，將其代入式(7.28)有(7.35)式(7.35)稱為變形的歐拉公式。7.3.3高階龍格－庫(kù)塔方法

1.三階龍格－庫(kù)塔方法

根據(jù)龍格－庫(kù)塔方法的思想，為了進(jìn)一步提高微分方程數(shù)值求解的精度，可以在區(qū)間［xi，xi+1］上，除了點(diǎn)xi、xi+p外，再增加節(jié)點(diǎn)。如再增加一個(gè)點(diǎn)xi+m，xi+m=xi+mh，其中p<m≤1，則利用該三點(diǎn)處的斜率K1、K2、K3的加權(quán)平均值作為K*近似值，即(7.36)可得計(jì)算公式為(7.37)式中K1和K2的取法與二階龍格－庫(kù)塔方法相同，分別為而K3=f(xi+m，yi+m)需要得到點(diǎn)xi+m處的函數(shù)值y(xi+m)，為此用二階龍格－庫(kù)塔公式對(duì)yi+m進(jìn)行估計(jì),即(7.38)因此可得(7.39)由此得計(jì)算公式為(7.40)式(7.40)中，λ1、λ2、λ3、p、m、μ1、μ2為待定的未知數(shù)，如果能夠確定它們的值，就可以得到微分方程的數(shù)值解。為了確定未知數(shù)的值，與二階龍格－庫(kù)塔公式的推導(dǎo)過程類似，利用式(7.40)具有三階精度的條件，借助泰勒展開，計(jì)算其局部截?cái)嗾`差，從而得到待定未知數(shù)滿足的條件(方程)。具體推導(dǎo)過程略，得到的待定未知數(shù)滿足的條件為(7.41)顯然,由式(7.41)構(gòu)成的方程組有無窮多個(gè)解，每一個(gè)解都可以得到一個(gè)具有三階精度的龍格－庫(kù)塔公式，因此稱滿足條件(7.41)的計(jì)算公式(7.40)為三階龍格－庫(kù)塔公式。常用的三階龍格－庫(kù)塔公式為(7.42)

2.四階龍格－庫(kù)塔方法

與推導(dǎo)三階龍格－庫(kù)塔公式的方法類似,可以得到四階龍格－庫(kù)塔公式，只需在區(qū)間［xi，xi+1］上用四個(gè)點(diǎn)處的斜率的加權(quán)平均作為K*的近似值即可，計(jì)算公式的推導(dǎo)過程略。常用的經(jīng)典四階龍格－庫(kù)塔公式如下：(7.43)基于龍格－庫(kù)塔公式的思想，也可推導(dǎo)出其更高階的公式，龍格－庫(kù)塔公式的一般表達(dá)式為(7.44)其中，需要說明的是，龍格－庫(kù)塔公式的推導(dǎo)基于泰勒展開方法，所以它要求解函數(shù)有很好的光滑性，即y(x)要具有所要求的導(dǎo)數(shù)。不然，高階龍格－庫(kù)塔公式可能不如低階龍格－庫(kù)塔公式的效果好。通常來說，對(duì)于一般實(shí)際問題，四階龍格－庫(kù)塔方法即可達(dá)到一般的精度要求。

例7－3分別用三階和四階龍格－庫(kù)塔公式計(jì)算例7－2的初值問題。

解由于h=0.1，區(qū)間為［0，1.5］，x0=0，所以需計(jì)算函數(shù)值的節(jié)點(diǎn)為xi=0+0.1i(i=1，2，…，15)。

(1)采用三階龍格－庫(kù)塔公式計(jì)算。

由公式(7.42)可得(7.43)基于龍格－庫(kù)塔公式的思想，也可推導(dǎo)出其更高階的公式，龍格－庫(kù)塔公式的一般表達(dá)式為(7.44)其中，需要說明的是，龍格－庫(kù)塔公式的推導(dǎo)基于泰勒展開方法，所以它要求解函數(shù)有很好的光滑性，即y(x)要具有所要求的導(dǎo)數(shù)。不然，高階龍格－庫(kù)塔公式可能不如低階龍格－庫(kù)塔公式的效果好。通常來說，對(duì)于一般實(shí)際問題，四階龍格－庫(kù)塔方法即可達(dá)到一般的精度要求。

例7－3分別用三階和四階龍格－庫(kù)塔公式計(jì)算例7－2的初值問題。

解由于h=0.1，區(qū)間為［0，1.5］，x0=0，所以需計(jì)算函數(shù)值的節(jié)點(diǎn)為xi=0+0.1i(i=1，2，…，15)。

(1)采用三階龍格－庫(kù)塔公式計(jì)算。

由公式(7.42)可得

(2)采用四階龍格－庫(kù)塔公式計(jì)算。

由公式(7.43)可得其中,i=0，1，2，…，14。將x0=0，y0=1代入進(jìn)行計(jì)算，計(jì)算過程略，計(jì)算結(jié)果見表7.3。同樣，根據(jù)該微分方程的解析解，將真實(shí)值y(xi)(i=0，1，2，…，15)也列于表7.3中。

從表7.3可以看出，四階龍格－庫(kù)塔公式的計(jì)算結(jié)果更接近于真實(shí)值。表7.3例7－3的計(jì)算結(jié)果 7.4線性多步法

7.4.1線性多步法的基本思想

前面所討論的求解一階微分方程初值問題的歐拉方法及龍格－庫(kù)塔方法，在計(jì)算節(jié)點(diǎn)xi+1上的函數(shù)值yi+1時(shí)，僅利用了節(jié)點(diǎn)xi上的函數(shù)值yi，而此時(shí)已經(jīng)計(jì)算得到了節(jié)點(diǎn)xi+1前面的一系列點(diǎn)上的函數(shù)值的近似值y0，y1，…，yi，但這些近似值沒有被利用，這樣的方法被稱為單步法。相反，在計(jì)算節(jié)點(diǎn)xi+1上的函數(shù)值yi+1時(shí)，不僅利用了節(jié)點(diǎn)xi上的函數(shù)值yi，還利用了節(jié)點(diǎn)xi前面的一系列節(jié)點(diǎn)上的函數(shù)值的近似值yi－1，yi－2，…，這樣的計(jì)算方法稱為多步法。例如r步法，即計(jì)算節(jié)點(diǎn)xi+1上的函數(shù)值yi+1時(shí)，用到了yi，yi－1，yi－2，…,yi－r+1。

多步法計(jì)算中，由于使用了前面多步的信息來計(jì)算當(dāng)前節(jié)點(diǎn)的函數(shù)值，可以獲得比單步法更高的精度。最常用的多步法是線性多步法。r步線性多步法的公式可表示為(7.45)(7.46)在已知節(jié)點(diǎn)x0，x1，…，xn及y0，y1，…，yi的情況下，選取不同的插值節(jié)點(diǎn)做插值多項(xiàng)式就會(huì)得到不同的數(shù)值解法，常用的阿達(dá)姆斯方法就是基于該思想構(gòu)造的。7.4.2阿達(dá)姆斯顯式公式

針對(duì)式(7.46)在已計(jì)算出y0，y1，…，yi的情況下，選取r個(gè)節(jié)點(diǎn)xi，xi－1，…，xi－r+1進(jìn)行拉格朗日插值，可得(7.47)其中，將式(7.47)作為插值函數(shù)P(x)代入式(7.46)，并記yk=y(xk)，可得(7.48)令則式（7.48）可表示為：(7.49)式(7.49)中，當(dāng)r=4時(shí)可以得到更加具體的四步法公式：(7.50)其中,由于是等間距節(jié)點(diǎn)，設(shè)x=xi+th(x∈［xi，xi+1］，t∈［0，1］)，將其代入上式可得同理可得則式(7.50)的具體表達(dá)式為(7.51)式(7.51)即為常用的四步阿達(dá)姆斯顯式公式，它具有四階精度，其局部截?cái)嗾`差可以根據(jù)拉格朗日插值余項(xiàng)計(jì)算得到,即(7.52)7.4.3阿達(dá)姆斯隱式公式

針對(duì)式(7.46)在已計(jì)算出y0，y1，…，yi的情況下，選取r個(gè)節(jié)點(diǎn)xi+1，xi，…，xi－r+2進(jìn)行拉格朗日插值，可得(7.53)其中，將式(7.53)作為插值函數(shù)P(x)代入式(7.46)，并記yk=f(xk)，可得(7.54)令，則式(7.54)可表示為(7.55)很顯然,式