常微分方程初值問題數(shù)值解法-課件

第9章常微分方程初值問題數(shù)值解法Euler公式改進的歐拉公式龍格—庫塔法亞當(dāng)斯法算法的穩(wěn)定性及收斂性2023/7/31§9.1引言包含自變量、未知函數(shù)及未知函數(shù)的導(dǎo)數(shù)或微分的方程稱為微分方程。在微分方程中,自變量的個數(shù)只有一個,稱為常微分方程。自變量的個數(shù)為兩個或兩個以上的微分方程叫偏微分方程。微分方程中出現(xiàn)的未知函數(shù)最高階導(dǎo)數(shù)的階數(shù)稱為微分方程的階數(shù)。如果未知函數(shù)y及其各階導(dǎo)數(shù)都是一次的,則稱它是線性的,否則稱為非線性的。2023/7/32

在高等數(shù)學(xué)中，對于常微分方程的求解，給出了一些典型方程求解析解的基本方法，如可分離變量法、常系數(shù)齊次線性方程的解法、常系數(shù)非齊次線性方程的解法等。但能求解的常微分方程仍然是有限的，大多數(shù)的常微分方程是不能給出解析解。譬如

這個一階微分方程就不能用初等函數(shù)來表達它的解。

2023/7/33再如，方程

的解,雖然有表可查,但對于表上沒有給出的值,仍需插值方法來計算2023/7/34從實際問題當(dāng)中歸納出來的微分方程，通常主要依靠數(shù)值解法來解決。本章主要討論一階常微分方程初值問題

(9.1)

在區(qū)間a≤x≤b上的數(shù)值解法。

可以證明,如果函數(shù)在帶形區(qū)域R:{a≤x≤b,-∞＜y＜∞}內(nèi)連續(xù)，且關(guān)于y滿足李普希茲(Lipschitz)條件，即存在常數(shù)L(它與x,y無關(guān))使

對R內(nèi)任意x及兩個都成立,則方程(9.1)的解在a,b上存在且唯一。

2023/7/35數(shù)值方法的基本思想對常微分方程初值問題(9.1)式的數(shù)值解法，就是要算出精確解y(x)在區(qū)間a,b上的一系列離散節(jié)點處的函數(shù)值的近似值相鄰兩個節(jié)點的間距稱為步長，步長可以相等，也可以不等。本章總是假定h為定數(shù)，稱為定步長，這時節(jié)點可表示為數(shù)值解法需要把連續(xù)性的問題加以離散化，從而求出離散節(jié)點的數(shù)值解。

2023/7/36

對常微分方程數(shù)值解法的基本出發(fā)點就是離散化。其數(shù)值解法有兩個基本特點，[1]它們都采用“步進式”，即求解過程順著節(jié)點排列的次序一步一步地向前推進，描述這類算法，要求給出用已知信息計算的遞推公式。[2]建立這類遞推公式的基本方法是在這些節(jié)點上用數(shù)值積分、數(shù)值微分、泰勒展開等離散化方法，對初值問題中的導(dǎo)數(shù)進行不同的離散化處理。

2023/7/37對于初值問題的數(shù)值解法，首先要解決的問題就是如何對微分方程進行離散化，建立求數(shù)值解的遞推公式。遞推公式通常有兩類，一類是計算yi+1時只用到xi+1,xi和yi，即前一步的值，因此有了初值以后就可以逐步往下計算，此類方法稱為單步法；其代表是龍格—庫塔法。另一類是計算yi+1時，除用到xi+1,xi和yi以外，還要用到，即前面k步的值，此類方法稱為多步法；其代表是亞當(dāng)斯法。

2023/7/38§9.2簡單的數(shù)值方法與基本概念9.2.1Euler公式歐拉（Euler）方法是解初值問題的最簡單的數(shù)值方法。初值問題的解y=y(x)為通過點的一條積分曲線。積分曲線上每一點的切線的斜率等于函數(shù)在這點的值。

2023/7/39Euler法的求解過程是:從初始點P0(即點(x0,y0))出發(fā),作積分曲線y=y(x)在P0點上切線(其斜率為

),與x=x1直線相交于P1點(即點(x1,y1),得到y(tǒng)1作為y(x1)的近似值,如上圖所示。過點(x0,y0),以f(x0,y0)為斜率的切線方程為

當(dāng)時,得這樣就獲得了P1點的坐標(biāo)。

2023/7/310同樣,過點P1(x1,y1),作積分曲線y=y(x)的切線交直線x=x2于P2點,切線的斜率直線方程為當(dāng)時,得2023/7/311當(dāng)時,得由此獲得了P2的坐標(biāo)。重復(fù)以上過程,就可獲得一系列的點:P1,P1,…,Pn。對已求得點以為斜率作直線

取2023/7/312從圖形上看,就獲得了一條近似于曲線y=y(x)的折線。這樣,從x0出發(fā)逐個算出2023/7/313通常取(常數(shù)),則

微分方程的初值問題Euler法的計算格式為：

用折線近似于曲線得到的計算公式稱為歐拉公式2023/7/314

Euler格式還可用數(shù)值微分、數(shù)值積分和泰勒展開等方法得到。

2023/7/315

Euler格式用泰勒展開方法得到。2023/7/316

Euler格式用數(shù)值微分方法得到。

選擇不同的計算方法計算上式的積分項就會得到不同的計算公式。

將方程的兩端在區(qū)間上積分得，2023/7/317

用左矩形方法計算積分項

代入(9.3)式,并用yi近似代替式中y(xi)即可得到向前歐拉（Euler）公式

由于數(shù)值積分的矩形方法精度很低，所以歐拉（Euler）公式當(dāng)然很粗糙。

2023/7/318例9.1用歐拉法解初值問題

取步長h=0.2,計算過程保留4位小數(shù)

解:當(dāng)k=0,x1=0.2時，已知x0=0，y0=1，有

y(0.2)y1=0.2×1(4－0×1)＝0.8歐拉迭代格式當(dāng)k=1,x2=0.4時，已知x1=0.2，y1=0.8，有

y(0.4)y2

=0.2×0.8×(4－0.2×0.8)＝0.6144當(dāng)k=2,x3=0.6時，已知x2=0.4，y2=0.6144，有

y(0.6)y3=0.2×0.6144×(4-0.4×0.6144)=0.46132023/7/319clear;y=1,x=0,%初始化forn=1:10y=1.1*y-0.2*x/y,x=x+0.1,endy=1x=0y=1.1000x=0.1000y=1.1918x=0.2000y=1.2774x=0.3000y=1.3582x=0.4000y=1.4351x=0.5000y=1.5090x=0.6000y=1.5803x=0.7000y=1.6498x=0.8000y=1.7178x=0.9000y=1.7848x=1.00002023/7/320對方程的兩端在區(qū)間上積分得，代入(9.4)式,即可得到梯形公式由于數(shù)值積分的梯形公式比矩形公式的精度高，因此梯形公式（9.5）是比歐拉公式(9.2)精度高的一個數(shù)值方法。為了提高精度,改用梯形方法計算其積分項，即

9.2.2梯形公式2023/7/321(9.5)

(9.2)和(9.5)式粗看差不多，

(9.2)

Euler公式梯形公式

但(9.5)式的右端含有未知的yi+1,它是一個關(guān)于yi+1的函數(shù)方程,這類數(shù)值方法稱為隱式方法。相反地,歐拉法是關(guān)于yi+1的一個直接的計算公式，這類數(shù)值方法稱為顯式方法。

2023/7/3229.2.3兩步歐拉公式

對方程的兩端在區(qū)間上積分得

(9.6)

改用中矩形公式計算其積分項，即

代入上式,并用yi近似代替式中y(xi)即可得到兩步歐拉公式2023/7/323前面介紹過的數(shù)值方法,無論是歐拉方法,還是梯形方法，它們都是單步法,其特點是在計算yi+1時只用到前一步的信息yi;可是公式(9.7)中除了yi外,還用到更前一步的信息yi-1,即調(diào)用了前兩步的信息,故稱其為兩步歐拉公式

2023/7/324定義9.1在yi準(zhǔn)確的前提下,即時,用數(shù)值方法計算yi+1的誤差,稱為該數(shù)值方法計算時yi+1的局部截斷誤差。9.2.4．歐拉法的局部截斷誤差衡量求解公式好壞的一個主要標(biāo)準(zhǔn)是求解公式的精度,因此引入局部截斷誤差和階數(shù)的概念。2023/7/325假定

由歐拉公式，則有因此有

局部截斷誤差2023/7/326定義9.2如果數(shù)值方法的局部截斷誤差為

,則稱這種數(shù)值方法的階數(shù)是P。顯然，步長(h<1)越小,P越高,則局部截斷誤差越小,計算精度越高。兩步歐拉公式是比歐拉公式精度高的一個數(shù)值方法,歐拉公式的局部截斷誤差為,歐拉方法僅為一階方法。2023/7/327設(shè)

前兩步準(zhǔn)確,則兩步歐拉公式

把y(xi-1)在xi處展開成Taylor級數(shù)，即

①把y(xi-1)代人①，即

(9.8)2023/7/328再將y(xi+1)在xi處進行泰勒展開(9.8)(9.9)所以,由(9.8)和(9.9)可得兩步歐拉公式的局部截斷誤差為：

故兩步歐拉公式為二階方法。2023/7/329顯式歐拉公式計算工作量小,但精度低。梯形公式雖提高了精度,但為隱式公式,需用迭代法求解,計算工作量非常大。將歐拉公式和梯形公式綜合，便可得到改進的歐拉公式。9.2.5改進的歐拉公式歐拉公式梯形公式2023/7/330先用歐拉公式(9.2)求出一個初步的近似值,稱為預(yù)測值,它的精度不高,再用梯形公式(9.5)對它校正一次,即再迭代一次,求得yi+1,稱為校正值,歐拉公式梯形公式這種預(yù)測-校正方法稱為改進的歐拉公式。2023/7/331(9.10)

預(yù)測

校正左矩形公式梯形公式改進的歐拉公式右矩形公式2023/7/332可以證明,改進的歐拉公式(9.10)的精度為二階。這是一種一步顯式格式,它可以表示為嵌套形式?；蛘弑硎境上铝衅骄问?023/7/333

改進的歐拉公式左矩形公式右矩形公式左、右矩形公式的平均值，即梯形公式.2023/7/3349.2.6改進歐拉法算法實現(xiàn)（1）計算步驟①輸入,h,N②使用以下改進歐拉法公式進行計算

③

輸出，并使轉(zhuǎn)到

②

直至n>N結(jié)束。

2023/7/335（2）改進歐拉法的流程圖

2023/7/336例9.2用改進歐拉法解初值問題區(qū)間為0,1,取步長h=0.1

解:改進歐拉法的具體形式本題的精確解為2023/7/337clearx=0,yn=1%初始化forn=1:10yp=yn+0.1*(yn-2*x/yn);%預(yù)測x=x+0.1;yc=yn+0.1*(yp-2*x/yp);yn=(yp+yc)/2%校正end2023/7/338例9.3對初值問題

證明用梯形公式求得的近似解為

并證明當(dāng)步長h0時,yn收斂于精確解整理成顯式

反復(fù)迭代,得到證明:解初值問題的梯形公式為2023/7/339由于，有

∴

證畢2023/7/340§9.3龍格-庫塔（Runge-Kutta）法只有對平均斜率提供一種算法便可得到一種計算格式9.3.1龍格-庫塔法的基本思想2023/7/341Euler公式可改寫成改進的Euler公式又可改寫成用左端點的斜率作為平均斜率得到的計算格式用左、右端點斜率的平均值作為平均斜率得到的計算格式2023/7/342上述兩組公式在形式上有一個共同點:都是用f(x,y)在某些點上值(斜率)的線性組合得出y(xi+1)的近似值yi+1,而且增加計算f(x,y)的次數(shù),可提高截斷誤差的階。如歐拉公式:每步計算一次f(x,y)的值,為一階方法。改進歐拉公式需計算兩次f(x,y)的值，它是二階方法。它的局部截斷誤差為。2023/7/343

于是可考慮在內(nèi)多預(yù)報幾個點的斜率值，然后將其加權(quán)平均作為平均斜率，構(gòu)造時要求近似公式在(xi,yi)處的Taylor展開式與解y(x)在xi處的Taylor展開式的前面幾項重合，從而使近似公式達到所需要的階數(shù)。則可構(gòu)造出更高精度的計算格式，這就是龍格—庫塔（Runge-Kutta）法的基本思想。

2023/7/344只有多計算幾個點的函數(shù)值，平均高度就更精確.龍格-庫塔法的基本思想多預(yù)報幾個點的斜率值，將其加權(quán)平均作為平均斜率.2023/7/3459.3.2二階龍格—庫塔法在上取兩點xi和,以該兩點處的斜率值k1和k2的加權(quán)平均(或稱為線性組合)來求取平均斜率k*的近似值K，即

式中:k1為xi點處的切線斜率值，

k2為點處的切線斜率值,比照改進的歐拉法,將視為，即可得對常微分方程初值問題(9.1)式的解y=y(x),根據(jù)微分中值定理，存在點，使得2023/7/346式中K可看作是y=y(x)在區(qū)間上的平均斜率。所以可得計算公式為：（9.14）

將在x=xi處進行二階Taylor展開：

（9.15）

也即（9.13）2023/7/347將以上結(jié)果代入2023/7/348對式(9.15)和(9.16)進行比較系數(shù)后可知,只要成立,格式(9.14)的局部截斷誤差就等于2023/7/349式(9.17)中具有三個未知量,但只有兩個方程,因而有無窮多解。若取,則P=1，這是無窮多解中的一個解，將以上所解的值代入式(9.14)并改寫可得

不難發(fā)現(xiàn)，上面的格式就是改進的歐拉格式。凡滿足條件式（9.17）有一簇形如上式的計算格式，這些格式統(tǒng)稱為二階龍格—庫塔格式。因此改進的歐拉格式是眾多的二階龍格—庫塔法中的一種特殊格式。2023/7/350若取,則，此時二階龍格-庫塔法的計算公式為

此計算公式稱為變形的二階龍格—庫塔法。式中為區(qū)間的中點。2023/7/3519.3.3三階龍格-庫塔法

為了進一步提高精度，設(shè)除外再增加一點并用三個點的斜率k1,k2,k3加權(quán)平均得出平均斜率k*的近似值，這時計算格式具有形式:（9.18）

為了預(yù)報點的斜率值k3,在區(qū)間內(nèi)有兩個斜率值k1和k2可以用,可將k1,k2加權(quán)平均得出上的平均斜率,從而得到的預(yù)報值2023/7/352于是可得

運用Taylor展開方法選擇參數(shù)

,可以使格式(9.18)的局部截斷誤差為

,即具有三階精度.2023/7/353局部截斷誤差為

,即具有三階精度，這類格式統(tǒng)稱為三階龍格—庫塔方法。（9.19）

是其中的一種，稱為庫塔（Kutta）公式。

2023/7/3549.3.4四階龍格—庫塔法

如果需要再提高精度，用類似上述的處理方法，只需在區(qū)間上用四個點處的斜率加權(quán)平均作為平均斜率k*的近似值，構(gòu)成一系列四階龍格—庫塔公式。具有四階精度，即局部截斷誤差是。由于推導(dǎo)復(fù)雜，這里從略，只介紹最常用的一種四階經(jīng)典R—K公式。

（9.20）

2023/7/3559.3.5四階龍格—庫塔法算法實現(xiàn)(1)

計算步驟①輸入,h,N②使用龍格—庫塔公式（9.20）計算出y1③輸出，并使轉(zhuǎn)到②直至n>N結(jié)束。2023/7/356(2)四階龍格—庫塔算法流程圖2023/7/357程序?qū)崿F(xiàn)(四階龍格-庫塔法計算常微分方程初值問題)

例9.4取步長h=0.2，用經(jīng)典格式求解初值問題解:

由四階龍格-庫塔公式可得

2023/7/358可同樣進行其余yi的計算。本例方程的解為，從表中看到所求的數(shù)值解具有4位有效數(shù)字。

龍格—庫塔法的推導(dǎo)基于Taylor展開方法，因而它要求所求的解具有較好的光滑性。如果解的光滑性差，那么，使用四階龍格—庫塔方法求得的數(shù)值解，其精度可能反而不如改進的歐拉方法。在實際計算時，應(yīng)當(dāng)針對問題的具體特點選擇合適的算法。

2023/7/3599.3.6變步長的龍格-庫塔法在微分方程的數(shù)值解中，選擇適當(dāng)?shù)牟介L是非常重要的。單從每一步看，步長越小，截斷誤差就越小；但隨著步長的縮小，在一定的求解區(qū)間內(nèi)所要完成的步數(shù)就增加了。這樣會引起計算量的增大，并且會引起舍入誤差的大量積累與傳播。因此微分方程數(shù)值解法也有選擇步長的問題。以經(jīng)典的四階龍格-庫塔法(9.20)為例。從節(jié)點xi出發(fā)，先以h為步長求出一個近似值，記為，由于局部截斷誤差為，故有

當(dāng)h值不大時，式中的系數(shù)c可近似地看作為常數(shù)。2023/7/360然后將步長折半,即以為步長,從節(jié)點xi出發(fā),跨兩步到節(jié)點xi+1,再求得一個近似值,每跨一步的截斷誤差是,因此有這樣由此可得這表明以作為的近似值，其誤差可用步長折半前后兩次計算結(jié)果的偏差來判斷所選步長是否適當(dāng)2023/7/361當(dāng)要求的數(shù)值精度為ε時：（1）如果Δ>ε，反復(fù)將步長折半進行計算，直至Δ<ε為止,并取其最后一次步長的計算結(jié)果作為（2）如果Δ<ε，反復(fù)將步長加倍，直到Δ>ε為止，并以上一次步長的計算結(jié)果作為。這種通過步長加倍或折半來處理步長的方法稱為變步長法。表面上看，為了選擇步長，每一步都要反復(fù)判斷Δ，增加了計算工作量，但在方程的解y(x)變化劇烈的情況下，總的計算工作量得到減少，結(jié)果還是合算的。2023/7/362

§9.4算法的穩(wěn)定性及收斂性9.4.1穩(wěn)定性穩(wěn)定性在微分方程的數(shù)值解法中是一個非常重要的問題。因為微分方程初值問題的數(shù)值方法是用差分格式進行計算的，而在差分方程的求解過程中，存在著各種計算誤差，這些計算誤差如舍入誤差等引起的擾動，在傳播過程中，可能會大量積累，對計算結(jié)果的準(zhǔn)確性將產(chǎn)生影響。這就涉及到算法穩(wěn)定性問題。

2023/7/363

當(dāng)在某節(jié)點上xi的yi值有大小為δ的擾動時，如果在其后的各節(jié)點上的值yi產(chǎn)生的偏差都不大于δ，則稱這種方法是穩(wěn)定的。

穩(wěn)定性不僅與算法有關(guān)，而且與方程中函數(shù)f(x,y)也有關(guān)，討論起來比較復(fù)雜。為簡單起見，通常只針對模型方程來討論。一般方程若局部線性化，也可化為上述形式。模型方程相對比較簡單，若一個數(shù)值方法對模型方程是穩(wěn)定的，并不能保證該方法對任何方程都穩(wěn)定，但若某方法對模型方程都不穩(wěn)定，也就很難用于其他方程的求解。2023/7/364先考察顯式Euler方法的穩(wěn)定性。模型方程的Euler公式為將上式反復(fù)遞推后，可得或式中2023/7/365

要使yi有界，其充要條件為即由于，故有（9.26）可見，如欲保證算法的穩(wěn)定，顯式Euler格式的步長h的選取要受到式（9.26）的限制。的絕對值越大，則限制的h值就越小。用隱式Euler格式，對模型方程的計算公式為，可化為2023/7/366由于,則恒有,故恒有因此，隱式Euler格式是絕對穩(wěn)定的（無條件穩(wěn)定）的（對任何h>0）。9.4.2收斂性常微分方程初值問題的求解，是將微分方程轉(zhuǎn)化為差分方程來求解，并用計算值yi來近似替代y(xi)，這種近似替代是否合理，還須看分割區(qū)間的長度h越來越小時，即時，是否成立。若成立，則稱該方法是收斂的，否則稱為不收斂。這里仍以Euler方法為例，來分析其收斂性。Euler格式2023/7/367設(shè)表示取時,按Euler公式的計算結(jié)果,即Euler方法局部截斷誤差為：設(shè)有常數(shù)，則（9.27）

總體截斷誤差（9.28）

又由于f(x,y)關(guān)于y滿足李普希茨條件,即2023/7/368代入(9.28)上式，有再利用式（9.27），式（9.28）即上式反復(fù)遞推后，可得（9.29）

2023/7/369設(shè)（T為常數(shù)）

因為

所以把上式代入式（9.29），得若不計初值誤差，即，則有（9.30）

式(9.30)說明,當(dāng)時,,即,所以Euler方法是收斂的，且其收斂速度為，即具有一階收斂速度。同時還說明Euler方法的整體截斷誤差為，因此算法的精度為一階。2023/7/370龍格-庫塔方法是一類重要算法，但這類算法在每一步都需要先預(yù)報幾個點上的斜率值，計算量比較大。考慮到計算yi+1之前已得出一系列節(jié)點上的斜率值，能否利用這些已知的斜率值來減少計算量呢？這就是亞當(dāng)姆斯（Adams）方法的設(shè)計思想。

§9.5亞當(dāng)姆斯方法9.5.1亞當(dāng)姆斯格式2023/7/371設(shè)用xi,xi+1兩點的斜率值加權(quán)平均作為區(qū)間上的平均斜率，有計算格式

（9.21）

選取參數(shù)λ，使格式（9.21）具有二階精度。2023/7/372將在xi處Taylor展開

代入計算格式化簡,因此有與y(xi+1)在xi處的Taylor展開式相比較,需取

才使格式(9.21)具有二階精度。這樣導(dǎo)出的計算格式稱之為二階亞當(dāng)姆斯格式。2023/7/373四階亞當(dāng)姆斯格式。

(9.22)上述幾種亞當(dāng)姆斯格式都是顯式的，算法比較簡單，但用節(jié)點的斜率值來預(yù)報區(qū)間上的平均斜率是個外推過程，效果不夠理想。為了進一步改善精度，變外推為內(nèi)插，即增加節(jié)點xi+1的斜率值來得出上的平均斜率。譬如考察形如類似地可以導(dǎo)出三階亞當(dāng)姆斯格式。

2023/7/374（9.23）

的隱式格式,設(shè)(9.23)右端的Taylor展開有

可見要使格式(9.23)具有二階精度,需令,這樣就可構(gòu)造二階隱式亞當(dāng)姆斯格式其實是梯形格式。類似可導(dǎo)出三階隱式亞當(dāng)姆斯格式2023/7/375和四階隱式亞當(dāng)姆斯格式

（9.24）

9.5.2亞當(dāng)姆斯預(yù)報-校正格式參照改進的歐拉格式的構(gòu)造方法，以四階亞當(dāng)姆斯為例，將顯式（9.22）和隱式（9.24）相結(jié)合，用顯式公式做預(yù)報，再用隱式公式做校正，可構(gòu)成亞當(dāng)姆斯預(yù)報-校正格式（9.25）

預(yù)報：

校正：

2023/7/376這種預(yù)報-校正格式是四步法，它在計算yi+1時不但用到前一步的信息，而且要用到再前面三步的信息，因此它不能自行啟動。在實際計算時，可借助于某種單步法，譬如四階龍格—庫塔法提供開始值。

2023/7/377例9.5取步長h=0.1,用亞當(dāng)姆斯預(yù)報-校正公式求解初值問題

的數(shù)值解。解:

用四階龍格-庫塔公式求出初值

,計算得：表中的,yi和y(xi)分別為預(yù)報值、校正值和準(zhǔn)確解(),以比較計算結(jié)果的精度。再使用亞當(dāng)姆斯預(yù)報-校正公式(9.25),2023/7/378§9.6一階常微分方程組與高階方程我們已介紹了一階常微分方程初值問題的各種數(shù)值解法，對于一階常微分方程組，可類似得到各種解法，而高階常微分方程可轉(zhuǎn)化為一階常微分方程組來求解。

9.6.1一階常微分方程組對于一階常微分方程組的初值問題

（9.31）

可以把單個方程中的f和y看作向量來處理，這樣就可把前面介紹的各種差分算法推廣到求一階方程組初值問題中來。2023/7/379設(shè)為節(jié)點上的近似解，則有改進的Euler格式為預(yù)報：校正：

（9.32）

又，相應(yīng)的四階龍格—庫塔格式（經(jīng)典格式）為（9.33）

2023/7/380式中

（9.34）

2023/7/381把節(jié)點xi上的yi和zi值代入式(7.34),依次算出

,然后把它們代入式(7.33),算出節(jié)點xi+1上的yi+1

和zi+1值。對于具有三個或三個以上方程的方程組的初