常微分方程的數(shù)值解法

1、第九章 常微分方程的數(shù)值解法,§1、引言 §2、初值問題的數(shù)值解法--單步法 §3、龍格-庫塔方法 §4、收斂性與穩(wěn)定性 §5、初值問題的數(shù)值解法―多步法 §6、方程組和剛性方程 §7、習題和總結(jié),,主 要 內(nèi) 容,主 要 內(nèi) 容,研究的問題數(shù)值解法的意義,§1、 引 言,現(xiàn)實世界中大多數(shù)事物,內(nèi)部聯(lián)系非常復雜,找出其狀態(tài)和狀態(tài)變化規(guī)律

2、之間的相互聯(lián)系，也即一個或一些函數(shù)與他們的導數(shù)之間的關系,此種關系的數(shù)學表達就為,,微分方程,1.什么是微分方程 ?,其狀態(tài)隨著時間、地點、條件的不同而不同,如何利用數(shù)值方法求解微分方程（組）的問題。,2.數(shù)值求解微分方程的意義,如何建立數(shù)學模型已在建模課程中得到討論，各類微分方程本身和他們的解所具有的特性已在常微分方程及數(shù)學物理方程中得以解釋，,本章專門討論,尋找解析解的過程稱為求解微分方程組。,一個或一組具有所要求階連續(xù)導數(shù)的

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

4、稱為微分方程在該范圍內(nèi)的,數(shù)值解,在大量的實際方程中出現(xiàn)的函數(shù)起碼的連續(xù)性都無法保證，更何況要求階的導數(shù),求解數(shù)值解,很多微分方程根本求不到問題的解析解！,重要手段。,,,常微分方程的數(shù)值解法常用來求近似解,根據(jù)提供的算法,通過計算機,便捷地實現(xiàn),5.常微分方程數(shù)值解法的特點,,,,數(shù)值解法得到的近似解(含誤差)是一個離散的函數(shù)表.,本章主要討論一階常微方程的初值問題,6.基本知識,其中f (x,y)是已知函數(shù)，(1.2

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

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

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

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

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

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

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

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

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

14、推法的基本思想是從 ( xi , yi ) 點出發(fā)，以某一斜率沿直線達到 ( xi+1 , yi+1 ) 點。歐拉法及其各種變形所能達到的最高精度為2階。,建立高精度的單步遞推格式。,首先希望能確定系數(shù) ?1、?2、p，使得到的算法格式有2階精度，即在 的前提假設下，使得,Step 1: 將 K2 在 ( xi , yi ) 點作 Taylor 展開,,Step 2: 將 K2 代入第1式，得到,

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

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

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

18、 在區(qū)間[0, 0.5]上的解。分別用歐拉顯、隱式格式和改進的歐拉格式計算數(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 ??!,一般分析時為簡單起見，只考慮試驗方程 /* test equation */,,常數(shù)，可以是復數(shù),,例：考察隱式歐拉法,,可見絕對穩(wěn)定區(qū)域為：,,注：一般來說，隱式歐拉法的絕對穩(wěn)定性比同階的顯式法的好。,例：隱式龍格-庫塔法,而顯式

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

21、梯形公式的預估校正法計算：,的數(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. 變步長方法,,龍格－

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

23、中有一個自由參數(shù)，方程組有無窮多解，從而得到一族2級2階R-K方法 。,3. 經(jīng)典Runge-Kutta方法,我們可以構(gòu)造出一族3級3階，一族4級4階和一族5級4階等R-K方法。最常用的4級4階是如下的經(jīng)典R-K方法：,,4．R-K-Fehhlberg 方法,R-K-Fehhlberg方法是在R-K方法的基礎上引進誤差和步長控制的辦法。即利用5階R-K法,,,估計4階R-K的局部誤差，其中,,注：R-K-Fehhlberg比4階R-

24、K方法具有更大的優(yōu)越性，他是計算穩(wěn)定，高精度的方法，他的顯著優(yōu)點是，每一步僅需計算f的6個值；若用4階R-K-L與5階R-K-L在一起使用，則每步需要計算f的10個值。推薦使用！,5. 隱式R-K方法,類似于顯式R-K公式，稍加改變，就得到隱式R-K方法。,,,它與顯式R-K公式的區(qū)別在于：顯式公式中對系數(shù)求和的上限是i-1，從而構(gòu)成的矩陣是一個嚴格下三角陣。而在隱式公式中對系數(shù)求和的上限是L，從而構(gòu)成的矩陣是方陣，需要用迭代法求出Ki

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

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

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,,,,,,,,,,,,,,表,（注：已指出過準確解,）,單步法的相容性、收斂性和穩(wěn)定性,1.單步法的相容性 2.單步法的收斂性 3.單步法的穩(wěn)定性 4.

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

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

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

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

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

34、誤差不增長則,則Euler 方法的絕對穩(wěn)定域是以 為半徑的圓的內(nèi)部。為了保證穩(wěn)定性步長有所控制。假如當 時h應滿足 ，當 時， h 應滿足 等等。同樣我們可以將試驗方程用到其它各種單步法當中，求出其絕對穩(wěn)定域。在實際應用中必須在這個范圍內(nèi)，

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

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

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

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

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

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

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

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

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

44、-校正公式。 它們的局部截斷誤差分別是,Adams預測-校正公式,利用外推原理，將上兩式作線性組合消去局部截斷誤差的主項。使計算精度至少提高一階，同時得到兩個修正公式，將它們和上兩式構(gòu)成了如下修正Adams公式：預測： 修正： 校正:,修正：

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

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

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

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

49、充要條件為： （1）相容性條件成立。 （2）特征多項式的零點皆不大于1，且等于1的零點是單重的。 （稱2為）特征根條件。,多步法的穩(wěn)定性,關于線性多步法成立以下定理： 定理四：若函數(shù)f(x,y)對變量y滿足Lipschtiz 條件在與h,x無關的常數(shù)L，對區(qū)域D= 任意兩點（x，y1）,(x,y2)成立 k步法

50、的相容性、收斂性、穩(wěn)定性有以下關系 對于常微分方程右端函數(shù)f(x,y)，在相容性條件成立情況下，k步法的收斂性和穩(wěn)定性是等價的。 事實上在相容性條件成立時，收斂性和穩(wěn)定性的充要條件都是特征根條件成立。,多步法的絕對穩(wěn)定性和絕對穩(wěn)定域,將線性多步法的公式應用到試驗方程 進行考察。這里假定Re <0,即試驗方程本身是穩(wěn)定的。 記 得

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

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

53、它根都為增根。 定義五：線性多步法的絕對穩(wěn)定區(qū)域?qū)o定的 ，如果特征方程的特征根 皆按模小于1，則線性多步法關于u是絕對穩(wěn)定的。使得 成立的 構(gòu)成絕對穩(wěn)定區(qū)域。注：從誤差角度來看絕對穩(wěn)定區(qū)域的方法是一個理想的方法。這樣，絕對穩(wěn)定區(qū)域越大，方法適用性越廣，因而越優(yōu)越。,§6 方程組和剛性方程,9.1 一階方程組 9.2 化高階方程為一階方程組 9.3 剛性方程

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

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

56、代入即可求得 節(jié)點上的 。,9.2化高階方程為一階方程組,關于高階微分方程的初值問題，原則上總可以歸結(jié)為一階方程組來求解，例如，考察下列m 階微分方程初始條件為只要引進新的變量即可將m 階方程化為如下的一階方程組,初始條件則相應的化為,,9.3 剛性方程組,在求解方程組時，經(jīng)常出現(xiàn)解的分量數(shù)量級差別很大的情形，這給數(shù)值求解帶來很大困難，這種問題稱為剛性問題。若線性系統(tǒng)

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

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

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

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

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

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

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