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

1、關(guān)于微分方程初值問(wèn)題的數(shù)值解法第1頁(yè)，共52頁(yè)，2022年，5月20日，14點(diǎn)53分，星期三 包含自變量、未知函數(shù)及未知函數(shù)的導(dǎo)數(shù)或微分的方程稱為微分方程。在微分方程中, 自變量的個(gè)數(shù)只有一個(gè), 稱為常微分方程。自變量的個(gè)數(shù)為兩個(gè)或兩個(gè)以上的微分方程叫偏微分方程。微分方程中出現(xiàn)的未知函數(shù)最高階導(dǎo)數(shù)的階數(shù)稱為微分方程的階數(shù)。如果未知函數(shù)y及其各階導(dǎo)數(shù)8.1 引 言8.1.1 微分方程知識(shí)回顧都是一次的,則稱它是線性的,否則稱為非線性的。 第2頁(yè)，共52頁(yè)，2022年，5月20日，14點(diǎn)53分，星期三 在高等數(shù)學(xué)中，對(duì)于常微分方程的求解，給出了一些典型方程求解析解的基本方法： 一階：可分離變量法、

2、齊次方程、一階線性方程以及伯努利方程 高階：可降階方程、 常系數(shù)齊次線性方程的解法、常系數(shù)非齊次線性方程的解法 但能求解的常微分方程仍然是有限的，大多數(shù)的常微分方程是求不出解析解的。 這個(gè)一階微分方程就不能用初等函數(shù)及其積分來(lái)表達(dá)它的解。 8.1.1 微分方程知識(shí)回顧 例如 第3頁(yè)，共52頁(yè)，2022年，5月20日，14點(diǎn)53分，星期三 從實(shí)際問(wèn)題當(dāng)中歸納出來(lái)的微分方程，通常主要依靠數(shù)值解法來(lái)解決。本章主要討論一階常微分方程初值問(wèn)題 ( 8.1 ) 在區(qū)間 可以證明,如果函數(shù)在帶形區(qū)域 R=axb,-y內(nèi)連續(xù)，且關(guān)于y滿足李普希茲(Lipschitz)條件，即存在常數(shù)L(它與x,y無(wú)關(guān))使 對(duì)

3、R內(nèi)任意兩個(gè) 都成立,則方程( 8.1 )的解 在a, b上存在且唯一。 上的數(shù)值解法。 第4頁(yè)，共52頁(yè)，2022年，5月20日，14點(diǎn)53分，星期三 常微分方程初值問(wèn)題(8.1)式的數(shù)值解法，首先要算出精確解y(x)在區(qū)間a,b上的一系列離散節(jié)點(diǎn) 處的函數(shù)值 的近似值: 相鄰兩個(gè)節(jié)點(diǎn)的間距 稱為步長(zhǎng)，本章總是假定h為定數(shù)，稱為定步長(zhǎng)，這時(shí)節(jié)點(diǎn)可表示為8.1.2 數(shù)值方法的基本思想 1、數(shù)值解法需要把連續(xù)性的問(wèn)題加以離散化，從而求出離散節(jié)點(diǎn)處的數(shù)值解。第5頁(yè)，共52頁(yè)，2022年，5月20日，14點(diǎn)53分，星期三 描述這類算法，要求給出用已知信息 計(jì)算 的遞推公式。 建立這類遞推公式的基本方

4、法是在這些節(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)。8.1.2 數(shù)值方法的基本思想第6頁(yè)，共52頁(yè)，2022年，5月20日，14點(diǎn)53分，星期三 遞推公式通常有兩類， 一類是計(jì)算yi+1時(shí)只用到xi+1, xi 和yi，即前一步的值，因此有了初值以后就可以逐步往下計(jì)算，此類方法稱為單步法，其代表是龍格庫(kù)塔法。 另一類是計(jì)算yi+1時(shí)，除用到xi+1,xi和yi以外，還要用到 ，即前面k步的值，此類方法稱為多步法，其代表是亞當(dāng)斯法。 8.1.2

5、數(shù)值方法的基本思想第7頁(yè)，共52頁(yè)，2022年，5月20日，14點(diǎn)53分，星期三一、Euler方法及其改進(jìn) 將a , b n 等分, 記 微分法: 積分法: 積分項(xiàng)利用矩形公式計(jì)算 1. 顯式Euler方法 ()第8頁(yè)，共52頁(yè)，2022年，5月20日，14點(diǎn)53分，星期三Taylor公式推導(dǎo): 第9頁(yè)，共52頁(yè)，2022年，5月20日，14點(diǎn)53分，星期三Oyxy=y(x)(x1,y1)p1p0 x0 x1x2xixi+1xn-1xnpipi+1pn-1pnp1p2pipi+1pn-1pn切線p0p1的斜率為f (x0,y0)p2(x2,y2) 歐拉公式的幾何意義：第10頁(yè)，共52頁(yè)，202

6、2年，5月20日，14點(diǎn)53分，星期三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,取第11頁(yè)，共52頁(yè)，2022年，5月20日，14點(diǎn)53分，星期三從圖形上看,就獲得了一條近似于曲線y=y(x)的折線通常取 (常數(shù)),則Euler法的計(jì)算格式 i=0,1,n ( 7.2 ) 第12頁(yè)，共52頁(yè)，2022年，5月20日，1

7、4點(diǎn)53分，星期三2. 梯形法 稱之為梯形公式.這是一個(gè)隱式公式,通常用迭代法求解.具體做法: 取 先用Euler法求出初值 ,即 ,將其代入梯形公式的右端,使之轉(zhuǎn)化為顯式公式,即 注: 當(dāng) f (x , y)關(guān)于y滿足Lipschitz條件且步長(zhǎng)h 滿足 直至滿足: 若采用梯形公式計(jì)算()中的積分項(xiàng),則有類似地,可得() 第13頁(yè)，共52頁(yè)，2022年，5月20日，14點(diǎn)53分，星期三時(shí),迭代格式 () 收斂 . 3. 改進(jìn)的Euler方法 把Euler法作為預(yù)報(bào)(稱為預(yù)估公式),把隱式的梯形公式作為校正(稱為校正公式 ),則得改進(jìn)的Euler方法:或也稱為預(yù)估-校正法.第14頁(yè)，共52頁(yè)，

8、2022年，5月20日，14點(diǎn)53分，星期三有時(shí)為了方便,預(yù)估-校正格式也寫成下面形式:第15頁(yè)，共52頁(yè)，2022年，5月20日，14點(diǎn)53分，星期三 改進(jìn)的歐拉公式比歐拉公式精度高的原因是：改進(jìn)歐拉公式用梯形面積代替曲邊梯形面積，而歐拉公式用矩形面積代替曲邊梯形面積。 數(shù)值積分的梯形公式比矩形公式的精度高。第16頁(yè)，共52頁(yè)，2022年，5月20日，14點(diǎn)53分，星期三二、單步法的局部截?cái)嗾`差及精度 Def 1: 先假設(shè) ,再估計(jì)誤差這種誤差稱為單步迭代法在 xk+1處的局部截?cái)嗾`差.Def 2: 若某種數(shù)值方法的局部截?cái)嗾`差為 ,則稱該數(shù)值方法的精度為P 階的.注: 通常情況下,P 越大

9、, h 越小,則截?cái)嗾`差越小,數(shù)值方法越精確. 第17頁(yè)，共52頁(yè)，2022年，5月20日，14點(diǎn)53分，星期三所以Euler方法為一階方法. 而 設(shè) 10 .Euler方法是一階方法.第18頁(yè)，共52頁(yè)，2022年，5月20日，14點(diǎn)53分，星期三20 . 梯形法是二階方法.Taylor展開(kāi) 第19頁(yè)，共52頁(yè)，2022年，5月20日，14點(diǎn)53分，星期三將 代入上式,得 而代入上式得:當(dāng)h充分小時(shí), 若 , 則可選取 h , 使得第20頁(yè)，共52頁(yè)，2022年，5月20日，14點(diǎn)53分，星期三故梯形法的精度為2 . 同樣可以證明改進(jìn)的Euler法也是二階方法. 梯形法的局部截?cái)嗾`差為: 從

10、而第21頁(yè)，共52頁(yè)，2022年，5月20日，14點(diǎn)53分，星期三例1: 取步長(zhǎng) h = 2/10, 2/20, 2/30, 2/40, 分別用歐拉法、改進(jìn)的歐拉法和梯形法求解.解: 記 f (x, y) = y x y2, xk= k h (k = 0, 1, 2, n )(1) . Euler法: yk+1 = yk + h( yk xk yk2 ) (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.02113

11、0.891690.783788第22頁(yè)，共52頁(yè)，2022年，5月20日，14點(diǎn)53分，星期三(2) . 改進(jìn)的Euler法: k01234yk+11.19121.343841.423481.419051.3473k56789yk+11.237261.114240.9941510.8847510.788666(3) . 梯形法(計(jì)算過(guò)程略) 第23頁(yè)，共52頁(yè)，2022年，5月20日，14點(diǎn)53分，星期三 n 10 20 30 40 h 0.2 0.1 0.0667 0.05誤差 0.1059 0.0521 0.0342 0.0256Euler法誤差:改進(jìn)的Euler法誤差: n 10 20

12、30 40 h 0.2 0.1 0.0667 0.05誤差 0.0123 0.0026 0.0011 5.9612e-004第24頁(yè)，共52頁(yè)，2022年，5月20日，14點(diǎn)53分，星期三預(yù)-校方法, h=0.2時(shí)誤差最大值: 0.0123歐拉方法, h=0.2時(shí)誤差最大值: 0.1059解析解:第25頁(yè)，共52頁(yè)，2022年，5月20日，14點(diǎn)53分，星期三三、Runge -Kutta 方法1、Taylor 級(jí)數(shù)法 設(shè)初值問(wèn)題 有解 y (x), 由Tayler公式得:令當(dāng) 時(shí), 有 . 此時(shí)為 p 階Taylor方法. p=1時(shí)即為Euler公式.稱之為Taylor級(jí)數(shù)法. 其中例2: 取

13、步長(zhǎng) h = 0.1, 用一階、二階和四階Taylor方法求解下列初值問(wèn)題第26頁(yè)，共52頁(yè)，2022年，5月20日，14點(diǎn)53分，星期三解: (1) 一階Taylor法k01234yk+11.11.2211.370081.557791.80046(2) 二階Taylor法k01234yk+11.111.246891.421751.652631.97088第27頁(yè)，共52頁(yè)，2022年，5月20日，14點(diǎn)53分，星期三(3) 四階Taylor法k01234yk+11.11111.249961.428481.666441.99942第28頁(yè)，共52頁(yè)，2022年，5月20日，14點(diǎn)53分，星期三

14、記由得稱為xk , xk+1上的平均斜率. 故2、Runge -Kutta方法只要對(duì)K*提供不同的算法, 就會(huì)得出不同的計(jì)算公式. 如取則得改進(jìn)的Euler公式, 它是利用xk , xk+1兩點(diǎn)的斜率值K1 , K2 的算術(shù)平均值作為K*, 精度比Euler法高.則得Euler公式; 取第29頁(yè)，共52頁(yè)，2022年，5月20日，14點(diǎn)53分，星期三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ì)算公式. 第30頁(yè)，

15、共52頁(yè)，2022年，5月20日，14點(diǎn)53分，星期三即可得 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 . 如要求: 注: 式中待定參數(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ù).第31頁(yè)，共52

16、頁(yè)，2022年，5月20日，14點(diǎn)53分，星期三Runge-Kutta方法的推導(dǎo)(以r =2為例): 當(dāng)r =2 時(shí)記第32頁(yè)，共52頁(yè)，2022年，5月20日，14點(diǎn)53分，星期三則又第33頁(yè)，共52頁(yè)，2022年，5月20日，14點(diǎn)53分，星期三這是一個(gè)四個(gè)參數(shù)三個(gè)方程的非線性方程組. 它有一個(gè)自由度. 稱滿足上述方程組的一族公式為二級(jí)二階Runge-Kutta方法.為使局部截?cái)嗾`差為 ,比較上述兩式右端同次冪系數(shù),應(yīng)取第34頁(yè)，共52頁(yè)，2022年，5月20日，14點(diǎn)53分，星期三(1) 常用的二階Runge-Kutta方法: 預(yù)估-校正算法(2)中間點(diǎn)方法 第35頁(yè)，共52頁(yè)，2022

17、年，5月20日，14點(diǎn)53分，星期三注: 二級(jí)Runge-Kutta方法的精度最高是二階的,不可能達(dá)到三階. 要提高計(jì)算方法的階, 就必須增加預(yù)報(bào)點(diǎn). 常用的三階Runge-Kutta方法(r =3): (1) Heun (休恩)方法 (3) 三階Kutta方法 第36頁(yè)，共52頁(yè)，2022年，5月20日，14點(diǎn)53分，星期三(1) 三階Heun方法 標(biāo)準(zhǔn)(經(jīng)典)四階Runge-Kutta方法 (2) 常用的四階Runge-Kutta方法(r =4): 第37頁(yè)，共52頁(yè)，2022年，5月20日，14點(diǎn)53分，星期三(2) 稱為Gill(吉爾)方法 注: 從理論上講,可以構(gòu)造任意高階的計(jì)算方法

18、. 但事實(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方法已可滿足精度要求. 第38頁(yè)，共52頁(yè)，2022年，5月20日，14點(diǎn)53分，星期三例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é)果. 第39頁(yè)，共52頁(yè)，2022年，5月20日，14

19、點(diǎn)53分，星期三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

20、.0112713.0744793.0766193.0766230.93.4976063.5761443.5788043.5788091.04.0551924.1515734.1548394.154845注: 用表中每種方法計(jì)算 yi 都需要計(jì)算四次 f 的值, 即它們的計(jì)算量基本相等. 第40頁(yè)，共52頁(yè)，2022年，5月20日，14點(diǎn)53分，星期三四、單步法的進(jìn)一步討論收斂性、相容性與穩(wěn)定性注: 由定義可知, 數(shù)值方法的收斂性并不涉及計(jì)算過(guò)程的舍入誤差, 只與方法的截?cái)嗾`差有關(guān). 若格式收斂, 則整體截?cái)嗾`差必趨于零. Def : (整體截?cái)嗾`差) 稱 為某一數(shù)值方法在點(diǎn) xk 處的整體截?cái)?/p>

21、誤差. 它不僅與 xk 有關(guān), 也與xk -1, xk -2 , x1, x0 有關(guān). 則稱該單步法收斂. Def : 對(duì)滿足解存在唯一性條件的初值問(wèn)題(1), 如果一個(gè)顯式單步法(3)產(chǎn)生的近似解對(duì)于任一固定的 , 均有 1. 收斂性第41頁(yè)，共52頁(yè)，2022年，5月20日，14點(diǎn)53分，星期三由于 , 且 關(guān)于 y 滿足Lipschitz條件, 得 則存在常數(shù) c 0 使得 且單步法中函數(shù) 關(guān)于 y 滿足Lipschitz條件, 則 定理1: 若初值問(wèn)題的一個(gè)單步法的局部截?cái)嗾`差為 記 證: 由局部截?cái)嗾`差的定義知第42頁(yè)，共52頁(yè)，2022年，5月20日，14點(diǎn)53分，星期三故 從而有

22、 故 若 y(x0) = y0, 則e0=0, 由不等式 得 第43頁(yè)，共52頁(yè)，2022年，5月20日，14點(diǎn)53分，星期三設(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條件, 故定理中的條件得到滿足, 解的收斂性得到保證.由于R n, h0(h0), 且 xn為任意點(diǎn), 故該式相當(dāng)于用近似方程 當(dāng)x = xn+1固定時(shí), , 所以有

23、2. 相容性第44頁(yè)，共52頁(yè)，2022年，5月20日，14點(diǎn)53分，星期三通過(guò)在 x = xn 處求解近似方程而獲得原方程的近似解. 因此,必須要求當(dāng)h0 時(shí), 近似方程應(yīng)逼近于原方程.來(lái)代替 因此, 要使 h0 時(shí), 近似方程的極限狀態(tài)為原微分方程, 需且只需下列極限成立:由于 由于假設(shè) 是連續(xù)函數(shù), 故上式可表示為 Def : 如果當(dāng)h0時(shí), 近似方程逼近微分方程, 則稱數(shù)值公式與原微分方程相容. 相容的充要條件: 第45頁(yè)，共52頁(yè)，2022年，5月20日，14點(diǎn)53分，星期三事實(shí)上: Remark: 可以證明若單步法的階大于或等于1, 則單步法與微分方程相容; 反之, 如果單步法與微

24、分方程相容, 且 關(guān)于 y 滿足Lipschitz條件, 則單步法至少為一階方法. (h0) (1) 若單步法的階大于或等于1, 由 知 即單步法與微分方程相容. 故有 (2) 如果單步法與微分方程相容, 且 關(guān)于y 滿足Lipschitz條件, 則第46頁(yè)，共52頁(yè)，2022年，5月20日，14點(diǎn)53分，星期三關(guān)于 單步法的收斂性以及收斂性定理都是在計(jì)算過(guò)程中無(wú)任何舍入誤差的前提條件下建立的, 但在實(shí)際計(jì)算時(shí)通常會(huì)有舍入誤差及其積累, 數(shù)值求解微分方程的過(guò)程是一個(gè)遞推公式, 必須考 即與微分方程相容的單步法至少為一階方法. Remark: 在定理?xiàng)l件下, Euler方法、預(yù)估-校正方法以及Runge-Kutta方法都與原微分方程相容.中連續(xù), 且關(guān)于變量y 滿足Lipschitz條件, 則單步法收斂的充要條件為相容性條件成立. Th1. 設(shè)增量函數(shù) 在區(qū)域 3. 穩(wěn)定性第47頁(yè)，共52頁(yè)，2022年，5月20日，14點(diǎn)53分，星期三 如果數(shù)值方法在計(jì)算過(guò)程中舍入誤差的積累越來(lái)越大, 得不到有效控制, 則稱其是不穩(wěn)定的; 反之如果計(jì)算結(jié)果對(duì)初始數(shù)據(jù)的誤差及計(jì)算過(guò)程中的誤差不敏感, 即舍入誤差不增長(zhǎng), 則稱相應(yīng)的算法是穩(wěn)定的. 數(shù)值方法的穩(wěn)定性有各種定義, 這里僅考慮