MATLAB在求解常微分方程中的应用.docx

1、MATLAB在求解常微分方程中的应用M ATLA语言课程论文MATLAE在求解常微分方程中的应用姓名：袁学婷学号：12010245278专业：通信工程班级:通信一班指导老师：汤全武学院：物理电气信息学院完成日期：2012年12月10日MATLAB在求解常微分方程中的应用（袁学婷 12010245278 2010 通信班）【摘要】 MATLAB成为许多学科的解题工具，将 MATLAB融入其它课程的学习中，可以大大提高运算效率和准确性。随着计算机的普及和国民整体素质的提高， 科学计算将会更加的普及。 MATLAB在矩阵及数值计算、多项式和线形代数、符号数学的基本方法等方面都有较好的应用。本文概括地

2、介绍了MATLAB在求解常微分方程中的应用。【关键词】MATLAB 求解 常微分方程 解析解 数值解一、问题的提出自20世纪80年代以来，出现了多种科学计算语言，亦称数学软件，比较流行的有 MATLAB Mathematica、Maple等。因为他们具有功能强、效率高、简单易学等特点，在在 许多领域等到广泛应用。MATLA便是一种影响大、流行广的科学计算语言。MATLA 的语法 规则简单，更加贴近人的思维方式。MATLA是英文MATrix LABoratory（矩阵实验室）的缩写。自1984年由美国MathWorks 公司推向市场以来，得到了广泛的应用和发展。在欧美各高等院校MATLA已经成为

3、线性代 数、自动控制理论、数字信号处理、时间序列分析、动态系统仿真、图像处理等诸多课程 的基本教学工具，成为大学生、硕士生以及博士生必须掌握的基本技能。在设计研究单位 和工业部门，MATLA已被广泛的应用于研究和解决各种具体的工程问题。近年来， MATLAB在我国也开始流行，应用 MATLAB勺单位和个人急剧增加。可以预见， MATLAB将在我国科学研究和工程应用中发挥越来越大的作用。 现在简单的介绍一下MATLABE解微分方程中的应用。众所周知，只有对一些典型的常微分方程，才能求出它们的一般解表达式并用初始条件 确定表达式中的任意常数。所以能够求解的微分方程是十分有限的，特别是高阶方程（组）

4、这就要求我们必须研究微分方程（组）的解法，既要研究微分方程（组）的解析解 法（精确解），更要研究微分方程（组）的数值解法（近似解）对微分方程（组）的解析解法（精确解），Matlab有专门的函数可以用，本文将作一定 的介绍.二、常微分方程的概述1、微分方程的概念未知的函数以及它的某些阶的导数连同自变量都由一已知方程联系在一起的方程称为 微分方程。如果未知函数是一元函数，称为常微分方程。常微分方程的一般形式为F （t,y,y,y, ,y（n）二 0（1）2、常微分方程的解析解有些微分方程可直接通过积分求解 .例如，一解常系数常微分方程 df =y 1可化为dy =dt t 4y 1 ,两边积分可得

5、通解为y =ce -1.其中c为任意常数.有些常微分方程可用一些技巧,如分离变量法，积分因子法，常数变异法，降阶法等可化为可积分的方程而求得解析解.线性常微分方程的解满足叠加原理，从而他们的求解可归结为求一个特解和相应齐次 微分方程的通解一阶变系数线性微分方程总可用这一思路求得显式解。高阶线性常系数微 分方程可用特征根法求得相应齐次微分方程的基本解，再用常数变异法求特解。一阶常微分方程与高阶微分方程可以互化，已给一个 n阶方程设 =y,y2 二 y, = y(y(n) =f(t,y,yr ,y(nJJ)，可将上式化为一阶方程组乂= y2 yn 严 ynyn= f (t,yi,y2, , yn)

6、反过来，在许多情况下，一阶微分方程组也可化为高阶方程。所以一阶微分方程组与高阶常微分方程的理论与方法在许多方面是相通的，一阶常系数线性微分方程组也可用特 征根法求解。3、微分方程的数值解法除常系数线性微分方程可用特征根法求解，少数特殊方程可用初等积分法求解外，大部分微分方程无限世界，应用中主要依靠数值解法。考虑一阶常微分方程初值问题V(t) = f (t,y(t) )t0 tfiy(to) = yo (4)其中 y 二(yi, y2,，ym)，f 二(fl, f2,fm),y0 二(y10, y20,ym0). 所谓数值解法，就是寻求 y(t)在一系列离散节点to :：t :：t-tf上的近似

7、值丫小二0，1，n称h/t-tk为步 长，通常取为常量h。最简单的数值解法是Euler法。Euler法的思路极其简单：在节点出用差商近似代替导数y(tk 1)- y(tk)y(tk) 一 -h这样导出计算公式(称为Euler格式)yk 1 二 yhf (tk, yQ,k 二 0,1,2,他能求解各种形式的微分方程。Euler法也称折线法。4、解微分方程的MATLAB命令MATLAB中主要用dsolve求符号解析解，0de45,0de23,ode15s求数值解。在MATLAB，由函数dsolve()解决常微分方程(组)的求解问题，其具体格式如下:r = dsolve(eq1,eq2,., con

8、d1,cond2,., v)eq1,eq2,. 为微分方程或微分方程组，con d1,c on d2,., 是初始条件或边界条件，v是独立变量，默认的独立变量是t。函数dsolve用来解符号常微分方程、方程组，如果没有初始条件，则求出通解，如果 有初始条件，则求出特解。基于龙格-库塔法，MATLAB!供了求常微分方程的数值解的函数，一般调用格式为： t,y=ode23(file name,tspa n,yO)t,y=ode45(file name,tspa n,y0)其中file name是定义f(t,y)的函数文件名，该函数文件必须返回一个列向量。 ,tspa n形式为t0,tf, 表示求解

9、区间，y0是初始状态列向量。t和y分别给出时间限量和相应的状态向量。ezplot(x,y,tmin,tmax):符号函数的作图命令.x,y为关于参数t的符号函数，tmin,tmax 为t的取值范围.inline():建立一个内联函数.格式:inline(expr, var1,var2,)，注意括号里的表达 式要加引号.三、实例应用1、直接用Matlab求微分方程的解析解2 2问题一：求解微分方程 矽/ 垄的通解.dx 2x2解：方程求解的Matlab程序为y=dsolve(Dy-(xA2+yA2)/xA2/2,x) ; % 求岀的微分方程的解运行结果为：y=x*(-log(x)+2+C1”(-

10、log(x)+C1)问题二：求X2鱼2x ex的通解。dx解：方程求解的Matlab程序为y=dsolve(Dy*xA2+2*x*y-exp(x),x) % 求岀微分方程的解运行结果为：y=1/xA2*exp(x)+1/xA2*C1问题三：求微分方程xyy-ex =0在初始条件y(1) = 2e下的特解，并画出解函数的图形.解：方程求解的Matlab程序为：syms x y; % 定义x,y 为符号变量y=dsolve(x*Dy+y-exp(x)=0, y( 1)=2*exp(1), x); % 求岀的微分方程在初始条件y(1) =2e下的特解ezplot(y); % 作岀解函数的图象微分方程

11、的特解为：y=1/x*exp(x)+1/x* exp (1)即yJF，解函解的图形如图1所示。x图1解函数的图2EMFikt f匪 Tgk H卡D 序 口闿 A罟x- y问题四：求q 的通解。先2x y I dt解：方程求解的Matlab程序为：% 求解方程组x,y=dsolve( Dx=4*x - 2*y , Dy=2*x - y)运行结果为：x=C1+C2*exp(3*t) y=1/2*C2*exp(3*t)+2*C1问题五：求下列微分方程的解析解(2)(3)方程(1)求解的MATLAB?序为:clear;s=dsolve(Dy=3*y+4);运行结果为清屏求岀微分方程的解-4/a+exp

12、(3*t)*C1方程(2)求解的MATLAB?序为：clear; %s=dsolve(D2y-si n(2*x)+y,y(0)=0,Dy(0)=1,x) simplify(s) %运行结果为s =清屏求岀微分方程的解 以最简形式显示 s5/3*si n(x)-1/3*si n(2*x)ans =-1/3*(-5+2*cos(x)*si n(x)方程(3)求解的MATLAB?序为：clear; %s=dsolve(Df-f-g,Dg-g+f,f(0)=1,g(0)=1); simplify(s.f); %ssimplify(s.g);运行结果为：清屏求岀微分方程的解是一个结构ans =exp(t

13、)*cos(t)+exp(t)*s in(t)ans =-exp(t)*s in( t)+exp(t)*cos(t)问题六：利用 MATLAB求常微分方程的初值问题（x2）yJ2xy, y解。解：MATLAB?序为：y=dsolve(D2y*(1+xA2)-2*x*Dy=0,y(0)=1,Dy(0)=3,x)运行结果：求解微分方程y =1+3*x+xA3问题七：利用MATLAB常微分方程y-2yy = 0的解。求解微分方程（组）的初值问题的数值解（近似解）解：MATLAB?序为：y=dsolve(D4y-2*D3y+D2y,x)运行结果为：y =C1+C2*x+C3*exp(x)+C4*exp

14、(x)*x2、用ode23、ode45等求解常微分方程问题八：求解微分方程初值问题 = _2y 2 2x的数值解，求解范围为区间0, 0.5。卜（0）=1解：MATLAB程序为：fun=inlin e(-2*y+2*xA2+2*x,x,y);%定义函数x,y=ode23(fu n,0,0.5,1);%在区间0,0.5求岀微分方程的数值解x；%输岀y(0)=1时x的值y;%输岀y(0)=1时y的值plot(x,y,o-)；%画岀函数的图像运行结果为：xans =0.0000 0.0400 0.0900 0.1400 0.1900 0.24000.2900 0.3400 0.3900 0.4400

15、 0.4900 0.5000 yans =1.0000 0.9247 0.8434 0.7754 0.7199 0.67640.6440 0.6222 0.6105 0.6084 0.6154 0.6179图形结果如图2所示。图2 函数图象问题九：用Euler折线法求解微分方程初值问题2x2 , y的数值解（步长h取0.4），求解范围为区间0,2. 解：本问题的差分方程为k = 0,1,2, , n -1(其中:f(x,y)二 y 乡).y= Xk +h,y“ = yk +hf(Xk, yQ I相应的Matlab程序为：图3函数图像问题十：求解微分方程yy t 1, y(0)= 1,清屏求微分

16、方程的解析解 以最简形式显示 s先求解析解，再求数值解，并进行比较。由clear; %s=dsolve(Dy=-y+t+1,y(0)=1,t) ; %simplify(s) %运行结果：ans =t+exp(-t)可得解析解为厂te0下面再求其数值解，先编与 M文件fun8.mM 函数 fun8.mfun cti on f=fu n8(t,y);%编写 M文件fun8.mf=-y+t+1;%f的表达式再用命令clear;%清屏close;t=0:0.1:1;%确定求解区间y=t+exp(-t);%解析解plot(t,y);%画解析解的图形hold on;%保留已经画好的图形，如果下面再画图，两

17、个图形和并在一起t,y=ode45(fun8,0,1,1); %调用龙格库塔求解函数求解数值解plot(t,y);%画数值解图形标岀t,y 轴xlabel(t),ylabel(y)结果如图4所示。图4解析解与数值解由图4可见，解析解和数值解吻合得很好。运行结果如图5所示问题十一：用欧拉方法和龙格库塔方法求下列微分方程初值问题的数值解，画出解的图形解：MATLAB?序为：fun cti on f=f(x,y)%函数文件f=y+2*X；%f的表达式clc;clear;%清屏a=0;%求解区间b=1;x1,y_r=ode45(f,a b,1);%调用龙格库塔求解函数求解数值解;y=1;N=100;h

18、=(b-a)/N;%禾U用Euler方法求解x=a:h:b;%x的取值for i=1:N%i的取值y(i+1)=y(i)+h*f(x(i),y(i);%y的表达式end%循环结束plot(x1,y_r,r*,x,y,b+,x,3*exp(x)-2*x-2,k-); %画出数值解与真解图title（ 数值解与真解图）； %标注图形名称lege nd(RK4,Euler,真解); %用红色*在折线上标出真解xlabel(x);ylabel(y);%标岀x,y轴0图5数值解与真解图%建立函数文件%y的表达式%确定t的范围%t=0时y的值%求数值解%求精确解%绘出图形并比较问题十二：设有初值问题:，y

19、2 -1-2y = 4(t 1) ,0 汨y(0) = 2试求其数值解，并与精确解相比较(精确解为 y(t), t 1 1 )。解：(1)建立函数文件funt.mfun cti on y=fu nt(t,y) y=(yA2-t-2)/4/(t+1);(2)求解微分方程tO=O；tf=1O；y0=2;t,y=ode23(fu nt,tO,tf,yO);y1= sqrt(t+1)+1;plot(t,y,b,t,y1,r-);运行结果如图6所示。图6微分方程精确解与数值解的比较问题十三：已知一个二阶线性系统的微分方程为:x(0) =0,x(0) =1其中a=2,绘制系统的时间响应曲线和相平面图解：令

20、X2 =X,Xi二X ,则得到系统的状态方程；X2 = X1Ix1 = -ax22(0)=0,X1(0)=1建立一个函数文件sys.m:fun cti on xdot=sys(t,x) xdot=-2*x (2) ;x(1);取t0=0,tf=20,求微分方程；t0=0；tf=20;t,x=ode45(sys,t0,tf,1,0);t,xsubplot(1,2,1);plot(t,x(:,2); subplot(1,2,2);plot(x(:,2),x(:,1); axis equal%建立函数文件% xdot 的表达式%确定t的值%求数值解%输出结果%以子图形式绘出解的曲线%以子图形式绘出相

21、平面曲线运行结果为:ans =01.000000.00011.00000.00010.00011.00000.00010.00021.00000.00020.00021.00000.00020.00051.00000.00050.00071.00000.00070.00101.00000.00100.00121.00000.001219.1332-0.34980.661819.2670-0.51960.603619.4007-0.67080.523819.5344-0.79800.425319.6681-0.89680.311619.7511-0.94220.235219.8340-0.974

22、70.155619.9170-0.99370.073820.0000-0.9991-0.0090方程的时间相应及相平面曲线如图 7所示。图7时间相应及相平面曲线四、 结论利用MATLA语言求解微分方程的分析我们不难得出以下结论：1利用MATLAB勺dsolve函数求微分方程的解析解非常方便快捷。2、 龙格-库塔法求常微分方程的数值解准确方便。3、 在求解常微分方程的过程中，函数文件的使用可是程序简单许多，用起来方便快 捷。4、 MATLABt很强大的数学计算和绘图功能，把 MATLAB勺的数学计算功能与绘图功 能充分的应用，可以方便快捷的解决很多数学方面的问题。5 、在应用MATLAB勺过程中

23、，不但可以普及计算机的利用，充分地利用资源，更可以 极大的提高学习者的兴趣，为其他学科的学习提供了较好的计算工具，并且对未来更进一 的学习打下了良好的基础。五、 课程体会1通过运用MATLA可以很轻松的计算出一些方程（组）的解，绘制图形，以及更多 的如对高等数学上的应用等。只要在经过合理的编排后运用 MATLAB软件运行即可得出结果，如果所得结果并非自己想要的，还可以随时修改，直到得出想要结果，如本文的例题， 输入程序后，经过检查及修改无误后，再运行程序既可在瞬间的到准确的图形，使用起来 也非常方便快捷。2、通过在本次写论文的过程中，我真的是受益匪浅。以前也写过不少的论文，但也 没有注意论文正

24、规的写作格式，所以在写本篇论文的过程中，我不仅清楚地知道了论文的 写作格式，还学会了使用 word软件中许多特殊的用法。在我选择 MATLA论文的题目的过 程中，我通过查找大量的资料，才发现 MATLA涉及的范围非常的大，例如在高等数学中的 应用，在大学物理中的应用，在线性代数中的应用等，我才感受到 MATLA对我们是如此的实用。其次，在刚学习的MATLA 的时候，我就觉得它和我们大一第一学期学的 C语言有很 多相似之处，所以在老师讲的过程中，我就争取好好的听每一节课，每次在实验课上我也 认真的做好每一个实验的题目，在程序结果运行出来的时候，我真的特别的兴奋，此时我 感觉到了 MATLA功能的

25、强大性。就如本篇论文例举的实例，在平时解决此类问题的时候， 非常的麻烦，但通过运用 MATLAB军决此类问题，就显得非常的简单了。但遗憾的就是在老 师讲MATLA绘图时，我几乎听的是稀里糊涂，做作业时我也只会编写一些简单的二维和三 维的绘图程序，对于图形的精细处理和光照处理等一些复杂一点的操作编程我几乎是一点 都不懂。看到老师绘出的那些经过处理的图形，我十分的羡慕，并且听老师说在我们以后 所学到的课程中很多地方都要用到 MATLAB所以我一定要好好的学习 MATLAB争取让它为 我以后的学习中提供更大的帮助。参考文献1 刘卫国.MATLA程序设计与应用（第二版）M.北京：高等教育出版社，2006.2 同济大学数学系.高等数学（第六版）M.北京：高等教育出版社，2007.