3第三讲：常微分方程数值解法

第三讲常微分方程数值解§5.1引言1、常微分方程与解为n阶常微分方程。0),,'',',,()(nyyyyxF如果函数在区间[a,b]内n阶可导，称方程)(xyy)(xyy满足方程的函数称为微分方程的解。则如为任意常数）xy2'CCxy(,2为方程的解,一般称为方程的通解。12xy如果则有10)(y为方程满足定解条件的解。第三讲常微分方程数值解00)(),('yxyyxfy如1)0(2'yxyCCxy1212xy第三讲常微分方程数值解微分关系（方程）方程的通解满足定解条件的解解的图示2、数值解的思想（1）将连续变量离散为],[baxbxxxxank10nkxyykk,,,)(21（2）用代数的方法求出解函数在点的近似值)(xyykx)(kxy*ky)(xyy数学界关注工程师关注如果找不到解函数数学界还关注：解的存在性解的唯一性解的光滑性解的振动性解的周期性解的稳定性解的混沌性解的分岔性……第三讲常微分方程数值解§5.2Euler方法kp0p1p1npnpkx0x1x1nxnx),()('1111212yxfxyxxyy第一步：连续变量离散化,,,,,nkxxxxx10第二步：用直线步进),()('0000101yxfxyxxyy),(),(111nnnnnnnnnnyxhfyyyxfxxyy·····Euler格式1、Euler格式第三讲常微分方程数值解第三讲常微分方程数值解1、Euler方法的误差分析及改进（1）Euler方法的误差)()()(''2)(2211hohOyhyxynnEuler格式具有局部一阶精度)(')())(,()(),(1nnnnnnnnnxhyxyxyxhfxyyxhfyy事实上，假设nnyxy)()(''21)(')()(21nnnnxyhxhyxyxy第三讲常微分方程数值解（2）后退Euler格式后退Euler格式也只具有局部一阶精度1nxnx1ny1~ny),(111nnnnyxhfyy))(,()()(1111nnnnnnxyxfxxxyxy由2、梯形格式)],(),([2111nnnnnnyxfyxfhyy同理可证梯形格式具有局部二阶精度第三讲常微分方程数值解3、改进的Euler格式),(1nnnnyxhfyy)],(),([2111nnnnnnyxfyxfhyy预测校正4、两步Euler格式),(211nnnnyxhfyy)],(),([2111nnnnnnyxfyxfhyy预测校正两步Euler格式具有局部二阶精度1nxnx1nx§5.3Lunge-Kutta方法二阶Lunge-Kutta方法第三讲常微分方程数值解令),(),()(12122111phKyxfKyxfKKKhyynpnnnnn因为)(')(')(211pnnnnxyhxyhxyy)(')(11nnnxyhxyy)(''')(21)('')('22nnnxyphxphyxyh···hxxnn1nxphxxnpn进而第三讲常微分方程数值解)(')(11nnnxyhxyy)('''21)('')('322222nnnxyhpxyphxyh与)('''61)(''21)(')()(321nnnnnxyhxyhxhyxyxy211221p比较得解得211021p1212121p得具有二阶精度Runge-Kutta公式2,2),(12121hKyhxfKyxfKhKyynnnnnn121211,),()(2hKyhxfKyxfKKKhyynnnnnn三阶Lunge-Kutta方法思想与方法同二阶（略）第三讲常微分方程数值解···hxxnn1nxphxxnpn·qhxxnqn令))(,(),(),()(2131213322111sKrKqhyxfKphKyxfKyxfKKKKhyynqnnpnnnnn四阶Lunge-Kutta经典格式),()2,2()2,2(),(6)22(342312133211hKyhxfKhKyhxfKhKyhxfKyxfKKKKKhyynnnnnnnnnn例xyyxyxyy211)0()10(2'§5.4几种方法的数值计算改进的Euler格式Euler格式第三讲常微分方程数值解四阶经典Lunge-Kutta方法第三讲常微分方程数值解x精确值Euler方法改进Euler方法四阶Lunge-KuttaEuler方法误差改进Euler误差四阶Lunge-Kutta误差01.00000001.00000001.00000001.00000000.00000000.00000000.00000000.11.09544511.10000001.09590911.09544550.00455490.00046400.00000040.21.18321601.19181821.18409661.18321670.00860220.00088060.00000080.31.26491111.27743781.26620141.26491220.01252680.00129030.00000120.41.34164081.35821261.34336021.34164240.01657180.00171940.00000160.51.41421361.43513291.41640191.41421560.02091940.00218840.00000200.61.48323971.50896631.48595561.48324220.02572660.00271590.00000250.71.54919331.58033821.55251411.54919650.03114490.00332080.00000310.81.61245151.64978341.61647481.61245530.03733190.00402320.00000380.91.67332011.71777931.67816641.67332470.04445930.00484630.000004611.73205081.78477081.73786741.73205640.05272000.00581660.0000056几种方法的结果与误差第三讲常微分方程数值解参考程序-Eulerx=0:0.01:1;y=sqrt(1+2.*x);a=0.0;b=1.0;n=10;h=(b-a)/n;x0=a:h:b;y0(1)=1.0;fork=1:10y0(k+1)=y0(k)+h*(y0(k)-2*x0(k)/y0(k));endfori=1:10y1(1)=1.0;y1(i+1)=y1(i)+h*(y1(i)-2*x0(i)/y1(i));y1(i+1)=y1(i)+h*(y1(i)-2*x0(i)/y1(i)+y1(i+1)-2*x0(i+1)/y1(i+1))/2;endplot(x,y,'b');holdon;plot(x0,y0,'or');holdon;plot(x0,y1,'*');第三讲常微分方程数值解参考程序-Lunge_Kuttax=0:0.01:1;y=sqrt(1+2.*x);a=0.0;b=1.0;n=10;h=(b-a)/n;x0=a:h:b;y0(1)=1.0;fork=1:10k1=y0(k)-2*x0(k)/y0(k);k2=y0(k)+h*k1/2-(2*x0(k)+h)/(y0(k)+h*k1/2);k3=y0(k)+h*k2/2-(2*x0(k)+h)/(y0(k)+h*k2/2);k4=y0(k)+h*k3-2*(x0(k)+h)/(y0(k)+h*k3);y0(k+1)=y0(k)+h*(k1+2*k2+2*k3+k4)/6;endholdon;plot(x,y,'b');holdon;plot(x0,y0,'or');第三讲常微分方程数值解),('yxfy§5.5线性多步方法1、Adams显式格式dxyxfdxynnnnxxxx11),('dxyxfxyxynnxxnn1),()()(1nnnnxxyxyxyxhfdxyxfnn)()),(,(),(1如果令),()(1nnnnyxhfyxy有1kxkx),(1nnnnyxhfyy)))(,(),((2)(111nnnnnnxyxfyxfhyxy)),(),((2111nnnnnnyxfyxfhyy第三讲常微分方程数值解dxxPdxyxfxyxynnnnxxrxxnn11)(),()()(1kxnx1nxrnx····P27(2.5.12))(thxNnrnrnnnfrrtttfttftf!)1()1(!2)1(2njrjjfjt0)1(knknkff11nnnnffff222112nnnnnnffffffjnjrjjfjt0)1(第三讲常微分方程数值解dxxPdxyxfxyxynnnnxxrxxnn11)(),()()(1jnjrjjnrfjtthxN0)1()(dththxNyynrnn101)(dtfjthyjnjrjjn100)1(jnjrjjnfdtjthy010)1(jnjrjjnfhy0第三讲常微分方程数值解10)1(dtjtjj1110100tdt2121)(1102101tdtt125213121!2)1(1023102ttdttt834161!3)2)(1(1034103tttdtttt第三讲常微分方程数值解jnjrjjnnfhyy01)(1101nnnnffhyy)3(21))(21(11nnnnnnnffhyfffhyr=3时见P125（5.5.6）第三讲常微分方程数值解（5.5.6）式的误差分析)9375955(243211nnnnnnffffhyy3,2,1,0)(kyxyknkn假设))('9)('37)('59)('55(24)(321nnnnnxyxyxyxyhxy)('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)(''24272474245922nnxyhxyh)(61)(4881481484859)3(3)3(