微分方程初值問題的數(shù)值解法

微分方程初值問題的數(shù)值解法第1頁(yè)，課件共41頁(yè)，創(chuàng)作于2023年2月引言初值問題的數(shù)值解法:求初值問題的解在一系列節(jié)點(diǎn)的值y(xn)的近似值yn的方法.本章數(shù)值解法的特點(diǎn):都是采用“步進(jìn)式”,即求解過程順著節(jié)點(diǎn)排列的次序一步步向前推進(jìn).基本知識(shí):(1)定理1:如果函數(shù)f(x,y)在區(qū)域上連續(xù),且關(guān)于y滿足Lipschitz條件常微分方程初值問題:求未知函數(shù)y=y(x).

第2頁(yè)，課件共41頁(yè)，創(chuàng)作于2023年2月此時(shí)Lipschitz條件顯然成立.故常用在D上連續(xù)有界來代替f(x,y)關(guān)于y滿足Lipschitz條件.注:如無特別說明,總假設(shè)(1)的解存在唯一且足夠光滑.在

f(x,y)對(duì)變量y可微的情形下,若偏導(dǎo)數(shù)連續(xù)有界,則可取L為除了要保證(1)有唯一解外,還需保證微分方程本身是穩(wěn)定的,即(1)的解連續(xù)依賴于初始值和函數(shù)f(x,y).也就是說,當(dāng)初始值y0及函數(shù)f(x,y)有微小變化時(shí),只能引起解的微小變化.(其中L稱為L(zhǎng)ipschitz常數(shù)),則對(duì)任何,初值問題(1)在[a,b]上存在唯一連續(xù)可微解y=y(x).定理2:如果函數(shù)f(x,y)在區(qū)域上關(guān)于y滿足Lipschitz條件,則(1)是穩(wěn)定的.第3頁(yè)，課件共41頁(yè)，創(chuàng)作于2023年2月單步迭代:計(jì)算yn+1時(shí)僅用yn;初值問題(1)與下列積分方程的解等價(jià):

初值問題的數(shù)值解就是求一系列節(jié)點(diǎn)上函數(shù)y=y(x)的近似值.

稱為步長(zhǎng).一般取等步長(zhǎng)h.多步迭代:計(jì)算yn+1時(shí)除用yn

外,還要用到y(tǒng)n-1,yn-2,…;k步迭代要用到y(tǒng)n-1,yn-2,…,yn-k+1.

顯式單步迭代:

隱式單步迭代:

(2)第4頁(yè)，課件共41頁(yè)，創(chuàng)作于2023年2月一、Euler方法及其改進(jìn)

將[a,b]n等分,記

微分法:

積分法:

積分項(xiàng)利用矩形公式計(jì)算

1.顯式Euler方法

(★)第5頁(yè)，課件共41頁(yè)，創(chuàng)作于2023年2月Taylor公式推導(dǎo):

Euler公式幾何意義:

P1P2Pk也稱折線法

P0xy第6頁(yè)，課件共41頁(yè)，創(chuàng)作于2023年2月2.梯形法

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

取

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

注:

當(dāng)f(x,y)關(guān)于y滿足Lipschitz條件且步長(zhǎng)h滿足

直至滿足:

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

第7頁(yè)，課件共41頁(yè)，創(chuàng)作于2023年2月時(shí),迭代格式

(☆)收斂

.

3.改進(jìn)的Euler方法

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

),則得改進(jìn)的Euler方法:或也稱為預(yù)估-校正法.有時(shí)為了方便,預(yù)估-校正格式也寫成下面形式:第8頁(yè)，課件共41頁(yè)，創(chuàng)作于2023年2月二、單步法的局部截?cái)嗾`差及精度

Def1:先假設(shè),再估計(jì)誤差這種誤差稱為單步迭代法在xk+1處的局部截?cái)嗾`差.Def2:若某種數(shù)值方法的局部截?cái)嗾`差為,則稱該數(shù)值方法的精度為P階的.注:通常情況下,P越大,h越小,則截?cái)嗾`差越小,數(shù)值方法越精確.設(shè)10.Euler方法是一階方法.第9頁(yè)，課件共41頁(yè)，創(chuàng)作于2023年2月所以Euler方法為一階方法.而20.梯形法是二階方法.Taylor展開

第10頁(yè)，課件共41頁(yè)，創(chuàng)作于2023年2月將代入上式,得而代入上式得:當(dāng)h充分小時(shí),若,則可選取h,使得第11頁(yè)，課件共41頁(yè)，創(chuàng)作于2023年2月故梯形法的精度為2.同樣可以證明改進(jìn)的Euler法也是二階方法.梯形法的局部截?cái)嗾`差為:從而第12頁(yè)，課件共41頁(yè)，創(chuàng)作于2023年2月例1:

取步長(zhǎng)h=2/10,2/20,2/30,2/40,分別用歐拉法、改進(jìn)的歐拉法和梯形法求解.解:記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當(dāng)h=2/10時(shí),n=10.由Euler公式可得:k01234yk+11.21.38241.5061.535041.46503k56789yk+11.328771.170771.021130.891690.783788第13頁(yè)，課件共41頁(yè)，創(chuàng)作于2023年2月(2).改進(jìn)的Euler法:

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

第14頁(yè)，課件共41頁(yè)，創(chuàng)作于2023年2月n10203040h0.20.10.06670.05誤差0.10590.05210.03420.0256Euler法誤差:改進(jìn)的Euler法誤差:n10203040h0.20.10.06670.05誤差0.01230.00260.00115.9612e-004第15頁(yè)，課件共41頁(yè)，創(chuàng)作于2023年2月預(yù)-校方法,h=0.2時(shí)誤差最大值:0.0123歐拉方法,h=0.2時(shí)誤差最大值:0.1059解析解:第16頁(yè)，課件共41頁(yè)，創(chuàng)作于2023年2月三、Runge-Kutta方法1、Taylor級(jí)數(shù)法設(shè)初值問題有解y(x),由Tayler公式得:令當(dāng)時(shí),有.此時(shí)①為p階Taylor方法.p=1時(shí)即為Euler公式.稱之為Taylor級(jí)數(shù)法.其中例2:取步長(zhǎng)h=0.1,用一階、二階和四階Taylor方法求解下列初值問題①第17頁(yè)，課件共41頁(yè)，創(chuàng)作于2023年2月解:(1)一階Taylor法k01234yk+11.11.2211.370081.557791.80046(2)二階Taylor法k01234yk+11.111.246891.421751.652631.97088第18頁(yè)，課件共41頁(yè)，創(chuàng)作于2023年2月(3)四階Taylor法k01234yk+11.11111.249961.428481.666441.99942第19頁(yè)，課件共41頁(yè)，創(chuàng)作于2023年2月記由得稱為[xk,xk+1]上的平均斜率.故2、Runge-Kutta方法只要對(duì)K*提供不同的算法,就會(huì)得出不同的計(jì)算公式.如取則得改進(jìn)的Euler公式,它是利用xk,xk+1兩點(diǎn)的斜率值K1,K2的算術(shù)平均值作為K*,精度比Euler法高.則得Euler公式;取第20頁(yè)，課件共41頁(yè)，創(chuàng)作于2023年2月Runge-Kutta法的基本思想:

設(shè)法在[xk,xk+1]內(nèi)多預(yù)報(bào)幾個(gè)點(diǎn)的斜率,再將它們的加權(quán)平均值作為平均斜率K*一般顯式Runge-Kutta公式為:其中為待定參數(shù),且.稱為r級(jí)Runge-Kutta方法計(jì)算公式.②注:

式中待定參數(shù)的確定:先將②式右端在(xk,yk)處展成h的冪級(jí)數(shù)(即將yk+1展成h的冪級(jí)數(shù));再將y(xk+1)作Taylor級(jí)數(shù)展開;最后比較兩式中hk(k=0,1,2,…)的系數(shù),以確定出所有待定參數(shù).第21頁(yè)，課件共41頁(yè)，創(chuàng)作于2023年2月即可得p個(gè)方程,從而確定出待定參數(shù).代入表達(dá)式即可得到計(jì)算公式.如果要求兩個(gè)表達(dá)式的前p+1項(xiàng)完全重合,即局部截?cái)嗾`差達(dá)到,則稱②式為p階r級(jí)的Runge-Kutta方法.常用的是r=2,3,4級(jí)的R-K方法,且適當(dāng)選取參數(shù)使得p=r

.如要求:Runge-Kutta方法的推導(dǎo)(以r=2為例):當(dāng)r=2時(shí)第22頁(yè)，課件共41頁(yè)，創(chuàng)作于2023年2月則記又第23頁(yè)，課件共41頁(yè)，創(chuàng)作于2023年2月(1)常用的二階Runge-Kutta方法:預(yù)估-校正算法(2)這是一個(gè)四個(gè)參數(shù)三個(gè)方程的非線性方程組.它有一個(gè)自由度.稱滿足上述方程組的一族公式為二級(jí)二階Runge-Kutta方法.為使局部截?cái)嗾`差為,比較上述兩式右端同次冪系數(shù),應(yīng)取第24頁(yè)，課件共41頁(yè)，創(chuàng)作于2023年2月注:二級(jí)Runge-Kutta方法的精度最高是二階的,不可能達(dá)到三階.要提高計(jì)算方法的階,就必須增加預(yù)報(bào)點(diǎn).常用的三階Runge-Kutta方法(r=3):

(1)Heun(休恩)方法

中間點(diǎn)方法

(3)三階Kutta方法

第25頁(yè)，課件共41頁(yè)，創(chuàng)作于2023年2月(1)三階Heun方法

標(biāo)準(zhǔn)(經(jīng)典)四階Runge-Kutta方法

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

第26頁(yè)，課件共41頁(yè)，創(chuàng)作于2023年2月(2)稱為Gill(吉爾)方法

注:從理論上講,可以構(gòu)造任意高階的計(jì)算方法.但事實(shí)上,精度的階數(shù)與預(yù)報(bào)點(diǎn)的個(gè)數(shù)之間并非等量關(guān)系.預(yù)報(bào)點(diǎn)的個(gè)數(shù)r123456789r≥10精度的階數(shù)123445667≤r-2一般情況下,四階Runge-Kutta方法已可滿足精度要求.第27頁(yè)，課件共41頁(yè)，創(chuàng)作于2023年2月例3:用經(jīng)典Runge-Kutta方法求解下列初值問題(取h=0.1)解:標(biāo)準(zhǔn)Runge-Kutta公式為:計(jì)算結(jié)果見下表.為比較在相同計(jì)算量條件下近似解的精度,表中列出了Euler法(h=0.025)和改進(jìn)的Euler法(h=0.05)在相應(yīng)節(jié)點(diǎn)上的計(jì)算結(jié)果.第28頁(yè)，課件共41頁(yè)，創(chuàng)作于2023年2月xiEuler法h=0.025改進(jìn)Euler法h=0.05經(jīng)典R-K法h=0.1準(zhǔn)確解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注:用表中每種方法計(jì)算yi都需要計(jì)算四次f的值,即它們的計(jì)算量基本相等.第29頁(yè)，課件共41頁(yè)，創(chuàng)作于2023年2月四、單步法的進(jìn)一步討論—收斂性、相容性與穩(wěn)定性注:由定義可知,數(shù)值方法的收斂性并不涉及計(jì)算過程的舍入誤差,只與方法的截?cái)嗾`差有關(guān).若格式收斂,則整體截?cái)嗾`差必趨于零.Def:

(整體截?cái)嗾`差)稱為某一數(shù)值方法在點(diǎn)xk處的整體截?cái)嗾`差.它不僅與xk有關(guān),也與xk-1,xk-2,…,x1,x0有關(guān).則稱該單步法收斂.Def:對(duì)滿足解存在唯一性條件的初值問題(1),如果一個(gè)顯式單步法(3)產(chǎn)生的近似解對(duì)于任一固定的,均有1.收斂性第30頁(yè)，課件共41頁(yè)，創(chuàng)作于2023年2月由于,且關(guān)于y滿足Lipschitz條件,得則存在常數(shù)c>0使得且單步法中函數(shù)關(guān)于y滿足Lipschitz條件,則定理1:若初值問題的一個(gè)單步法的局部截?cái)嗾`差為記證:由局部截?cái)嗾`差的定義知第31頁(yè)，課件共41頁(yè)，創(chuàng)作于2023年2月故從而有故若y(x0)=y0,則e0=0,由不等式得第32頁(yè)，課件共41頁(yè)，創(chuàng)作于2023年2月設(shè)單步法為注:定理表明,數(shù)值方法的整體截?cái)嗾`差比局部截?cái)嗾`差低一階.收斂的方法至少是一階方法.在該定義條件下,Euler方法是一階的,預(yù)估-校正方法是二階.當(dāng)f(x,y)關(guān)于y也滿足Lipschitz條件,r級(jí)Runge-Kutta方法中的φ關(guān)于y也滿足Lipschitz條件,故定理中的條件得到滿足,解的收斂性得到保證.由于Rn,h→0(h→0),且xn為任意點(diǎn),故該式相當(dāng)于用近似方程當(dāng)x=xn+1固定時(shí),,所以有2.相容性第33頁(yè)，課件共41頁(yè)，創(chuàng)作于2023年2月通過在x=xn處求解近似方程而獲得原方程的近似解.因此,必須要求當(dāng)h→0時(shí),近似方程應(yīng)逼近于原方程.來代替因此,要使h→0時(shí),近似方程的極限狀態(tài)為原微分方程,需且只需下列極限成立:由于由于假設(shè)是連續(xù)函數(shù),故上式可表示為Def:如果當(dāng)h→0時(shí),近似方程逼近微分方程,則稱數(shù)值公式與原微分方程相容.相容的充要條件:

第34頁(yè)，課件共41頁(yè)，創(chuàng)作于2023年2月事實(shí)上:Remark:可以證明若單步法的階大于或等于1,則單步法與微分方程相容;反之,如果單步法與微分方程相容,且關(guān)于y滿足Lipschitz條件,則單步法至少為一階方法.(h→0)(1)若單步法的階大于或等于1,由知即單步法與微分方程相容.故有(2)如果單步法與微分方程相容,且關(guān)于y滿足Lipschitz條件,則第35頁(yè)，課件共41頁(yè)，創(chuàng)作于2023年2月關(guān)于單步法的收斂性以及收斂性定理都是在計(jì)算過程中無任何舍入誤差的前提條件下建立的,但在實(shí)際計(jì)算時(shí)通常會(huì)有舍入誤差及其積累,數(shù)值求解微分方程的過程是一個(gè)遞推公式,必須考即與微分方程相容的單步法至少為一階方法.Remark:在定理?xiàng)l件下,Euler方法、預(yù)估-校正方法以及Runge-Kutta方法都與原微分方程相容.中連續(xù),且關(guān)于變量y滿足Lipschitz條件,則單步法收斂的充要條件為相容性條件成立.Th1.設(shè)增量函數(shù)在區(qū)域3.穩(wěn)定性第36頁(yè)，課件共41頁(yè)，創(chuàng)作于2023年2月如果數(shù)值方法在計(jì)算過程中舍入誤差的積累越來越大,得不到有效控制,則稱其是不穩(wěn)定的;反之如果計(jì)算結(jié)果對(duì)初始數(shù)據(jù)的誤差及計(jì)算過程中的誤差不敏感,即舍入誤差不增長(zhǎng),則稱相應(yīng)的算法是穩(wěn)定的.數(shù)值方法的穩(wěn)定性有各種定義,這里僅考慮絕對(duì)穩(wěn)定性概念.慮誤差積累能否得到控制.Remark:從上面的定義可知,單步法是絕對(duì)穩(wěn)定的,與模型方程設(shè)某數(shù)值方法在節(jié)點(diǎn)xn處對(duì)初值問題的數(shù)值解為yn,實(shí)際計(jì)算得到的近似解為,稱為第n步數(shù)值解得擾動(dòng).=-Def:若某種數(shù)值方法在計(jì)算yn

時(shí)有擾動(dòng)但在計(jì)算后面的ym(m>n)由引起的誤差滿足