微分方程初值問題的數值解法

微分方程初值問題的數值解法第一頁，共四十九頁，2022年，8月28日

包含自變量、未知函數及未知函數的導數或微分的方程稱為微分方程。在微分方程中,自變量的個數只有一個,稱為常微分方程。自變量的個數為兩個或兩個以上的微分方程叫偏微分方程。微分方程中出現的未知函數最高階導數的階數稱為微分方程的階數。如果未知函數y及其各階導數8.1

引言

微分方程知識回顧都是一次的,則稱它是線性的,否則稱為非線性的。

第二頁，共四十九頁，2022年，8月28日

在高等數學中，對于常微分方程的求解，給出了一些典型方程求解析解的基本方法：

一階：可分離變量法、齊次方程、一階線性方程以及伯努利方程

高階：可降階方程、常系數齊次線性方程的解法、常系數非齊次線性方程的解法但能求解的常微分方程仍然是有限的，大多數的常微分方程是求不出解析解的。

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

微分方程知識回顧

例如第三頁，共四十九頁，2022年，8月28日

從實際問題當中歸納出來的微分方程，通常主要依靠數值解法來解決。本章主要討論一階常微分方程初值問題

(8.1)

在區(qū)間

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

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

上的數值解法。第四頁，共四十九頁，2022年，8月28日

常微分方程初值問題(8.1)式的數值解法，首先要算出精確解y(x)在區(qū)間a,b上的一系列離散節(jié)點處的函數值的近似值:

相鄰兩個節(jié)點的間距稱為步長，本章總是假定h為定數，稱為定步長，這時節(jié)點可表示為

數值方法的基本思想1、數值解法需要把連續(xù)性的問題加以離散化，從而求出離散節(jié)點處的數值解。第五頁，共四十九頁，2022年，8月28日

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

數值方法的基本思想第六頁，共四十九頁，2022年，8月28日

遞推公式通常有兩類，一類是計算yi+1時只用到xi+1,xi和yi，即前一步的值，因此有了初值以后就可以逐步往下計算，此類方法稱為單步法，其代表是龍格—庫塔法。另一類是計算yi+1時，除用到xi+1,xi和yi以外，還要用到，即前面k步的值，此類方法稱為多步法，其代表是亞當斯法。

數值方法的基本思想第七頁，共四十九頁，2022年，8月28日一、Euler方法及其改進

將[a,b]n等分,記

微分法:

積分法:

積分項利用矩形公式計算

1.顯式Euler方法

(★)第八頁，共四十九頁，2022年，8月28日Taylor公式推導:

第九頁，共四十九頁，2022年，8月28日Oyxy=y(x)(x1,y1)p1p0x0x1x2xixi+1xn-1xnpipi+1pn-1pnp’1p’2p’ip’i+1p’n-1p’n切線p0p1的斜率為f(x0,y0)p2(x2,y2)

歐拉公式的幾何意義：第十頁，共四十九頁，2022年，8月28日Euler法的求解過程是:

從初始點P0(即點(x0,y0))出發(fā),作積分曲線y=y(x)在P0點上切線(其斜率為),與直線x=x1相交于P1點(即點(x1,y1),得到y1作為y(x1)的近似值）這樣就獲得了P1點的坐標。當時,得

重復以上過程,就可獲得一系列的點:p1,p2,…,pn，相應的可求出y1,y2,…,yn,取第十一頁，共四十九頁，2022年，8月28日從圖形上看,就獲得了一條近似于曲線y=y(x)的折線通常取(常數),則Euler法的計算格式

i=0,1,…,n(7.2)第十二頁，共四十九頁，2022年，8月28日2.梯形法

稱之為梯形公式.這是一個隱式公式,通常用迭代法求解.具體做法:

取

先用Euler法求出初值,即,將其代入梯形公式的右端,使之轉化為顯式公式,即

注:

當f(x,y)關于y滿足Lipschitz條件且步長h滿足

直至滿足:

若采用梯形公式計算(★)中的積分項,則有類似地,可得(☆)

第十三頁，共四十九頁，2022年，8月28日時,迭代格式

(☆)收斂

.

3.改進的Euler方法

把Euler法作為預報(稱為預估公式),把隱式的梯形公式作為校正(稱為校正公式

),則得改進的Euler方法:或也稱為預估-校正法.第十四頁，共四十九頁，2022年，8月28日有時為了方便,預估-校正格式也寫成下面形式:第十五頁，共四十九頁，2022年，8月28日

改進的歐拉公式比歐拉公式精度高的原因是：改進歐拉公式用梯形面積代替曲邊梯形面積，而歐拉公式用矩形面積代替曲邊梯形面積。數值積分的梯形公式比矩形公式的精度高。第十六頁，共四十九頁，2022年，8月28日二、單步法的局部截斷誤差及精度

Def1:

先假設,再估計誤差這種誤差稱為單步迭代法在xk+1處的局部截斷誤差.Def2:

若某種數值方法的局部截斷誤差為,則稱該數值方法的精度為P階的.注:

通常情況下,P越大,h越小,則截斷誤差越小,數值方法越精確.第十七頁，共四十九頁，2022年，8月28日所以Euler方法為一階方法.而設10.Euler方法是一階方法.第十八頁，共四十九頁，2022年，8月28日20.梯形法是二階方法.Taylor展開

第十九頁，共四十九頁，2022年，8月28日將代入上式,得而代入上式得:當h充分小時,若,則可選取h,使得第二十頁，共四十九頁，2022年，8月28日故梯形法的精度為2.同樣可以證明改進的Euler法也是二階方法.梯形法的局部截斷誤差為:從而第二十一頁，共四十九頁，2022年，8月28日例1:

取步長

h=2/10,2/20,2/30,2/40,分別用歐拉法、改進的歐拉法和梯形法求解.解:

記

f(x,y)=y－xy2,xk=kh(k=0,1,2,···,n)(1).Euler法:yk+1=yk+h(yk－xkyk2)(k=0,1,···,n)

y0=1當

h=2/10時,n=10.由Euler公式可得:k01234yk+11.21.38241.5061.535041.46503k56789yk+11.328771.170771.021130.891690.783788第二十二頁，共四十九頁，2022年，8月28日(2).改進的Euler法:

k01234yk+11.19121.343841.423481.419051.3473k56789yk+11.237261.114240.9941510.8847510.788666(3).梯形法(計算過程略)

第二十三頁，共四十九頁，2022年，8月28日n10203040h0.20.10.06670.05誤差

0.10590.05210.03420.0256Euler法誤差:改進的Euler法誤差:n10203040h0.20.10.06670.05誤差

0.01230.00260.00115.9612e-004第二十四頁，共四十九頁，2022年，8月28日預-校方法,h=0.2時誤差最大值:0.0123歐拉方法,h=0.2時誤差最大值:0.1059解析解:第二十五頁，共四十九頁，2022年，8月28日三、Runge-Kutta方法1、Taylor級數法

設初值問題有解y(x),由Tayler公式得:令當時,有.此時①為p階Taylor方法.p=1時即為Euler公式.稱之為Taylor級數法.其中例2:

取步長h=0.1,用一階、二階和四階Taylor方法求解下列初值問題①第二十六頁，共四十九頁，2022年，8月28日解:(1)一階Taylor法k01234yk+11.11.2211.370081.557791.80046(2)二階Taylor法k01234yk+11.111.246891.421751.652631.97088第二十七頁，共四十九頁，2022年，8月28日(3)四階Taylor法k01234yk+11.11111.249961.428481.666441.99942第二十八頁，共四十九頁，2022年，8月28日記由得稱為[xk,xk+1]上的平均斜率.故2、Runge-Kutta方法只要對K*提供不同的算法,就會得出不同的計算公式.如取則得改進的Euler公式,它是利用xk,xk+1兩點的斜率值K1,K2的算術平均值作為K*,精度比Euler法高.則得Euler公式;取第二十九頁，共四十九頁，2022年，8月28日Runge-Kutta法的基本思想:

設法在[xk,xk+1]內多預報幾個點的斜率,再將它們的加權平均值作為平均斜率K*一般顯式Runge-Kutta公式為:其中為待定參數,且.稱為r級Runge-Kutta方法計算公式.②第三十頁，共四十九頁，2022年，8月28日即可得p個方程,從而確定出待定參數.代入表達式即可得到計算公式.如果要求兩個表達式的前p+1項完全重合,即局部截斷誤差達到,則稱②式為p階r級的Runge-Kutta方法.常用的是r=2,3,4

級的R-K方法,且適當選取參數使得p=r

.如要求:注:

式中待定參數的確定:

先將②式右端在(xk,yk)處展成h的冪級數(即將yk+1展成h的冪級數);再將y(xk+1)作Taylor級數展開;最后比較兩式中hk(k=0,1,2,…)的系數,以確定出所有待定參數.第三十一頁，共四十九頁，2022年，8月28日Runge-Kutta方法的推導(以r=2為例):當r=2時記第三十二頁，共四十九頁，2022年，8月28日則又第三十三頁，共四十九頁，2022年，8月28日這是一個四個參數三個方程的非線性方程組.它有一個自由度.稱滿足上述方程組的一族公式為二級二階Runge-Kutta方法.為使局部截斷誤差為,比較上述兩式右端同次冪系數,應取第三十四頁，共四十九頁，2022年，8月28日(1)常用的二階Runge-Kutta方法:預估-校正算法(2)中間點方法

第三十五頁，共四十九頁，2022年，8月28日注:

二級Runge-Kutta方法的精度最高是二階的,不可能達到三階.要提高計算方法的階,就必須增加預報點.常用的三階Runge-Kutta方法(r=3):

(1)Heun(休恩)方法

(3)三階Kutta方法

第三十六頁，共四十九頁，2022年，8月28日(1)三階Heun方法

標準(經典)四階Runge-Kutta方法

(2)常用的四階Runge-Kutta方法(r=4):

第三十七頁，共四十九頁，2022年，8月28日(2)稱為Gill(吉爾)方法

注:

從理論上講,可以構造任意高階的計算方法.但事實上,精度的階數與預報點的個數之間并非等量關系.預報點的個數r123456789r≥10精度的階數123445667≤r-2一般情況下,四階Runge-Kutta方法已可滿足精度要求.第三十八頁，共四十九頁，2022年，8月28日例3:

用經典Runge-Kutta方法求解下列初值問題(取h=0.1)解:標準Runge-Kutta公式為:計算結果見下表.為比較在相同計算量條件下近似解的精度,表中列出了Euler法(h=0.025)和改進的Euler法(h=0.05)在相應節(jié)點上的計算結果.第三十九頁，共四十九頁，2022年，8月28日xiEuler法h=0.025改進Euler法h=0.05經典R-K法h=0.1準確解0.11.1114391.1153801.1155121.1155130.21.2552091.2639141.2642081.2642080.31.4346671.4490891.4495761.4495760.41.6535171.6747561.6754731.6754740.51.9158491.9451711.9461621.9461640.62.2261782.2650402.2663542.2663560.72.5894852.6395612.6412552.6412580.83.0112713.0744793.0766193.0766230.93.4976063.5761443.5788043.5788091.04.0551924.1515734.1548394.154845注:

用表中每種方法計算yi都需要計算四次f的值,即它們的計算量基本相等.第四十頁，共四十九頁，2022年，8月28日四、單步法的進一步討論—收斂性、相容性與穩(wěn)定性注:

由定義可知,數值方法的收斂性并不涉及計算過程的舍入誤差,只與方法的截斷誤差有關.若格式收斂,則整體截斷誤差必趨于零.Def:

(整體截斷誤差)

稱為某一數值方法在點xk處的整體截斷誤差.它不僅與xk有關,也與xk-1,xk-2,…,x1,x0有關.則稱該單步法收斂.Def:

對滿足解存在唯一性條件的初值問題(1),如果一個顯式單步法(3)產生的近似解對于任一固定的,均有1.收斂性第四十一頁，共四十九頁，2022年，8月28日由于,且關于y滿足Lipschitz條件,得則存在常數c>0使得且單步法中函數關于y滿足Lipschitz條件,則定理1:

若初值問題的一個單步法的局部截斷誤差為記證:由局部截斷誤差的定義知第四十二頁，共四十九頁，2022年，8月28日故從而有故若y(x0)=y0,則e0=0,由不等式得第四十三頁，共四十九頁，2022年，8月28日設單步法為注:定理表明,數值方法的整體截斷誤差比局部截斷誤差低一階.收斂的方法至少是一階方法.在該定義條件下,Euler方法是一階的,預估-校正方法是二階.當f(x,y)關于y也滿足Lipschitz條件,r級Runge-Kutta方法中的φ

關于y也滿足Lipschitz條件,故定理中的條件得到滿足,解的收斂性得到保證.由于Rn,h→0(h→0),且xn為任意點,故該式相當于用近似方程當x=xn+1固定時,,所以有2.相容性第四十四頁，共四十九頁，2022年，8月28日通過在x=xn處求解近似方程而獲得原方程的近似解.因此,必須要求當h→0

時,近似方程應逼近于原方程.來代替因此,要使h→0

時,近似方程的極限狀態(tài)為原微分方程,需且只需下列極限成立:由于由于假設是連續(xù)函數,故上式可表示為Def:

如果當h→0時,

近似方程逼近微分方程,則稱數值公式與原微分方程相容.相容的充要條件:

第四十五頁，共四十九頁，2022年，8月28日事實上:Remark:

可以證明若單步法的階大于或等于1,則單步法與微分方程相容;反之,如果單步法與微分方程相容,且關于y滿足Lipschitz條件,則單步法至少為一階方法.(h→0)(1)若單步法的階大于或等于1,由知即單步法與微分方程相容.故有(2)如果單步法與微分方程相容,且關于y滿足Lipschitz條件,則第四十六頁，共四十九頁，2022年，8月28日關于單步法的收斂性以及收斂性定理都是在計算過程中無任何舍入誤差的前提條件下建立的,但在實際計算時通常會有舍入誤差及其積累,數值求解微分方程的過程是一個遞推公式,必須考

即與微分方程相容的單步法至少為一階方法.Remark:

在定理條件下,Euler方法、預估-校正方法以及Runge-Kutta方法都與原微分方程相容.中連續(xù),且關于變量y滿足Lipschitz條件,則單步法收斂的充要條件為相容性條件成立.Th1.

設增量函數