第四章 常微分方程數(shù)值解法.ppt

1、1,四常微分方程數(shù)值解法,2,常微分方程數(shù)值解法,引言（常微分方程數(shù)值解法概述） 顯式歐拉法、隱式歐拉法、二步歐拉法 局部截?cái)嗾`差與精度 改進(jìn)的歐拉方法 龍格-庫(kù)塔方法 收斂性與穩(wěn)定性簡(jiǎn)述 一階常微分方程組與高階常微分方程,3,引言,一階常微分方程初值問(wèn)題：,定理：若 f (x, y) 在某閉區(qū)域 R ：,上連續(xù)，且在 R 域內(nèi)滿(mǎn)足李普希茲 (Lipschitz) 條件，即存在正數(shù) L，使得對(duì)于 R 域內(nèi)的任意兩值 y1, y2，下列不等式成立：,則上述初值問(wèn)題的連續(xù)可微的解 y(x) 存在并且唯一。,4,引言（續(xù)）,實(shí)際生產(chǎn)與科研中，除少數(shù)簡(jiǎn)單情況能獲得初值問(wèn)題的初等解（用初等函數(shù)表示的解）

2、外，絕大多數(shù)情況下是求不出初等解的。 有些初值問(wèn)題即便有初等解，也往往由于形式過(guò)于復(fù)雜而不便處理。 實(shí)用的方法是在計(jì)算機(jī)上進(jìn)行數(shù)值求解：即不直接求 y(x) 的顯式解，而是在解所存在的區(qū)間上，求得一系列點(diǎn) xn (n = 0, 1, 2, ) 上解的近似值。,5,歐拉（Euler）方法,方法一化導(dǎo)數(shù)為差商的方法,由于在逐步求解的過(guò)程中，y(xn) 的準(zhǔn)確值無(wú)法求解出來(lái)，因此用其近似值代替。 為避免混淆，以下學(xué)習(xí)簡(jiǎn)記：,y(xn)：待求函數(shù) y(x) 在 xn 處的精確函數(shù)值 yn ：待求函數(shù) y(x) 在 xn 處的近似函數(shù)值,6,代入初值問(wèn)題表達(dá)式可得：,根據(jù) y0 可以一步步計(jì)算出函數(shù) y

3、 = y(x) 在 x1, x2, x3 x4, 上的近似值 y1, y2, y3, y4 , 常微分方程數(shù)值解是一組離散的函數(shù)值數(shù)據(jù)，它的精確表達(dá)式很難求解得到，但可以進(jìn)行插值計(jì)算后用插值函數(shù)逼近 y(x),7,歐拉方法（續(xù)）,方法二數(shù)值積分法,同樣以近似值 yn 代替精確值 y(xn) 可得：,將微分方程 y = f (x, y) 在區(qū)間 xn, xn+1 上積分：,8,歐拉方法的幾何意義,x,y,0,9,隱式歐拉法,在數(shù)值積分法推導(dǎo)中，積分的近似值取為積分區(qū)間寬度與右端點(diǎn)處的函數(shù)值乘積，即：,這樣便得到了隱式歐拉法：,隱式歐拉法沒(méi)有顯式歐拉法方便,10,二步歐拉法,在數(shù)值積分法推導(dǎo)中，積

4、分區(qū)間寬度選為兩步步長(zhǎng)，即積分區(qū)間為：xn-1, xn+1，則：,以 y(x) 在 xn -1, xn 上的近似值代替精確值可得：,中矩形公式,11,梯形公式歐拉法,在數(shù)值積分法中，如果用梯形公式近似計(jì)算 f (x, y) 在區(qū)間 xn, xn+1 上的積分，即：,用近似值代替精確值可得梯形公式歐拉法：,上式右端出現(xiàn)了未知項(xiàng)，可見(jiàn)梯形法是隱式歐拉法的一種；實(shí)際上，梯形公式歐拉法是顯式歐拉法與隱式歐拉法的算術(shù)平均。,12,例,用顯式歐拉法、隱式歐拉法、梯形法求解初值問(wèn)題：,取 h = 0.1，計(jì)算到 x = 0.5，并與精確解進(jìn)行比較,解：由已知條件可得：h = 0.1，x0 = 0， y0 =

5、 1， f (x, y) = - y + x + 1,顯式歐拉法：,13,例：（續(xù)）,隱式歐拉法：,化簡(jiǎn)得：,梯形公式歐拉法：,14,計(jì)算結(jié)果：,本題的精確解為：,15,局部截?cái)嗾`差,為了簡(jiǎn)化分析某常微分方程數(shù)值算法的誤差，現(xiàn)假設(shè) yn = y(xn)，即在前一步 yn 準(zhǔn)確的前提下，估計(jì)：,稱(chēng)上述誤差 Tn+1 為該常微分方程數(shù)值算法的局部截?cái)嗾`差,如果某個(gè)常微分方程數(shù)值算法的局部截?cái)嗾`差可表示為 O(h p+1)，則稱(chēng)該數(shù)值算法的精度是 p 階,歐拉法的精度為一階；二步歐拉法的精度為二階；梯形公式歐拉法的精度為二階。,16,泰勒展開(kāi)法,如果初值問(wèn)題中的 f (x, y) 充分可微，則可將

6、y(xn+1) 在點(diǎn) xn 處展開(kāi)：,如果只保留線(xiàn)性項(xiàng)，忽略 h2 及以上各項(xiàng)，則：,顯式歐拉公式,17,局部截?cái)嗾`差的分析,利用泰勒公式展開(kāi)，比較各算法與展開(kāi)式的前幾項(xiàng),將 y(xn+1) 在 xn 點(diǎn)處用泰勒公式展開(kāi)：,顯式歐拉法的局部截?cái)嗾`差：,歐拉法,1 階精度,18,補(bǔ)充：二元函數(shù)微分中值定理,19,y(xn+1) 在 xn 點(diǎn)處展開(kāi)：,隱式歐拉法：,1 階精度,20,分別將 y(xn+1), y(xn-1) 在 xn 點(diǎn)處用泰勒公式展開(kāi)：,二步歐拉法的局部截?cái)嗾`差：,二步歐拉法：,2 階精度,21,梯形公式歐拉法：,y(xn+1) 在 xn 點(diǎn)處展開(kāi)：,2 階精度,22,各種歐拉法

7、的比較,23,改進(jìn)的歐拉法,從上述例子可以看到，梯形法由于具有二階精度，其局部截?cái)嗾`差比顯式歐拉法和隱式歐拉法小，但梯形法實(shí)質(zhì)上是一種隱式算法 顯式歐拉法是一個(gè)顯式算法，雖然計(jì)算量較小，但是精度不高 綜合兩種方法的長(zhǎng)處，可以先用顯式歐拉法求出 y(xn+1) 的一個(gè)粗略近似值，然后用它代入梯形法公式的右端，用梯形法計(jì)算 y(xn+1) 的較為精確的近似值。,24,改進(jìn)的歐拉法（續(xù)）,按照上述思想，可以建立如下預(yù)報(bào)-校正系統(tǒng)：,按以上兩式求解常微分方程的算法稱(chēng)為改進(jìn)的歐拉法，它還可以表示為：,嵌套形式,平均化形式,2 階精度,25,用改進(jìn)歐拉法求上例所述的初值問(wèn)題并與歐拉法和梯形法比較誤差的大小

8、。,解：采用改進(jìn)歐拉法的嵌套形式：,26,計(jì)算結(jié)果,可見(jiàn)，改進(jìn)歐拉法的誤差數(shù)量級(jí)與梯形法大致相同，而比歐拉法小得多。,27,改進(jìn)的歐拉法的意義,改進(jìn)的歐拉法的平均化形式,y (xn+1) 在點(diǎn) xn 處的一階展開(kāi)式為：,28,改進(jìn)的歐拉法的幾何意義,0,x,y,29,龍格-庫(kù)塔(Runge-Kutta)方法,改進(jìn)的歐拉法（2 階精度）,y (xn+1) 在點(diǎn) xn 處的一階泰勒展開(kāi)式為：,顯式歐拉法（1 階精度）,30,龍格-庫(kù)塔方法（續(xù)）,顯式歐拉法用一個(gè)點(diǎn)的值 k1 作為 k* 的近似值 改進(jìn)的歐拉公式用二個(gè)點(diǎn)的值 k1 和 k2 的平均值作為 k* 近似值； 改進(jìn)的歐拉法比顯式歐拉法精度

9、高； 在 xn, xn+1 內(nèi)多預(yù)報(bào)幾個(gè)點(diǎn)的 ki 值，并用其加權(quán)平均值作為 k* 的近似值從而構(gòu)造出具有更高精度的計(jì)算公式，這就是龍格-庫(kù)塔方法的基本思想。,31,二階龍格-庫(kù)塔方法,以 k1 和 k2 的加權(quán)平均來(lái)近似取代 k*,為分析局部截?cái)嗾`差，令 yn = y(xn)，由泰勒公式得：,32,補(bǔ)充：二元泰勒展開(kāi)式,33,用二元泰勒公式展開(kāi),將 k1, k2 代入 中可得：,34,二階龍格-庫(kù)塔方法（續(xù)）,2 階精度,35,四個(gè)未知變量，只有三個(gè)方程，有無(wú)窮多組解 每組解的構(gòu)成的龍格-庫(kù)塔方法均為二階,二階龍格-庫(kù)塔方法即為改進(jìn)的歐拉方法,變形的歐拉法 中 點(diǎn) 方 法,36,三階龍格-庫(kù)

10、塔方法,三階龍格-庫(kù)塔方法是用三個(gè)值 k1, k2, k3 的加權(quán)平均來(lái)近似取代 k*,要使三階龍格-庫(kù)塔方法具有三階精度，必須使其局部截?cái)嗾`差為 O(h4) 將 k1, k2, k3 代入 yn+1 的表達(dá)式中，在 (xn, yn) 處用二元泰勒公式展開(kāi)，與 y(xn+1) 在 xn 處的泰勒展開(kāi)式比較,37,三階龍格-庫(kù)塔方法（續(xù)）,類(lèi)似二階龍格-庫(kù)塔方法的推導(dǎo)過(guò)程，8 個(gè)待定系數(shù) c1, c2, c3, a2, a3, b21, b31, b32 應(yīng)滿(mǎn)足：,8 個(gè)未知參數(shù)，6 個(gè)方程，有無(wú)窮多組解,庫(kù)塔公式,38,四階龍格-庫(kù)塔方法,類(lèi)似可以推出四階龍格-庫(kù)塔公式，常用的有：,標(biāo)準(zhǔn)四階龍

11、格-庫(kù)塔公式,39,四階龍格-庫(kù)塔方法（續(xù)）,吉爾（Gill）公式,4 階以上龍格-庫(kù)塔公式的計(jì)算量太大，并且精度不一定提高，有時(shí)反而會(huì)降低，因此實(shí)際應(yīng)用中一般選用四階龍格-庫(kù)塔已足可滿(mǎn)足精度要求。,40,用經(jīng)典四階龍格-庫(kù)塔方法求解前例的初值問(wèn)題，并與改進(jìn) 歐拉 法、梯形法在 x5 = 0.5 處比較其誤差大小,解：采用經(jīng)典四階龍格-庫(kù)塔公式：,41,四階R-K方法的精度比二階方法高得多,精確解為：,R-K方法的誤差：,改進(jìn)歐拉法的誤差：,梯形法的誤差：,42,變步長(zhǎng)的龍格-庫(kù)塔方法,設(shè) y (xn) 在 xn 處的值 yn = y (xn)，當(dāng) xn+1 = xn+ h 時(shí) y (xn+1

12、) 的近似值為 ，由于四階 R-K 方法的精度為 4 階，故局部截?cái)嗾`差為：,用四階R-K方法求解初值問(wèn)題精度較高，但要從理論上給出誤差 | y (xn) - yn | 的估計(jì)式則比較困難；那么應(yīng)如何判斷計(jì)算結(jié)果的精度以及如何選擇合適的步長(zhǎng) h？ 通常是通過(guò)不同步長(zhǎng)在計(jì)算機(jī)上的計(jì)算結(jié)果進(jìn)行近似估計(jì)。,43,若以 h/2 為步長(zhǎng)，從 xn 出發(fā)，經(jīng)過(guò)兩步計(jì)算，得到 y(xn+1) 的近似值,變步長(zhǎng)的龍格-庫(kù)塔方法（續(xù)）,以上每步的截?cái)嗾`差約為 cn(h/2)5，于是兩步的局部截?cái)嗾`差為：,于是：,整理得：,44,變步長(zhǎng)的龍格-庫(kù)塔方法（續(xù)）,記： ，給定的精度要求為 e,D e，反復(fù)將步長(zhǎng)折半計(jì)

13、算，直至 D e，取最終得到的 作為 y(xn+1) 的近似值。,D e，再將步長(zhǎng)折半一次計(jì)算，最終得到符合精度要求的 y(xn+1) 的近似值。,45,單步法的收斂性,顯式單步法可統(tǒng)一寫(xiě)成：,增量函數(shù)，僅依賴(lài)于函數(shù) f，且僅僅是 xn, yn, h 的函數(shù),求 y = y(x),h 0 時(shí)，近似解 是否收斂到精確解,，它應(yīng)當(dāng) 是一個(gè)固定節(jié)點(diǎn)，因 此 h 0 時(shí)應(yīng)同時(shí)附 帶 n ,46,單步法的收斂性（續(xù)）,對(duì)于 p 階的常微分方程數(shù)值算法，當(dāng) h 0， n 時(shí)，是否 yn+1 y(xn+1)？,p 階算法的局部截?cái)嗾`差為：,顯然：,局部截?cái)嗾`差的前提假設(shè)是：,局部截?cái)嗾`差 0 并不能保證算法

14、收斂,47,單步法的收斂性（續(xù)）,定義：若求解某初值問(wèn)題的單步數(shù)值法，對(duì)于固定的 當(dāng) h 0 且 n 時(shí)，它的近似 解趨向于精確解 y(xn)，即：,則稱(chēng)該單步法是收斂的。,定義：稱(chēng) y(xn) - yn 為單步法的近似解 yn 的整體截?cái)?誤差。,單步法收斂,48,單步法的收斂性（續(xù)）,收斂性定理,若某單步法滿(mǎn)足以上條件，則該方法是收斂的,則該單步法的整體截?cái)嗾`差為：,若單步法 具有 p 階精度，且增量函數(shù) 關(guān)于 y 滿(mǎn)足：,Lipschitz 條件：,初值 y0 是準(zhǔn)確的,49,假設(shè)在前一步 yn 準(zhǔn)確的前提下求得的近似值為：,算法精度為 p 階，局部截?cái)嗾`差：,50,51,即：,若初值是

15、準(zhǔn)確的，則 e 0 = 0 ，從而整體截?cái)嗾`差為：,y = e x 為單調(diào)增函數(shù)，當(dāng) 時(shí),當(dāng) h 0 且 n 時(shí)，則,52,單步法的穩(wěn)定性,在討論單步法收斂性時(shí)一般認(rèn)為數(shù)值方法本身的計(jì)算過(guò)程是準(zhǔn)確的，實(shí)際上并非如此：,初始值 y0 有誤差 d = y0 - y(x0) 后續(xù)的每一步計(jì)算均有舍入誤差,這些初始和舍入誤差在計(jì)算過(guò)程的傳播中是逐步衰減的還是惡性增長(zhǎng)就是數(shù)值方法的穩(wěn)定性問(wèn)題,53,定義：若一種數(shù)值方法在節(jié)點(diǎn) xn 處的數(shù)值解 yn 的擾動(dòng) ，而在以后各節(jié)點(diǎn) ym (m n) 上產(chǎn)生的擾動(dòng)為 ，如果：,單步法的穩(wěn)定性（續(xù)）,定義：設(shè)在節(jié)點(diǎn) xn 處用數(shù)值算法得到的理想數(shù)值解為 yn，而實(shí)

16、際計(jì)算得到的近似解為 ，稱(chēng)差值：,為第 n 步的數(shù)值解的擾動(dòng)。,則稱(chēng)該數(shù)值方法是穩(wěn)定的。,54,單步法的穩(wěn)定性（續(xù)）,歐拉法：,由于函數(shù) f (x, y) 的多樣性，數(shù)值穩(wěn)定性的分析相當(dāng)復(fù)雜，通常只研究模型方程,考察模型方程：,即：,假設(shè)在節(jié)點(diǎn)值 yn 上有擾動(dòng) n，在節(jié)點(diǎn)值 yn+1 上有擾動(dòng) n+1，且 n+1 僅由 n 引起（即：計(jì)算過(guò)程中不再引起新的誤差）,55,歐拉法穩(wěn)定,即：,歐拉法穩(wěn)定的條件：,針對(duì)模型方程： 的顯式歐拉法：,化簡(jiǎn)得：,56,隱式歐拉法：,考察模型方程：,即：,化簡(jiǎn)為：,假設(shè) yn 上有擾動(dòng) ，則 yn+1 的擾動(dòng)為：,隱式歐拉法穩(wěn)定,，上式均成立，所以：,隱式歐拉法穩(wěn)定是恒穩(wěn)