第九章常微分方程的數(shù)值解法課件

1、 第九章 常微分方程的數(shù)值解法 1、引言 2、初值問題的數(shù)值解法-單步法 3、龍格-庫塔方法 4、收斂性與穩(wěn)定性 5、初值問題的數(shù)值解法多步法 6、方程組和剛性方程 7、習(xí)題和總結(jié)主 要 內(nèi) 容主 要 內(nèi) 容研究的問題數(shù)值解法的意義1、 引 言現(xiàn)實(shí)世界中大多數(shù)事物內(nèi)部聯(lián)系非常復(fù)雜找出其狀態(tài)和狀態(tài)變化規(guī)律之間的相互聯(lián)系，也即一個(gè)或一些函數(shù)與他們的導(dǎo)數(shù)之間的關(guān)系此種關(guān)系的數(shù)學(xué)表達(dá)就為微分方程1.什么是微分方程 ?其狀態(tài)隨著時(shí)間、地點(diǎn)、條件的不同而不同如何利用數(shù)值方法求解微分方程（組）的問題。2.數(shù)值求解微分方程的意義如何建立數(shù)學(xué)模型已在建模課程中得到討論，各類微分方程本身和他們的解所具有的特性已在

2、常微分方程及數(shù)學(xué)物理方程中得以解釋，本章專門討論尋找解析解的過程稱為求解微分方程組。 一個(gè)或一組具有所要求階連續(xù)導(dǎo)數(shù)的解析函數(shù)，將它代入微分方程（組），恰使其所有條件都得到滿足的解稱為解析解(或古典解),稱為真解或解。3.什么是微分方程(組)的解析解?3.什么是微分方程(組)的解析解?4.什么是微分方程的數(shù)值解?雖然求解微分方程有許多解析方法,但解析方法只能夠求解一些特殊類型的方程,從實(shí)際意義上來講我們更關(guān)心的是某些 特定的自變量在某一個(gè)定義范圍內(nèi)的一系列離散點(diǎn)上的近似值.尋找數(shù)值解的過程稱為數(shù)值求解微分方程。把這樣一組近似解稱為微分方程在該范圍內(nèi)的數(shù)值解在大量的實(shí)際方程中出現(xiàn)的函數(shù)起碼的連續(xù)

3、性都無法保證，更何況要求階的導(dǎo)數(shù)求解數(shù)值解很多微分方程根本求不到問題的解析解！重要手段。常微分方程的數(shù)值解法常用來求近似解根據(jù)提供的算法通過計(jì)算機(jī)便捷地實(shí)現(xiàn)5.常微分方程數(shù)值解法的特點(diǎn)數(shù)值解法得到的近似解(含誤差)是一個(gè)離散的函數(shù)表.本章主要討論一階常微方程的初值問題6.基本知識(shí)其中f (x,y)是已知函數(shù)，(1.2)是定解條件也稱為初值條件。各種數(shù)值解法則稱 f (x,y) 對(duì)y 滿足李普希茲條件，L 稱為Lipschitz常數(shù).常微分方程的理論指出：當(dāng) f (x,y) 定義在區(qū)域 G=(axb,y)若存在正的常數(shù) L 使：就可保證方程解的存在唯一性(Lipschitz)條件我們以下的討論，

4、都在滿足上述條件下進(jìn)行.若 f (x,y) 在區(qū)域 G連續(xù)，關(guān)于y一階常微分方程的初值問題的解存在且唯一.滿足李普希茲條件一階常微分方程組常表述為：方程組初值條件寫成向量形式為高階常微分方程定解問題如二階定解問題：也就解決了高階方程的定解問題.這些解法都可以寫成向量形式用于一階常微分方程組的初值問題.簡單的數(shù)值方法與基本概念 1. 簡單歐拉法（Euler） 2后退的歐拉法 3梯形法 4改進(jìn)Euler法 2、初值問題的數(shù)值解法單步法1. 簡單的歐拉(Euler)方法考慮模型：在精度要求不高時(shí)通過歐拉方法的討論弄清常微方程初值問題數(shù)值解法的一些基本概念和構(gòu)造方法的思路.歐拉方法最簡單而直觀實(shí)用方法

5、2. 歐拉方法的導(dǎo)出把區(qū)間a,b分為n個(gè)小區(qū)間步長為要計(jì)算出解函數(shù) y(x) 在一系列節(jié)點(diǎn)節(jié)點(diǎn)處的近似值N等分對(duì)微分方程(1.1)兩端從進(jìn)行積分右端積分用左矩形數(shù)值求積公式得x0 x1亦稱為歐拉折線法 /* Eulers polygonal arc method*/ 每步計(jì)算只用到或用向前差商近似導(dǎo)數(shù)依上述公式逐次計(jì)算可得：例題3.歐拉公式有明顯的幾何意義依此類推得到一折線故也稱Euler為單步法。公式右端只含有已知項(xiàng)所以又稱為顯格式的單步法。也稱歐拉折線法.就是用這條折線近似地代替曲線歐拉方法從上述幾何意義上得知，由Euler法所得的折線明顯偏離了積分曲線，可見此方法非常粗糙。4.歐拉法的局

6、部截?cái)嗾`差：在假設(shè)第 i 步計(jì)算是精確的前提下，考慮截?cái)嗾`差稱為局部截?cái)嗾`差/* local truncation error */。定義若某算法的局部截?cái)嗾`差為O(hp+1)，則稱該算法有p 階精度。Ri 的主項(xiàng)/* leading term */ 歐拉法的局部截?cái)嗾`差：歐拉法具有 1 階精度。定義在假設(shè) yi = y(xi)，即第 i 步計(jì)算是精確的前提下，考慮的截?cái)嗾`差 Ri = y(xi+1) yi+1 稱為局部截?cái)嗾`差 /* local truncation error */。如果單步差分公式的局部截?cái)嗾`差為O(hp+1),則稱該公式為p階方法.這里p為非負(fù)整數(shù).顯然,階數(shù)越高,方法

7、的精度越高.Taylor展開式,一元函數(shù)的Taylor展開式為:若某算法的局部截?cái)嗾`差為O(hp+1)，則稱該算法有p 階精度。Ri 的主項(xiàng)/* leading term */5. 歐拉公式的改進(jìn)： 隱式歐拉法 /* implicit Euler method */向后差商近似導(dǎo)數(shù)x0 x1)(,()(1101xyxfhyxy+)1,.,0(),(111-=+=+niyxfhyyiiii由于未知數(shù) yi+1 同時(shí)出現(xiàn)在等式的兩邊，不能直接得到，故稱為隱式 /* implicit */ 歐拉公式，而前者稱為顯式 /* explicit */ 歐拉公式。一般先用顯式計(jì)算一個(gè)初值，再迭代求解。 隱式

8、歐拉法的局部截?cái)嗾`差：即隱式歐拉公式具有 1 階精度。6.梯形公式 /* trapezoid formula */ 顯、隱式兩種算法的平均注：的確有局部截?cái)嗾`差 ， 即梯形公式具有2 階精度，比歐拉方法有了進(jìn)步。但注意到該公式是隱式公式，計(jì)算時(shí)不得不用到迭代法，其迭代收斂性與歐拉公式相似。 中點(diǎn)歐拉公式 /* midpoint formula */中心差商近似導(dǎo)數(shù)x0 x2x1假設(shè) ，則可以導(dǎo)出即中點(diǎn)公式具有 2 階精度。需要2個(gè)初值 y0和 y1來啟動(dòng)遞推過程，這樣的算法稱為雙步法 /* double-step method */，而前面的三種算法都是單步法 /* single-step m

9、ethod */。方 法顯式歐拉隱式歐拉梯形公式中點(diǎn)公式簡單精度低穩(wěn)定性最好精度低, 計(jì)算量大精度提高計(jì)算量大精度提高, 顯式多一個(gè)初值, 可能影響精度 改進(jìn)歐拉法 /* modified Eulers method */Step 1: 先用顯式歐拉公式作預(yù)測，算出),(1iiiiyxfhyy+=+Step 2: 再將 代入隱式梯形公式的右邊作校正，得到1+iy),(),(2111+=iiiiiiyxfyxfhyy注：此法亦稱為預(yù)測-校正法 /* predictor-corrector method */?？梢宰C明該算法具有 2 階精度，同時(shí)可以看到它是個(gè)單步遞推格式，比隱式公式的迭代求解過程

10、簡單。后面將看到，它的穩(wěn)定性高于顯式歐拉法。3 龍格 - 庫塔法 /* Runge-Kutta Method */ 考察改進(jìn)的歐拉法，可以將其改寫為：斜率一定取K1 K2 的平均值嗎？步長一定是一個(gè)h 嗎？單步遞推法的基本思想是從 ( xi , yi ) 點(diǎn)出發(fā)，以某一斜率沿直線達(dá)到 ( xi+1 , yi+1 ) 點(diǎn)。歐拉法及其各種變形所能達(dá)到的最高精度為2階。建立高精度的單步遞推格式。首先希望能確定系數(shù) 1、2、p，使得到的算法格式有2階精度，即在 的前提假設(shè)下，使得 Step 1: 將 K2 在 ( xi , yi ) 點(diǎn)作 Taylor 展開將改進(jìn)歐拉法推廣為：),(),(121221

11、11phKyphxfKyxfKKKhyyiiiiii+=+=+llStep 2: 將 K2 代入第1式，得到Step 3: 將 yi+1 與 y( xi+1 ) 在 xi 點(diǎn)的泰勒展開作比較要求 ，則必須有：這里有 個(gè)未知數(shù)， 個(gè)方程。32存在無窮多個(gè)解。所有滿足上式的格式統(tǒng)稱為2階龍格 - 庫塔格式。注意到， 就是改進(jìn)的歐拉法。 Q: 為獲得更高的精度，應(yīng)該如何進(jìn)一步推廣？其中i ( i = 1, , m )，i ( i = 2, , m ) 和 ij ( i = 2, , m; j = 1, , i1 ) 均為待定系數(shù)，確定這些系數(shù)的步驟與前面相似。 ).,(.),(),(),(.1122

12、112321313312122122111-+=+=+=+=mm mmmmimiiiiiimmiihKhKhKyhxfKhKhKyhxfKhKyhxfKyxfKKKKhyybbbabbaballl 最常用為四級(jí)4階經(jīng)典龍格-庫塔法 /* Classical Runge-Kutta Method */ ：注： 龍格-庫塔法的主要運(yùn)算在于計(jì)算 Ki 的值，即計(jì)算 f 的值。Butcher 于1965年給出了計(jì)算量與可達(dá)到的最高精度階數(shù)的關(guān)系：753可達(dá)到的最高精度642每步須算Ki 的個(gè)數(shù) 由于龍格-庫塔法的導(dǎo)出基于泰勒展開，故精度主要受解函數(shù)的光滑性影響。對(duì)于光滑性不太好的解，最好采用低階算法而

13、將步長h 取小。深入研究龍格-庫塔法請(qǐng)看此處！4 收斂性與穩(wěn)定性 /* Convergency and Stability */ 1.收斂性 /* Convergency */定義 若某算法對(duì)于任意固定的 x = xi = x0 + i h，當(dāng) h0 ( 同時(shí) i ) 時(shí)有 yi y( xi )，則稱該算法是收斂的。 例：就初值問題 考察歐拉顯式格式的收斂性。解：該問題的精確解為 歐拉公式為對(duì)任意固定的 x = xi = i h ，有 2.穩(wěn)定性 /* Stability */例：考察初值問題 在區(qū)間0, 0.5上的解。分別用歐拉顯、隱式格式和改進(jìn)的歐拉格式計(jì)算數(shù)值解。0.00.10.20.3

14、0.40.5精確解改進(jìn)歐拉法 歐拉隱式歐拉顯式 節(jié)點(diǎn) xi 1.00002.0000 4.00008.0000 1.6000101 3.2000101 1.00002.5000101 6.25001021.56251023.90631039.76561041.00002.50006.25001.56261013.90631019.76561011.00004.97871022.47881031.23411046.14421063.0590107What is wrong ?!定義若某算法在計(jì)算過程中任一步產(chǎn)生的誤差在以后的計(jì)算中都逐步衰減，則稱該算法是絕對(duì)穩(wěn)定的 /*absolutely st

15、able */。一般分析時(shí)為簡單起見，只考慮試驗(yàn)方程 /* test equation */常數(shù)，可以是復(fù)數(shù)當(dāng)步長取為 h 時(shí)，將某算法應(yīng)用于上式，并假設(shè)只在初值產(chǎn)生誤差 ，則若此誤差以后逐步衰減，就稱該算法相對(duì)于 絕對(duì)穩(wěn)定， 的全體構(gòu)成絕對(duì)穩(wěn)定區(qū)域。我們稱算法A 比算法B 穩(wěn)定，就是指 A 的絕對(duì)穩(wěn)定區(qū)域比 B 的大。hl h=h例：考察顯式歐拉法由此可見，要保證初始誤差0 以后逐步衰減，必須滿足：0-1-2ReImg例：考察隱式歐拉法可見絕對(duì)穩(wěn)定區(qū)域?yàn)椋?10ReImg注：一般來說，隱式歐拉法的絕對(duì)穩(wěn)定性比同階的顯式法的好。例：隱式龍格-庫塔法而顯式 1 4 階方法的絕對(duì)穩(wěn)定區(qū)域?yàn)槠渲?階

16、方法 的絕對(duì)穩(wěn)定區(qū)域?yàn)?ReImgk=1k=2k=3k=4-1-2-3-123ReImg無條件穩(wěn)定例1 用歐拉方法，隱式歐拉方法和歐拉中點(diǎn)公式計(jì)算的近似解，取步長h=0.1，并與精確值比較解:歐拉方法的算式為：歐拉隱式方法在本題可解出方程，不必迭代，公式為：歐拉中點(diǎn)公式是兩步法，第一步y(tǒng)1用歐拉公式，以后用公式本題精確解為y=e-x,計(jì)算結(jié)果見表9-1例2 用歐拉公式和梯形公式的預(yù)估校正法計(jì)算：的數(shù)值解，取h=0.1，梯形公式只迭代一次，并與精確值比較.方程的解析解為:解: 本例中歐拉公式為：梯形公式只校正一次的格式為結(jié)果列入表9.21. Runge-Kutta 法的一般形式 2. 2階Run

17、ge-Kutta 方法 3. 經(jīng)典Runge-Kutta 方法 4R-K-Fehhlberg 方法 5. 隱式R-K方法 6. 變步長方法 龍格庫塔法深入研究1.Runge-Kutta 法的一般形式Runge-Kutta 法是用區(qū)間上若干個(gè)點(diǎn)上的導(dǎo)數(shù)的線性組合得到平均斜率，以提高方法的階。其一般形式為 ：式(9.11) 稱L級(jí)p階Runge-Kutta方法(簡稱R-K法)。當(dāng)L1就是歐拉法，當(dāng)L2時(shí)為改進(jìn)的歐拉法。 其中它的局部截?cái)嗾`差是(9.11)2. 2級(jí)2階Runge-Kutta方法令 L=2，則 其局部截?cái)嗾`差是： 這是3個(gè)未知量的兩個(gè)方程，其中有一個(gè)自由參數(shù)，方程組有無窮多解，從而得

18、到一族2級(jí)2階R-K方法 。3. 經(jīng)典Runge-Kutta方法我們可以構(gòu)造出一族3級(jí)3階，一族4級(jí)4階和一族5級(jí)4階等R-K方法。最常用的4級(jí)4階是如下的經(jīng)典R-K方法：4R-K-Fehhlberg 方法R-K-Fehhlberg方法是在R-K方法的基礎(chǔ)上引進(jìn)誤差和步長控制的辦法。即利用5階R-K法 估計(jì)4階R-K的局部誤差，其中 注：R-K-Fehhlberg比4階R-K方法具有更大的優(yōu)越性，他是計(jì)算穩(wěn)定，高精度的方法，他的顯著優(yōu)點(diǎn)是，每一步僅需計(jì)算f的6個(gè)值；若用4階R-K-L與5階R-K-L在一起使用，則每步需要計(jì)算f的10個(gè)值。推薦使用！ 5. 隱式R-K方法類似于顯式R-K公式，稍

19、加改變，就得到隱式R-K方法。它與顯式R-K公式的區(qū)別在于：顯式公式中對(duì)系數(shù)求和的上限是i-1，從而構(gòu)成的矩陣是一個(gè)嚴(yán)格下三角陣。而在隱式公式中對(duì)系數(shù)求和的上限是L，從而構(gòu)成的矩陣是方陣，需要用迭代法求出Ki，。推導(dǎo)隱式公式的思路和方法與顯式R-K法類似。通常，同級(jí)的隱式公式獲得比顯式公式更高的階。通常，同級(jí)的隱式公式獲得比顯式公式更高的階。常用的隱式R-K法有:1級(jí)2階中點(diǎn)公式 :2級(jí)2階梯形公式： 2級(jí)4階R-K公式： 6.變步長方法在單步法中每一積分步步長實(shí)際上是相互獨(dú)立的，步長的選擇具有了靈活性。根據(jù)合理地選擇每一積分步的步長，既保證精度的要求，又可以減少計(jì)算量，從而減小舍入誤差。其方

20、便的控制手段是基于誤差的事后估計(jì)式。對(duì)于給定的精度 ,如果 ,反復(fù)將步長折半進(jìn)行計(jì)算，直至 為止,這時(shí)取最終得到的作為結(jié)果；如果 為止，這時(shí)再將步長折半一次，就得到所要的結(jié)果。這種通過加倍或折半處理步長的計(jì)算方法稱為變步長方法。 注：推薦使用精度好計(jì)算量低的變步長方法。用四階顯式R-K方法做變步長方法是實(shí)踐中較好的方法！ 例 分別用改進(jìn)的歐拉格式和四階龍格庫塔格式解初值問題（取步長h=0.2）：節(jié)點(diǎn) 改進(jìn)歐拉法 四階龍格庫塔法 準(zhǔn)確解 0 1 1 1 0.2 1.186667 1.183229 1.1832160.4 1.348312 1.341667 1.3416410.6 1.493704

21、 1.483281 1.4832400.8 1.627861 1.612514 1.612452 1 1.754205 1.732142 1.732051表（注：已指出過準(zhǔn)確解）單步法的相容性、收斂性和穩(wěn)定性 1.單步法的相容性 2.單步法的收斂性 3.單步法的穩(wěn)定性 4.相容性和收斂性的關(guān)系 5.相容性和方法階的關(guān)系 6.穩(wěn)定性和收斂性 7.絕對(duì)穩(wěn)定性和絕對(duì)穩(wěn)定域 單步法的相容性定義一：對(duì)于（9.1.1）常微分方程初值問題單步法的形式可以變表示為 （9.2.19）其中 h為步長若對(duì)求解區(qū)間中任一固定的x，當(dāng) 時(shí)皆成立 則稱由（9.2.19）確定的單步法與微分方程初值問題是相容的注意到上式左邊

22、極限為由（9.1.1）知它應(yīng)等于從而由相容性定義得稱相容性條件。 單步法的收斂性定義二： 設(shè) y(x) 是（9.1.1）的解， 是單步法（9.2.19）產(chǎn)生的數(shù)值解，對(duì)于每一個(gè)固定的 , ， 當(dāng) 即 。若成立 ， 則稱該方法是收斂的。單步法的穩(wěn)定性定義三： 若一個(gè)數(shù)值方法在基點(diǎn) 處的值有 的擾動(dòng)，在 此后各基點(diǎn) （mn）處的值產(chǎn)生的偏差均不超 過 ，則稱該方法是穩(wěn)定的。 單步法的穩(wěn)定性有以下定理定理二： 若單步法的增量函數(shù)對(duì)變量y滿足 Lipschtiz 條件， 則單步方法是穩(wěn)定的。相容性和收斂性的關(guān)系 定理一： 若單步法的增量函數(shù)對(duì)變量y滿足Lipschtiz 條件，即存在與 h , x 無

23、關(guān)的常數(shù) L，對(duì)區(qū)域D= 任意兩點(diǎn)（x，y1）,(x,y2)成立，則單步法收斂的充分必要條件是相容性條件成立。（讀者自證）相容性和方法階的關(guān)系 若單步法是p階方法則成立 若單步法滿足相容性條件，得 所以 =0也就是說單步法的階數(shù)一定要是正數(shù)。由于我們考慮的單步法皆為正整數(shù)，p至少為一。因此我們考慮的單步法都滿足相容性條件。 穩(wěn)定性和收斂性的關(guān)系若單步法的增量函數(shù)滿足定理二的條件即單步法是穩(wěn)定的則單步法收斂的充分必要條件是相容性條件成立。 絕對(duì)穩(wěn)定性和絕對(duì)穩(wěn)定域 穩(wěn)定性問題是一個(gè)比較復(fù)雜的問題。為了簡化討論一般僅對(duì)試驗(yàn)方程 進(jìn)行考察。這里假定Re0和的允許范圍，稱為該方法的 絕對(duì)穩(wěn)定域。 絕對(duì)穩(wěn)

24、定性和絕對(duì)穩(wěn)定域2將Euler方法應(yīng)用到試驗(yàn)方程得誤差方程是要求誤差不增長則則Euler 方法的絕對(duì)穩(wěn)定域是以 為半徑的圓的內(nèi)部。為了保證穩(wěn)定性步長有所控制。假如當(dāng) 時(shí)h應(yīng)滿足 ，當(dāng) 時(shí)， h 應(yīng)滿足 等等。同樣我們可以將試驗(yàn)方程用到其它各種單步法當(dāng)中，求出其絕對(duì)穩(wěn)定域。在實(shí)際應(yīng)用中必須在這個(gè)范圍內(nèi)，否則誤差傳播相當(dāng)大，即不穩(wěn)定。 1.Adams方法 2.米爾尼方法、漢明方法及辛普森方法 3.預(yù)測校正方法 4.多步法的相容性、穩(wěn)定性和收斂性 5 初值問題的數(shù)值解法多步法 考慮型如 的k步法， 稱為阿當(dāng)姆斯（Adams）方法1.Adams方法拿其中一種為例，推導(dǎo)其公式。對(duì)上面所說的法一般形式若取

25、 ，有方程組同樣解之，得到3步4階隱式Adams公式和它的余項(xiàng)。其它請(qǐng)讀者自證。我們僅將結(jié)論列于下表。 Adams顯式公式kp公式11223344Adams隱式公式kp公式12233445注：在這些方法中，4階的Adams預(yù)測校正方法具有方法簡潔、使用方便、易排序、高精度等優(yōu)點(diǎn)。尤其當(dāng)函數(shù)f比較復(fù)雜時(shí)更能顯出它的優(yōu)越性。同理得到5個(gè)待定參數(shù)方程組。解之得 ， ， ， ， 。構(gòu)成著名的Miline 4步4階顯式公式和它的余項(xiàng)。 ， 2.米爾尼方法、漢明方法及辛普森方法 同理得到5個(gè)參數(shù)方程組。求解后就構(gòu)成著名的3步4階隱式Hamming公式和它的余項(xiàng)。 若取，也得到5個(gè)參數(shù)方程組。求解后就構(gòu)成S

26、impson隱式公式和其余項(xiàng)。 米爾尼方法、漢明方法及辛普森方法 不論單步法或多步法，隱式公式比顯式公式穩(wěn)定性好，但在實(shí)際使用隱式公式時(shí)，都會(huì)遇到兩個(gè)問題：一個(gè)是隱式公式如何能方便地進(jìn)行計(jì)算；另一個(gè)是實(shí)際計(jì)算步長取多大。如隱式梯形公式，每往前推進(jìn)一步，不必進(jìn)行多次迭代，而是采用一階顯式Euler公式預(yù)測，二階隱式梯形公式校正一次，構(gòu)成顯式改進(jìn)Euler公式，能達(dá)到與梯形公式同階的精度，即二階精度。 3. 預(yù)測校正方法 對(duì)于線性多步公式，一般地，不難驗(yàn)證，如果 預(yù)測公式是階或階精度，校正公式為階精度， 則用預(yù)測公式提供初值，校正公式迭代一次的 效果也能達(dá)到階精度，再迭代下去，效果就不 明顯了。預(yù)

27、測-校正技術(shù)即保證了計(jì)算精度，又 使隱式計(jì)算顯式化，克服了隱式公式要反復(fù)迭 代的困難。至于實(shí)際計(jì)算步長的選取，我們對(duì) 預(yù)測-校正公式使用外推原理，得到誤差估計(jì)式 ，用來調(diào)整計(jì)算步長，使達(dá)到控制誤差要求。對(duì)于線性多步方法常用的預(yù)測-校正公式有Miline-Hamming方法和4階Adams顯隱式預(yù)測-校正公式 將Miline 公式和Hamming公式結(jié)合，構(gòu)成預(yù)測-校正公式 Miline公式和Hamming公式的局部截?cái)嗾`差分別為修正Miline-Hamming公式利用外推原理，即加上后消去局部截?cái)嗾`差的主項(xiàng)。使 說明經(jīng)過外推后的算法其精度提高一階。忽略誤差項(xiàng)，上式可改寫為 由于 和是 在計(jì)算過

28、程中獲得的數(shù)據(jù)，稱為Miline 公式和Hamming公式的事后誤差估計(jì)式。我們可以用它們來調(diào)節(jié)計(jì)算步長的大小，即選擇一個(gè)合適的的步長，使 ，其中的是要求達(dá)到的計(jì)算精度。 又可得到 Miline公式和Hamming公式的修正公式，它們分別是 從而構(gòu)成如下的修正Hamming預(yù)測-校正公式，簡稱修正Hamming公式： 預(yù)測：修正：校正： 修正： 在應(yīng)用這套公式時(shí)，先由同階單步法提供初值 ， ， ， 。計(jì)算 時(shí)可取。 將4步4階顯式Adams公式作為預(yù)測公式和3步4階式Adams公式作為校正公式，構(gòu)成Adams預(yù)測-校正公式。 它們的局部截?cái)嗾`差分別是 Adams預(yù)測-校正公式利用外推原理，將上

29、兩式作線性組合消去局部截?cái)嗾`差的主項(xiàng)。使計(jì)算精度至少提高一階，同時(shí)得到兩個(gè)修正公式，將它們和上兩式構(gòu)成了如下修正Adams公式：預(yù)測： 修正： 校正: 修正： 同理，在計(jì)算時(shí)，調(diào)節(jié)計(jì)算步長 使 , 由同階單步法提供初值 , ， ， 。計(jì)算 時(shí)，可取 。 理論分析得出它們的絕對(duì)穩(wěn)定區(qū)域，修正Hamming法： ; 4階Adams預(yù)測-校正法： 其中 我們也可以在教學(xué)演示上看出修正的預(yù)測校正格式的誤差非常的小。 多步法的相容性條件 k步法的一般形式為 其中 由對(duì)單步法的討論可知,當(dāng) 時(shí),方法階數(shù)至少為1。即對(duì) y1,y=x 應(yīng)精確成立。令 y 1，所以令y=x得 4.多步法的相容性、穩(wěn)定性和收斂性

30、所以我們稱為線性多步法的相容性條件。 k步法需要k個(gè)出發(fā)值，而初值問題只提供了一個(gè)初值，在使用k步法時(shí)尚需要其它方法作補(bǔ)充k-1個(gè)出發(fā)值，今假定它們 是 ，當(dāng) 這k-1個(gè)值都 應(yīng)收斂于共同的極限y（a）即 在討論多步法收斂性時(shí)我們總假定（9.3.12）成立 定義四： 多步法的收斂問題 的解 收斂于初值問題的解 y(x)。這里 xa+nh. 定義五：稱多項(xiàng)式 為k步法的特征多項(xiàng)式。如果特征多項(xiàng)式的零點(diǎn)皆不大于1，且等 于1的零點(diǎn)是單重的，稱根條件成立。稱多項(xiàng)式 為第二特征多項(xiàng)式。 顯然根條件可以表示 定理三：k步法收斂的充要條件為： （1）相容性條件成立。 （2）特征多項(xiàng)式的零點(diǎn)皆不大于1，且等

31、于1的零點(diǎn)是單重的。 （稱2為）特征根條件。 多步法的穩(wěn)定性 關(guān)于線性多步法成立以下定理： 定理四：若函數(shù)f(x,y)對(duì)變量y滿足Lipschtiz 條件在與h,x無關(guān)的常數(shù)L，對(duì)區(qū)域D= 任意兩點(diǎn)（x，y1）,(x,y2)成立 k步法的相容性、收斂性、穩(wěn)定性有以下關(guān)系 對(duì)于常微分方程右端函數(shù)f(x,y)，在相容性條件成立情況下，k步法的收斂性和穩(wěn)定性是等價(jià)的。 事實(shí)上在相容性條件成立時(shí)，收斂性和穩(wěn)定性的充要條件都是特征根條件成立。多步法的絕對(duì)穩(wěn)定性和絕對(duì)穩(wěn)定域 將線性多步法的公式應(yīng)用到試驗(yàn)方程 進(jìn)行考察。這里假定Re 0,即試驗(yàn)方程本身是穩(wěn)定的。 記 得 是常系數(shù)差分方程，其特征方程為 記它

32、的 k個(gè)特征根為 并設(shè)它們是互異的。顯然根與 有關(guān)，不妨記為 注意到當(dāng) 時(shí) 這時(shí)由特征方程得 由線性多步法的相容性條件得 是一個(gè)根。不妨設(shè), 差分方程的解為 其中系數(shù)由線性多步法的出發(fā)值確定。另一方面，y(0)=1的試驗(yàn)方程的精確解為 ， 設(shè)多步法截?cái)嗾`差為 ，由此可得 ，我們稱 為主根，其它根都為增根。 定義五：線性多步法的絕對(duì)穩(wěn)定區(qū)域?qū)o定的 ，如果特征方程的特征根 皆按模小于1，則線性多步法關(guān)于u是絕對(duì)穩(wěn)定的。使得 成立的 構(gòu)成絕對(duì)穩(wěn)定區(qū)域。注：從誤差角度來看絕對(duì)穩(wěn)定區(qū)域的方法是一個(gè)理想的方法。這樣，絕對(duì)穩(wěn)定區(qū)域越大，方法適用性越廣，因而越優(yōu)越。 6 方程組和剛性方程 9.1 一階方程組

33、 9.2 化高階方程為一階方程組 9.3 剛性方程組9.1 一階方程組 前面我們研究了單個(gè)方程y=f的數(shù)值解法，只要把y和f理解為向量，那么，所提供的各種計(jì)算公式即可應(yīng)用到一階方程組的情形。 考慮一階方程組的出值問題 若采用向量的記號(hào)，記 則常微分方程組的出值問題可以表示為 前幾節(jié)我們主要討論了常微分方程的出值問題的數(shù)值解法。只要將 y 和f 改寫為向量，那么前面提供的各種計(jì)算公式即可適用于一階常微分方程組的出值問題。 Runge-Kutta方法 對(duì)于方程組 四級(jí)四階顯示Runge- Kutta公式 若寫成分量形式就是 i=1,2,m 為了幫助理解這一公式的計(jì)算過程，我們?cè)倏紤]兩個(gè)方程的特殊情

34、形： 這時(shí)四階龍格庫塔公式具有形式 其中 這是一步法，利用節(jié)點(diǎn) 上的值 ， ，由上式順序計(jì) 算 ，然后代入即可求得 節(jié)點(diǎn)上的 。 9.2化高階方程為一階方程組關(guān)于高階微分方程的初值問題，原則上總可以歸結(jié)為一階方程組來求解，例如，考察下列m 階微分方程初始條件為只要引進(jìn)新的變量即可將m 階方程化為如下的一階方程組 初始條件則相應(yīng)的化為9.3 剛性方程組在求解方程組時(shí)，經(jīng)常出現(xiàn)解的分量數(shù)量級(jí)差別很大的情形，這給數(shù)值求解帶來很大困難，這種問題稱為剛性問題。若線性系統(tǒng) ，其中，中A的特征值 滿足條件0(j=1,N), 且 則稱 系統(tǒng)為剛性方程，稱s為剛性比 7 習(xí)題與總結(jié)1數(shù)值例題2數(shù)值練習(xí) 3總結(jié)1

35、、數(shù)值例題 我們已經(jīng)學(xué)習(xí)了很多數(shù)值算法，他們的效果到底如何呢？下面我們來分析一道例題，看看那些方法，就這個(gè)問題，最能接近真實(shí)值 求初值問題的數(shù)值解，取h=0.1并同精確解 比較 （1）用歐拉法來計(jì)算這個(gè)初值問題（2）用各階的RungeKutta方法來計(jì)算這個(gè)初值問題（3）用四階的Adams定步長，Adams定步長預(yù)測校正方 法來 計(jì)算這個(gè)初值問題。然后將數(shù)值結(jié)果列成表格：自變量 解析解 歐拉法 后退歐拉法 梯形法 改進(jìn)的歐拉法 0.11.0048 1.0000 1.0091 1.0048 1.0050 0.21.0187 1.0100 1.0264 1.0186 1.0190 0.31.040

36、8 1.0290 1.0513 1.0406 1.0412 0.41.0703 1.0561 1.0830 1.0701 1.0708 0.51.1065 1.0905 1.1209 1.1063 1.1071 0.61.1488 1.1314 1.1645 1.1485 1.1494 0.71.1966 1.1783 1.2132 1.1963 1.1972 0.81.2493 1.2305 1.2665 1.2490 1.2500 0.91.3066 1.2874 1.3241 1.3063 1.3072 1.01.3679 1.3488 1.3855 1.3676 1.3685 自變量

37、解析解 4階顯RK 4階Adams顯法 預(yù)測校正法 0.11.0048 1.0048 1.0048 1.0048 0.21.0187 1.0187 1.0187 1.0187 0.31.0408 1.0408 1.0408 1.0408 0.41.0703 1.0703 1.0703 1.0703 0.51.1065 1.1065 1.1065 1.1065 0.61.1488 1.1488 1.1488 1.1488 0.71.1966 1.1966 1.1966 1.1966 0.81.2493 1.2493 1.2493 1.2493 0.91.3066 1.3066 1.3066 1.3066 1.01.3679 1.3679 1.3679 1.3679 2、數(shù)值練習(xí)計(jì)算實(shí)習(xí)題：1.使用改進(jìn)的Euler算法求初值問題 的數(shù)值解，取步長h=0.1,并同精確解 比較 2. 利用4階Runge-Kutta算法求初值問題 的數(shù)值解，取步長：(1)h=0.1;(2)h=0.01.求數(shù)值解及精確解 3.使用4階定步長Adams預(yù)測校正算法求初值問題 的數(shù)值解，取h=0.1,并同精確解 比較 4.使用4階定步長Adams預(yù)測校正算法求初值問題的數(shù)值解，使用h=0.5 5綜合題請(qǐng)用已學(xué)的數(shù)值方法來計(jì)算下面初值問題，看看那