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

關(guān)于微分方程初值問題的數(shù)值解法

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

引言8.1.1

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

第2頁,共52頁，2024年2月25日，星期天

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

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

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

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

微分方程知識回顧

例如第3頁,共52頁，2024年2月25日，星期天

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

(8.1)

在區(qū)間

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

對R內(nèi)任意兩個都成立,則方程(8.1)的解在

a,b

上存在且唯一。

上的數(shù)值解法。第4頁,共52頁，2024年2月25日，星期天

常微分方程初值問題(8.1)式的數(shù)值解法，首先要算出精確解y(x)在區(qū)間

a,b

上的一系列離散節(jié)點(diǎn)處的函數(shù)值的近似值:

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

數(shù)值方法的基本思想1、數(shù)值解法需要把連續(xù)性的問題加以離散化，從而求出離散節(jié)點(diǎn)處的數(shù)值解。第5頁,共52頁，2024年2月25日，星期天

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

數(shù)值方法的基本思想第6頁,共52頁，2024年2月25日，星期天

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

數(shù)值方法的基本思想第7頁,共52頁，2024年2月25日，星期天一、Euler方法及其改進(jìn)

將[a,b]n等分,記

微分法:

積分法:

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

1.顯式Euler方法

(★)第8頁,共52頁，2024年2月25日，星期天Taylor公式推導(dǎo):

第9頁,共52頁，2024年2月25日，星期天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)

歐拉公式的幾何意義：第10頁,共52頁，2024年2月25日，星期天Euler法的求解過程是:

從初始點(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)時,得

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

i=0,1,…,n(7.2)第12頁,共52頁，2024年2月25日，星期天2.梯形法

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

取

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

注:

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

直至滿足:

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

第13頁,共52頁，2024年2月25日，星期天時,迭代格式

(☆)收斂

.

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

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

),則得改進(jìn)的Euler方法:或也稱為預(yù)估-校正法.第14頁,共52頁，2024年2月25日，星期天有時為了方便,預(yù)估-校正格式也寫成下面形式:第15頁,共52頁，2024年2月25日，星期天

改進(jìn)的歐拉公式比歐拉公式精度高的原因是：改進(jìn)歐拉公式用梯形面積代替曲邊梯形面積，而歐拉公式用矩形面積代替曲邊梯形面積。數(shù)值積分的梯形公式比矩形公式的精度高。第16頁,共52頁，2024年2月25日，星期天二、單步法的局部截?cái)嗾`差及精度

Def1:

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

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

通常情況下,P越大,h越小,則截?cái)嗾`差越小,數(shù)值方法越精確.第17頁,共52頁，2024年2月25日，星期天所以Euler方法為一階方法.而設(shè)10.Euler方法是一階方法.第18頁,共52頁，2024年2月25日，星期天20.梯形法是二階方法.Taylor展開

第19頁,共52頁，2024年2月25日，星期天將代入上式,得而代入上式得:當(dāng)h充分小時,若,則可選取h,使得第20頁,共52頁，2024年2月25日，星期天故梯形法的精度為2.同樣可以證明改進(jìn)的Euler法也是二階方法.梯形法的局部截?cái)嗾`差為:從而第21頁,共52頁，2024年2月25日，星期天例1:

取步長

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時,n=10.由Euler公式可得:k01234yk+11.21.38241.5061.535041.46503k56789yk+11.328771.170771.021130.891690.783788第22頁,共52頁，2024年2月25日，星期天(2).改進(jìn)的Euler法:

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

第23頁,共52頁，2024年2月25日，星期天n10203040h0.20.10.06670.05誤差

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

0.01230.00260.00115.9612e-004第24頁,共52頁，2024年2月25日，星期天預(yù)-校方法,h=0.2時誤差最大值:0.0123歐拉方法,h=0.2時誤差最大值:0.1059解析解:第25頁,共52頁，2024年2月25日，星期天三、Runge-Kutta方法1、Taylor級數(shù)法

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

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

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

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

.如要求:注:

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

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

第35頁,共52頁，2024年2月25日，星期天注:

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

(1)Heun(休恩)方法

(3)三階Kutta方法

第36頁,共52頁，2024年2月25日，星期天(1)三階Heun方法

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

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

第37頁,共52頁，2024年2月25日，星期天(2)稱為Gill(吉爾)方法

注:

從理論上講,可以構(gòu)造任意高階的計(jì)算方法.但事實(shí)上,精度的階數(shù)與預(yù)報(bào)點(diǎn)的個數(shù)之間并非等量關(guān)系.預(yù)報(bào)點(diǎn)的個數(shù)r123456789r≥10精度的階數(shù)123445667≤r-2一般情況下,四階Runge-Kutta方法已可滿足精度要求.第38頁,共52頁，2024年2月25日，星期天例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é)果.第39頁,共52頁，2024年2月25日，星期天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ì)算量基本相等.第40頁,共52頁，2024年2月25日，星期天四、單步法的進(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:

對滿足解存在唯一性條件的初值問題(1),如果一個顯式單步法(3)產(chǎn)生的近似解對于任一固定的,均有1.收斂性第41頁,共52頁，2024年2月25日，星期天由于,且關(guān)于y滿足Lipschitz條件,得則存在常數(shù)c>0使得且單步法中函數(shù)關(guān)于y滿足Lipschitz條件,則定理1:

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

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

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

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

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

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

第45頁,共52頁，2024年2月25日，星期天事實(shí)上:Remark:

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

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

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

設(shè)增量函數(shù)在區(qū)域3.穩(wěn)定性第47頁,共52頁，2024年2月25日，星期天

如果數(shù)值方法在計(jì)算過程中舍入誤差的積累越來越大,得不到有效控制,則稱其是不穩(wěn)定的;反之如果計(jì)算結(jié)果對初始數(shù)據(jù)的誤差及計(jì)算過程中的誤差不敏感,即舍入誤差不增長,則稱相應(yīng)的算法是穩(wěn)定的.數(shù)值方法的穩(wěn)定性有各種定義,這里僅考慮絕對穩(wěn)定性概念.慮誤差積累能否得到控制.Remark:從上面的定義可知,單步法是絕對穩(wěn)定的,與模型方程

設(shè)某數(shù)值方法在節(jié)點(diǎn)xn處對初值問題的數(shù)值解為yn,實(shí)際計(jì)算得到的近似解為,