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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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