計(jì)算方法-常微分方程

1、第八章第八章 常微分方程數(shù)值解常微分方程數(shù)值解/* Numerical Methods for Ordinary Differential Equations */y 為人口數(shù)量，為人口數(shù)量，r為人口增長(zhǎng)率，為人口增長(zhǎng)率， y0為時(shí)刻為時(shí)刻x0的人口數(shù)量。的人口數(shù)量。 000)()(yxyxxryy 求常微分方程求常微分方程解析解解析解的方法有多種多樣，但是利用這些方法，的方法有多種多樣，但是利用這些方法，我們只能對(duì)極少數(shù)特殊類(lèi)型的常微分方程求解；科學(xué)研究和我們只能對(duì)極少數(shù)特殊類(lèi)型的常微分方程求解；科學(xué)研究和工程技術(shù)上的大量常微分方程的求解需借助于工程技術(shù)上的大量常微分方程的求解需借助于數(shù)值計(jì)

2、算方法數(shù)值計(jì)算方法。 比如，描述人口增長(zhǎng)的著名比如，描述人口增長(zhǎng)的著名人口模型：人口模型：常微分方程常微分方程是常用的數(shù)學(xué)模型。是常用的數(shù)學(xué)模型。該方程有該方程有解析解解析解y(x)=y0er(x-xo)。 考慮考慮一階一階常微分方程的常微分方程的初值問(wèn)題初值問(wèn)題 /* Initial-Value Problem */: 0)(,),(yaybaxyxfdxdy只要只要 f (x, y) 在在a, b R1 上連續(xù)，且關(guān)于上連續(xù)，且關(guān)于 y 滿(mǎn)足滿(mǎn)足 Lipschitz 條條件件，即存在與，即存在與 x, y 無(wú)關(guān)的常數(shù)無(wú)關(guān)的常數(shù) L 使使對(duì)任意定義在對(duì)任意定義在 a, b 上的上的 y1(x

3、) 和和 y2(x) 都成立，則上述都成立，則上述IVP存存在唯一解在唯一解。| ),(),(|2121yyLyxfyxf 要計(jì)算出解函數(shù)要計(jì)算出解函數(shù) y(x) 在一系列節(jié)點(diǎn)在一系列節(jié)點(diǎn) a = x0 x1 xn= b 處的近似值處的近似值),., 1()(nixyyii 節(jié)點(diǎn)間距節(jié)點(diǎn)間距hi=xi+1-xi稱(chēng)為稱(chēng)為步長(zhǎng)步長(zhǎng)，當(dāng)，當(dāng) hi = h 為常數(shù)時(shí)稱(chēng)為為常數(shù)時(shí)稱(chēng)為等步長(zhǎng)等步長(zhǎng)。利用利用y0求節(jié)點(diǎn)求節(jié)點(diǎn)x1處的近似值處的近似值y1, 再?gòu)脑購(gòu)膟1來(lái)求出來(lái)求出y2, 直至直至求出所有的求出所有的yn. 稱(chēng)之為稱(chēng)之為步進(jìn)法步進(jìn)法. 1 歐拉方法歐拉方法 /* Eulers Method *

4、/ 歐拉公式：歐拉公式：xixi+1向前差商近似導(dǎo)數(shù)向前差商近似導(dǎo)數(shù)hxyxyxyiii)()()(1 )1,., 0(),(1 niyxfhyyiiii對(duì)于對(duì)于xia,b, 有有 )(,()(iiixyxfxy )(,()()(1iiiixyxhfxyxy 則則 有有 將上式中的函數(shù)值將上式中的函數(shù)值y(xi)都用近似值都用近似值yi來(lái)表示來(lái)表示, 則有數(shù)值計(jì)算則有數(shù)值計(jì)算格式格式7.1 Eulers Method例例7.1 求解初值問(wèn)題求解初值問(wèn)題 取步長(zhǎng)分別為取步長(zhǎng)分別為h=0.1和和h=0.05, 進(jìn)行計(jì)算進(jìn)行計(jì)算, 結(jié)果有結(jié)果有 1010 ,2)(yxyxyy解解: 該方程的解析解是

5、該方程的解析解是y=(1+2x)1/2. 歐拉格式是歐拉格式是)2(1iiiiiyxyhyy k xi y(xi) yi (h=0.1) yi (h=0.05) 012345.0.00.10.20.30.40.5.1.01.10001.19181.27741.35821.4351.1.01.09771.18761.27131.35021.4174.1.01.09541.18321.26491.34161.4142.定義定義在假設(shè)在假設(shè) yi = y(xi)，即第，即第 i 步計(jì)算是精確的前提下，考慮步計(jì)算是精確的前提下，考慮的截?cái)嗾`差的截?cái)嗾`差 Ri = y(xi+1) yi+1 稱(chēng)為稱(chēng)為局部

6、截?cái)嗾`差局部截?cái)嗾`差 /* local truncation error */。定義定義若某算法的局部截?cái)嗾`差為若某算法的局部截?cái)嗾`差為O(hp+1)，則稱(chēng)該算法有，則稱(chēng)該算法有p 階階精度。精度。 歐拉法的局部截?cái)嗾`差：歐拉法的局部截?cái)嗾`差：),()()()()()(32112iiiihiiiiiyxhfyhOxyxyhxyyxyR )()(322hOxyih 歐拉法具有歐拉法具有 1 階精度。階精度。Ri 的的主項(xiàng)主項(xiàng)/* leading term */ 局部截?cái)嗾`差：局部截?cái)嗾`差：7.1 Eulers Method 歐拉公式的改進(jìn)：歐拉公式的改進(jìn)： 隱式歐拉法隱式歐拉法 /* impli

7、cit Euler method */向后差商近似導(dǎo)數(shù)向后差商近似導(dǎo)數(shù)hxyxyxyiii)()()(11 xixi+1y(xi+1)yi+hf(xi+1, y(xi+1)1,., 0(),(111 niyxfhyyiiii由于未知數(shù)由于未知數(shù) yi+1 同時(shí)出現(xiàn)在等式的兩邊，不能直接得到，故同時(shí)出現(xiàn)在等式的兩邊，不能直接得到，故稱(chēng)為稱(chēng)為隱式隱式 /* implicit */ 歐拉公式，而前者稱(chēng)為歐拉公式，而前者稱(chēng)為顯式顯式 /* explicit */ 歐拉公式。歐拉公式。 隱式隱式歐拉法的局部截?cái)嗾`差：歐拉法的局部截?cái)嗾`差：11)(iiiyxyR)()(322hOxyih 即隱式歐拉公式具

8、有即隱式歐拉公式具有 1 階精度。階精度。7.1 Eulers Method 梯形公式梯形公式 /* trapezoid formula */ 顯、隱式兩種算法的顯、隱式兩種算法的平均平均)1,., 0(),(),(2111 niyxfyxfhyyiiiiii注：注：的確有局部截?cái)嗾`差的確有局部截?cái)嗾`差 ， 即梯形公式具有即梯形公式具有2 階精度，比歐拉方法有了進(jìn)步。階精度，比歐拉方法有了進(jìn)步。但注意到該公式是但注意到該公式是隱式隱式公式，不便于實(shí)際計(jì)算。公式，不便于實(shí)際計(jì)算。)()(311hOyxyRiii)1,., 0(),(1 niyxfhyyiiii)1,., 0(),(111 niy

9、xfhyyiiii)()(322hOxyRihi )()(322hOxyRihi 7.1 Eulers Method)(3hORi 改進(jìn)歐拉法改進(jìn)歐拉法 /* modified Eulers method */Step 1: 先用先用顯式顯式歐拉公式作歐拉公式作預(yù)測(cè)預(yù)測(cè)，算出，算出),(1iiiiyxfhyy Step 2: 再將再將 代入代入隱式隱式梯形公式的右邊作梯形公式的右邊作校正校正，得到，得到1 iy),(),(2111 iiiiiiyxfyxfhyy )1,., 0(),(,),(211 niyxfhyxfyxfhyyiiiiiiii7.1 Eulers Method 2/ )()

10、,(),(11cpipiiciiipyyyyxhfyyyxhfyy或者或者 例例7.2 求解初值問(wèn)題求解初值問(wèn)題 取步長(zhǎng)分別為取步長(zhǎng)分別為h=0.1和和h=0.05, 進(jìn)行計(jì)算進(jìn)行計(jì)算, 結(jié)果有結(jié)果有 1010 ,2)(yxyxyy解解: 該方程的解析解是該方程的解析解是y=(1+2x)1/2. 改進(jìn)歐拉格式改進(jìn)歐拉格式是是)2( ),22(211111iiiiiiiiiiiiiyxyhyyyxyyxyhyy k xi y(xi) yi(h=0.1) yi(h=0.05) 012345.0.00.10.20.30.40.5.1.01.09591.18411.26621.34321.4164.1

11、.01.09561.18351.26531.34211.4148.1.01.09541.18321.26491.34161.4142.7.1 Eulers Method2 龍格龍格 - 庫(kù)塔法庫(kù)塔法 /* Runge-Kutta Method */建立高精度的單步遞推格式。建立高精度的單步遞推格式。hxyxyii)()(1)(hxyi )(,(hxyhxfii *)()(,()()(hKxyhxyhxhfxyxyiiiii 1xnxn+1xn+hK*)(,(iiiixyxhfyy1考察歐拉公式考察歐拉公式 K* K1=f(xi, yi)改進(jìn)歐拉公式改進(jìn)歐拉公式 ),(),(11212112hK

12、yxfKyxfKKKhyyiiiiiiK* (K1+K2)/2斜率斜率一定取一定取K1 K2 的的平均值平均值嗎嗎？步長(zhǎng)一定是步長(zhǎng)一定是一個(gè)一個(gè)h 嗎嗎？這里稱(chēng)為函這里稱(chēng)為函數(shù)數(shù)y(x)在區(qū)間在區(qū)間xn, xn+1上的上的平均斜率平均斜率. 幾何意義是幾何意義是y(x)在區(qū)間在區(qū)間xn, xn+1上點(diǎn)上點(diǎn)xn+h處的處的斜率斜率. )(,(*hxyhxfKii 根據(jù)微分中值定根據(jù)微分中值定理理, 存在存在01 2 Runge-Kutta Method首先希望能確定系數(shù)首先希望能確定系數(shù) 1、 2、p，使得到的算法格式有，使得到的算法格式有2階階精度，即在精度，即在 的前提假設(shè)下，使得的前提假設(shè)

13、下，使得 )(iixyy )()(311hOyxyRiii Step 1: 將將 K2 在在 ( xi , yi ) 點(diǎn)作點(diǎn)作 Taylor 展開(kāi)展開(kāi))(),(),(),(),(2112hOyxfphKyxphfyxfphKyphxfKiiyiixiiii )()()(2hOxyphxyii 將改進(jìn)歐拉法推廣為：將改進(jìn)歐拉法推廣為：),(),(12122111phKyphxfKyxfKKKhyyiiiiii ),(),(),(),(),(),()(yxfyxfyxfdxdyyxfyxfyxfdxdxyyxyx Step 2: 將將 K2 代入公式，得到代入公式，得到 )()()()()()()

14、()(322212211hOxyphxyhyhOxyphxyxyhyyiiiiiiii 2 Runge-Kutta MethodStep 3: 將將 yi+1 與與 y( xi+1 ) 在在 xi 點(diǎn)的點(diǎn)的泰勒泰勒展開(kāi)作比較展開(kāi)作比較)()()()(322211hOxyphxyhyyiiii )()(2)()()(321hOxyhxyhxyxyiiii 要求要求 ，則必須有：，則必須有：)()(311hOyxyRiii21,1221 p 這里有這里有 個(gè)未知個(gè)未知數(shù)，數(shù)， 個(gè)方程。個(gè)方程。32存在存在無(wú)窮多個(gè)解無(wú)窮多個(gè)解。所有滿(mǎn)足上式的格式統(tǒng)稱(chēng)為。所有滿(mǎn)足上式的格式統(tǒng)稱(chēng)為2階龍格階龍格 - 庫(kù)

15、庫(kù)塔格式塔格式。21, 121 p注意到，注意到， 就是改進(jìn)的歐拉法。就是改進(jìn)的歐拉法。 Q: 為獲得更高的精度，應(yīng)該如何進(jìn)一步推廣？為獲得更高的精度，應(yīng)該如何進(jìn)一步推廣？其中其中 i ( i = 1, , m )， i ( i = 2, , m ) 和和 ij ( i = 2, , m; j = 1, , i 1 ) 均為待定均為待定系數(shù)，確定這些系數(shù)的系數(shù)，確定這些系數(shù)的步驟與前面相似。步驟與前面相似。 2 Runge-Kutta Method).,(.),(),(),(.1122112321313312122122111 mm mmmmimiiiiiimmiihKhKhKyhxfKhKh

16、KyhxfKhKyhxfKyxfKKKKhyy 常用的三階常用的三階龍格龍格-庫(kù)塔法庫(kù)塔法 /* Classical Runge-Kutta Method */ ：),(),(),()(213121321122246hKhKyhxfKKhyhxfKyxfKKKKhyyiiiiiiii2 Runge-Kutta Method 常用的四階常用的四階龍格龍格-庫(kù)塔法庫(kù)塔法 ：),(),(),(),()22(34222312221432161hKyhxfKKyxfKKyxfKyxfKKKKKyyiihihihihiiihii 例例7.3 用四階龍格庫(kù)塔公式求解例題用四階龍格庫(kù)塔公式求解例題7.1 取取

17、h=0.1計(jì)算計(jì)算, 結(jié)果有結(jié)果有解解: 對(duì)于例題對(duì)于例題7.1有以下四階龍格庫(kù)塔計(jì)算格式有以下四階龍格庫(kù)塔計(jì)算格式k xi y(xi) yi K1 K2 K3 K4 334223112143211)(22222222)22(6hKyhxhKyKKhyhxKhyKKhyhxKhyKyxyKKKKKhyynnnnnnnnnnnnnn , ,012345.0.00.10.20.30.40.5.1.00.912870.845120.790570.74536.0.954760.878180.817470.767810.72622.1.01.09541.18321.26491.34161.4142.0.

18、952290.876050.815620.766200.72481.0.912620.844940.790390.745200.70697.1.01.09541.18321.26491.34161.4142.2 Runge-Kutta Method3 微分方程組與高階方程微分方程組與高階方程 /* Systems of Differential Equations and Higher-Order Equations */ 一階微分方程組一階微分方程組IVP的一般形式為：的一般形式為： )(,.),(,()(.)(,.),(,()(1111xyxyxfxyxyxyxfxymmmm初值初值0002020101)(,.,)(,)(mmyxyyxyyxy 將問(wèn)題記作將問(wèn)題記作向量形式向量形式，令：，令： 001011.,.,.mmmyyyfffyyy 00)(),()(yxyyxfxy前述所有公式皆適前述所有公式皆適用于向量形式。用于向量形式。 改進(jìn)歐拉公式的形式改進(jìn)歐拉公式的形式 3 Systems of DEs and