常微分方程初值問題數值解法-課件

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

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

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

2023/7/33再如，方程

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

(9.1)

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

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

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

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

2023/7/36

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

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

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

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

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

當時,得這樣就獲得了P1點的坐標。

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

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

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

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

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

2023/7/315

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

Euler格式用數值微分方法得到。

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

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

用左矩形方法計算積分項

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

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

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

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

解:當k=0,x1=0.2時，已知x0=0，y0=1，有

y(0.2)y1=0.2×1(4－0×1)＝0.8歐拉迭代格式當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當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)式,即可得到梯形公式由于數值積分的梯形公式比矩形公式的精度高，因此梯形公式（9.5）是比歐拉公式(9.2)精度高的一個數值方法。為了提高精度,改用梯形方法計算其積分項，即

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

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

(9.2)

Euler公式梯形公式

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

2023/7/3229.2.3兩步歐拉公式

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

(9.6)

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

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

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

由歐拉公式，則有因此有

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

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

前兩步準確,則兩步歐拉公式

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

①把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)求出一個初步的近似值,稱為預測值,它的精度不高,再用梯形公式(9.5)對它校正一次,即再迭代一次,求得yi+1,稱為校正值,歐拉公式梯形公式這種預測-校正方法稱為改進的歐拉公式。2023/7/331(9.10)

預測

校正左矩形公式梯形公式改進的歐拉公式右矩形公式2023/7/332可以證明,改進的歐拉公式(9.10)的精度為二階。這是一種一步顯式格式,它可以表示為嵌套形式。或者表示成下列平均化形式2023/7/333

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

③

輸出，并使轉到

②

直至n>N結束。

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);%預測x=x+0.1;yc=yn+0.1*(yp-2*x/yp);yn=(yp+yc)/2%校正end2023/7/338例9.3對初值問題

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

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

反復迭代,得到證明:解初值問題的梯形公式為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)的次數,可提高截斷誤差的階。如歐拉公式:每步計算一次f(x,y)的值,為一階方法。改進歐拉公式需計算兩次f(x,y)的值，它是二階方法。它的局部截斷誤差為。2023/7/343

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

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

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

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

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

（9.15）

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

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

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

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

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

運用Taylor展開方法選擇參數

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

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

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

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

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

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

（9.20）

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

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

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

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

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

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

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

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

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

2023/7/363

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

穩(wěn)定性不僅與算法有關，而且與方程中函數f(x,y)也有關，討論起來比較復雜。為簡單起見，通常只針對模型方程來討論。一般方程若局部線性化，也可化為上述形式。模型方程相對比較簡單，若一個數值方法對模型方程是穩(wěn)定的，并不能保證該方法對任何方程都穩(wěn)定，但若某方法對模型方程都不穩(wěn)定，也就很難用于其他方程的求解。2023/7/364先考察顯式Euler方法的穩(wěn)定性。模型方程的Euler公式為將上式反復遞推后，可得或式中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收斂性常微分方程初值問題的求解，是將微分方程轉化為差分方程來求解，并用計算值yi來近似替代y(xi)，這種近似替代是否合理，還須看分割區(qū)間的長度h越來越小時，即時，是否成立。若成立，則稱該方法是收斂的，否則稱為不收斂。這里仍以Euler方法為例，來分析其收斂性。Euler格式2023/7/367設表示取時,按Euler公式的計算結果,即Euler方法局部截斷誤差為：設有常數，則（9.27）

總體截斷誤差（9.28）

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

2023/7/369設（T為常數）

因為

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

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

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

（9.21）

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

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

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

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

2023/7/374（9.23）

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

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

（9.24）

9.5.2亞當姆斯預報-校正格式參照改進的歐拉格式的構造方法，以四階亞當姆斯為例，將顯式（9.22）和隱式（9.24）相結合，用顯式公式做預報，再用隱式公式做校正，可構成亞當姆斯預報-校正格式（9.25）

預報：

校正：

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

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

的數值解。解:

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

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

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

（9.31）

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

（9.32）

又，相應的四階龍格—庫塔格式（經典格式）為（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值。對于具有三個或三個以上方程的方程組的初