第7章_常微分方程初值問(wèn)題的數(shù)值解法

1、第第7章章 常微分方程初值問(wèn)題的數(shù)值解法常微分方程初值問(wèn)題的數(shù)值解法7.1 引言引言7.2 尤拉方法尤拉方法7.3 龍格龍格庫(kù)塔法庫(kù)塔法7.4 收斂性和穩(wěn)定性收斂性和穩(wěn)定性7.1 引引 言言常微分方程（ordinary differential equation ODE)的求解： 分離變量法、齊次方程的求解、可降階高階微分方程求解特殊類型的微分方程。微分方程的近似解法：(1)近似解析法：逐次逼近法、級(jí)數(shù)解法(2)數(shù)值解法：求離散點(diǎn)上的近似值。定解問(wèn)題：微分方程定解條件(初值條件、邊界條件)分別稱為初值問(wèn)題和邊值問(wèn)題。00( , )(7 1)()(72)yf x yaxby xy 一階常微分方程

2、的初值問(wèn)題：1212( , )( , ),()( ,)( ,)(7-1)(7-2)f x yx y axbyx yyLipschitzLf x yf x yL yy 函數(shù)在帶形區(qū)域 為連續(xù)函數(shù)，且對(duì)任意的 滿足李普希茲條件，即存在常數(shù) ，使得由常微分方程理論知，初值問(wèn)題的解必存在且唯一。012(7-1)(7-2)naxbaxxxxb：在初值問(wèn)題解的區(qū)間上插入一系列節(jié)點(diǎn)在節(jié)點(diǎn)上用離散化方法將連續(xù)型微分方程轉(zhuǎn)化成離散型代數(shù)方程即數(shù)值解差分方法的程基本思想來(lái)求解。0,1,2,iy i ：尋求唯一存在的解在離散節(jié)點(diǎn)數(shù)值方法的任務(wù)上的近似值 ()。01112()()ny xy xyyyy：利用求出的近似

3、值 ，再由求出 ， ，直到求出 為止。該具體作法步進(jìn)式遞推算法稱為或式算法。微分方程離散化常用方法微分方程離散化常用方法11: ( , )(0,)11, nnnnxxxxdydxf x y dxndx用數(shù)值積分方法離散（數(shù)值積分化）1111,(), (),( , )(,) (,)(0,1,)nnnnnnxnnxnnnnyyy xy xf x y dxhf xyyyhf xyn用代替對(duì)右端積分采用取左端點(diǎn)的 則有矩形公式1(,)1()2()(, ()nnnnnnxynny xy xdyf xy xdxxx（ ）數(shù)值微分11111 ,(), () (,)(,)0, 1, 2,nnnnnnnnnnn

4、nnnhxxyy xyy xyyf xyhyyhf xyn用代替，則： 22 ()3 ()()()2 ()(, ()() ) (2nnnnnnnnnhy xhy xhxy xTayy xy xhy xhf xy xyrxlo（ ） 在 附近的展開：1 () ()()0,1,2, nnnnnnnhyy xy xyyhf x y xnTaylor取的線性部分，且得的近似值： 展開法不僅可得到求數(shù)值解的公式，且容易估計(jì)截?cái)嗾`差。7.2 尤拉方法尤拉方法1. Euler 公式1111( )( , ( )()()( , ( )(73)nnnnnnxxxxxnnxy t dtf t y t dty xy

5、xf t y t dt對(duì)右端積分采用取左端點(diǎn)矩形公式的,則有：00 ( , )(7 1)()(72)yf x yaxby xy 初值問(wèn)題：12( , ( )(, ()()2nnxnnnxhf t y t dth f xy xy21()2ihhRy當(dāng) 充分小，舍去誤差項(xiàng)，則：1(0,1,),nnnnxx1()()(, ()nnnny xy xh f xy x11,(), (),7-17-2nnnnyyy xy x用代替得出初值問(wèn)題（）（）的離散化形式：001() ( (,() (0 1 274)nnnnyy xyyh f xyn, , ,000000(,)0000 (,) ( )(,)(,) (

6、,) xyxy yy xdyf xydxf xyxy由出發(fā)取曲線的切線（存在)，則斜率由于及已知，(1) 幾何必有切意義線方程。00000000(,) ()() (,)xydyyyxxyxxf xydx由點(diǎn)斜式寫出切線方程：01hxxhy等步長(zhǎng)為 ，則，可由切線算出1000(,)yyhf xy11 (,) 0 1 2 nnnnnyy xxyyhf xyn逐步計(jì)算出（ ） 在點(diǎn)的值：，11()()( , ( )(73)nnxnnxy xy xf t y t dt（2）梯形公式11111(),()( (,)(,)2nnnnnnnnnnyy xyy xhyyf xyf xy積分用梯形公式，且令：則得

7、：31111()( )12nnnnnhRy xyyxx (0)1(1)( )111Euler012(,) (75) ( (,)(,) 0 1 22nnnnkknnnnnnnyyhf xyhyyf xyf xyk, , ,與法結(jié)合，形成迭代算法，對(duì)， ， ，隱式方法1ny2. Euler方法的截?cái)嗾`差(1) 一步局部截?cái)嗾`差在中產(chǎn)生的誤差而非累積誤差：1() Taylor nny xx將在 點(diǎn)展開：111()nnnRy xy1()nnxx21()()()(, ()( )2nnnnnhy xy xhy xhf xy xy1(,)nnnnyyhf xy22111()( )()2nnnhRy xyyO

8、 h11()pnRO h：數(shù)值方法的局部截?cái)嗾`差為：定義pp稱該方法具有或。顯然，尤拉階精度階方法一階精度公式具有，梯形公式具有二階精度。3. 改進(jìn)的尤拉方法梯形公式雖然提高了精度，但算法復(fù)雜。而在實(shí)際計(jì)算中只迭代一次，這樣建立的預(yù)測(cè)校正系統(tǒng)稱作改進(jìn)的尤拉公式。1111 (,) (,)(,)2nnnnnnnnnnyyhf xyhyyf xyf xy預(yù)測(cè)校正1 (,)(,(,)2nnnnnnnnhyyf xyf xh yhf xy也可表示為11(,)(,)(78)()/2pnnncnnpnpcyyhf xyyyhf xyyyy或者表示為下列平均化形式1112(01)(0)10.12 ()2()2

9、 ()1()2nnnnnnpnnnncnppnpcxyyxyyhxyyh yyxyyh yyxyyh yyyyy ：用尤拉公式和改進(jìn)的尤拉公式解初值問(wèn)題解：取步長(zhǎng)，尤拉公式為：改進(jìn)的尤拉公式為：例17.3 龍格龍格-庫(kù)塔（庫(kù)塔（R-K）法）法 考察改進(jìn)的尤拉法，可以將其改寫為：考察改進(jìn)的尤拉法，可以將其改寫為：),(),(2121121211hKyhxfKyxfKKKhyyiiiiii 尤拉法尤拉法iiiiyxfkhkyy,111局部截?cái)嗾`差局部截?cái)嗾`差O(hO(h2 2) )局部截?cái)嗾`差局部截?cái)嗾`差O(hO(h3 3) )增加計(jì)算增加計(jì)算f(x,y)f(x,y)在不同點(diǎn)的在不同點(diǎn)的 值，能否

10、提高局部截?cái)嗾`差的階？值，能否提高局部截?cái)嗾`差的階？所以）使得（有微分中值定理，存在)(,()()()(1 , 01hxyhxfhxyhxyxyiiiii）（ 9-7)(,()()(1hxyhxhfxyxyiiii 不同的計(jì)算公式?？梢缘玫讲捎貌煌品椒ǎ鸵虼?，平均斜率是一般無(wú)法求得）式是精確表達(dá)式，但（中值定理知上的平均斜率，有微分在可以看作這里kkhxyxxxhxyhxfhxyiiiiii_1_,)(9-7,y)(,()(啟示如果在區(qū)間內(nèi)多取幾個(gè)點(diǎn)的斜率值，然后把它們的線性組合作為平均斜率的近似值，則有可能構(gòu)造出更高精度的計(jì)算公式，這就是龍格-庫(kù)塔法的基本思想二階龍格-庫(kù)塔公式)(22

11、111kkhyyii),(1iiyxfk ),(12phkyphxfkii（7-10）選取適當(dāng)?shù)?.32,1hop使其截?cái)嗾`差為，121),(,yyxfkkkii做泰勒展開假定對(duì)）得代入（10-7)()()()(),(),(),(),(),(2 21112hxphyxyhyxfphkyxphfyxfphkyphxfphkyxfkiiiiyiixiiiiipi )()()()()(3 22211hxyphxyhxyyiiii 之間與介于！的二階泰勒展開式為在又iiiiiiiixxxxyxxxyxxxyxyxyxxy3 2 )(3)()(2)()()()(選取過(guò)程，有！時(shí)，當(dāng))()(3)(33 1

12、hxxyxxii成立二階精度只要下列條件使公式具有）的系數(shù)即可發(fā)現(xiàn)，要）和（比較（）（！)()(2)()()(3 21hxyhxhyxyxyiiii將滿足條件（7-11）式到一簇公式（7-10）統(tǒng)稱為二階龍格-庫(kù)塔公式這里有這里有3個(gè)未知數(shù)，個(gè)未知數(shù)，2 個(gè)方程。個(gè)方程。存在存在無(wú)窮多個(gè)解無(wú)窮多個(gè)解121212p7-1121hkyyii),(1iiyxfk )2,2(12khyhxfkii中點(diǎn)公式變形尤拉公式庫(kù)塔公式二階龍格改進(jìn)尤拉公式庫(kù)塔公式時(shí)，二階龍格-, 1, 0,21-21, 12121pp11 12233121312,iiiiiiiiyyhkkkkf x ykf xph yphkkf

13、 xqh yqh rksk同理可得三階龍格-庫(kù)塔公式112312131246,22,2iiiiiiiihyykkkkf x yhhkfxykkf xh yhkhk局部截?cái)嗾`差o(h4)標(biāo)準(zhǔn)四階龍格-庫(kù)塔公式112341213243226,22,22,iiiiiiiiiihyykkkkkfxyhhkfxykhhkfxykkfxh yhk局部截?cái)嗾`差o(h5) -20101xyxyyyxx用四階龍格 庫(kù)塔公式求解初值問(wèn)題從到（取h=0.2）各節(jié)點(diǎn)上的數(shù)值解100200130024003,1,0.91818222,0.90863722,0.843239kfxyhhkfxykhhkfxykkfxh y

14、hk得2,0.2xf x yyhy由000,1xy101234221.1832296hyykkkk所以例題2111211131124113,22,22,ykfxyhhkfxykhhkfxykkfxhyh k在 計(jì) 算211234221.3416676hyykkkk所以基點(diǎn)改進(jìn)尤拉法龍格-庫(kù)塔法精確解01110.21.1840961.1832291.1832160.41.3433601.3416671.341641這樣繼續(xù)下去，計(jì)算結(jié)果列于表步長(zhǎng)的自動(dòng)選擇單從第一步看，步長(zhǎng)h越小，局部截?cái)嗾`差就越??； 但隨著h的縮小，不但引起計(jì)算量的增加， 而且也引起舍入誤差的嚴(yán)重積累。所以選取合適的步長(zhǎng)h在實(shí)

15、際計(jì)算中是很重要的。 兩種情況處理適，分以下來(lái)判斷選取步長(zhǎng)是否合記得近似值計(jì)算兩步到，從然后折半取步長(zhǎng)，求出近似值記為出發(fā)，步長(zhǎng)以四階為例，從節(jié)點(diǎn))(12/1111,2/,hihihiiihiiyyyxxhyhx前的老值作為結(jié)果為止，這時(shí)取步長(zhǎng)加倍，直到長(zhǎng)加倍反復(fù)將步長(zhǎng)加倍，將步）如果（作為結(jié)果。的新值為止，取步長(zhǎng)折半后折半計(jì)算，直到反復(fù)將步長(zhǎng)，如果）對(duì)于給定的精度（,2,012/1hiyihxxi0對(duì)于任意固定的), 2 , 1(i0h當(dāng)）（iiy如果數(shù)值解0)(iixyy0)0(yyyy初值問(wèn)題尤拉法的收斂性), 0(且為常數(shù)相應(yīng)的尤拉公式為精確解為xeyxy0)(iiyhy)1 (1ii

16、hyy)1(0逐步遞推得，由于這里00 x，于是有ihxiihhihyy10)1(7.4收斂性與穩(wěn)定性時(shí)當(dāng)0h，ehh1)1 (iy所以尤拉公式的解0h當(dāng))(16-7ixy）式的精確解收斂到（控制。穩(wěn)定的，這是擾動(dòng)得以則稱該數(shù)值方法是絕對(duì)若及數(shù)值方法本身，僅依賴于擾動(dòng)）（值解的擾動(dòng)步有擾動(dòng)，此后各步數(shù)即第假設(shè)步數(shù)值解得擾動(dòng)。稱為第差值值為而實(shí)際計(jì)算得到的近似）的數(shù)值解（）處給出初值問(wèn)題（某一數(shù)值方法在節(jié)點(diǎn),.2, 1,.2, 1,0,2-71-7yiikiikiiyiiyxikikiiiiiy穩(wěn)定性滿足下述方程，于是擾動(dòng)其中尤拉公式尤拉法的穩(wěn)定性分析iiiiiihySyxhfyy)(),(1稱上式為擾動(dòng)方程，）（18-7)(1iiS時(shí)，有當(dāng)1)(Sii1ik由此歸納可得, 2, 1 iik區(qū)域定區(qū)域稱為尤拉公式的絕對(duì)穩(wěn)1)(Sh可見，適當(dāng)選擇的步長(zhǎng)值方法絕對(duì)穩(wěn)定成立，從而就可保證數(shù)可使1)(S/78. 20-,0/20hhh庫(kù)塔方法對(duì)于標(biāo)準(zhǔn)四階龍格無(wú)條件穩(wěn)定的此時(shí)稱梯形法是對(duì)于梯形公式對(duì)于尤拉公式例題應(yīng)滿足長(zhǎng)按絕對(duì)穩(wěn)定性要求，步這里解有何限制？考慮，對(duì)步長(zhǎng)穩(wěn)定性庫(kù)塔方法求解，從絕對(duì)用標(biāo)準(zhǔn)四階龍格對(duì)于初值問(wèn)題hyyy20-h-1t01)0(20算法不穩(wěn)定。則超出了上述范圍，固的，而取是滿足絕對(duì)穩(wěn)定性要求所以選取