第七章 常微分方程數(shù)值解法.

1、 第七章第七章 常微分方程初值問題的數(shù)值解法常微分方程初值問題的數(shù)值解法如何求解一階常微分方程的初值問題如何求解一階常微分方程的初值問題0( , ), ( )yf x yaxby ay 本章主要內(nèi)容本章主要內(nèi)容一、歐拉方法及其改進(jìn)形式一、歐拉方法及其改進(jìn)形式二、龍格二、龍格- -庫塔方法庫塔方法三三 、線性多步法、線性多步法7.17.1 引引 言言一、求常微分方程數(shù)值解的必要性一、求常微分方程數(shù)值解的必要性1、方程本身很復(fù)雜，難以求出解析解；、方程本身很復(fù)雜，難以求出解析解；2、即使可以獲得解析解，計(jì)算量太大或者計(jì)算、即使可以獲得解析解，計(jì)算量太大或者計(jì)算 過程太復(fù)雜而不實(shí)用；過程太復(fù)雜而不實(shí)

2、用；3、在實(shí)際應(yīng)用中，只需要求得解在某些特殊點(diǎn)、在實(shí)際應(yīng)用中，只需要求得解在某些特殊點(diǎn) 上的近似值。上的近似值。0)(),(yaybxayxfyyyLyxfyxf),(),(考慮一階常微分方程的初值問題考慮一階常微分方程的初值問題 (1) 則（則（1）的解存在且唯一）的解存在且唯一。 ，對任意，對任意),(yxfL,bax 滿足滿足Lipschitz條件，即存在常數(shù)條件，即存在常數(shù)，均有，均有若若二、初值問題及其二、初值問題及其數(shù)值解數(shù)值解的概念的概念常用的一些常用的一些解析解法解析解法：常數(shù)變易等常數(shù)變易等分離變量法、變量代換、分離變量法、變量代換、一階常微分方程初值問題：一階常微分方程初值

3、問題：000( , )();dyf x ydxy xyxx 數(shù)值解數(shù)值解是指：在解的是指：在解的存在區(qū)間存在區(qū)間上取一系列點(diǎn)上取一系列點(diǎn)012.nxxxx1 2 3(, , ,.)iy i 逐個(gè)求出的解逐個(gè)求出的解 的近似值的近似值()iy x0;ixxih等距等距節(jié)點(diǎn)：節(jié)點(diǎn)：:h步長步長定義：定義：(1) 初值問題初值問題(1)(1)的解析解及其數(shù)值解的幾何意義：的解析解及其數(shù)值解的幾何意義：oxy初值問題初值問題(1)(1)的解析解表示過點(diǎn)的解析解表示過點(diǎn) 的一條曲線的一條曲線00(,)xynx (,)nnxy ),(00yx( )yy x 0 x),(00yx 2x),(22yx 1x)

4、,(11yx 初值問題初值問題(1)(1)的數(shù)值解表示一組離散點(diǎn)列的數(shù)值解表示一組離散點(diǎn)列(,)iixy可用擬合方法求該組數(shù)據(jù)可用擬合方法求該組數(shù)據(jù) 的近似曲線的近似曲線(,)iixy積分曲線積分曲線 7.2 歐拉歐拉(Euler)方法方法一、歐拉方法的推導(dǎo)一、歐拉方法的推導(dǎo)在常微分方程中直接用在常微分方程中直接用向前差商向前差商代替代替微商微商.設(shè)在點(diǎn)設(shè)在點(diǎn) xi 列出方程列出方程(1):)(,()(iiixyxfxy并用并用向前差商向前差商,)()()()(111hxyxyxxxyxyiiiiii則有則有若用若用y(xi)的近似值的近似值yi代入上式右端，并記所得結(jié)果為代入上式右端，并記所

5、得結(jié)果為yi+1, 這樣得到計(jì)算公式：這樣得到計(jì)算公式：)(,()()(1iiiixyxhfxyxyyi+1= yi+h f (xi, yi) (2)i=0,1,2,當(dāng)初值當(dāng)初值 y0已知時(shí)已知時(shí),利用利用(2)可以遞推求出可以遞推求出12,.,.nyyy稱稱(2)式為式為歐拉方法歐拉方法.歐拉方法也稱為歐拉方法也稱為歐拉折線法歐拉折線法。oxy),(00yx( )yy x 0 x),(00yx nx (,)nnxy 2x),(22yx 1x),(11yx 00(,)xy現(xiàn)在從現(xiàn)在從 出發(fā)出發(fā)作解曲線的切線，作解曲線的切線，切線方程為：切線方程為：0000(,)()yyf x yx x 切線切

6、線斜率斜率為為00( ,)f x y該切線與直線該切線與直線 的的交點(diǎn)的縱坐標(biāo)交點(diǎn)的縱坐標(biāo)1xx 11( )yy x 100010(,)()yyf xyxx 11(,)x y再從再從 出發(fā)，以出發(fā)，以 為斜率作直線為斜率作直線11(,)f x y1111(,)()yyf x yxx 二、歐拉法的幾何解釋二、歐拉法的幾何解釋該直線與該直線與直線直線 的交點(diǎn)的的交點(diǎn)的縱坐標(biāo)縱坐標(biāo)2xx 22()yy x 211121(,)()yyf x yxx 依次類推，依次類推， 的近似值的近似值 的計(jì)算公式：的計(jì)算公式：1()ny x 1ny 1(,)nnnnyyhf xy 為為等距節(jié)點(diǎn)等距節(jié)點(diǎn)的步長的步長h

7、稱上述公式為稱上述公式為歐拉公式歐拉公式。歐拉方法的幾何意義：歐拉方法的幾何意義：用一條折線近似代替積分曲線用一條折線近似代替積分曲線所以，歐拉方法又稱所以，歐拉方法又稱折線法折線法歐拉方法歐拉方法的缺陷：誤差比較大！的缺陷：誤差比較大！三、歐拉方法的積分學(xué)解釋三、歐拉方法的積分學(xué)解釋11( ) d( , ( ) dkkkkxxxxy xxf x y xx 即即11( , ( ) dkkxkkxyyf x y xx 利用左矩形公式近似積分項(xiàng)，有利用左矩形公式近似積分項(xiàng)，有21(, ()()kkkkyyhf xy xO h 在區(qū)間在區(qū)間 上對微分方程積分得上對微分方程積分得1,kkxx 略去高階

8、項(xiàng)，便得歐拉公式略去高階項(xiàng)，便得歐拉公式1(,)kkkkyyhf xy 四、顯式格式和隱式格式四、顯式格式和隱式格式若在微分方程若在微分方程 中用向后差商代替微商，則中用向后差商代替微商，則( , )yf x y 1111(,),iiiiiiyyf xyxx 或或1111(,), ()iiiiiiyyhf xyhxx此公式此公式(3)稱為稱為向后的歐拉公式。向后的歐拉公式。注：注：向后的歐拉公式是關(guān)于向后的歐拉公式是關(guān)于 的隱式表達(dá)式，這類的隱式表達(dá)式，這類1iy 格式稱為是格式稱為是隱式格式隱式格式；而歐拉公式是關(guān)于；而歐拉公式是關(guān)于 的顯式的顯式1iy 表達(dá)式，稱為是表達(dá)式，稱為是顯式格式

9、顯式格式。(3)利用右矩形公式近似積分項(xiàng)也可導(dǎo)出利用右矩形公式近似積分項(xiàng)也可導(dǎo)出隱式歐拉格式。隱式歐拉格式。隱式隱式格式格式五、梯形法五、梯形法11()()( , )nnxnnxy xy xf x y dx 112 (, ()(, ()nnnnhf xy xf xy x在區(qū)間在區(qū)間 上對方程上對方程(1) 兩端積分，并用兩端積分，并用梯形梯形1,nnxx 然后用然后用 代替代替 ，得，得ny()ny x1112 (,)(,)nnnnnnhyyf xyf xy 公式公式近似積分：近似積分：(4)此公式此公式(4)稱為稱為梯形法梯形法。將將EulerEuler公式與隱式公式與隱式EulerEule

10、r公式做算術(shù)平均，也可得出梯形公式公式做算術(shù)平均，也可得出梯形公式。),(),(211nnnnnyxfyxfhy 梯形公式與梯形公式與EulerEuler公式相比要精確的多，但公式相比要精確的多，但是梯形公式的計(jì)算量要大一些。每步計(jì)算要解一是梯形公式的計(jì)算量要大一些。每步計(jì)算要解一個(gè)關(guān)于個(gè)關(guān)于 的非線性方程，從而要用如下迭代的非線性方程，從而要用如下迭代公式：公式：1ny ),(),(21111knnnnnknyxfyxfhyy, 2, 1, 0k取初值為取初值為 ，反復(fù)迭代，即反復(fù)迭代，即 nnyy011ny011223 01ny 11ny 21ny， ， 若序列若序列 收斂于收斂于 ，當(dāng)，

11、當(dāng) 時(shí)，得到：時(shí)，得到： 01kkny *1nyk ),(),(2*11*1nnnnnnyxfyxfhyy則取則取 為第為第 個(gè)近似值。個(gè)近似值。 1ny *1ny1n如此迭代下去如此迭代下去得到迭代序列得到迭代序列： kny1， knknyy111在實(shí)際計(jì)算中，通常要求滿足在實(shí)際計(jì)算中，通常要求滿足為終止條件，此時(shí)取為終止條件，此時(shí)取 作為作為 的近似值的近似值 。 11kny)(1nty1ny六、歐拉預(yù)估六、歐拉預(yù)估校正公式校正公式0101112( )( )(,)(,)(,)nnnnnnnnnnyyhf xyhyyf xyf xy 將梯形公式和歐拉公式聯(lián)合使用：將梯形公式和歐拉公式聯(lián)合使用

12、：即先通過歐拉格式即先通過歐拉格式預(yù)估預(yù)估一個(gè)粗糙的近似值，再利用一個(gè)粗糙的近似值，再利用梯形公式對其進(jìn)行梯形公式對其進(jìn)行校正。校正。這種這種預(yù)估預(yù)估- -校正格式校正格式稱為稱為改進(jìn)的歐拉方法。改進(jìn)的歐拉方法。它可表為如下嵌套形式：它可表為如下嵌套形式：(5)或表為下列平均化形式：或表為下列平均化形式：(6)121211,),(2hKyhxfKyxfKKKhyynnnnn),(,(),(211nnnnnnnnyxhfyxfyxfhyy解：解：201( )dyxydxyy 01( , )x 例例1：分別利用歐拉公式、改進(jìn)的分別利用歐拉公式、改進(jìn)的歐拉公式求解下列初值問題歐拉公式求解下列初值問題

13、的數(shù)值解（取步長為的數(shù)值解（取步長為 ）0 1 .h 歐拉公式：歐拉公式：1(,)nnnnyyhf xy 0 21 1.nnnxyy 0.1 0.2 0.3 0.4 0.5 1.1000 1.1918 1.2774 1.3582 1.4351nxny 0.6 0.7 0.8 0.9 1.0 1.5090 1.5803 1.6498 1.7178 1.7848nxny改進(jìn)的歐拉公式：改進(jìn)的歐拉公式： 0.1 0.2 0.3 0.4 0.5 1.0959 1.1841 1.2662 1.3434 1.4164nxny 0.6 0.7 0.8 0.9 1.0 1.4860 1.5525 1.6165

14、 1.6782 1.7379010111010 21 1220 12( )( )( ).()()nnnnnnnnnnnnxyyyxxyyyyyy nxny 0.1 1.1000 1.0959 1.0954 0.2 1.1918 1.1841 1.1832 0.3 1.2774 1.2662 1.2649 0.4 1.3582 1.3434 1.3416 0.5 1.4351 1.4164 1.4142 0.6 1.5090 1.4860 1.4832 0.7 1.5803 1.5525 1.5492 0.8 1.6498 1.6165 1.6125 0.9 1.7178 1.6782 1.67

15、33 1.0 1.7848 1.7379 1.7321 nx)()()(nnnxyEuleryEulery改進(jìn)注：本例有解析解注：本例有解析解xy21 六、六、 誤差分析誤差分析 在分析一個(gè)數(shù)值方法的誤差時(shí)在分析一個(gè)數(shù)值方法的誤差時(shí),為了簡化分析為了簡化分析, 一一 般是假定般是假定 yn = y(xn)精確成立的前提下，估計(jì)誤差精確成立的前提下，估計(jì)誤差111)(nnnyxye這種誤差稱為這種誤差稱為局部截?cái)嗾`差局部截?cái)嗾`差。定義定義如果一種數(shù)值方法的局部截?cái)嗾`差為如果一種數(shù)值方法的局部截?cái)嗾`差為O( h p+1 ),則稱該數(shù)值方法的精度是則稱該數(shù)值方法的精度是 p 階階的。的。首先估計(jì)首先

16、估計(jì)Euler方法方法的局部截?cái)嗾`差。的局部截?cái)嗾`差。假定假定 yn = y(xn)，則，則Euler公式可以寫成公式可以寫成)( )()(,()(1nnnnnnxyhxyxyxhfxyy假定解假定解 y(x)是足夠光滑，從而可將是足夠光滑，從而可將 y(xn+1)在點(diǎn)在點(diǎn) xn作作Taylor展開：展開：)()(21)( )()()(321hOxyhxyhxyhxyxynnnnn于是有于是有)()()( 21)(23211hOhOxyhyxynnn由定義知，由定義知，Euler方法的精度是一階方法的精度是一階)(21)( )()()(11yxhyxyhhxyxynnnn由由)(21)( )(

17、)(11yxhyxyxynnn)(21)(,()(211yhxyxhfxynnn(將將y(xn)在點(diǎn)在點(diǎn)xn+1處展開處展開)據(jù)此可知據(jù)此可知隱式隱式Euler方法的精度也是一階。方法的精度也是一階?？梢宰C明可以證明改進(jìn)的改進(jìn)的Euler方法方法是二階精度的方法是二階精度的方法。121211,),(2hKyhxfKyxfKKKhyynnnnn（6） 假定解函數(shù)假定解函數(shù)y(x)和函數(shù)和函數(shù)f(x, y(x)足夠光滑，并假定足夠光滑，并假定yn=y(xn), 則有則有),()(,()(nnnnnyxfxyxfxy)()(,)(,)(nnnnnnxyxyxyfxyxxfxy)(,()(,()(,(nnynnnnxxyxfxyxfxyxf)()(,(),(1nnnnnxyxyxfyxfK)(,(12hKxyhxfKnn)()(,()(,()(,(