1、1,第九章 常微分方程數(shù)值解法,第一節(jié) Euler方法,第三節(jié) 單步法的收斂性和穩(wěn)定性,第二節(jié) Runge-Kutta方法,上一頁(yè) 下一頁(yè) 返回,2,上一頁(yè) 下一頁(yè) 返回,本章介紹求解微分方程數(shù)值解的基本思想和方法.,含有自變量、未知函數(shù)和它的一階導(dǎo)數(shù)和高階導(dǎo)數(shù)的方程.,常微分方程,它是描述運(yùn)動(dòng)、變化規(guī)律的重要數(shù)學(xué)方法之一，分為兩類：,1. 初值問(wèn)題 即給出未知函數(shù)及導(dǎo)數(shù)在初始點(diǎn)的值；,2. 邊值問(wèn)題 即給出未知函數(shù)及（或）它的某些導(dǎo)數(shù)在區(qū)間兩個(gè)端點(diǎn)的值 。,3, 考慮一階常微分方程的初值問(wèn)題 :,只要 f (x, y) 在a, b R1 上連續(xù)，且關(guān)于 y 滿足 Lipschitz 條件，即

2、存在與 x, y 無(wú)關(guān)的常數(shù) L 使 對(duì)任意定義在 a, b 上的 y1(x) 和 y2(x) 都成立，則上述問(wèn)題解存在唯一解。,所謂數(shù)值解法就是要計(jì)算出初值問(wèn)題的解函數(shù) y(x) 在一系列離散點(diǎn) a = x0 x1 xN= b上的近似值： y0 , y1 ,yN .,節(jié)點(diǎn)間距 為步長(zhǎng)，,通常采用等距節(jié)點(diǎn)，即取 hi = h (常數(shù))。, yn 稱為問(wèn)題的數(shù)值解.,數(shù)值解所滿足的離散方程統(tǒng)稱為差分格式.,上一頁(yè) 下一頁(yè) 返回,4,第一節(jié) 歐拉方法,一、 歐拉公式,令yn為y(xn)的近似值，將上式代入(*)式可得,此式稱為歐拉(Euler)公式.,為Euler方法的局部截?cái)嗾`差.,上一頁(yè) 下一

3、頁(yè) 返回,5,例1 用歐拉公式解初值問(wèn)題,解： 取步長(zhǎng) h=0.1, 歐拉公式的具體形式為：,依次計(jì)算可得, ,y3，y4，y5，y6，y7，y8，y9，y10,上一頁(yè) 下一頁(yè) 返回,6,其部分結(jié)果見(jiàn)下表,可見(jiàn)Euler方法的計(jì)算結(jié)果精度不太高。,上一頁(yè) 下一頁(yè) 返回,7, 歐拉公式的幾何意義：,幾何意義：,用折線近似代替方程的解曲線,因而也稱Euler方法為折線法.,上一頁(yè) 下一頁(yè) 返回,8,二、 后退的歐拉公式,也用一階差商逼近導(dǎo)數(shù),令yn+1為y(xn+1)的近似值，則可得,稱為后退Euler公式,已知 yn時(shí),必須通過(guò)解方程才能求出yn1 ,這樣的公式稱為隱式公式, 而Euler公式為

4、顯式公式.,Euler公式和后退Euler公式都是由yn去計(jì)算yn+1 ，因此，稱 它們?yōu)閱尾椒ā?上一頁(yè) 下一頁(yè) 返回,9,顯然, p越大, 精度越高.,三、 局部截?cái)嗾`差與方法的階,Euler方法的精度,其中：,上一頁(yè) 下一頁(yè) 返回,10,所以， Euler方法具有 1 階精度。,將,在點(diǎn),處一階Taylor展開(kāi),上一頁(yè) 下一頁(yè) 返回,11,所以， 后退的Euler方法也具有 1 階精度。,隱式Euler方法的精度,上一頁(yè) 下一頁(yè) 返回,12, 顯、隱式兩種算法的平均, 歐拉公式的改進(jìn),其局部誤差為：,此公式具有2階精度.,稱平均公式或梯形公式,梯形公式可由下迭代式計(jì)算：其中迭代初值是Eu

5、ler公式提供.,上一頁(yè) 下一頁(yè) 返回,13,四、改進(jìn)的歐拉公式,注：此法亦稱為預(yù)測(cè)-校正法 。 可以證明該算法具有 2 階精度，同時(shí)可以看到它是個(gè)單步遞推格式，比隱式公式的迭代求解過(guò)程簡(jiǎn)單。它的精度高于顯式歐拉法。,上一頁(yè) 下一頁(yè) 返回,14,為了便于編程, 常將改進(jìn)的歐拉公式寫為：,上一頁(yè) 下一頁(yè) 返回,15,例2 用改進(jìn)的歐拉法解例1中的初值問(wèn)題.,解：取步長(zhǎng) h=0.1, 改進(jìn)歐拉法的具體 形式為,具體計(jì)算過(guò)程如下,上一頁(yè) 下一頁(yè) 返回,16,依次計(jì)算可得, ,y3，y4，y5，y6，y7，y8，y9，y10,其部分結(jié)果見(jiàn)下表,上一頁(yè) 下一頁(yè) 返回,17,例3 對(duì)下面的初值問(wèn)題,解 （

6、1）取步長(zhǎng) h=0.1, 歐拉方法的具體公式為,（2）取步長(zhǎng) h=0.1, 改進(jìn)的歐拉方法的具體公式為,取步長(zhǎng)h=0.1，分別用Euler方法、改進(jìn)的Euler方法求數(shù)值解。,上一頁(yè) 下一頁(yè) 返回,18,計(jì)算結(jié)果見(jiàn)下表,上一頁(yè) 下一頁(yè) 返回,19,第二節(jié) 龍格 - 庫(kù)塔法,基本思想,考察改進(jìn)的歐拉法，可以將其改寫為：,斜率 一定取k1 k2 的平均值嗎？,步長(zhǎng)一定是一個(gè)h 嗎？,只要能對(duì)平均斜率提供一種近似算法,就能得到一種對(duì)應(yīng)的差分格式.,上一頁(yè) 下一頁(yè) 返回,20,例如取 m 個(gè)點(diǎn)的斜率構(gòu)造如下形式的公式,該公式稱為m級(jí)龍格庫(kù)塔(Runge-Kutta)公式,簡(jiǎn)稱R-K公式.,求解：只需將

7、公式的局部截?cái)嗾`差在xn點(diǎn)進(jìn)行Taylor展開(kāi),令其前面盡可能多的項(xiàng)為0, 便可導(dǎo)出ai，bij，ci所滿足的方程組,即可從中求出這些系數(shù).,上一頁(yè) 下一頁(yè) 返回,21,以 m=2 的情形為例說(shuō)明建立R-K公式的方法.,其局部截?cái)嗾`差為：,上一頁(yè) 下一頁(yè) 返回,22,因此有：,而對(duì)于h3，若將k2的Taylor展開(kāi)式多取一項(xiàng),會(huì)發(fā)現(xiàn)h3項(xiàng)的系數(shù)不可能為0.,而對(duì)于上式有無(wú)窮多個(gè)解,它的每一組解都給出了一個(gè)局部截?cái)嗾`差為 的二級(jí)R-K公式,即二階R-K公式.,這里有 個(gè)未知數(shù)， 個(gè)方程。,3,2,上一頁(yè) 下一頁(yè) 返回,23, 常用的標(biāo)準(zhǔn)四階RK公式（經(jīng)典R-K方法）,最常用的四階標(biāo)準(zhǔn)RK公式（經(jīng)

8、典R-K方法）為：,上一頁(yè) 下一頁(yè) 返回,24,例 用四階標(biāo)準(zhǔn)R-K公式解初值問(wèn)題,解：取 h=0.2 ，四階標(biāo)準(zhǔn)R-K法的具體格式如下:,上一頁(yè) 下一頁(yè) 返回,25,已知,上一頁(yè) 下一頁(yè) 返回,26,至少具有四位有效數(shù)字.,比較:上節(jié)用改進(jìn)的Euler公式計(jì)算,取h=0.1，最多具有四位有效數(shù)字 。,上一頁(yè) 下一頁(yè) 返回,27,改進(jìn)的Euler公式每前進(jìn)一步只要計(jì)算兩次f 值,而4階R-K公式每前進(jìn)一步要計(jì)算四次f 值,但改進(jìn)的Euler法的步長(zhǎng)比4階R-K法的小一半,兩者計(jì)算總量差不多.,而4階R-K法的效果要比改進(jìn)的Euler法好., 由于龍格-庫(kù)塔法的導(dǎo)出基于泰勒展開(kāi)，故精度主要受解函

9、數(shù)的光滑性影響。對(duì)于光滑性不太好的解，最好采用低階算法而將步長(zhǎng)h 取小。,上一頁(yè) 下一頁(yè) 返回,28,第三節(jié) 單步法的收斂性與穩(wěn)定性, 收斂性 /* Convergency */,例：就初值問(wèn)題 考察歐拉顯式格式的收斂性。,解：該問(wèn)題的精確解為,歐拉公式為,對(duì)任意固定的 x = xi = i h ，有,上一頁(yè) 下一頁(yè) 返回,29, 穩(wěn)定性 /* Stability */,例：考察初值問(wèn)題 在區(qū)間0, 0.5上的解。 分別用歐拉法、隱式歐拉法和改進(jìn)的歐拉格式計(jì)算數(shù)值解。,1.0000 2.0000 4.0000 8.0000 1.6000101 3.2000101,1.0000 2.500010

10、1 6.2500102 1.5625102 3.9063103 9.7656104,1.0000 2.5000 6.2500 1.5626101 3.9063101 9.7656101,1.0000 4.9787102 2.4788103 1.2341104 6.1442106 3.0590107,What is wrong ?!,上一頁(yè) 下一頁(yè) 返回,30,一般分析時(shí)為簡(jiǎn)單起見(jiàn)，只考慮試驗(yàn)方程,常數(shù)l0，可以是復(fù)數(shù),當(dāng)步長(zhǎng)取為 h 時(shí)，將某算法應(yīng)用于上式，并假設(shè)在初值產(chǎn)生誤差 ，則若此誤差以后逐步衰減，就稱該算法相對(duì)于z= l h 絕對(duì)穩(wěn)定， z 的全體構(gòu)成絕對(duì)穩(wěn)定區(qū)域。,我們稱算法A 比算法B 穩(wěn)定，就是指 A 的絕對(duì)穩(wěn)定區(qū)域比 B 的大。,上一頁(yè) 下一頁(yè) 返回,31,例：考察隱式歐拉法,可見(jiàn)絕對(duì)穩(wěn)定區(qū)域?yàn)椋?注：一般來(lái)說(shuō)，隱式歐拉法的絕對(duì)穩(wěn)定性比同階的顯式法的好。,上一頁(yè) 下一頁(yè) 返回,32,一些單步法的絕對(duì)穩(wěn)定區(qū)間見(jiàn)下表,Euler方法 改進(jìn)的Euler方法 三階R-K法 四階R-K法 隱式Euler法 梯形法,-2 z 0