常微分方程的數(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ī)律

2、之間的相互聯(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é)模型已在建模課程中得到討論，各類微分方程本身和他們的解所具有的特性已在常微分方程及數(shù)學(xué)物理方程中得以解釋，,本章專門討論,尋找解析解的過程稱為求解微分方程組。,一個(gè)或一組具有所要求階連續(xù)導(dǎo)數(shù)的

3、解析函數(shù)，將它代入微分方程（組），恰使其所有條件都得到滿足的解稱為解析解(或古典解),稱為真解或解。,3.什么是微分方程(組)的解析解?,3.什么是微分方程(組)的解析解?,4.什么是微分方程的數(shù)值解?,雖然求解微分方程有許多解析方法,但解析方法只能夠求解一些特殊類型的方程,從實(shí)際意義上來講,我們更關(guān)心的是某些 特定的自變量在某一個(gè)定義范圍內(nèi)的一系列離散點(diǎn)上的近似值.,尋找數(shù)值解的過程稱為數(shù)值求解微分方程。,把這樣一組近似解

4、稱為微分方程在該范圍內(nèi)的,數(shù)值解,在大量的實(shí)際方程中出現(xiàn)的函數(shù)起碼的連續(xù)性都無法保證，更何況要求階的導(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

5、)是定解條件也稱為初值條件。,各種數(shù)值解法,,則稱 f (x,y) 對(duì)y 滿足李普希茲條件，L 稱為L(zhǎng)ipschitz常數(shù).,常微分方程的理論指出：,當(dāng) f (x,y) 定義在區(qū)域 G=(a≤x≤b,｜y｜＜∞),若存在正的常數(shù) L 使：,,就可保證方程解的存在唯一性,(Lipschitz)條件,我們以下的討論，都在滿足上述條件下進(jìn)行.,若 f (x,y) 在區(qū)域 G連續(xù)，關(guān)于y,一階常微分方程的初值問題的解存在且唯一.,滿足李普希

6、茲條件,一階常微分方程組常表述為：,方程組,初值條件,寫成向量形式為,,高階常微分方程定解問題如二階定解問題：,也就解決了高階方程的定解問題.,這些解法都可以寫成向量形式,用于一階常微分方程組的初值問題.,,,簡(jiǎn)單的數(shù)值方法與基本概念,1. 簡(jiǎn)單歐拉法（Euler） 2．后退的歐拉法 3．梯形法 4．改進(jìn)Euler法,,§2、初值問題的數(shù)值解法―單步法,1. 簡(jiǎn)單的歐拉(Euler)方法,考慮模型：,在精度要求

7、不高時(shí),通過歐拉方法的討論,弄清常微方程初值問題數(shù)值解法的一些基本概念和構(gòu)造方法的思路.,,歐拉方法,最簡(jiǎn)單而直觀實(shí)用方法,2. 歐拉方法的導(dǎo)出,把區(qū)間［a,b］,,分為n個(gè)小區(qū)間,步長(zhǎng)為,要計(jì)算出解函數(shù) y(x) 在一系列節(jié)點(diǎn),,節(jié)點(diǎn),處的近似值,N等分,對(duì)微分方程(1.1)兩端從,進(jìn)行積分,,,右端積分用左矩形數(shù)值求積公式,得,,,,,,,,,,,,亦稱為歐拉折線法 /* Euler’s polygonal arc method

8、*/,每步計(jì)算,只用到,或用向前差商近似導(dǎo)數(shù),,依上述公式逐次計(jì)算可得：,例題,3.歐拉公式有明顯的幾何意義,,依此類推得到一折線,故也稱Euler為單步法。,公式右端只含有已知項(xiàng),所以又稱為顯格式的單步法。,,也稱歐拉折線法.,就是用這條折線近似地代替曲線,歐拉方法,,,,,,,,從上述幾何意義上得知，由Euler法所得的折線明顯偏離了積分曲線，可見此方法非常粗糙。,4.歐拉法的局部截?cái)嗾`差：,在假設(shè)第 i 步計(jì)算是精確的前提

9、下，考慮截?cái)嗾`差,稱為局部截?cái)嗾`差,/* local truncation error */。,Ri 的主項(xiàng)/* leading term */,? 歐拉法的局部截?cái)嗾`差：,歐拉法具有 1 階精度。,如果單步差分公式的局部截?cái)嗾`差為O(hp+1),則稱該公式為p階方法.這里p為非負(fù)整數(shù).顯然,階數(shù)越高,方法的精度越高.,Taylor展開式,一元函數(shù)的Taylor展開式為:,若某算法的局部截?cái)嗾`差為O(hp+1)，則稱該算法有p

10、階精度。,Ri 的主項(xiàng)/* leading term */,5. 歐拉公式的改進(jìn)：,? 隱式歐拉法 /* implicit Euler method */,,,由于未知數(shù) yi+1 同時(shí)出現(xiàn)在等式的兩邊，不能直接得到，故稱為隱式 /* implicit */ 歐拉公式，而前者稱為顯式 /* explicit */ 歐拉公式。,,,,,,,,,,一般先用顯式計(jì)算一個(gè)初值，再迭代求解。,? 隱式歐拉法的局部截?cái)嗾`差：,即隱式歐拉公式具有

11、 1 階精度。,6.梯形公式 /* trapezoid formula */,— 顯、隱式兩種算法的平均,,注：的確有局部截?cái)嗾`差 ， 即梯形公式具有2 階精度，比歐拉方法有了進(jìn)步。但注意到該公式是隱式公式，計(jì)算時(shí)不得不用到迭代法，其迭代收斂性與歐拉公式相似。,? 中點(diǎn)歐拉公式 /* midpoint formula */,,,,,假設(shè)

12、 ，則可以導(dǎo)出即中點(diǎn)公式具有 2 階精度。,需要2個(gè)初值 y0和 y1來啟動(dòng)遞推過程，這樣的算法稱為雙步法 /* double-step method */，而前面的三種算法都是單步法 /* single-step method */。,簡(jiǎn)單,精度低,穩(wěn)定性最好,精度低, 計(jì)算量大,精度提高,計(jì)算量大,精度提高, 顯式,多一個(gè)初值, 可能影響精度,? 改進(jìn)歐拉法 /*

13、modified Euler’s method */,注：此法亦稱為預(yù)測(cè)-校正法 /* predictor-corrector method */。可以證明該算法具有 2 階精度，同時(shí)可以看到它是個(gè)單步遞推格式，比隱式公式的迭代求解過程簡(jiǎn)單。后面將看到，它的穩(wěn)定性高于顯式歐拉法。,,§3 龍格 - 庫塔法 /* Runge-Kutta Method */,斜率一定取K1 K2 的平均值嗎？,步長(zhǎng)一定是一個(gè)h 嗎？,單步遞

14、推法的基本思想是從 ( 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 展開,,Step 2: 將 K2 代入第1式，得到,

15、Step 3: 將 yi+1 與 y( xi+1 ) 在 xi 點(diǎn)的泰勒展開作比較,要求 ，則必須有：,,這里有 個(gè)未知數(shù)， 個(gè)方程。,3,2,存在無窮多個(gè)解。所有滿足上式的格式統(tǒng)稱為2階龍格 - 庫塔格式。,注意到， 就是改進(jìn)的歐拉法。,Q: 為獲得更高的精度，應(yīng)該如何進(jìn)一步推廣？,其中?i ( i = 1

16、, …, m )，?i ( i = 2, …, m ) 和 ?ij ( i = 2, …, m; j = 1, …, i?1 ) 均為待定系數(shù)，確定這些系數(shù)的步驟與前面相似。,? 最常用為四級(jí)4階經(jīng)典龍格-庫塔法 /* Classical Runge-Kutta Method */ ：,? 由于龍格-庫塔法的導(dǎo)出基于泰勒展開，故精度主要受解函數(shù)的光滑性影響。對(duì)于光滑性不太好的解，最好采用低階算法而將步長(zhǎng)h 取小。,深入研究龍格-庫塔

17、法請(qǐng)看此處！,§4 收斂性與穩(wěn)定性 /* Convergency and Stability */,? 1.收斂性 /* Convergency */,例：就初值問題 考察歐拉顯式格式的收斂性。,解：該問題的精確解為,歐拉公式為,,對(duì)任意固定的 x = xi = i h ，有,,?,?,? 2.穩(wěn)定性 /* Stability */,例：考察初值問題

18、 在區(qū)間[0, 0.5]上的解。分別用歐拉顯、隱式格式和改進(jìn)的歐拉格式計(jì)算數(shù)值解。,1.0000?2.0000 4.0000?8.0000 1.6000?101 ?3.2000?101,1.00002.5000?10?1 6.2500?10?21.5625?10?23.9063?10?39.7656?10?4,1.00002.50006.25001.5626?1013.90

19、63?1019.7656?101,1.00004.9787?10?22.4788?10?31.2341?10?46.1442?10?63.0590?10?7,What is wrong ??!,一般分析時(shí)為簡(jiǎn)單起見，只考慮試驗(yàn)方程 /* test equation */,,常數(shù)，可以是復(fù)數(shù),,例：考察隱式歐拉法,,可見絕對(duì)穩(wěn)定區(qū)域?yàn)椋?,注：一般來說，隱式歐拉法的絕對(duì)穩(wěn)定性比同階的顯式法的好。,例：隱式龍格-庫塔法,而顯式

20、1~ 4 階方法的絕對(duì)穩(wěn)定區(qū)域?yàn)?其中2階方法 的絕對(duì)穩(wěn)定區(qū)域?yàn)?無條件穩(wěn)定,例1 用歐拉方法，隱式歐拉方法和歐拉中點(diǎn)公式計(jì)算,的近似解，取步長(zhǎng)h=0.1，并與精確值比較,解:歐拉方法的算式為：,歐拉隱式方法在本題可解出方程，不必迭代，公式為：,歐拉中點(diǎn)公式是兩步法，第一步y(tǒng)1用歐拉公式，以后用公式,本題精確解為y=e-x,計(jì)算結(jié)果見表9-1,例2 用歐拉公式和

21、梯形公式的預(yù)估校正法計(jì)算：,的數(shù)值解，取h=0.1，梯形公式只迭代一次，并與精確值比較.方程的解析解為:,解: 本例中歐拉公式為：,梯形公式只校正一次的格式為,結(jié)果列入表9.2,1. Runge-Kutta 法的一般形式 2. 2階Runge-Kutta 方法 3. 經(jīng)典Runge-Kutta 方法 4．R-K-Fehhlberg 方法 5. 隱式R-K方法 6. 變步長(zhǎng)方法,,龍格－

22、庫塔法深入研究,1.Runge-Kutta 法的一般形式,Runge-Kutta 法是用區(qū)間上若干個(gè)點(diǎn)上的導(dǎo)數(shù)的線性組合得到平均斜率，以提高方法的階。其一般形式為 ：,,,式(9.11) 稱L級(jí)p階Runge-Kutta方法(簡(jiǎn)稱R-K法)。當(dāng)L＝1就是歐拉法，當(dāng)L＝2時(shí)為改進(jìn)的歐拉法。,其中,它的局部截?cái)嗾`差是,(9.11),2. 2級(jí)2階Runge-Kutta方法,令 L=2，則,,其局部截?cái)嗾`差是：,這是3個(gè)未知量的兩個(gè)方程，其

23、中有一個(gè)自由參數(shù)，方程組有無窮多解，從而得到一族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方法：,,4．R-K-Fehhlberg 方法,R-K-Fehhlberg方法是在R-K方法的基礎(chǔ)上引進(jìn)誤差和步長(zhǎng)控制的辦法。即利用5階R-K法,,,估計(jì)4階R-K的局部誤差，其中,,注：R-K-Fehhlberg比4階R-

24、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公式，稍加改變，就得到隱式R-K方法。,,,它與顯式R-K公式的區(qū)別在于：顯式公式中對(duì)系數(shù)求和的上限是i-1，從而構(gòu)成的矩陣是一個(gè)嚴(yán)格下三角陣。而在隱式公式中對(duì)系數(shù)求和的上限是L，從而構(gòu)成的矩陣是方陣，需要用迭代法求出Ki

25、，。推導(dǎo)隱式公式的思路和方法與顯式R-K法類似。通常，同級(jí)的隱式公式獲得比顯式公式更高的階。,通常，同級(jí)的隱式公式獲得比顯式公式更高的階。常用的隱式R-K法有:,,1級(jí)2階中點(diǎn)公式 :,2級(jí)2階梯形公式：,2級(jí)4階R-K公式：,6.變步長(zhǎng)方法,在單步法中每一積分步步長(zhǎng)實(shí)際上是相互獨(dú)立的，步長(zhǎng)的選擇具有了靈活性。根據(jù)合理地選擇每一積分步的步長(zhǎng)，既保證精度的要求，又可以減少計(jì)算量，從而減小舍入誤差。其方便的控制手段是基于誤差的事后估

26、計(jì)式。,,對(duì)于給定的精度 ,如果 > ,反復(fù)將步長(zhǎng)折半進(jìn)行計(jì)算，直至 為止，這時(shí)再將步長(zhǎng)折半一次，就得到所要的結(jié)果。,這種通過加倍或折半處理步長(zhǎng)的計(jì)算方法稱為變步長(zhǎng)方法。,注：推薦使用精度好計(jì)算量低的變步長(zhǎng)方法。用四階顯式R-K方法做變步長(zhǎng)方法是實(shí)踐中較好的方法！,例 分別用改進(jìn)的歐拉格式和四階龍格—庫塔格式解初值問題（取步長(zhǎng)h=0.2）：,節(jié)點(diǎn) 改進(jìn)歐拉法 四階龍格—庫塔法 準(zhǔn)

27、確解 0 1 1 1 0.2 1.186667 1.183229 1.1832160.4 1.348312 1.341667 1.3416410.6 1.493704 1.483281

28、 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.

29、相容性和收斂性的關(guān)系 5.相容性和方法階的關(guān)系 6.穩(wěn)定性和收斂性 7.絕對(duì)穩(wěn)定性和絕對(duì)穩(wěn)定域,,單步法的相容性,定義一：對(duì)于（9.1.1）常微分方程初值問題單步法的形式可以變表示為 （9.2.19）其中 h為步長(zhǎng)若對(duì)求解區(qū)間中任一固定的x，當(dāng) 時(shí)皆成立,,則稱由（9.2.19）確定的單步法與微分方程初值問題是相容的注意到上式左邊極限為由（9.1.1）知它應(yīng)等于從而由相容性定義得稱相容性

30、條件。,單步法的收斂性,定義二： 設(shè) y(x) 是（9.1.1）的解， 是單步法（9.2.19）產(chǎn)生的數(shù)值解，對(duì)于每一個(gè)固定的 , ， 當(dāng) 即 。若成立 ， 則稱該方法是收斂的。,單步法的穩(wěn)定性,定義三： 若一個(gè)數(shù)值方法在基點(diǎn) 處的值有 的擾動(dòng)，在 此后各

31、基點(diǎn) （m>n）處的值產(chǎn)生的偏差均不超 過 ，則稱該方法是穩(wěn)定的。 單步法的穩(wěn)定性有以下定理,,定理二： 若單步法的增量函數(shù)對(duì)變量y滿足 Lipschtiz 條件， 則單步方法是穩(wěn)定的。,相容性和收斂性的關(guān)系,定理一： 若單步法的增量函數(shù)對(duì)變量y滿足Lipschtiz 條件，即存在與 h , x

32、 無關(guān)的常數(shù) L，對(duì)區(qū)域D= 任意兩點(diǎn)（x，y1）,(x,y2)成立，則單步法收斂的充分必要條件是相容性條件成立。（讀者自證）,相容性和方法階的關(guān)系,若單步法是p階方法則成立 若單步法滿足相容性條件，得 所以 =0也就是說單步法的階數(shù)一定要是正數(shù)。由于我們考慮的單步法皆為正整數(shù)，p至少為一。因

33、此我們考慮的單步法都滿足相容性條件。,穩(wěn)定性和收斂性的關(guān)系,若單步法的增量函數(shù)滿足定理二的條件即單步法是穩(wěn)定的則單步法收斂的充分必要條件是相容性條件成立。,絕對(duì)穩(wěn)定性和絕對(duì)穩(wěn)定域,穩(wěn)定性問題是一個(gè)比較復(fù)雜的問題。為了簡(jiǎn)化討論一般僅對(duì)試驗(yàn)方程 進(jìn)行考察。這里假定Re0和的允許范圍，稱為該方法的 絕對(duì)穩(wěn)定域。,絕對(duì)穩(wěn)定性和絕對(duì)穩(wěn)定域2,將Euler方法應(yīng)用到試驗(yàn)方程得誤差方程是要求

34、誤差不增長(zhǎng)則,則Euler 方法的絕對(duì)穩(wěn)定域是以 為半徑的圓的內(nèi)部。為了保證穩(wěn)定性步長(zhǎng)有所控制。假如當(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)，

35、否則誤差傳播相當(dāng)大，即不穩(wěn)定。,,1.Adams方法 2.米爾尼方法、漢明方法及辛普森方法 3.預(yù)測(cè)校正方法 4.多步法的相容性、穩(wěn)定性和收斂性,§5 初值問題的數(shù)值解法―多步法,考慮型如 的k步法， 稱為阿當(dāng)姆斯（Adams）方法,1.Adams方法,拿其中一種為例，推導(dǎo)其公式。對(duì)上面所說的法一

36、般形式若取 ，有方程組同樣解之，得到3步4階隱式Adams公式和它的余項(xiàng)。,,其它請(qǐng)讀者自證。我們僅將結(jié)論列于下表。 Adams顯式公式,Adams隱式公式,注：在這些方法中，4階的Adams預(yù)測(cè)校正方法具有方法簡(jiǎn)潔、使用方便、易排序、高精度等優(yōu)點(diǎn)。尤其當(dāng)函數(shù)f比較復(fù)雜時(shí)更能顯出它的優(yōu)越性。,同理得到5個(gè)待定參數(shù)

37、方程組。解之得 ， ， ， ， 。構(gòu)成著名的Miline 4步4階顯式公式和它的余項(xiàng)。,，,2.米爾尼方法、漢明方法及辛普森方法,同理得到5個(gè)參數(shù)方程組。求解后就構(gòu)成著名的3步4階隱式Hamming公式和它的余項(xiàng)。

38、 若取，也得到5個(gè)參數(shù)方程組。求解后就構(gòu)成Simpson隱式公式和其余項(xiàng)。,米爾尼

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

40、公式是階或階精度，校正公式為階精度， 則用預(yù)測(cè)公式提供初值，校正公式迭代一次的 效果也能達(dá)到階精度，再迭代下去，效果就不 明顯了。預(yù)測(cè)-校正技術(shù)即保證了計(jì)算精度，又 使隱式計(jì)算顯式化，克服了隱式公式要反復(fù)迭 代的困難。至于實(shí)際計(jì)算步長(zhǎng)的選取，我們對(duì) 預(yù)測(cè)-校正公式使用外推原理，得到誤差估計(jì)式 ，用來調(diào)整計(jì)算步長(zhǎng)，使達(dá)到控制誤差要求。,對(duì)于線性多步方法常用的預(yù)測(cè)-校正公式有Miline-H

41、amming方法和4階Adams顯隱式預(yù)測(cè)-校正公式,將Miline 公式和Hamming公式結(jié)合，構(gòu)成預(yù)測(cè)-校正公式 Miline公式和Hamming公式的局部截?cái)嗾`差分別為,修正Miline-Hamming公式,,,利用外推原理，即加上后消去局部截?cái)嗾`差的主項(xiàng)。使 說明經(jīng)過外推后的算法其精度提高一階。忽略誤差項(xiàng)，上式可改

42、寫為,,由于 和是 在計(jì)算過程中獲得的數(shù)據(jù)，稱為Miline 公式和Hamming公式的事后誤差估計(jì)式。我們可以用它們來調(diào)節(jié)計(jì)算步長(zhǎng)的大小，即選擇一個(gè)合適的的步長(zhǎng)，使 ，其中的是要求達(dá)到的計(jì)算精度。,又可得到 Miline公式和Hamming公式的修正公式，它們分別是

43、 從而構(gòu)成如下的修正Hamming預(yù)測(cè)-校正公式，簡(jiǎn)稱修正Hamming公式：,預(yù)測(cè)：修正：校正： 修正： 在應(yīng)用這套公式時(shí)，先由同階單步法提供初值 ， ， ， 。計(jì)算 時(shí)可取。,將4步4階顯式Adams公式作為預(yù)測(cè)公式和3步4階式Adams公式作為校正公式，構(gòu)成Adams預(yù)測(cè)

44、-校正公式。 它們的局部截?cái)嗾`差分別是,Adams預(yù)測(cè)-校正公式,利用外推原理，將上兩式作線性組合消去局部截?cái)嗾`差的主項(xiàng)。使計(jì)算精度至少提高一階，同時(shí)得到兩個(gè)修正公式，將它們和上兩式構(gòu)成了如下修正Adams公式：預(yù)測(cè)： 修正： 校正:,修正：

45、 同理，在計(jì)算時(shí)，調(diào)節(jié)計(jì)算步長(zhǎng) 使 , 由同階單步法提供初值 , ， ， 。計(jì)算 時(shí)，可取 。 理論分析得出它們的絕對(duì)穩(wěn)定區(qū)域，修正Hamming法： ; 4階Adams預(yù)測(cè)-校正法： 其中,我們也

46、可以在教學(xué)演示上看出修正的預(yù)測(cè)校正格式的誤差非常的小。,多步法的相容性條件 k步法的一般形式為 其中 由對(duì)單步法的討論可知,當(dāng) 時(shí),方法階數(shù)至少為1。即對(duì) y＝1,y=x 應(yīng)精確成立。令 y ＝1，所以令y=x得,4.多步法的相容性、穩(wěn)定性和收斂性,所以我們稱為線性多步法的相容性條件。,

47、k步法需要k個(gè)出發(fā)值，而初值問題只提供了一個(gè)初值，在使用k步法時(shí)尚需要其它方法作補(bǔ)充k-1個(gè)出發(fā)值，今假定它們 是 ，當(dāng) 這k-1個(gè)值都 應(yīng)收斂于共同的極限y（a）即

48、 在討論多步法收斂性時(shí)我們總假定（9.3.12）成立 定義四：,多步法的收斂問題,的解 收斂于初值問題的解 y(x)。這里 ＝x＝a+nh. 定義五：稱多項(xiàng)式 為k步法的特征多項(xiàng)式。如果特征多項(xiàng)式的零點(diǎn)皆不大于1，且等 于1的零點(diǎn)是單重的，稱根條件成立。稱多項(xiàng)式 為第二特征多項(xiàng)式。 顯然根條件可以表示 定理三：k步法收斂的

49、充要條件為： （1）相容性條件成立。 （2）特征多項(xiàng)式的零點(diǎn)皆不大于1，且等于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步法

50、的相容性、收斂性、穩(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)定的。 記 得

51、 是常系數(shù)差分方程，其特征方程為 記它的 k個(gè)特征根為 并設(shè)它們是互異的。顯然根與 有關(guān)，不妨記為 注意到當(dāng) 時(shí)

52、 這時(shí)由特征方程得 由線性多步法的相容性條件得 是一個(gè)根。不妨設(shè), 差分方程的解為 其中系數(shù)由線性多步法的出發(fā)值確定。,另一方面，y(0)=1的試驗(yàn)方程的精確解為 ， 設(shè)多步法截?cái)嗾`差為 ，由此可得 ，我們稱 為主根，其

53、它根都為增根。 定義五：線性多步法的絕對(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 一階方程組 9.2 化高階方程為一階方程組 9.3 剛性方程

54、組,9.1 一階方程組,前面我們研究了單個(gè)方程y’=f的數(shù)值解法，只要把y和f理解為向量，那么，所提供的各種計(jì)算公式即可應(yīng)用到一階方程組的情形。 考慮一階方程組的出值問題 若采用向量的記號(hào)，記,則常微分方程組的出值問題可以表示為 前幾節(jié)我們主要討論了常微分方程的出值問題的數(shù)值解法。只要將 y 和f 改寫為向量，那么前面提供的各種計(jì)算公式即可適用于一階常微分方程組的出值問題。 Runge-Kut

55、ta方法 對(duì)于方程組 四級(jí)四階顯示Runge- Kutta公式,若寫成分量形式就是 i=1,2,…,m 為了幫助理解這一公式的計(jì)算過程，我們?cè)倏紤]兩個(gè)方程的特殊情形：,這時(shí)四階龍格－庫塔公式具有形式,其中,這是一步法，利用節(jié)點(diǎn) 上的值 ， ，由上式順序計(jì) 算 ，然后

56、代入即可求得 節(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)

57、 ，其中，中A的特征值 滿足條件<0(j=1,…,N),,,,,且 則稱 系統(tǒng)為剛性方程，稱s為剛性比,§7 習(xí)題與總結(jié),,1數(shù)值例題2數(shù)值練習(xí) 3總結(jié),1、數(shù)值例題,,我們已經(jīng)學(xué)習(xí)了很多數(shù)值算法，他們的效果到底如何呢？下面我們來分析一道例題，看看那些方法，就這個(gè)問題，最能接近真實(shí)值 求初值問題

58、的數(shù)值解，取h=0.1并同精確解 比較 （1）用歐拉法來計(jì)算這個(gè)初值問題（2）用各階的Runge—Kutta方法來計(jì)算這個(gè)初值問題（3）用四階的Adams定步長(zhǎng)，Adams定步長(zhǎng)預(yù)測(cè)校正方 法來 計(jì)算這個(gè)初值問題。然后將數(shù)值結(jié)果列成表格：,2、數(shù)值練習(xí),計(jì)算實(shí)習(xí)題：1.使用改進(jìn)的Euler算法求初值問題

59、 的數(shù)值解，取步長(zhǎng)h=0.1,并同精確解 比較 2. 利用4階Runge-Kutta算法求初值問題,,的數(shù)值解，取步長(zhǎng)：(1)h=0.1;(2)h=0.01.求數(shù)值解及精確解 3.使用4階定步長(zhǎng)Adams預(yù)測(cè)－校正算法求初值問題 的數(shù)值解，取h=0.1,并同精確解 比較,4.使用4階定步長(zhǎng)Adams預(yù)測(cè)－校正

60、算法求初值問題的數(shù)值解，使用h=0.5,,5．綜合題請(qǐng)用已學(xué)的數(shù)值方法來計(jì)算下面初值問題，看看那種精度最好， 的數(shù)值解。,總 結(jié),,,本章研究的是常微分方程初值問題的數(shù)值解法，我們從單步法、多步法兩大類來介紹經(jīng)典算法，并且對(duì)相容性、收斂性、穩(wěn)定性這些理論進(jìn)行系統(tǒng)的、全面介紹，為同學(xué)們學(xué)習(xí)偏微分方程數(shù)值解法作了較好的鋪墊。梯形法、

61、R-K方法、變步長(zhǎng)法以及Adams方法，都是一些高精度的方法，但是它們也有缺點(diǎn)，四階R-K方法必須要求f(x,y)有較好的光滑性，若光滑性差，還不如歐拉法或改進(jìn)歐拉公式。它的另一個(gè)缺點(diǎn)是計(jì)算量較大，需要好較多的機(jī)器時(shí)間（每一步需要計(jì)算4次f(x,y)），相比之下漢明方法節(jié)省計(jì)算量，但它是四步方法不是自開始的需要借助四階R-K方法提供初始值。 對(duì)數(shù)值方法的分析還涉及到局部截?cái)嗾`差，整體誤差，相容性、收斂性、穩(wěn)定性等概念，他們涉及

62、到步長(zhǎng)的選取，如步長(zhǎng)不合適，舍入誤差惡性增長(zhǎng)，結(jié)果會(huì)完全錯(cuò)誤。大家在實(shí)際應(yīng)用時(shí)一定要注意這些問題。,,剛性方程組具有很重要的應(yīng)用價(jià)值，它的理論和解法很多，并且還在進(jìn)一步發(fā)展，我們僅是簡(jiǎn)單介紹一下，詳細(xì)介紹請(qǐng)查看參考書目。,抽取客觀事物的本質(zhì)特征及其內(nèi)部聯(lián)系，用適當(dāng)?shù)臄?shù)學(xué)語言和工具所建立的一個(gè)數(shù)學(xué)結(jié)構(gòu),什么是數(shù)學(xué)模型,例如著名的Malthus人口模型,,,,又如弱肉強(qiáng)食模型,例1 用歐拉方法，隱式歐拉方法和歐拉中點(diǎn)公式計(jì)算,的近似解，取步

63、長(zhǎng)h=0.1，并與精確值比較,解:歐拉方法的算式為：,歐拉隱式方法在本題可解出方程，不必迭代，公式為：,歐拉中點(diǎn)公式是兩步法，第一步y(tǒng)1用歐拉公式，以后用公式,利用 Euler方法求初值問題,解 此時(shí)的 Euler公式為,例1,的數(shù)值解.此問題的精確解是,分別取步長(zhǎng)h=0.2 ,0.1 ,0.05,計(jì)算結(jié)果如下,,例2 .用歐拉公式和梯形公式的預(yù)估校正法計(jì)算：,的數(shù)值解，取h=0.1，梯形公式只迭代一次，并與精確值比較.方程的解析解為