数值分析第五讲：常微分方程数值解

1、第五章：常微分方程数值解第五章：常微分方程数值解5.1 5.1 引言引言1 1、常微分方程与解、常微分方程与解为为n n阶常微分方程阶常微分方程。0 ), , ,()(nyyyyxf如果函数如果函数 在区间在区间a,ba,b内内n n阶可导，称方程阶可导，称方程)(xyy )(xyy 满足方程的函数满足方程的函数称为微分方程的称为微分方程的解解。则则如如为任意常数）为任意常数）xy2 ccxy(, 2一般称为方程的一般称为方程的通解通解。为方程的解。为方程的解。12 xy如果如果则有则有10 )(y为方程满足定解条件的解。为方程满足定解条件的解。第五章：常微分方程数值解第五章：常微分方程数值解

2、 102)(yxy ccxy1212 xy方程的通解方程的通解满足定解条件的解满足定解条件的解微分关系（方程）微分关系（方程）解的图示解的图示第五章：常微分方程数值解第五章：常微分方程数值解本教材重点讨论定解问题本教材重点讨论定解问题( (初值问题）初值问题）定解条件（初始条件）定解条件（初始条件） 00yxyyxfy)(),(),(yxf是否能够找到定解问题的解取决于是否能够找到定解问题的解取决于仅有极少数的方程可以通过仅有极少数的方程可以通过“常数变易法常数变易法”、“可分可分离变量法离变量法”等特殊方法求得初等函数形式的解，绝大等特殊方法求得初等函数形式的解，绝大部分方程至今无法理论求解

3、。部分方程至今无法理论求解。如如xyxyeyxyxyy 212,),sin(sin等等等等2 2、数值解的思想、数值解的思想第五章：常微分方程数值解第五章：常微分方程数值解（1 1）将连续变量）将连续变量 离散为离散为,bax bxxxxank 10nkxyykk,)(21 （2 2）用代数的方法求出解函数）用代数的方法求出解函数 在在 点的近似值点的近似值)(xyy kx)(kxy* *ky)(xyy 数学界关注数学界关注工程师关注工程师关注如果找不到解函数如果找不到解函数数学界还关注：数学界还关注：解的存在性解的存在性解的唯一性解的唯一性解的光滑性解的光滑性解的振动性解的振动性解的周期性解

4、的周期性解的稳定性解的稳定性解的混沌性解的混沌性5.2 euler5.2 euler方法方法kp0p1p1npnpkx0 x1x1nxnx),(111212yxfxxyy 第五章：常微分方程数值解第五章：常微分方程数值解第一步：连续变量离散化第一步：连续变量离散化,nkxxxxx10第二步：用直线步进第二步：用直线步进),(000101yxfxxyy ),(),(nnnnnnnnnnyxhfyyyxfxxyy 111eulereuler格式格式1 1、eulereuler格式格式l 18 18世纪最杰出的数学家之一，世纪最杰出的数学家之一，1313岁岁时入读巴塞尔大学，时入读巴塞尔大学，151

5、5岁大学毕业，岁大学毕业，1616岁获得硕士学位。岁获得硕士学位。l 17271727年年-1741-1741年（年（2020岁岁-34-34岁）在彼岁）在彼得堡科学院从事研究工作，在分析学、得堡科学院从事研究工作，在分析学、数论、力学方面均有出色成就，并应数论、力学方面均有出色成就，并应俄国政府要求，解决了不少地图学、俄国政府要求，解决了不少地图学、造船业等实际问题。造船业等实际问题。l 2424岁晋升物理学教授。岁晋升物理学教授。l 17351735年（年（2828岁）右眼失明。岁）右眼失明。第五章：常微分方程数值解第五章：常微分方程数值解l 1741 1741年年 - 1766- 176

6、6（3434岁岁-59-59岁）任德国科学院物理数学所所岁）任德国科学院物理数学所所长，任职长，任职2525年。在行星运动、刚体运动、热力学、弹道学、人年。在行星运动、刚体运动、热力学、弹道学、人口学、微分方程、曲面微分几何等研究领域均有开创性的工作。口学、微分方程、曲面微分几何等研究领域均有开创性的工作。l 17661766年应沙皇礼聘重回彼得堡，在年应沙皇礼聘重回彼得堡，在17711771年（年（6464岁）左眼失岁）左眼失明。明。l eulereuler是数学史上最多产的数学家，平均以每年是数学史上最多产的数学家，平均以每年800800页的速页的速度写出创造性论文。他去世后，人们用度写出

7、创造性论文。他去世后，人们用3535年整理出他的研究成年整理出他的研究成果果7474卷。卷。 第五章：常微分方程数值解第五章：常微分方程数值解例例 p106p106第五章：常微分方程数值解第五章：常微分方程数值解 10102)(/yxyxyy初值问题初值问题bernoullibernoulli方程方程由由bernoullibernoulli方程的求解方法可得解析解方程的求解方法可得解析解xy21 ),(nnnnyxhfyy 1eulereuler格式为格式为 nnnnnyxyhyy21令令10. h将将1000 yx,代入代入eulereuler格式格式步进计算结果见步进计算结果见p106p1

8、06表表5.15.1第五章：常微分方程数值解第五章：常微分方程数值解eulereuler值值xy21 eulereuler格式的误差分析格式的误差分析事实上事实上eulereuler格式的每一步都存在误差，为了方便讨论算格式的每一步都存在误差，为了方便讨论算法的好坏，假定第法的好坏，假定第 n n 步准确的前提下分析第步准确的前提下分析第n+1n+1步的误步的误差，称为差，称为局部截断误差。局部截断误差。第五章：常微分方程数值解第五章：常微分方程数值解即即11 nnyxy)()(nnxyy 讨论讨论221hyhxyxynn)( )( )( 211121)( )( )()(nnnnnnnxxyx

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

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预测预测校正校正第五章：常微分方程数值解第五章：常微分方程数值解),(),(1112 nnnnnnyxf

11、yxfhyy隐式格式隐式格式为方便计算，一般用以下改进格式计算为方便计算，一般用以下改进格式计算用改进格式计算例用改进格式计算例5.15.1的结果见的结果见p110p110表表5.25.2第五章：常微分方程数值解第五章：常微分方程数值解5 5、两步、两步eulereuler格式格式第五章：常微分方程数值解第五章：常微分方程数值解)()(111 pnnhoyxy一般，如果一般，如果称计算格式具有称计算格式具有 阶精度。阶精度。p)( )(nnnxyhyxy2211 ),(nnnnyxhfyy 1已知已知eulereuler格式格式即即eulereuler格式具有一阶精度格式具有一阶精度),()(

12、 )()(nnnnnyxfxyhxyxy 211如果令如果令并假定并假定nnnnyxyyxy )(,)(11第五章：常微分方程数值解第五章：常微分方程数值解),()(nnnnyxhfyxy211 则有则有),(nnnnyxhfyy211 记记),()( )( )( )(nnnnnnnnyxhfyhyhxyhxyyyxy2612113211 则则将将 在在 点点taylortaylor展开展开)(1 nxynx1 ny的计算格式的计算格式其中其中),(nnnnyxhfyy211 )( )(nnxhyxy21 假设假设第五章：常微分方程数值解第五章：常微分方程数值解),()(,()( nnnnny

13、xfxyxfxy 特别要注意的是：一般特别要注意的是：一般331hy)( 两点预测格式具有二阶精度。两点预测格式具有二阶精度。当当 时时nnyxy )(),()( nnnyxfxy )( )( )( )( nnnnxhyhyhxyhxyy2612132 所以所以32116121hyhxyhxyyyxynnnnn)( )( )( )( 第五章：常微分方程数值解第五章：常微分方程数值解),(),(1112 nnnnnnyxfyxfhyy考察两点校正格式的精度考察两点校正格式的精度为便于处理，通常假定为便于处理，通常假定),()(,()( 11111 nnnnnyxfxyxfxy否则见否则见p108

14、p108),()(,()( 11111 nnnnnyxfxyxfxy一般情况下一般情况下112121 nnnnhyhyyy则则又又 ，nnyxy )(),()( nnnyxfxy 并记并记)( nnyxy 第五章：常微分方程数值解第五章：常微分方程数值解 )( 32212121hohyhyyhyynnnnn112121 nnnnhyhyyy)( 3324121hoyhyhhyynnnn 比较比较)( )(43216121hoyhyhhyyxynnnnn 得得)()(311hoyxynn 即梯形格式具有二阶精度，因此两步格式从预测到校即梯形格式具有二阶精度，因此两步格式从预测到校正均达到二阶精度

15、。正均达到二阶精度。因此得具有二阶精度的两步因此得具有二阶精度的两步eulereuler格式格式),(nnnnyxhfyy211 ),(),(1112 nnnnnnyxfyxfhyy预测预测校正校正第五章：常微分方程数值解第五章：常微分方程数值解nx1 nx1 nx)( nxyhyynn 1hyynn211 5.3 lunge-kutta5.3 lunge-kutta方法方法1 1、二阶、二阶lunge-kuttalunge-kutta方法（方法（p113-p115)p113-p115)第五章：常微分方程数值解第五章：常微分方程数值解依据精度要求的待定系数法依据精度要求的待定系数法令令101

16、pxxphxxnnnpn),(),(nnyxfk 1)(22111kkhyynn ),(pnpnyxfk 2确定确定 使使,2121kkp )(22111kkhyynn 具有二阶精度具有二阶精度),(1phkyxfnpn 平均斜率平均斜率 在点在点 的斜率的斜率),(nnyx)(xy 在点在点 的斜率的斜率),(pnpnyx )(xy用用eulereuler格式预格式预测测pny nx1 nxpnx 第五章：常微分方程数值解第五章：常微分方程数值解)(22111kkhyynn )( )( pnnnxyxyhy 21 phxyxyhxyhynnnn)( )( )( 1121 )()( 32221

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

18、nge-kutta方法方法第五章：常微分方程数值解第五章：常微分方程数值解补补充充 110qpqhxxpphxxnqnnpn,),(nnyxfk 1)(3322111kkkhyynn ),(),(12phkyphxfyxfknnpnpn 确定参数确定参数 使使srkkkqp,321321 )(22111kkhyynn 具有二阶精度具有二阶精度nx1 nxpnx qnx )(,(213skrkqhyqhxfknn 第五章：常微分方程数值解第五章：常微分方程数值解)(3322111kkkhyynn 使使具有三阶精度具有三阶精度)()(411hoyxynn 即即分别将分别将),(nnyxfk 1),

19、(12phkyphxfknn )(,(213skrkqhyqhxfknn )(3322111kkkhyynn 在点在点 taylortaylor展开，代入（展开，代入（p116) p116) ),(nnyx与与 的的taylortaylor展开比较展开比较 )(1 nxy第五章：常微分方程数值解第五章：常微分方程数值解 613121113232232321/pqsqpqpsr 得得可解得可解得 21616461121321srqp,/,/,/,/ 第五章：常微分方程数值解第五章：常微分方程数值解得一个三阶精度的得一个三阶精度的runge-kuttarunge-kutta格式格式 )(,()/,

20、/(),()(123121321122246kkhyhxfkhkyhxfkyxfkkkkhyynnnnnnnn4 4、四阶、四阶lunge-kuttalunge-kutta方法见方法见p117p1175.4 5.4 几种方法的数值计算几种方法的数值计算例例5.1 p1065.1 p106改进的改进的eulereuler格式格式eulereuler格式格式第五章：常微分方程数值解第五章：常微分方程数值解四阶经典四阶经典lunge-kuttalunge-kutta方法方法例例5.1 p1065.1 p106第五章：常微分方程数值解第五章：常微分方程数值解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几种方法的结果与误差几种方法的结果与误差第五章：常微分方程数值解第五章：常微分方程数值解参考程序参考程序-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,* *););第五章：常微分方程数值解第五章：常微分方程数值解参考程序参考程序-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);第五章：常微分方程数值解第五章：常微分方程数值解),(yxfy 5.5 5.5

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

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

33、jtjj)( 1110100 tdt 21211102101 tdtt)( 125213121211023102 ttdttt!)( 8341613211034103 tttdtttt!)( 第五章：常微分方程数值解第五章：常微分方程数值解jnjrjjnnfhyy 01 )(1101 nnnnffhyy )()(1132121 nnnnnnnffhyfffhyr=3 r=3 时见时见p125p125（5.5.65.5.6）第五章：常微分方程数值解第五章：常微分方程数值解)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第五章：常微分方程数值解第五章：常微分方程数值解)( )( 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第五章：常微分方程数值解第五章：常微分方程数值解)()()(4511hohoyxynn 例例5.1 p1065.1 p106第五章：常微分方程数值解第五章：常微分方程数值解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);第五章：常微分方程数值解第五章：常微分方程数值解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