欧拉方法及其改进的欧拉方法的Matlab实现

-1-常微分方程初值问题的欧拉方法及其改进的欧拉方法的Matlab实现纪秀浩辽宁工程技术大学理学院，辽宁阜新（123000）E-mail：[e-mail]摘要：欧拉(Euler)方法及改进的欧拉方法是解决常微分方程初值问题常用的数值解法，但Matlab的工具箱中没有Euler方法的功能函数。本文在简要介绍Euler方法及其改进的Euler方法的基础上，通过编写Matlab程序实现两种数值解法，并通过作图形式对比其精度，加深对两种方法的认识。关键词：欧拉方法；改进的欧拉方法；matlab实现1．引言常微分方程是解决工程实例的常用的工具[1]，建立微分方程只是解决问题的第一步，通常需要求出方程的解来说明实际现象，并加以检验。如果能得到解析形式的解固然是便于分析和应用的，但是我们知道，虽然求解常微分方程有各种各样的解析方法，但解析方法只能用来求解一些特殊类型的方程，实际问题终归结出来的微分方程主要靠数值解法[2]。数值解法就是一个十分重要的手段，而欧拉方法又是数值解法昀基础昀常用的方法。2.欧拉方法、改进的欧拉方法及Matlab实现下面主要讨论一阶常微分方程的初值问题，其一般形式为：'00(,)()yfxyyxy⎧=⎨=⎩(1)我们知道，只要函数(,)fxy适当光滑——譬如关于y满足利普希茨(Lipschitz)条件(,)(,)fxyfxyLyy−≤−理论上就可以保证初值问题(1)的解()yyx=存在并且唯一[3]。所谓数值解法，就是求问题(1)在某些离散点01Naxxxb==的近似值012,,,Nyyyy的方法。012,,,Nyyyy就称为问题(1)的数值解。1nnnhxx+=−成为nx到nx的步长，我们为了方便取为常量h。2.1.欧拉方法2.1.1欧拉方法将微分方程离散化，用向前差商1()()nnyxyxh+−代替微分'()nyx，代入(1)中的微分方程，可得1()()(,())(1,2,3,)nnnnyxyxfxyxnh+−==化简可得1()()(,())(1,2,3,)nnnnyxyxfxyxhn+=+=近似()nyx代入上式便可得到1()nyx+的近似值1ny+，计算式为：1(,)(1,2,3,)nnnnyyfxyhn+=+=(2)这样问题(1)的近似解可通过求解下面的差分初值问题：10(,)(1,2,3,)()nnnnyyfxyhnyya+=+=⎧⎨=⎩(3)得到，按(3)式由初值0y可逐次求出12,,yy。Eule方法就是用差分方程初值问题(3)的解来近似微分方程初值问题(1)的解。即由公式(3)算出()nyx的近似值。这组公式求问题(1)的数值解就是著名的欧拉(Euler)公式。2.1.2欧拉方法的误差估计对于Euler公式(3)我们看到，当1,2,n=时公式右端的ny都是近似的，所以用它计算的1ny+会有累积误差，分析累积误差比较复杂，这里先讨论比较简单的所谓局部截断误差。假定用(3)式时右端的ny没有误差，即()nnyyx=那么由此算出1()(,())nnnnyyxfxyxh+=+(4)局部截断误差指的是，按(4)式计算由nx到1nx+这一步的计算值1ny+与精确值1()nyx+之差11()nnyxy++−。为了估计它，由Taylor展开得到的精确值1()nyx+是2'''31()()()()()2nnnnhyxyxhyxyxOh+=+++(5)(4)、(5)两式相减（注意到'(,)yfxy=）得2''3211()()()()2nnnhyxyyxOhOh++−=+≈(6)即局部截断误差是2h阶的，而数值算法的精度定义为：若一种算法的局部截断误差为1()pOh+，则称该算法具有p阶精度。显然p越大，方法的精度越高。式(6)说明，Euler方法是一阶方法，因此它的精度不高。2.2改进的欧拉方法2.2.1改进的欧拉方法用数值积分方法离散化问题(1)，两端积分可得11()()(,())(0,1,2,)nnxnnxyxyxfxyxdxn++−==∫对右端积分使用梯形公式可得,111(,())[(,())(,())]2nnxnnnnxhfxyxdxfxyxfxyx+++≈+∫再用1,nnyy+代替1(),()nnyxyx+，则得计算公式111[(,)(,)]2nnnnnnhyyfxyfxy+++=++(7)很明显可以注意到(7)式为隐式形式。改进的欧拉方法是先用欧拉公式求1()nyx+的一个近似值1ny+，称为预测值，然后用梯形公式进行矫正并求得近似值1ny+。即1111(,)[(,)(,)]2nnnnnnnnnnyyfxyhhyyfxyfxy++++⎧=+⎪⎨=++⎪⎩(8)2.2.2改进的欧拉方法的误差估计由(7)式可知，11232434()()['()'()]2'()()'()23![()'()''()'''()]()22'''()()12nnnnnnnnnnnnhyxyxyxyxhhhyxyxyxhhyxyxhyxyxOhhyxOh++−−+=++−++++=−+所以改进的欧拉方法是二阶的，精度较欧拉方法要高，实用性更加广泛。3．实验及结果分析3.1欧拉方法的matlab实现及实例欧拉方法的Matlab实现程序如下：function[x,y]=euler(fun,x0,xfinal,y0,n);ifnargin5,n=50;endh=(xfinal-x0)/n;x(1)=x0;y(1)=y0;fori=1:nx(i+1)=x(i)+h;y(i+1)=y(i)+h*feval(fun,x(i),y(i));end例[4]求解初值问题'2(01)(0)1xyyxyy⎧=−⎪⎨⎪=⎩(9)解编写函数文件doty.m如下：functionf=doty(x,y);f=y-2*x/y;在Matlab命令窗口输入：[x,y]=euler('doty',0,1,1,10)便可得到结果：x=00.10000.20000.30000.40000.50000.60000.70000.80000.90001.0000y=11.10001.19181.27741.35821.43511.50901.58031.64981.71781.7848初值问题(6)的精确解为：12yx=+，此时对应nx的精确解()nyx可编程：y1=sqrt(1+2.*x)得到y1=11.09541.18321.26491.34161.41421.48321.54921.61251.67331.7321下面通过作图比较欧拉数值解与精确解之间的误差，编程：plot(x,abs(y-y1),'*')，可得到图2-1：00.10.20.30.40.50.60.70.80.9100.010.020.030.040.050.06图2-1欧拉数值解与精确解的误差图从图中可以看出欧拉方法精度很差。3.1改进欧拉方法的matlab实现为了编程方便，常把(5)式改写为：11(,)[(,)(,)]21()2pnnnqnnnnpnpqyyfxyhhyyfxyfxyyyy++⎧=+⎪⎪⎪=++⎨⎪⎪=+⎪⎩则Matlab实现程序为：function[x,y]=eulerpro(fun,x0,xfinal,y0,n);ifnargin5,n=50;endh=(xfinal-x0)/n;x(1)=x0;y(1)=y0;fori=1:nx(i+1)=x(i)+h;y1=y(i)+h*feval(fun,x(i),y(i));y2=y(i)+h*feval(fun,x(i+1),y1);(i+1)=(y1+y2)/2;end仍然选用上面的实例(6)，函数doty.m文件相同，在Matlab命令窗口中输入：[x,y]=eulerpro('doty',0,1,1,10)便可得到x=00.10000.20000.30000.40000.50000.60000.70000.80000.90001.0000y=11.09591.18411.26621.34341.41641.48601.55251.61651.67821.7379y1=11.09541.18321.26491.34161.41421.48321.54921.61251.67331.7321y为改进的欧拉数值解，y1为精确解。从上面的数据可看出，改进的欧拉方法有较好的精确度。作图plot(x,abs(y-y1),'*')可得到误差如图2-2：00.10.20.30.40.50.60.70.80.910123456x10-3图2-2改进的欧拉方法数值解与精确解的误差图4．总结本文通过编程实现了常微分方程初值问题数值解法中的欧拉方法及其改进的算法，并比较了其数值解与精确解之间的误差。可以看出欧拉方法得到的数值解精确度较差，而改进的欧拉方法得到的结果则相对较好。参考文献[1]王高熊，周之铭，朱思铭，等．数值计算原理［M］．北京：清华大学出版社，2000．[2]李庆扬，王能超，易大义．数值分析［M］．北京：清华大学出版社，施普林格出版社，2001．[3]李庆扬，关治，白峰杉．数值计算原理［M］．北京：清华大学出版社，2000．[4]李庆扬，王能超，易大义．数值分析［M］．武汉：华中科技大学出版社，2001．MatlabProgramingForEuler’sMethordandItsimprovedMethordJiXiuhaoCollegeofScience,LiaoningTechnicalUniversity,Fuxin,Liaoning(123000)AbstractThispapermakesmatlabprogramstosolvetheinitialvalueproblemofordinarydifferentialequationsusingEuler'smethordandEuler'simprovedmethord.Itshowstheaccurationsofthosemethordsusingfigures.ItcangiveusbetterunderstandingofEuler’smethordanditsimprovemethordforinitialvalueproblemofordinarydifferentialequations.Keywords:Euler'smethod;Euler'simprovedmethod;MatlabProgramming中国科技论文在线