第9章 微分方程初值問(wèn)題的數(shù)值解法-1

1、第九章 微分方程初值問(wèn)題的數(shù)值解法,內(nèi)容提綱,引言 Euler法及其改進(jìn) Runge-Kutta方法 線性多步法 誤差分析 數(shù)值解法的收斂性、相容性和穩(wěn)定性 邊值問(wèn)題數(shù)值解法簡(jiǎn)介,引言,初值問(wèn)題的數(shù)值解法:求初值問(wèn)題的解在一系列節(jié)點(diǎn)的值 y ( xn )的近似值 yn 的方法.本章數(shù)值解法的特點(diǎn):都是采用“步進(jìn)式”,即求解過(guò)程順著節(jié)點(diǎn)排列的次序一步步向前推進(jìn).,基本知識(shí):,(1),定理1: 如果函數(shù) f (x , y)在區(qū)域 上連續(xù),且關(guān)于 y 滿(mǎn)足Lipschitz條件,常微分方程初值問(wèn)題:,求未知函數(shù) y= y (x) .,此時(shí)Lipschitz條件顯然成立. 故常用 在D上連續(xù)有界來(lái)代替

2、 f (x , y)關(guān)于 y 滿(mǎn)足Lipschitz條件.,注: 如無(wú)特別說(shuō)明,總假設(shè)(1)的解存在唯一且足夠光滑. 在f (x , y)對(duì)變量 y 可微的情形下, 若偏導(dǎo)數(shù) 連續(xù)有界, 則可取L為,除了要保證(1)有唯一解外,還需保證微分方程本身是穩(wěn)定的,即(1)的解連續(xù)依賴(lài)于初始值和函數(shù) f (x , y). 也就是說(shuō), 當(dāng)初始值 y0 及函數(shù) f (x , y)有微小變化時(shí), 只能引起解的微小變化.,(其中L 稱(chēng)為L(zhǎng)ipschitz常數(shù)),則對(duì)任何 ,初值問(wèn)題(1)在a ,b上存在唯一連續(xù)可微解 y = y (x).,定理2: 如果函數(shù) f (x , y)在區(qū)域 上關(guān)于 y 滿(mǎn)足Lips

3、chitz條件, 則(1)是穩(wěn)定的.,單步迭代: 計(jì)算 yn+1時(shí)僅用 yn ;,初值問(wèn)題(1)與下列積分方程的解等價(jià):,初值問(wèn)題的數(shù)值解就是求一系列節(jié)點(diǎn) 上函數(shù) y= y (x)的近似值 . 稱(chēng)為步長(zhǎng). 一般取等步長(zhǎng) h .,多步迭代: 計(jì)算 yn+1時(shí)除用 yn 外, 還要用到 yn-1, yn-2 ,; k 步迭代要用到 yn-1, yn-2 , yn-k+1 .,顯式單步迭代:,隱式單步迭代:,(2),一、Euler方法及其改進(jìn),將a , b n 等分, 記,微分法:,積分法:,積分項(xiàng)利用矩形公式計(jì)算,1. 顯式Euler方法,(),Taylor公式推導(dǎo):,Euler公式幾何意義:,2

4、. 梯形法,稱(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)似地,可得,(),時(shí),迭代格式 () 收斂 .,3. 改進(jìn)的Euler方法,把Euler法作為預(yù)報(bào)(稱(chēng)為預(yù)估公式),把隱式的梯形公式作為校正(稱(chēng)為校正公式 ),則得改進(jìn)的Euler方法:,或,也稱(chēng)為預(yù)估-校正法.,有時(shí)為了方便,預(yù)估-校正格式也寫(xiě)成下面形式:,二、單步法的局部截?cái)嗾`差及精度,Def 1: 先

5、假設(shè) ,再估計(jì)誤差 這種誤差稱(chēng)為單步迭代法在 xk+1處的局部截?cái)嗾`差.,Def 2: 若某種數(shù)值方法的局部截?cái)嗾`差為 ,則稱(chēng)該數(shù)值方法的精度為P 階的.,注: 通常情況下,P 越大, h 越小,則截?cái)嗾`差越小,數(shù)值方法越精確.,設(shè),10 .Euler方法是一階方法.,所以Euler方法為一階方法.,而,20 . 梯形法是二階方法.,Taylor展開(kāi),將 代入上式,得,而,代入上式得:,當(dāng)h充分小時(shí), 若 , 則可選取 h , 使得,故梯形法的精度為2 .,同樣可以證明改進(jìn)的Euler法也是二階方法.,梯形法的局部截?cái)嗾`差為:,從而,例1:,取步長(zhǎng) h = 2/10, 2/20, 2/30,

6、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公式可得:,(2) . 改進(jìn)的Euler法:,(3) . 梯形法(計(jì)算過(guò)程略),Euler法誤差:,改進(jìn)的Euler法誤差:,預(yù)-校方法, h=0.2時(shí) 誤差最大值: 0.0123,歐拉方法, h=0.2時(shí) 誤差最大值: 0.1059,解析解:,三、Runge -Kut

7、ta 方法,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)題,解: (1) 一階Taylor法,(2) 二階Taylor法,(3) 四階Taylor法,記,由,得,稱(chēng)為xk , xk+1上的平均斜率. 故,2、Runge -Kutta方法,只要對(duì)K*提供不同的算法, 就會(huì)得出不同的計(jì)算公式. 如取,則得改進(jìn)的Euler公式, 它是利用xk , xk+1兩點(diǎn)的斜

8、率值K1 , K2 的算術(shù)平均值作為K*, 精度比Euler法高.,則得Euler公式; 取,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ì)算公式.,注: 式中待定參數(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ù).,即可得 p 個(gè)方程,

9、從而確定出待定參數(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 .,如要求:,Runge-Kutta方法的推導(dǎo)(以r =2為例):,當(dāng)r =2 時(shí),則,記,又,(1),常用的二階Runge-Kutta方法:,預(yù)估-校正算法,(2),這是一個(gè)四個(gè)參數(shù)三個(gè)方程的非線性方程組. 它有一個(gè)自由度. 稱(chēng)滿(mǎn)足上述方程組的一族公式為二級(jí)二階Runge-Kutta方法.,為使局部截?cái)嗾`差為 ,比較上述兩式右端同次冪系數(shù),應(yīng)

10、取,注: 二級(jí)Runge-Kutta方法的精度最高是二階的,不可能達(dá)到三階. 要提高計(jì)算方法的階, 就必須增加預(yù)報(bào)點(diǎn).,常用的三階Runge-Kutta方法(r =3):,(1),Heun (休恩)方法,中間點(diǎn)方法,(3),三階Kutta方法,(1),三階Heun方法,標(biāo)準(zhǔn)(經(jīng)典)四階Runge-Kutta方法,(2),常用的四階Runge-Kutta方法(r =4):,(2),稱(chēng)為Gill(吉爾)方法,注: 從理論上講,可以構(gòu)造任意高階的計(jì)算方法. 但事實(shí)上, 精度的階數(shù)與預(yù)報(bào)點(diǎn)的個(gè)數(shù)之間并非等量關(guān)系.,一般情況下, 四階Runge-Kutta方法已可滿(mǎn)足精度要求.,例3: 用經(jīng)典Runge

11、-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é)果.,注: 用表中每種方法計(jì)算 yi 都需要計(jì)算四次 f 的值, 即它們的計(jì)算量基本相等.,四、單步法的進(jìn)一步討論收斂性、相容性與穩(wěn)定性,注: 由定義可知, 數(shù)值方法的收斂性并不涉及計(jì)算過(guò)程的舍入誤差, 只與方法的截?cái)嗾`差有關(guān). 若格式收斂, 則整體截?cái)嗾`差必趨于零.,1. 收斂性,由于 , 且 關(guān)于 y 滿(mǎn)足Lipschitz條件, 得,則存

12、在常數(shù) c 0 使得,且單步法中函數(shù) 關(guān)于 y 滿(mǎn)足Lipschitz條件, 則,定理1: 若初值問(wèn)題的一個(gè)單步法的局部截?cái)嗾`差為,記,證: 由局部截?cái)嗾`差的定義知,故,從而有,故,若 y(x0) = y0, 則e0=0, 由不等式,得,設(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)足, 解的收斂性得到保證.,

13、由于R n, h0(h0), 且 xn為任意點(diǎn), 故該式相當(dāng)于用近似方程,當(dāng)x = xn+1固定時(shí), , 所以有,2. 相容性,通過(guò)在 x = xn 處求解近似方程而獲得原方程的近似解. 因此,必須要求當(dāng)h0 時(shí), 近似方程應(yīng)逼近于原方程.,來(lái)代替,因此, 要使 h0 時(shí), 近似方程的極限狀態(tài)為原微分方程, 需且只需下列極限成立:,由于,由于假設(shè) 是連續(xù)函數(shù), 故上式可表示為,相容的充要條件:,事實(shí)上:,Remark: 可以證明若單步法的階大于或等于1, 則單步法與微分方程相容; 反之, 如果單步法與微分方程相容, 且 關(guān)于 y 滿(mǎn)足Lipschitz條件, 則單步法至少為一階方法.,(h0)

14、,(1) 若單步法的階大于或等于1, 由,知,即單步法與微分方程相容.,故有,(2) 如果單步法與微分方程相容, 且 關(guān)于y 滿(mǎn)足Lipschitz條件, 則,關(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è)

15、增量函數(shù) 在區(qū)域,3. 穩(wěn)定性,如果數(shù)值方法在計(jì)算過(guò)程中舍入誤差的積累越來(lái)越大, 得不到有效控制, 則稱(chēng)其是不穩(wěn)定的; 反之如果計(jì)算結(jié)果對(duì)初始數(shù)據(jù)的誤差及計(jì)算過(guò)程中的誤差不敏感, 即舍入誤差不增長(zhǎng), 則稱(chēng)相應(yīng)的算法是穩(wěn)定的. 數(shù)值方法的穩(wěn)定性有各種定義, 這里僅考慮絕對(duì)穩(wěn)定性概念.,慮誤差積累能否得到控制.,Remark: 從上面的定義可知, 單步法是絕對(duì)穩(wěn)定的, 與模型方程,Def : 若某種數(shù)值方法在計(jì)算 yn 時(shí)有擾動(dòng) 但在計(jì)算后面的 ym (m n)由 引起的誤差 滿(mǎn)足 則稱(chēng)該數(shù)值方法是絕對(duì)穩(wěn)定的.,設(shè) f (x , y)關(guān)于 y 滿(mǎn)足Lipschitz條件, 下面僅對(duì)典型方程(模型方程) 進(jìn)行討論. (其中為復(fù)常數(shù), 且 ),(1) Euler顯式公式:,幾個(gè)常用公式的穩(wěn)定性,要求誤差不增加, 即 必須,中的復(fù)數(shù)以及所用步長(zhǎng) h 有關(guān). 若對(duì)復(fù)平面上的某個(gè)區(qū)域G, 當(dāng) 時(shí)單步法絕對(duì)收斂, 則稱(chēng)G為單步法的絕對(duì)收斂域, G與實(shí)軸的交集稱(chēng)為絕對(duì)穩(wěn)定區(qū)間. 顯然絕對(duì)穩(wěn)定域越大, 數(shù)值方法的絕對(duì)穩(wěn)定性越好.,將Euler顯式公式用于模型方程 , 則有,設(shè) yn 有誤差 , 參與運(yùn)算的量為 由此引起的 yn+1 有誤差 ,則實(shí)際得