欧拉法与Matlab数值求解

《欧拉法与Matlab数值求解》由会员上传分享，免费在线阅读，更多相关内容在行业资料-天天文库。

1、第8卷第1期武汉交通职业学院学报2006年3月Vol.8No.1JournalofWuhanJiaotongPolytechnicMar.20063欧拉法与Matlab数值求解陈誌敏(武汉工交职业学院,湖北武汉430205)摘要:欧拉法在常微分方程近似数值求解,特别是理论分析中占有很重要的地位,本文通过Matlab的平台,显示执行欧拉法各步骤的细节。由于数学软件Matlab自身不带有欧拉法执行程序,笔者编写了三个欧拉法计算程序,以满足常微分方程近似数值求解和理论分析的需要。关键词:欧拉法;Matlab;数值法;初值问题中图分类号:O175文献标识

2、码:A文章编号:1672-9846(2006)01-0067-041引言法,它特别适合于快速编程,尽管精度不是很高。由于大多数常微分方程初值问题不能求得解欧拉法包含以下三种模式:前进欧拉法、改进欧拉析解,所以采用其它近似方法求解是十分重要的。法、后退欧拉法。近似方法,如级数解法、逐步逼近法等,可以给出2.1前进欧拉法解的近似表达式,通常称为近似解析方法。还有前进欧拉法是基于向前微分近似[1]一类近似方法称为数值方法,它可以给出解在(yn+1-yn)/h≈y′n一些离散点上的近似值。本文介绍一种最古老变形为:[2]的,也是最简单的方法,称为切线法或

3、欧拉法。yi+1=yi+hf(yi,xi)(2)欧拉法以固定步长h生成近似解。由于欧拉法在y(x0)=y0(i=0,1,2,⋯,n-1)常微分方程近似数值求解,特别是理论分析中占其中y′i=f(yi,xi),以n循环,得到:y1=y0+hy′0有很重要的地位,本文通过Matlab的平台,显示y1=y0+hf(y0,x0)执行欧拉法各步骤的细节。y2=y1+hf(y1,x1)数学软件Matlab自身不带有欧拉法执行程y3=y2+hf(y2,x2)序,笔者编写了三个欧拉法计算程序,即创造三个⋯⋯⋯用户自定义函数文件,以满足常微分方程近似数yn=yn-

4、1+hf(yn-1,xn-1)2.2改进欧拉法值求解和理论分析的需要。改进欧拉法比前进欧拉法的精度高,而且更2欧拉法稳定。它是对y′=f(y,x)的解利用梯形法则而我们讨论一阶常微分方程,而高阶常微分方程可化为一阶常微分方程(组)。一阶初值问题得到的。即一般可以写成h)+f(y)]yi+1=yi+[f(yi+1,xi+1i,xi2(3)dydx=y(x,y),x∈[x0,X];(1)y(x0)=y0(i=0,1,2,⋯,n-1)y(x0)=y02.3后退欧拉法欧拉法是数值求解常微分方程最简单的方后退欧拉法是基于后退微分近似并表示成3收稿日期:20

5、06-02-02作者简介:陈誌敏(1952-),男,广东潮阳人,武汉工交职业学院副教授。—67—陈誌敏:欧拉法与Matlab数值求解yi+1=yi+hf(yi+1,xi+1)[x,y]=euler_forward(f,xinit,yinit,(4)xfinal,n)y(x0)=y0(i=0,1,2,⋯,n-1)h=(xfinal-xinit)/n;后退欧拉法的精度和前进欧拉法的基本相同。x=[xinitzeros(1,n)];3基本原理[3,6]y=[yinitzeros(1,n)];3.1唯一性定理fori=1∶nx(i+1)=x(i)+h;设

6、初值问题(1)满足:y(i+1)=y(i)+h3f(x(i),y(i));a)f在区域R,x0≤x≤X,y<∞内连续endb)f关于y满足Lipschitz条件:常数Lend>0,对②euler_modified.m(改进欧拉法)(8)function(x,y1),(x,y2)∈R[x,y]=euler_modified(f,xinit,yinit,│f(x,y1)-f(x,y2)│≤L│y1-y2│(5)xfinal,n)那么(1)就存在唯一的解。h=(xfinal-xinit)/n;3.2截断误差估计x=[xinitzeros(1,n)];估

7、计近似替代使数值解与真解之间产生误y=[yinitzeros(1,n)];fori=1∶n差。截断误差x(i+1)=x(i)+h;∈m=ym-y(xm)(6)ynew=y(i)+h3f(x(i),y(i));其中ym是式(2)的解,y(x)是式(1)的解。y(i+1)=y(i)+(h/2)3(f(x(i),y(i))+f(x(i+1),ynew));3.3收敛性end当步长h取得充分小时,数值解能否足够精end确地逼近初值问题的真解。若limmax│ym-③euler_backward.m(后退欧拉法)(9)h→01≤m≤ny(xm)│=0成立,

8、则Euler方法是收敛的。function[x,y]=euler_backward(f,xinit,yinit,3.4稳定性xfinal