资源描述
免费查阅标准与论文：常微分方程初值问题的欧拉方法及其改进的欧拉方法的 Matlab 实现纪秀浩 辽宁工程技术大学理学院，辽宁阜新（123000） E-mail：jixiuhao摘要：欧拉(Euler)方法及改进的欧拉方法是解决常微分方程初值问题常用的数值解法， 但 Matlab 的工具箱中没有 Euler 方法的功能函数。本文在简要介绍 Euler 方法及其改进的 Euler 方法的基础上，通过编写 Matlab 程序实现两种数值解法，并通过作图形式对比其精度，加深对两种方法的认识。关键词：欧拉方法；改进的欧拉方法；matlab 实现1引言常微分方程是解决工程实例的常用的工具1，建立微分方程只是解决问题的第一步，通 常需要求出方程的解来说明实际现象，并加以检验。如果能得到解析形式的解固然是便于分 析和应用的，但是我们知道，虽然求解常微分方程有各种各样的解析方法，但解析方法只能 用来求解一些特殊类型的方程，实际问题终归结出来的微分方程主要靠数值解法2。数值解 法就是一个十分重要的手段，而欧拉方法又是数值解法最基础最常用的方法。2. 欧拉方法、改进的欧拉方法及 Matlab 实现下面主要讨论一阶常微分方程的初值问题，其一般形式为：-6- y =f ( x, y)(1) y( x0 ) = y0我们知道，只要函数 f ( x, y) 适当光滑譬如关于 y 满足利普希茨(Lipschitz)条件f ( x, y) f ( x, y) L y y理论上就可以保证初值问题(1)的解 y = y( x) 存在并且唯一3。所谓数值解 法，就 是求 问题 (1) 在某些离散 点 a = x0 x1 xN = b 的近似值 y0 , y1 , y2 , yN 的方法。 y0 , y1 , y2 , yN 就称为问题(1)的数值解。 hn = xn +1 xn 成为 xn 到xn 的步长，我们为了方便取为常量 h 。2.1. 欧拉方法2.1.1 欧拉方法y( xn+1 ) y( xn )将微分方程离散化，用向前差商代替微分 y ( x ) ，代入(1)中的微分方hn程，可得化简可得y( xn+1 ) y( xn ) =hf ( xn , y( xn )(n = 1, 2, 3,)y( xn +1 ) = y( xn ) + f ( xn , y( xn )h(n = 1, 2, 3,)如果用 yn 近似 y( xn ) 代入上式便可得到 y( xn+1 ) 的近似值 yn+1 ，计算式为：yn+1 = yn + f ( xn , yn )h(n = 1, 2, 3,)(2)这样问题(1)的近似解可通过求解下面的差分初值问题： y0 yn+1 = yn + f ( xn , yn )h= y(a)得到，按(3)式由初值 y0 可逐次求出 y1 , y2 , 。(n = 1, 2, 3,)(3)Eule 方法就是用差分方程初值问题(3)的解来近似微分方程初值问题(1)的解。即由公式(3)算出 y( xn ) 的近似值。这组公式求问题(1)的数值解就是著名的欧拉(Euler)公式。2.1.2 欧拉方法的误差估计对于Euler公式(3)我们看到，当 n = 1, 2, 时公式右端的 yn 都是近似的，所以用它计算的 yn+1 会有累积误差，分析累积误差比较复杂，这里先讨论比较简单的所谓局部截断误差。 假定用(3)式时右端的 yn 没有误差，即 yn = y( xn ) 那么由此算出yn+1 = y( xn ) + f ( xn , y( xn )h(4)局部截断误差指的是，按(4)式计算由 xn 到 xn+1 这一步的计算值 yn+1 与精确值 y( xn+1 ) 之差y( xn+1 ) yn+1 。为了估计它，由Taylor展开得到的精确值 y( xn+1 ) 是n2n+1 n ny( x) = y( x ) + hy ( x ) + h2y ( x ) + O(h3 )(5)(4)、(5)两式相减（注意到 y =f ( x, y) ）得2h 3 2y( xn+1 ) yn+1 =2 y ( xn ) + O(h) O(h )(6)即局部截断误差是 h2 阶的，而数值算法的精度定义为：若一种算法的局部截断误差为O(h p +1 ) ，则称该算法具有 p 阶精度。显然 p 越大，方法的精度越高。式(6)说明，Euler方法是一阶方法，因此它的精度不高。2.2 改进的欧拉方法2.2.1 改进的欧拉方法 用数值积分方法离散化问题(1)，两端积分可得xn+1ny( xn +1 ) y( xn ) = xf ( x, y( x)dx(n = 0,1, 2,)对右端积分使用梯形公式可得,xn+1hnx2f ( x, y( x)dx f ( xn , y( xn ) + f ( xn+1 , y( xn+1 )再用 yn , yn+1 代替 y( xn ), y( xn +1 ) ，则得计算公式y= y + hf xy+ f xy(7)n+1n ( n ,2n )(n+1 ,n+1 )很明显可以注意到(7)式为隐式形式。改进的欧拉方法是先用欧拉公式求 y( xn+1 ) 的一个近似值 yn +1 ，称为预测值，然后用梯 形公式进行矫正并求得近似值 yn+1 。即 yn+1 = yn + f ( xn , yn )hhy= y +f xy+ f xy(8) n+1n ( n ,2n )(n+1 ,n+1 )2.2.2 改进的欧拉方法的误差估计由(7)式可知，y( x) y( x ) h y ( x ) + y ( x)n+1 n 2n n+1h2= hy ( xn ) + 2y ( xn ) +h323! y ( xn ) h y( x) + y ( x) + hy ( x) + hy ( x) + O(h4 )2n n n 2n3= h y ( x) + O(h4 )12n所以改进的欧拉方法是二阶的，精度较欧拉方法要高，实用性更加广泛。3实验及结果分析3.1 欧拉方法的 matlab 实现及实例 欧拉方法的 Matlab 实现程序如下： function x,y=euler(fun,x0,xfinal,y0,n);if nargin5,n=50;endh=(xfinal-x0)/n; x(1)=x0;y(1)=y0; for i=1:n x(i+1)=x(i)+h;y(i+1)=y(i)+h*feval(fun,x(i),y(i);end例4 求解初值问题 y = y 2x(0 x 1)y y(0) = 1解 编写函数文件 doty.m 如下：function f=doty(x,y);f=y-2*x/y;在 Matlab 命令窗口输入：x,y=euler(doty,0,1,1,10)(9)便可得到结果：x=00.10000.20000.30000.40000.50000.60000.70000.80000.90001.0000 y=11.10001.19181.27741.35821.43511.50901.58031.64981.71781.7848初值问题(6)的精确解为：y =得到1 + 2 x ，此时对应 xn 的精确解 y( xn ) 可编程：y1=sqrt(1+2.*x)y1=1 1.09541.18321.26491.34161.41421.48321.54921.61251.67331.7321下面通过作图比较欧拉数值解与精确解之间的误差，编程：plot(x,abs(y-y1),*)，可得到 图 2-1：0.060.050.040.030.020.0100 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1图 2-1 欧拉数值解与精确解的误差图从图中可以看出欧拉方法精度很差。3.1 改进欧拉方法的 matlab 实现为了编程方便，常把(5)式改写为： y p = yn + f ( xn , yn )hh2 yq = yn + f ( xn , yn ) + f ( xn +1 , y p ) yn+1= 1 ( y2p + yq )则 Matlab 实现程序为：function x,y=eulerpro(fun,x0,xfinal,y0,n);if nargin5,n=50;endh=(xfinal-x0)/n; x(1)=x0;y(1)=y0; for i=1:n x(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);y(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=1 1.09541.18321.26491.34161.41421.48321.54921.61251.67331.7321y 为改进的欧拉数值解，y1 为精确解。从上面的数据可看出，改进的欧拉方法有较好的精确度。作图 plot(x,abs(y-y1),*)可得到误差如图 2-2：-3x 10654321 00 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 14总结图 2-2 改进的欧拉方法数值解与精确解的误差图本文通过编程实现了常微分方程初值问题数值解法中的欧拉方法及其改进的算法，并比 较了其数值解与精确解之间的误差。可以看出欧拉方法得到的数值解精确度较差，而改进的 欧拉方法得到的结果则相对较好。参考文献1 王高熊，周之铭，朱思铭，等数值计算原理M北京：清华大学出版社，20002 李庆扬，王能超，易大义 数值分析M北京：清华大学出版社，施普林格出版社，20013 李庆扬，关治，白峰杉数值计算原理M北京：清华大学出版社，20004 李庆扬，王能超，易大义数值分析M武汉：华中科技大学出版社，2001Matlab Programing For Eulers Methord and Its improvedMethordJi XiuhaoCollege of Science, Liaoning Technical University, Fuxin, Liaoning (123000)AbstractThis paper makes matlab programs to solve the initial value problem of ordinary differential equationsusing Eulers methord and Eulers improved methord.It shows the accurations of those methords using figures.It can give us better understanding of Eulers methord and its improve methord for initial value problem of ordinary differential equations.Keywords: Eulersmethod; Eulers improved method; Matlab Programming
展开阅读全文