數(shù)值計算方法51

1、120001)2)12(2)1(21)(kjkkabjafabkTkTnkknknxfCabI0)()()()()(2)0(0bfafabTbadxxfA)(, 3 ,2 , 1k第五章 常微分方程數(shù)值解2第五章 常微分方程數(shù)值解 5.1 引言引言(基本求解公式基本求解公式) 5.2 Runge-Kutta法法 5.3 微分方程組和高階方程解法簡介微分方程組和高階方程解法簡介3本章要點：本章作業(yè)本章主要研究基于微積分數(shù)值解法的常微分方程數(shù)值解，主要方法有線性單步法中的Euler方法、Simpson方法、 Runge-Kutta方法高階微分方程和微分方程組的數(shù)值解法P208. 1. 3. 4.

2、7. 8. 10. 11. 12.4本章應(yīng)用題:驅(qū)逐艦在濃霧中搜索潛艇，其時發(fā)現(xiàn)潛艇在3英里的海面上，但潛艇立即下潛，驅(qū)逐艦速度兩倍于潛艇，且已知潛艇下潛后即以全速朝某一未知方向直線前進，問驅(qū)逐艦應(yīng)采取什么路線才能保證它會開過潛艇的上方以投放深水炸彈？提示取極坐標(biāo)，并以發(fā)現(xiàn)潛艇時潛艇的位置為原點反潛5 5.1 引言引言(基本求解公式基本求解公式)在工程和科學(xué)技術(shù)的實際問題中，常需要求解微分方程只有簡單的和典型的微分方程可以求出解析解而在實際問題中的微分方程往往無法求出解析解在高等數(shù)學(xué)中我們見過以下常微分方程：0)(),(yaybxayxfy )(,)(),(0ayyaybxayyxfy-(1)

3、-(2)6 nybyyaybxayyxfy)(,)(),(0-(3)(1),(2)式稱為初值問題,(3)式稱為邊值問題2002212210012111)(),()(),(yxyyyxfyyxyyyxfy-(4)另外,在實際應(yīng)用中還經(jīng)常需要求解常微分方程組:本課程主要研究問題(1)的數(shù)值解法,對(2)(4)只作簡單介紹我們首先介紹初值問題(1)的解存在的條件7定理1. 條件，即滿足如果Lipschitzyxf),(，均有使得正數(shù),baxL|),(),(|2121yyLyxfyxf的解存在且唯一則初值問題)1(對于問題(1),要求它的數(shù)值解)(,)(節(jié)點上的一系列離散點在區(qū)間就是求未知函數(shù)baxy

4、bxxxxan210),2 , 1()(nkyxykk的近似值上函數(shù)值的數(shù)值解就是問題而)1(),2 , 1(nkyk80)(),(yaybxayxfy-(1)從(1)的表達式可以看出,求它的數(shù)值解的關(guān)鍵在于數(shù)值計算問題)(xy中程或者它的等價的積分方xadttytfyxy)(,()(0的數(shù)值計算問題積分xadttytf)(,(而數(shù)值微分或數(shù)值積分問題我們都已經(jīng)學(xué)習(xí)過9一、基于數(shù)值微分的常微分方程數(shù)值解法0)(),(yaybxayxfy-(1)對于初值問題(1)在下列子區(qū)間上分別應(yīng)用兩點數(shù)值微分公式bxxxxan210為了討論方便,假設(shè)以下節(jié)點為等距節(jié)點khaxnabhk,10)()(1)(0

5、1xyxyhay)(20yh ,1xa)()(1)(011xyxyhxy)(20yh )()(1)(121xyxyhxy)(21yh ,21xx)()(1)(122xyxyhxy)(21yh )()(1)(1jjjxyxyhxy)(2jyh ,1jjxx)()(1)(11jjjxyxyhxy)(2jyh -(5)(一) Euler公式11由(5)式每組的前一半可得)()()(01ayhxyxy)(202yh )()()(112xyhxyxy)(212yh )()()(1jjjxyhxyxy)(22jyh -(6)1, 1 ,0nj),(1jjjjyxhfyy)(2)(21jjyhhe )(22

6、jxyh -(7)記)(jjxyy 1, 1 ,0nj其中11)(jjyxy),()(jjjyxfxy(6)和(7)式稱為求解初值問題(1)的(前進)Euler公式和誤差項12由(5)式每組的后一半可得)()()(11jjjxyhxyxy)(22jyh 記),(111jjjjyxhfyy)(2)(21jjyhhe )(212 jxyh)(jjxyy 其中11)(jjyxy-(8)-(9)(8)和(9)式稱為求解初值問題(1)的后退Euler公式和誤差項式形公此類公式稱為隱式右端含有注意)(,)8(1jy1, 1 ,0nj),()(111jjjyxfxy式形公此類公式稱為顯公式右端不含而前進)(

7、,1jyEuler13從(6)或(8)式不難看出,jjyy時只要用到前一個值在計算1這種類型的方法稱為單步格式或單步法Euler方法的幾何體現(xiàn):)(jjxyhy),(1jjjjyxhfyy前進Euler公式)(1jjxyhy),(111jjjjyxhfyy后退Euler公式0123456-1.5-1-0.500.511.50123456-1.5-1-0.500.511.50123456-1.5-1-0.500.511.50123456-1.5-1-0.500.511.50123456-1.5-1-0.500.511.514Euler1.m例1.公式求解初值問題用前進Euler1)0(102yx

8、yxyy解:yxyyxf2),(顯然1, 1,10,000ybnax由前進Euler公式),(1jjjjyxhfyynj,2 , 1)2(jjjjyxyhy1 . 0h取1500.10.20.30.40.50.60.70.80.9111.11.21.31.41.51.61.71.8得)2(00001yxyhyy1 . 1)1021( 1 . 01)2(11112yxyhyy1918. 1)1 . 11 . 021 . 1( 1 . 01 . 1依此類推,有,yx00.10.20.30.40.50.60.70.80.9111.11.21.31.41.51.61.71.8前 進 Eul er公 式

9、 精確 解 0 1.0000 0.1000 1.1000 0.2000 1.1918 0.3000 1.2774 0.4000 1.3582 0.5000 1.4351 0.6000 1.5090 0.7000 1.5803 0.8000 1.6498 0.9000 1.7178 1.0000 1.784816由于后退Euler公式是隱形公式,計算例1將很麻煩事實上大多數(shù)情況下用后退Euler公式都較困難),(111jjjjyxhfyy),(1jjjjyxhfyy11jjyyEuler的預(yù)測值公式得到如果用前進11jjyEulery公式計算代入后退然后將就可得到新的Euler公式0)(yay-

10、(10)1, 1 ,0nj此方法稱為預(yù)測校正系統(tǒng)17用Euler公式的預(yù)測校正系統(tǒng)求解例1.例2.解:由(10)式,有)2(00001yxyhyy1 . 1)1021( 1 . 01)2(11101yxyhyy0918. 1)1 . 11 . 021 . 1( 1 . 01)2(11112yxyhyy1827. 1)0918. 11 . 020918. 1( 1 . 00918. 1)2(22212yxyhyy1763. 1)1827. 12 . 021827. 1( 1 . 00918. 1Euler1.m18依此類推,得 0 1.0000 0.1000 1.0918 0.2000 1.17

11、63 0.3000 1.2546 0.4000 1.3278 0.5000 1.3964 0.6000 1.4609 0.7000 1.5216 0.8000 1.5786 0.9000 1.6321 1.0000 1.6819,yx比較不同的結(jié)果00.10.20.30.40.50.60.70.80.9111.11.21.31.41.51.61.71.8前 進 Eul er公 式改 進 Eul er公 式精 確值校正系統(tǒng)預(yù)測19(二) 常微分方程數(shù)值解的截斷誤差評價一個微分方程求解公式的標(biāo)準(zhǔn)當(dāng)然是其精度:)(11的差與計算值也就是精確值jjyxy)()(111heyxyjjj而在求解公式 中)

12、,(1jjjjyxhfyy誤差項1, 1 ,0nj)(jjjxyyy都是近似值，即一般步的誤差只能表示求解公式第1)()(111jyxyhejjj20定義1. 部截斷誤差步的局的求解公式第為計算稱jyyxyhejjjj)()(因為一般情況下,求解公式的每一步都存在誤差,因此有定義2. 且步的截斷誤差的求解公式第為計算設(shè),)(jyhejjkjjkhehE1)()(步的累計截斷誤差為該求解公式第則稱khEk)(點上的總體截斷誤差即該求解公式在kx定義3. )()(1pjhOhe誤差為若求解公式的局部截斷階精度則稱該求積公式具有p21-2-1.8-1.6-1.4-1.2-1-0.8-0.6-0.4-

13、0.2000.511.522.53210 xxxx)()(1122xyyxyy)(1he)(2he)(3he)(1hE)(2hE)(3hE22Euler公式的局部截斷誤差為)(2)(21jjyhhe 具有1階精度后退Euler公式的局部截斷誤差為)(2)(21jjyhhe 也具有1階精度顯然一個求解公式的精度越高,計算解的精確性也就越好從前面的分析可知,Euler法的精度并不算高因此有必要找尋精度更高的求解公式23二、基于數(shù)值積分的常微分方程數(shù)值解法0)(),(yaybxayxfy-(1)對于初值問題上積分對上式在區(qū)間,1kkxxkkxxkkdxyxfxyxy1),()()(1-(11)kkx

14、xkkdxyxfxyxy1),()(,)(1就要求計算積分則計算已知假設(shè)kkkkxxxxdxyxfdxy11),(24:),(1的計算kkxxdxyxf),(),(111kkxxyxhfdxyxfkk矩形求積公式)(,(),(2),(111kkkkxxxyxfyxfhdxyxfkk已知假設(shè)11)(kkyxy梯形求積公式,誤差為)(,()2(,2(4),(6),(11111kkkkkkxxxyxfhxyhxfyxfhdxyxfkkSimpson求積公式,誤差為將以上求積公式代入(11)式,并加以處理就可得到相對應(yīng)的求解公式)(,(12)(3kkkyfhTR )(,(16180)()4(5kkky

15、fhSR25(一) 矩形求解公式kkxxkkdxyxfxyxy1),()()(1),(),(111kkxxyxhfdxyxfkk由可得),()(111kkkkyxhfyxy令),(111kkkkyxhfyy-(12)(12)式稱為矩形公式(矩形法)實際上就是Euler求解公式26(二) 梯形求解公式kkxxkkdxyxfxyxy1),()()(1由)(,(),(2),(111kkkkxxxyxfyxfhdxyxfkk可得)(,(),(2)(111kkkkkkxyxfyxfhyxy),(),(2111kkkkkkyxfyxfhyy令-(13)稱(13)式為梯形求解公式(梯形法)注意:(13)式是

16、隱形公式27)(,(12)()(3kkkkyfhTRhe )()13(11kkxyy式中如果在則梯形公式第k步的截斷誤差為)(123kyh ,1kkkxx)(3hO顯然梯形法具有二階精度由于梯形公式為隱形公式,一般情況下不易顯化kkyyEuler的預(yù)測值求出矩形法公式可以先使用)12)(28即進行校正代入梯形公式然后將,)13(ky),(111kkkkyxhfyy),(),(2111kkkkkkyxfyxfhyy-(14)以上公式稱為改進的Euler求解公式(改進Euler法),即),(),(2111kkkkkkyxfyxfhyy),(,(),(2111111kkkkkkkyxhfyxfyxf

17、hy-(15)29例3.用Euler公式、梯形公式和改進Euler公式求解初值問題,并比較結(jié)果的精度1)0(5 . 0 ,0,yxydxdy解:1 . 0h取步長(1)Euler公式11jjjhyyy)(2)(12 jjxyhhe)(212jxyh122jyh001hyyy9 . 00212|)(|yhhe2105 . 030112hyyy81. 01222|)(|yhhe21045. 0223hyyy729. 02232|)(|yhhe210405. 0334hyyy6561. 03242|)(|yhhe2103645. 0445hyyy59049. 04252|)(|yhhe2103280

18、5. 0(2)梯形公式211kkkkyyhyyhyhykk2)2(1即312)()(hyhekk 312hyk311 . 29 . 101yy 31112|)(|hyhe904752. 05105397. 71 . 29 . 112yy 32212|)(|hyhe818585. 05108215. 61 . 29 . 123yy 33312|)(|hyhe740625. 05101719. 61 . 29 . 134yy 34412|)(|hyhe670089. 05105841. 51 . 29 . 145yy 35512|)(|hyhe606271. 05100523. 532(3)改進E

19、uler公式211kkkkyyhyy x y 0 1.00000.1 0.90500.2 0.81900.3 0.74120.4 0.67080.5 0.6071使用MATLAB軟件Euler2.m結(jié)果為11kkkhyyy339 . 02105 . 081. 021045. 0729. 0210405. 06561. 02103645. 05905. 021032805. 0904752. 05105397. 7818585. 05108215. 6740625. 05101719. 6670089. 05105841. 5606271. 05100523. 50.90500.81900.74120.67080.6071Euler公式j(luò)y|)(|hej梯形公式j(luò)y|)(|hej改進Euler公式j(luò)y結(jié)果比較Euler法的精度不如梯形公式3400.050.10.150.20.250.30.350.40.450.50.550.60.650.70.750.80.850.90.951Euler公 式改 進 Eul er公 式梯