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

1、常微分方程初值问题的欧拉方法及其改进的欧拉方法的 Matlab 实现纪秀浩辽宁工程技术大学理学院,辽宁阜新(123000E-mail :摘 要 :欧拉 (Euler方法及改进的欧拉方法是解决常微分方程初值问题常用的数值解法, 但 Matlab 的工具箱中没有 Euler 方法的功能函数。本文在简要介绍 Euler 方法及其改进的 Euler 方法的基础上, 通过编写 Matlab 程序实现两种数值解法, 并通过作图形式对比其精度, 加深对两种方法的认识。关键词:欧拉方法;改进的欧拉方法; matlab 实现1.引言常微分方程是解决工程实例的常用的工具 1,建立微分方程只是解决问题的第一步,通

2、常需要求出方程的解来说明实际现象, 并加以检验。 如果能得到解析形式的解固然是便于分 析和应用的, 但是我们知道, 虽然求解常微分方程有各种各样的解析方法, 但解析方法只能 用来求解一些特殊类型的方程,实际问题终归结出来的微分方程主要靠数值解法 2。数值解 法就是一个十分重要的手段,而欧拉方法又是数值解法最基础最常用的方法。2. 欧拉方法、改进的欧拉方法及 Matlab 实现下面主要讨论一阶常微分方程的初值问题,其一般形式为:' 00(, ( y f x y y x y = (1 我们知道,只要函数 (, f x y 适当光滑 譬如关于 y 满足利普希茨 (Lipschitz条件(,

3、(, f x y f x y L y y 理论上就可以保证初值问题 (1的解 ( y y x =存在并且唯一 3。所 谓 数 值 解 法 , 就 是 求 问 题 (1 在 某 些 离 散 点 01N a x x x b =<<=" 的 近 似 值012, , , N y y y y " 的方法。 012, , , N y y y y " 就称为问题 (1的数值解。 1n n n h x x +=成为 n x 到 n x 的步长,我们为了方便取为常量 h 。 2.1. 欧拉方法将微分方程离散化,用向前差商 1( ( n n y x y x h+代替微分

4、'( n y x ,代入 (1中的微分方程,可得1( (, (1, 2,3, n n n n y x y x f x y x n h+="化简可得1( ( (, ( (1, 2,3, n n n n y x y x f x y x hn +=+="中国 科技论文在线如果用 n y 近似 ( n y x 代入上式便可得到 1( n y x +的近似值 1n y +,计算式为:1(, (1, 2,3, n n n n y y f x y h n +=+=" (2 这样问题 (1的近似解可通过求解下面的差分初值问题:10(, (1,2,3, (n n n n y

5、 y f x y h n y y a +=+=" (3得到,按 (3式由初值 0y 可逐次求出 12, , y y " 。对于 Euler 公式 (3我们看到,当 1, 2, n =" 时公式右端的 n y 都是近似的,所以用它计算 的 1n y +会有累积误差,分析累积误差比较复杂,这里先讨论比较简单的所谓局部截断误差。假定用 (3式时右端的 n y 没有误差,即 ( nn y y x =那么由此算出1( (, ( n n n n y y x f x y x h +=+ (4局部截断误差指的是, 按 (4式计算由 n x 到 1n x +这一步的计算值 1n y

6、 +与精确值 1( n y x +之差11( n n y x y +。为了估计它,由 Taylor 展开得到的精确值 1( n y x +是2'''31( ( ( ( ( 2n n n n h y x y x hy x y x O h +=+ (5(4、 (5两式相减(注意到 '(, y f x y =得2''3211( ( ( ( 2n n n h y x y y x O h O h +=+ (6即局部截断误差是 2h 阶的,而数值算法的精度定义为:若一种算法的局部截断误差为1( p O h +,则称该算法具有 p 阶精度。显然 p 越大,方法

7、的精度越高。式 (6说明, Euler方法是一阶方法,因此它的精度不高。2.2 改进的欧拉方法用数值积分方法离散化问题 (1,两端积分可得11( ( (, ( (0,1, 2, n nx n n x y x y x f x y x dx n +="对右端积分使用梯形公式可得 ,111(, ( (, ( (, (2n nx n n n n x hf x y x dx f x y x f x y x + 再用 1, n n y y +代替 1(, ( n n y x y x +,则得计算公式111(, (, 2n n n n n n hy y f x y f x y +=+ (7很明显可

8、以注意到 (7式为隐式形式。改进的欧拉方法是先用欧拉公式求 1( n y x +的一个近似值 1n y +,称为预测值,然后用梯 形公式进行矫正并求得近似值 1n y +。即1111(, (, (, 2n n n n n n n n n n y y f x y h hy y f x y f x y +=+=+由 (7式可知,11232434( ( '( '(2'( "( "'(23!( '( ''( '''( 22'''( ( 12n n n n n n n n n n

9、n n hy x y x y x y x h h hy x y x y x h hy x y x hy x y x O h h y x O h +=+=+所以改进的欧拉方法是二阶的,精度较欧拉方法要高,实用性更加广泛。3.实验及结果分析3.1 欧拉方法的 matlab 实现及实例欧拉方法的 Matlab 实现程序如下: function x,y=euler(fun,x0,xfinal,y0,n; if nargin<5,n=50; endh=(xfinal-x0/n; x(1=x0;y(1=y0; for i=1:nx(i+1=x(i+h;y(i+1=y(i+h*feval(fun,x(

10、i,y(i; end例 4 求解初值问题'2(01(01x y y x y y =<<=(9解 编写函数文件 doty.m 如下: function f=doty(x,y; f=y-2*x/y;在 Matlab 命令窗口输入: x,y=euler('doty',0,1,1,10便可得到结果:x=0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.0000 y=1 1.1000 1.1918 1.2774 1.3582 1.435 11.5090 1.5803 1.6498 1.

11、7178 1.7848 初值问题 (6 的精确解为:y =, 此时对应 n x 的精确解 ( n y x 可编程:y1=sqrt(1+2.*x得到 y1=1 1.0954 1.1832 1.2649 1.3416 1.4142 1.4832 1.5492 1.6125 1.6733 1.7321下面通过作图比较欧拉数值解与精确解之间的误差, 编程:plot(x,abs(y-y1,'*', 可得到 图 2-1: 图 2-1 欧拉数值解与精确解的误差图从图中可以看出欧拉方法精度很差。3.1 改进欧拉方法的 matlab 实现为了编程方便,常把 (5式改写为:11(, (, (, 2

12、1( 2p n n n q nn n n p n p q y y f x y h h y y f x y f x y y y y +=+=+=+则 Matlab 实现程序为:function x,y=eulerpro(fun,x0,xfinal,y0,n; if nargin<5,n=50; endh=(xfinal-x0/n; x(1=x0;y(1=y0; for i=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;中国 科技论文在线y(i+1=(y1+y2/2; end仍然选用上面的实例 (

13、6,函数 doty.m 文件相同,在 Matlab 命令窗口中输入: x,y=eulerpro('doty',0,1,1,10便可得到x=0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.0000 y=1 1.0959 1.1841 1.2662 1.3434 1.4164 1.4860 1.5525 1.6165 1.6782 1.7379 y1=1 1.0954 1.1832 1.2649 1.3416 1.4142 1.4832 1.5492 1.6125 1.6733 1.7321 y 为改进的欧拉数值解, y1为精确解。 从上面的数