數(shù)值分析第五講：常微分方程數(shù)值解

1、第五章：常微分方程數(shù)值解第五章：常微分方程數(shù)值解5.1 5.1 引言引言1 1、常微分方程與解、常微分方程與解為為n n階常微分方程階常微分方程。0 ), , ,()(nyyyyxf如果函數(shù)如果函數(shù) 在區(qū)間在區(qū)間a,ba,b內內n n階可導，稱方程階可導，稱方程)(xyy )(xyy 滿足方程的函數(shù)滿足方程的函數(shù)稱為微分方程的稱為微分方程的解解。則則如如為任意常數(shù)）為任意常數(shù)）xy2 ccxy(, 2一般稱為方程的一般稱為方程的通解通解。為方程的解。為方程的解。12 xy如果如果則有則有10 )(y為方程滿足定解條件的解。為方程滿足定解條件的解。第五章：常微分方程數(shù)值解第五章：常微分方程數(shù)值解

2、 102)(yxy ccxy1212 xy方程的通解方程的通解滿足定解條件的解滿足定解條件的解微分關系（方程）微分關系（方程）解的圖示解的圖示第五章：常微分方程數(shù)值解第五章：常微分方程數(shù)值解本教材重點討論定解問題本教材重點討論定解問題( (初值問題）初值問題）定解條件（初始條件）定解條件（初始條件） 00yxyyxfy)(),(),(yxf是否能夠找到定解問題的解取決于是否能夠找到定解問題的解取決于僅有極少數(shù)的方程可以通過僅有極少數(shù)的方程可以通過“常數(shù)變易法常數(shù)變易法”、“可分可分離變量法離變量法”等特殊方法求得初等函數(shù)形式的解，絕大等特殊方法求得初等函數(shù)形式的解，絕大部分方程至今無法理論求解

3、。部分方程至今無法理論求解。如如xyxyeyxyxyy 212,),sin(sin等等等等2 2、數(shù)值解的思想、數(shù)值解的思想第五章：常微分方程數(shù)值解第五章：常微分方程數(shù)值解（1 1）將連續(xù)變量）將連續(xù)變量 離散為離散為,bax bxxxxank 10nkxyykk,)(21 （2 2）用代數(shù)的方法求出解函數(shù)）用代數(shù)的方法求出解函數(shù) 在在 點的近似值點的近似值)(xyy kx)(kxy* *ky)(xyy 數(shù)學界關注數(shù)學界關注工程師關注工程師關注如果找不到解函數(shù)如果找不到解函數(shù)數(shù)學界還關注：數(shù)學界還關注：解的存在性解的存在性解的唯一性解的唯一性解的光滑性解的光滑性解的振動性解的振動性解的周期性解

4、的周期性解的穩(wěn)定性解的穩(wěn)定性解的混沌性解的混沌性5.2 euler5.2 euler方法方法kp0p1p1npnpkx0 x1x1nxnx),(111212yxfxxyy 第五章：常微分方程數(shù)值解第五章：常微分方程數(shù)值解第一步：連續(xù)變量離散化第一步：連續(xù)變量離散化,nkxxxxx10第二步：用直線步進第二步：用直線步進),(000101yxfxxyy ),(),(nnnnnnnnnnyxhfyyyxfxxyy 111eulereuler格式格式1 1、eulereuler格式格式l 18 18世紀最杰出的數(shù)學家之一，世紀最杰出的數(shù)學家之一，1313歲歲時入讀巴塞爾大學，時入讀巴塞爾大學，151

5、5歲大學畢業(yè)，歲大學畢業(yè)，1616歲獲得碩士學位。歲獲得碩士學位。l 17271727年年-1741-1741年（年（2020歲歲-34-34歲）在彼歲）在彼得堡科學院從事研究工作，在分析學、得堡科學院從事研究工作，在分析學、數(shù)論、力學方面均有出色成就，并應數(shù)論、力學方面均有出色成就，并應俄國政府要求，解決了不少地圖學、俄國政府要求，解決了不少地圖學、造船業(yè)等實際問題。造船業(yè)等實際問題。l 2424歲晉升物理學教授。歲晉升物理學教授。l 17351735年（年（2828歲）右眼失明。歲）右眼失明。第五章：常微分方程數(shù)值解第五章：常微分方程數(shù)值解l 1741 1741年年 - 1766- 176

6、6（3434歲歲-59-59歲）任德國科學院物理數(shù)學所所歲）任德國科學院物理數(shù)學所所長，任職長，任職2525年。在行星運動、剛體運動、熱力學、彈道學、人年。在行星運動、剛體運動、熱力學、彈道學、人口學、微分方程、曲面微分幾何等研究領域均有開創(chuàng)性的工作。口學、微分方程、曲面微分幾何等研究領域均有開創(chuàng)性的工作。l 17661766年應沙皇禮聘重回彼得堡，在年應沙皇禮聘重回彼得堡，在17711771年（年（6464歲）左眼失歲）左眼失明。明。l eulereuler是數(shù)學史上最多產的數(shù)學家，平均以每年是數(shù)學史上最多產的數(shù)學家，平均以每年800800頁的速頁的速度寫出創(chuàng)造性論文。他去世后，人們用度寫出

7、創(chuàng)造性論文。他去世后，人們用3535年整理出他的研究成年整理出他的研究成果果7474卷。卷。 第五章：常微分方程數(shù)值解第五章：常微分方程數(shù)值解例例 p106p106第五章：常微分方程數(shù)值解第五章：常微分方程數(shù)值解 10102)(/yxyxyy初值問題初值問題bernoullibernoulli方程方程由由bernoullibernoulli方程的求解方法可得解析解方程的求解方法可得解析解xy21 ),(nnnnyxhfyy 1eulereuler格式為格式為 nnnnnyxyhyy21令令10. h將將1000 yx,代入代入eulereuler格式格式步進計算結果見步進計算結果見p106p1

8、06表表5.15.1第五章：常微分方程數(shù)值解第五章：常微分方程數(shù)值解eulereuler值值xy21 eulereuler格式的誤差分析格式的誤差分析事實上事實上eulereuler格式的每一步都存在誤差，為了方便討論算格式的每一步都存在誤差，為了方便討論算法的好壞，假定第法的好壞，假定第 n n 步準確的前提下分析第步準確的前提下分析第n+1n+1步的誤步的誤差，稱為差，稱為局部截斷誤差。局部截斷誤差。第五章：常微分方程數(shù)值解第五章：常微分方程數(shù)值解即即11 nnyxy)()(nnxyy 討論討論221hyhxyxynn)( )( )( 211121)( )( )()(nnnnnnnxxyx

9、xxyxyxy 由由taylortaylor公式公式nxnp1 np1 nx)(xy第五章：常微分方程數(shù)值解第五章：常微分方程數(shù)值解),()( )( )(nnnnnnnyxhfyhyhxyyyxy 21121 eulereuler格式的誤差為格式的誤差為)( )( )( nnxhyhyhxy 221 )( )( nxyhyh2222 2 2、后退、后退eulereuler格式格式令令),()( )()(nnnnnyxfxyhxyxy 1得得),()()(nnnnyxhfxyxy 1令令),(nnnnyxhfyy 1eulereuler格式格式第五章：常微分方程數(shù)值解第五章：常微分方程數(shù)值解后退

10、后退eulereuler格式格式令令),()( )()(1111 nnnnnyxfxyhxyxy得得),()()(111 nnnnyxhfxyxy隱式格式隱式格式令令),(111 nnnnyxhfyynx1 nx1 np1 np后退后退eulereuler格式的值格式的值eulereuler格式的值格式的值211 nnpp)(1 nxy3 3、梯形格式、梯形格式4 4、改進的、改進的eulereuler格式格式),(nnnnyxhfyy 1),(),(1112 nnnnnnyxfyxfhyy預測預測校正校正第五章：常微分方程數(shù)值解第五章：常微分方程數(shù)值解),(),(1112 nnnnnnyxf

11、yxfhyy隱式格式隱式格式為方便計算，一般用以下改進格式計算為方便計算，一般用以下改進格式計算用改進格式計算例用改進格式計算例5.15.1的結果見的結果見p110p110表表5.25.2第五章：常微分方程數(shù)值解第五章：常微分方程數(shù)值解5 5、兩步、兩步eulereuler格式格式第五章：常微分方程數(shù)值解第五章：常微分方程數(shù)值解)()(111 pnnhoyxy一般，如果一般，如果稱計算格式具有稱計算格式具有 階精度。階精度。p)( )(nnnxyhyxy2211 ),(nnnnyxhfyy 1已知已知eulereuler格式格式即即eulereuler格式具有一階精度格式具有一階精度),()(

12、 )()(nnnnnyxfxyhxyxy 211如果令如果令并假定并假定nnnnyxyyxy )(,)(11第五章：常微分方程數(shù)值解第五章：常微分方程數(shù)值解),()(nnnnyxhfyxy211 則有則有),(nnnnyxhfyy211 記記),()( )( )( )(nnnnnnnnyxhfyhyhxyhxyyyxy2612113211 則則將將 在在 點點taylortaylor展開展開)(1 nxynx1 ny的計算格式的計算格式其中其中),(nnnnyxhfyy211 )( )(nnxhyxy21 假設假設第五章：常微分方程數(shù)值解第五章：常微分方程數(shù)值解),()(,()( nnnnny

13、xfxyxfxy 特別要注意的是：一般特別要注意的是：一般331hy)( 兩點預測格式具有二階精度。兩點預測格式具有二階精度。當當 時時nnyxy )(),()( nnnyxfxy )( )( )( )( nnnnxhyhyhxyhxyy2612132 所以所以32116121hyhxyhxyyyxynnnnn)( )( )( )( 第五章：常微分方程數(shù)值解第五章：常微分方程數(shù)值解),(),(1112 nnnnnnyxfyxfhyy考察兩點校正格式的精度考察兩點校正格式的精度為便于處理，通常假定為便于處理，通常假定),()(,()( 11111 nnnnnyxfxyxfxy否則見否則見p108

14、p108),()(,()( 11111 nnnnnyxfxyxfxy一般情況下一般情況下112121 nnnnhyhyyy則則又又 ，nnyxy )(),()( nnnyxfxy 并記并記)( nnyxy 第五章：常微分方程數(shù)值解第五章：常微分方程數(shù)值解 )( 32212121hohyhyyhyynnnnn112121 nnnnhyhyyy)( 3324121hoyhyhhyynnnn 比較比較)( )(43216121hoyhyhhyyxynnnnn 得得)()(311hoyxynn 即梯形格式具有二階精度，因此兩步格式從預測到校即梯形格式具有二階精度，因此兩步格式從預測到校正均達到二階精度

15、。正均達到二階精度。因此得具有二階精度的兩步因此得具有二階精度的兩步eulereuler格式格式),(nnnnyxhfyy211 ),(),(1112 nnnnnnyxfyxfhyy預測預測校正校正第五章：常微分方程數(shù)值解第五章：常微分方程數(shù)值解nx1 nx1 nx)( nxyhyynn 1hyynn211 5.3 lunge-kutta5.3 lunge-kutta方法方法1 1、二階、二階lunge-kuttalunge-kutta方法（方法（p113-p115)p113-p115)第五章：常微分方程數(shù)值解第五章：常微分方程數(shù)值解依據(jù)精度要求的待定系數(shù)法依據(jù)精度要求的待定系數(shù)法令令101

16、pxxphxxnnnpn),(),(nnyxfk 1)(22111kkhyynn ),(pnpnyxfk 2確定確定 使使,2121kkp )(22111kkhyynn 具有二階精度具有二階精度),(1phkyxfnpn 平均斜率平均斜率 在點在點 的斜率的斜率),(nnyx)(xy 在點在點 的斜率的斜率),(pnpnyx )(xy用用eulereuler格式預格式預測測pny nx1 nxpnx 第五章：常微分方程數(shù)值解第五章：常微分方程數(shù)值解)(22111kkhyynn )( )( pnnnxyxyhy 21 phxyxyhxyhynnnn)( )( )( 1121 )()( 32221

17、hohpxyn)()()( )( )()(43223222121hohpxyphxyxyhynnnn 2121hxyhxyyxynnnn)( )( )(對照對照第五章：常微分方程數(shù)值解第五章：常微分方程數(shù)值解 121221p 可解得可解得 211021/p 1212121p/ ),(),()(12122111phkyxfkyxfkkkhyynpnnnnn 得得 )/,/(),(2212121hkyhxfkyxfkhkyynnnnnn ),(),(/ )(11212112hkyxfkyxfkkkhyynnnnnn改進的改進的eulereuler格式格式3 3、三階、三階lunge-kuttalu

18、nge-kutta方法方法第五章：常微分方程數(shù)值解第五章：常微分方程數(shù)值解補補充充 110qpqhxxpphxxnqnnpn,),(nnyxfk 1)(3322111kkkhyynn ),(),(12phkyphxfyxfknnpnpn 確定參數(shù)確定參數(shù) 使使srkkkqp,321321 )(22111kkhyynn 具有二階精度具有二階精度nx1 nxpnx qnx )(,(213skrkqhyqhxfknn 第五章：常微分方程數(shù)值解第五章：常微分方程數(shù)值解)(3322111kkkhyynn 使使具有三階精度具有三階精度)()(411hoyxynn 即即分別將分別將),(nnyxfk 1),

19、(12phkyphxfknn )(,(213skrkqhyqhxfknn )(3322111kkkhyynn 在點在點 taylortaylor展開，代入（展開，代入（p116) p116) ),(nnyx與與 的的taylortaylor展開比較展開比較 )(1 nxy第五章：常微分方程數(shù)值解第五章：常微分方程數(shù)值解 613121113232232321/pqsqpqpsr 得得可解得可解得 21616461121321srqp,/,/,/,/ 第五章：常微分方程數(shù)值解第五章：常微分方程數(shù)值解得一個三階精度的得一個三階精度的runge-kuttarunge-kutta格式格式 )(,()/,

20、/(),()(123121321122246kkhyhxfkhkyhxfkyxfkkkkhyynnnnnnnn4 4、四階、四階lunge-kuttalunge-kutta方法見方法見p117p1175.4 5.4 幾種方法的數(shù)值計算幾種方法的數(shù)值計算例例5.1 p1065.1 p106改進的改進的eulereuler格式格式eulereuler格式格式第五章：常微分方程數(shù)值解第五章：常微分方程數(shù)值解四階經典四階經典lunge-kuttalunge-kutta方法方法例例5.1 p1065.1 p106第五章：常微分方程數(shù)值解第五章：常微分方程數(shù)值解x x精確值精確值eulereuler方法方

21、法改進改進eulereuler方法方法四階四階lunge-lunge-kuttakuttaeulereuler方法方法誤差誤差改進改進eulereuler誤差誤差四階四階lunge-lunge-kuttakutta誤差誤差0 01.00000001.00000001.00000001.00000001.00000001.00000001.00000001.00000000.00000000.00000000.00000000.00000000.00000000.00000000.10.11.09544511.09544511.10000001.10000001.09590911.0959091

22、1.09544551.09544550.00455490.00455490.00046400.00046400.00000040.00000040.20.21.18321601.18321601.19181821.19181821.18409661.18409661.18321671.18321670.00860220.00860220.00088060.00088060.00000080.00000080.30.31.26491111.26491111.27743781.27743781.26620141.26620141.26491221.26491220.01252680.0125268

23、0.00129030.00129030.00000120.00000120.40.41.34164081.34164081.35821261.35821261.34336021.34336021.34164241.34164240.01657180.01657180.00171940.00171940.00000160.00000160.50.51.41421361.41421361.43513291.43513291.41640191.41640191.41421561.41421560.02091940.02091940.00218840.00218840.00000200.0000020

24、0.60.61.48323971.48323971.50896631.50896631.48595561.48595561.48324221.48324220.02572660.02572660.00271590.00271590.00000250.00000250.70.71.54919331.54919331.58033821.58033821.55251411.55251411.54919651.54919650.03114490.03114490.00332080.00332080.00000310.00000310.80.81.61245151.61245151.64978341.6

25、4978341.61647481.61647481.61245531.61245530.03733190.03733190.00402320.00402320.00000380.00000380.90.91.67332011.67332011.71777931.71777931.67816641.67816641.67332471.67332470.04445930.04445930.00484630.00484630.00000460.00000461 11.73205081.73205081.78477081.78477081.73786741.73786741.73205641.7320

26、5640.05272000.05272000.00581660.00581660.00000560.0000056幾種方法的結果與誤差幾種方法的結果與誤差第五章：常微分方程數(shù)值解第五章：常微分方程數(shù)值解參考程序參考程序-euler-eulerx=0:0.01:1;x=0:0.01:1;y=sqrt(1+2.y=sqrt(1+2.* *x);x);a=0.0;b=1.0;n=10;a=0.0;b=1.0;n=10;h=(b-a)/n;h=(b-a)/n;x0=a:h:b;y0(1)=1.0;x0=a:h:b;y0(1)=1.0;forfor k=1:10 k=1:10 y0(k+1)=y0(k)

27、+h y0(k+1)=y0(k)+h* *(y0(k)-2(y0(k)-2* *x0(k)/y0(k);x0(k)/y0(k);endendforfor i=1:10 i=1:10 y1(1)=1.0; y1(1)=1.0; y1(i+1)=y1(i)+h y1(i+1)=y1(i)+h* *(y1(i)-2(y1(i)-2* *x0(i)/y1(i);x0(i)/y1(i); y1(i+1)=y1(i)+h y1(i+1)=y1(i)+h* *(y1(i)-2(y1(i)-2* *x0(i)/y1(i)+y1(i+1)- x0(i)/y1(i)+y1(i+1)- 2 2* *x0(i+1)/

28、y1(i+1)/2;x0(i+1)/y1(i+1)/2;endendplot(x,y,b);plot(x,y,b);hold on; hold on; plot(x0,y0,or);plot(x0,y0,or);hold on; hold on; plot(x0,y1,plot(x0,y1,* *););第五章：常微分方程數(shù)值解第五章：常微分方程數(shù)值解參考程序參考程序-lunge_kutta-lunge_kuttax=0:0.01:1;x=0:0.01:1;y=sqrt(1+2.y=sqrt(1+2.* *x);x);a=0.0;b=1.0;n=10;a=0.0;b=1.0;n=10;h=(b

29、-a)/n;h=(b-a)/n;x0=a:h:b;x0=a:h:b;y0(1)=1.0;y0(1)=1.0;forfor k=1:10 k=1:10 k1=y0(k)-2 k1=y0(k)-2* *x0(k)/y0(k);x0(k)/y0(k); k2=y0(k)+h k2=y0(k)+h* *k1/2-(2k1/2-(2* *x0(k)+h)/(y0(k)+hx0(k)+h)/(y0(k)+h* *k1/2);k1/2); k3=y0(k)+h k3=y0(k)+h* *k2/2-(2k2/2-(2* *x0(k)+h)/(y0(k)+hx0(k)+h)/(y0(k)+h* *k2/2);k

30、2/2); k4=y0(k)+h k4=y0(k)+h* *k3-2k3-2* *(x0(k)+h)/(y0(k)+h(x0(k)+h)/(y0(k)+h* *k3);k3); y0(k+1)=y0(k)+h y0(k+1)=y0(k)+h* *(k1+2(k1+2* *k2+2k2+2* *k3+k4)/6;k3+k4)/6;endendhold on; hold on; plot(x,y,b);plot(x,y,b);hold on; hold on; plot(x0,y0,or);plot(x0,y0,or);第五章：常微分方程數(shù)值解第五章：常微分方程數(shù)值解),(yxfy 5.5 5.5

31、 線性多步方法線性多步方法1 1、adamsadams顯式格式顯式格式dxyxfdxynnnnxxxx 11),(dxyxfxyxynnxxnn 11),()()(nnnnxxyxyxyxhfdxyxfnn )(),(,(),(1如果令),()(nnnnyxhfyxy 1有1kxkx),(nnnnyxhfyy 1)(,(),()(1112 nnnnnnxyxfyxfhyxy),(),(1112 nnnnnnyxfyxfhyy第五章：常微分方程數(shù)值解第五章：常微分方程數(shù)值解dxxpdxyxfxyxynnnnxxrxxnn 111)(),()()(kxnx1nxrnxp27 (2.5.12)p27

32、 (2.5.12)(thxnnr nrnnnfrrtttfttftf !)()(!)(11212njrjjfjt 01)(knknkff 11 nnnnffff222112 nnnnnnffffffjnjrjjfjt 01)(第五章：常微分方程數(shù)值解第五章：常微分方程數(shù)值解dxxpdxyxfxyxynnnnxxrxxnn 111)(),()()(jnjrjjnrfjtthxn 01)()(dththxnyynrnn 101)(dtfjthyjnjrjjn 1001)(jnjrjjnfdtjthy 0101)(jnjrjjnfhy 0 第五章：常微分方程數(shù)值解第五章：常微分方程數(shù)值解 101dt

33、jtjj)( 1110100 tdt 21211102101 tdtt)( 125213121211023102 ttdttt!)( 8341613211034103 tttdtttt!)( 第五章：常微分方程數(shù)值解第五章：常微分方程數(shù)值解jnjrjjnnfhyy 01 )(1101 nnnnffhyy )()(1132121 nnnnnnnffhyfffhyr=3 r=3 時見時見p125p125（5.5.65.5.6）第五章：常微分方程數(shù)值解第五章：常微分方程數(shù)值解)9375955(243211nnnnnnffffhyy（5.5.65.5.6）式的誤差分析）式的誤差分析)9375955(2

34、43211nnnnnnffffhyy3 , 2 , 1 , 0)(kyxyknkn假設)( 9)( 37)( 59)( 55(24)(321nnnnnxyxyxyxyhxy)( 2455)(nnxhyxy)(2462 24595)5(4)4(3)3(2hoyhyhyhhyyhnnnnn)(241668242 24375)5(4)4(3)3(2hoyhyhyhhyyhnnnnn第五章：常微分方程數(shù)值解第五章：常微分方程數(shù)值解)( )( 249243724592455nnxhyxhy)(2481627293 2495)5(4)4(3)3(2hoyhyhyhhyyhnnnnn)( 21)( 2427

35、2474245922nnxyhxyh)(61)(4881481484859)3(3)3(3nnxyhxyh)(241)(2424324296245961)4(4)4(4nnxyhxyh)(14449)(24729245922459241)5(5)5(5nnxyhxyh第五章：常微分方程數(shù)值解第五章：常微分方程數(shù)值解)()()(4511hohoyxynn 例例5.1 p1065.1 p106第五章：常微分方程數(shù)值解第五章：常微分方程數(shù)值解adams adams 程序程序x=0:0.01:1;x=0:0.01:1;y=sqrt(1+2.y=sqrt(1+2.* *x);x);a=0.0;b=1.0

36、;n=10;a=0.0;b=1.0;n=10;h=(b-a)/n;h=(b-a)/n;x0=a:h:b;x0=a:h:b;y0(1)=1.0;y0(2)=1.0954;y0(3)=1.1832;y0(4)=1.2649;y0(1)=1.0;y0(2)=1.0954;y0(3)=1.1832;y0(4)=1.2649;forfor k=4:10 k=4:10 y0(k+1)=y0(k)+h y0(k+1)=y0(k)+h* *(55(55* *(y0(k)-2(y0(k)-2* *x0(k)/y0(k)x0(k)/y0(k)- -5959* *(y0(k-1)-2(y0(k-1)-2* *x0(

37、k-1)/y0(k-1)x0(k-1)/y0(k-1)+37+37* *(y0(k-2)-2(y0(k-2)-2* *x0(k-x0(k-2)/y0(k-2)2)/y0(k-2)-9-9* *(y0(k-3)-2(y0(k-3)-2* *x0(k-3)/y0(k-3)x0(k-3)/y0(k-3)/24;)/24;endendhold on; hold on; plot(x,y,b);plot(x,y,b);hold on; hold on; plot(x0,y0,or);plot(x0,y0,or);第五章：常微分方程數(shù)值解第五章：常微分方程數(shù)值解x x精確值精確值四階四階lunge-lun

38、ge-kuttakuttaadamsadams四階四階lunge-lunge-kuttakutta誤差誤差adamsadams誤差誤差0 01.00000001.00000001.00000001.00000001.00000001.00000000.00000000.00000000.00000000.00000000.10.11.09544511.09544511.09544551.09544551.09540001.09540000.00000040.00000040.00004510.00004510.20.21.18321601.18321601.18321671.18321671.

39、18320001.18320000.00000080.00000080.00001600.00001600.30.31.26491111.26491111.26491221.26491221.26490001.26490000.00000120.00000120.00001110.00001110.40.41.34164081.34164081.34164241.34164241.34153281.34153280.00000160.00000160.00010800.00010800.50.51.41421361.41421361.41421561.41421561.41402401.41402400.00000200.00000200.000