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

微分方程初值問(wèn)題的數(shù)值解法第一頁(yè)，共四十九頁(yè)，2022年，8月28日

包含自變量、未知函數(shù)及未知函數(shù)的導(dǎo)數(shù)或微分的方程稱(chēng)為微分方程。在微分方程中,自變量的個(gè)數(shù)只有一個(gè),稱(chēng)為常微分方程。自變量的個(gè)數(shù)為兩個(gè)或兩個(gè)以上的微分方程叫偏微分方程。微分方程中出現(xiàn)的未知函數(shù)最高階導(dǎo)數(shù)的階數(shù)稱(chēng)為微分方程的階數(shù)。如果未知函數(shù)y及其各階導(dǎo)數(shù)8.1

引言

微分方程知識(shí)回顧都是一次的,則稱(chēng)它是線性的,否則稱(chēng)為非線性的。

第二頁(yè)，共四十九頁(yè)，2022年，8月28日

在高等數(shù)學(xué)中，對(duì)于常微分方程的求解，給出了一些典型方程求解析解的基本方法：

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

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

這個(gè)一階微分方程就不能用初等函數(shù)及其積分來(lái)表達(dá)它的解。

微分方程知識(shí)回顧

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

從實(shí)際問(wèn)題當(dāng)中歸納出來(lái)的微分方程，通常主要依靠數(shù)值解法來(lái)解決。本章主要討論一階常微分方程初值問(wèn)題

(8.1)

在區(qū)間

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

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

上的數(shù)值解法。第四頁(yè)，共四十九頁(yè)，2022年，8月28日

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

相鄰兩個(gè)節(jié)點(diǎn)的間距稱(chēng)為步長(zhǎng)，本章總是假定h為定數(shù)，稱(chēng)為定步長(zhǎng)，這時(shí)節(jié)點(diǎn)可表示為

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

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

數(shù)值方法的基本思想第六頁(yè)，共四十九頁(yè)，2022年，8月28日

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

數(shù)值方法的基本思想第七頁(yè)，共四十九頁(yè)，2022年，8月28日一、Euler方法及其改進(jìn)

將[a,b]n等分,記

微分法:

積分法:

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

1.顯式Euler方法

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

第九頁(yè)，共四十九頁(yè)，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)

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

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

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

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

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

取

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

注:

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

直至滿(mǎn)足:

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

第十三頁(yè)，共四十九頁(yè)，2022年，8月28日時(shí),迭代格式

(☆)收斂

.

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

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

),則得改進(jìn)的Euler方法:或也稱(chēng)為預(yù)估-校正法.第十四頁(yè)，共四十九頁(yè)，2022年，8月28日有時(shí)為了方便,預(yù)估-校正格式也寫(xiě)成下面形式:第十五頁(yè)，共四十九頁(yè)，2022年，8月28日

改進(jìn)的歐拉公式比歐拉公式精度高的原因是：改進(jìn)歐拉公式用梯形面積代替曲邊梯形面積，而歐拉公式用矩形面積代替曲邊梯形面積。數(shù)值積分的梯形公式比矩形公式的精度高。第十六頁(yè)，共四十九頁(yè)，2022年，8月28日二、單步法的局部截?cái)嗾`差及精度

Def1:

先假設(shè),再估計(jì)誤差這種誤差稱(chēng)為單步迭代法在xk+1處的局部截?cái)嗾`差.Def2:

若某種數(shù)值方法的局部截?cái)嗾`差為,則稱(chēng)該數(shù)值方法的精度為P階的.注:

通常情況下,P越大,h越小,則截?cái)嗾`差越小,數(shù)值方法越精確.第十七頁(yè)，共四十九頁(yè)，2022年，8月28日所以Euler方法為一階方法.而設(shè)10.Euler方法是一階方法.第十八頁(yè)，共四十九頁(yè)，2022年，8月28日20.梯形法是二階方法.Taylor展開(kāi)

第十九頁(yè)，共四十九頁(yè)，2022年，8月28日將代入上式,得而代入上式得:當(dāng)h充分小時(shí),若,則可選取h,使得第二十頁(yè)，共四十九頁(yè)，2022年，8月28日故梯形法的精度為2.同樣可以證明改進(jìn)的Euler法也是二階方法.梯形法的局部截?cái)嗾`差為:從而第二十一頁(yè)，共四十九頁(yè)，2022年，8月28日例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第二十二頁(yè)，共四十九頁(yè)，2022年，8月28日(2).改進(jìn)的Euler法:

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

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

0.10590.05210.03420.0256Euler法誤差:改進(jìn)的Euler法誤差:n10203040h0.20.10.06670.05誤差

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

設(shè)初值問(wèn)題有解y(x),由Tayler公式得:令當(dāng)時(shí),有.此時(shí)①為p階Taylor方法.p=1時(shí)即為Euler公式.稱(chēng)之為T(mén)aylor級(jí)數(shù)法.其中例2:

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

設(shè)法在[xk,xk+1]內(nèi)多預(yù)報(bào)幾個(gè)點(diǎn)的斜率,再將它們的加權(quán)平均值作為平均斜率K*一般顯式Runge-Kutta公式為:其中為待定參數(shù),且.稱(chēng)為r級(jí)Runge-Kutta方法計(jì)算公式.②第三十頁(yè)，共四十九頁(yè)，2022年，8月28日即可得p個(gè)方程,從而確定出待定參數(shù).代入表達(dá)式即可得到計(jì)算公式.如果要求兩個(gè)表達(dá)式的前p+1項(xiàng)完全重合,即局部截?cái)嗾`差達(dá)到,則稱(chēng)②式為p階r級(jí)的Runge-Kutta方法.常用的是r=2,3,4

級(jí)的R-K方法,且適當(dāng)選取參數(shù)使得p=r

.如要求:注:

式中待定參數(shù)的確定:

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

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

二級(jí)Runge-Kutta方法的精度最高是二階的,不可能達(dá)到三階.要提高計(jì)算方法的階,就必須增加預(yù)報(bào)點(diǎn).常用的三階Runge-Kutta方法(r=3):

(1)Heun(休恩)方法

(3)三階Kutta方法

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

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

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

第三十七頁(yè)，共四十九頁(yè)，2022年，8月28日(2)稱(chēng)為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方法已可滿(mǎn)足精度要求.第三十八頁(yè)，共四十九頁(yè)，2022年，8月28日例3:

用經(jīng)典Runge-Kutta方法求解下列初值問(wèn)題(取h=0.1)解:標(biāo)準(zhǔn)Runge-Kutta公式為:計(jì)算結(jié)果見(jiàn)下表.為比較在相同計(jì)算量條件下近似解的精度,表中列出了Euler法(h=0.025)和改進(jìn)的Euler法(h=0.05)在相應(yīng)節(jié)點(diǎn)上的計(jì)算結(jié)果.第三十九頁(yè)，共四十九頁(yè)，2022年，8月28日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ì)算量基本相等.第四十頁(yè)，共四十九頁(yè)，2022年，8月28日四、單步法的進(jìn)一步討論—收斂性、相容性與穩(wěn)定性注:

由定義可知,數(shù)值方法的收斂性并不涉及計(jì)算過(guò)程的舍入誤差,只與方法的截?cái)嗾`差有關(guān).若格式收斂,則整體截?cái)嗾`差必趨于零.Def:

(整體截?cái)嗾`差)

稱(chēng)為某一數(shù)值方法在點(diǎn)xk處的整體截?cái)嗾`差.它不僅與xk有關(guān),也與xk-1,xk-2,…,x1,x0有關(guān).則稱(chēng)該單步法收斂.Def:

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

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

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

時(shí),近似方程應(yīng)逼近于原方程.來(lái)代替因此,要使h→0

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

如果當(dāng)h→0時(shí),

近似方程逼近微分方程,則稱(chēng)數(shù)值公式與原微分方程相容.相容的充要條件:

第四十五頁(yè)，共四十九頁(yè)，2022年，8月28日事實(shí)上:Remark:

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

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

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

設(shè)增量函數(shù)