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

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

2、提出自20世纪80年代以来，出现了多种科学计算语言，亦称数学软件，比较流行的有MATLAB、Mathematica、Maple等。因为他们具有功能强、效率高、简单易学等特点，在在许多领域等到广泛应用。MATLAB便是一种影响大、流行广的科学计算语言。MATLAB的语法规则简单，更加贴近人的思维方式。MATLAB是英文MATrix LABoratory(矩阵实验室)的缩写。自1984年由美国MathWorks公司推向市场以来，得到了广泛的应用和发展。在欧美各高等院校MATLAB已经成为线性代数、自动控制理论、数字信号处理、时间序列分析、动态系统仿真、图像处理等诸多课程的基本教学工具，成为大学生、

3、硕士生以及博士生必须掌握的基本技能。在设计研究单位和工业部门，MATLAB已被广泛的应用于研究和解决各种具体的工程问题。近年来，MATLAB在我国也开始流行，应用MATLAB的单位和个人急剧增加。可以预见，MATLAB将在我国科学研究和工程应用中发挥越来越大的作用。现在简单的介绍一下MATLAB在解微分方程中的应用。众所周知，只有对一些典型的常微分方程，才能求出它们的一般解表达式并用初始条件确定表达式中的任意常数。所以能够求解的微分方程是十分有限的，特别是高阶方程（组）这就要求我们必须研究微分方程（组）的解法，既要研究微分方程（组）的解析解法（精确解），更要研究微分方程（组）的数值解法（近似解

4、）对微分方程（组）的解析解法(精确解)，Matlab 有专门的函数可以用，本文将作一定的介绍二、常微分方程的概述1、微分方程的概念未知的函数以及它的某些阶的导数连同自变量都由一已知方程联系在一起的方程称为微分方程。如果未知函数是一元函数，称为常微分方程。常微分方程的一般形式为 0 ),",',()(=nyyyytF. (1) 2、常微分方程的解析解有些微分方程可直接通过积分求解.例如,一解常系数常微分方程可化为,两边积分可得通解为.其中为任意常数.有些常微分方程可用一些技巧,如分离变量法,积分因子法,常数变异法,降阶法等可化为可积分的方程而求得解析解.线性常微分方程的解满足叠

5、加原理,从而他们的求解可归结为求一个特解和相应齐次微分方程的通解.一阶变系数线性微分方程总可用这一思路求得显式解。高阶线性常系数微分方程可用特征根法求得相应齐次微分方程的基本解，再用常数变异法求特解。一阶常微分方程与高阶微分方程可以互化，已给一个阶方程 (2)设，可将上式化为一阶方程组 (3) 反过来，在许多情况下，一阶微分方程组也可化为高阶方程。所以一阶微分方程组与高阶常微分方程的理论与方法在许多方面是相通的，一阶常系数线性微分方程组也可用特征根法求解。3、微分方程的数值解法 除常系数线性微分方程可用特征根法求解，少数特殊方程可用初等积分法求解外，大部分微分方程无限世界，应用中主要依靠数值解

6、法。考虑一阶常微分方程初值问题 (4)其中所谓数值解法，就是寻求在一系列离散节点上的近似值称为步长，通常取为常量。最简单的数值解法是Euler法。Euler法的思路极其简单：在节点出用差商近似代替导数 (5)这样导出计算公式（称为Euler格式） (6)他能求解各种形式的微分方程。Euler法也称折线法。4、解微分方程的MATLAB命令 MATLAB中主要用dsolve求符号解析解,0de45,0de23,ode15s求数值解 。 在MATLAB中，由函数dsolve()解决常微分方程（组）的求解问题，其具体格式如下：r = dsolve('eq1,eq2,.', 'c

7、ond1,cond2,.', 'v')'eq1,eq2,.'为微分方程或微分方程组，'cond1,cond2,.',是初始条件或边界条件，'v'是独立变量，默认的独立变量是't'。函数dsolve用来解符号常微分方程、方程组，如果没有初始条件，则求出通解，如果有初始条件，则求出特解。基于龙格-库塔法，MATLAB提供了求常微分方程的数值解的函数，一般调用格式为： t,y=ode23(filename,tspan,y0) t,y=ode45(filename,tspan,y0)其中filename是定义f(t

8、,y)的函数文件名，该函数文件必须返回一个列向量。,tspan形式为t0,tf,表示求解区间，y0是初始状态列向量。t和y分别给出时间限量和相应的状态向量。ezplot(x,y,tmin,tmax)：符号函数的作图命令x,y 为关于参数t 的符号函数，tmin,tmax 为 t 的取值范围inline()：建立一个内联函数格式：inline('expr', 'var1', 'var2',) ，注意括号里的表达式要加引号三、实例应用1、直接用 Matlab 求微分方程的解析解问题一：求解微分方程解：方程求解的 Matlab 程序为y=dsolve(

9、'Dy-(x2+y2)/x2/2','x') ; %求出的微分方程的解运行结果为：y=x*(-log(x)+2+C1)/(-log(x)+C1)问题二：求的通解。解：方程求解的 Matlab 程序为y=dsolve('Dy*x2+2*x*y-exp(x)','x') %求出微分方程的解运行结果为：y=1/x2*exp(x)+1/x2*C1问题三：求微分方程在初始条件下的特解，并画出解函数的图形解：方程求解的 Matlab 程序为：syms x y; %定义x,y为符号变量y=dsolve('x*Dy+y-exp(x)=0&

10、#39;,'y(1)=2*exp(1)','x'); %求出的微分方程在初始条 件 下的特解 ezplot(y); %作出解函数的图象微分方程的特解为：y= 1/x*exp(x)+1/x* exp (1) 即，解函解的图形如图1所示。图1 解函数的图问题四：求的通解。解：方程求解的 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问题五：求下列微分方程的解析解（1）（2）（3）方程（1）求解的MATLAB程序为：clear

11、; %清屏s=dsolve('Dy=3*y+4'); %求出微分方程的解运行结果为s = -4/a+exp(3*t)*C1方程（2）求解的MATLAB程序为：clear; %清屏s=dsolve('D2y-sin(2*x)+y','y(0)=0','Dy(0)=1','x') %求出微分方程的解simplify(s)  %以最简形式显示s运行结果为s = 5/3*sin(x)-1/3*sin(2*x)ans = -1/3*(-5+2*cos(x)*sin(x)方程（3）求解的MATLAB程序为：clear

12、; %清屏s=dsolve('Df-f-g','Dg-g+f','f(0)=1','g(0)=1'); %求出微分方程的解simplify(s.f); %s是一个结构simplify(s.g);运行结果为:ans = exp(t)*cos(t)+exp(t)*sin(t)ans = -exp(t)*sin(t)+exp(t)*cos(t)问题六：利用MATLAB求常微分方程的初值问题，的解。解：MATLAB程序为：y=dsolve('D2y*(1+x2)-2*x*Dy=0','y(0)=1,Dy(0)=3&

13、#39;,'x') %求解微分方程运行结果：y = 1+3*x+x3 问题七：利用MATLAB求常微分方程的解。解：MATLAB程序为：y=dsolve('D4y-2*D3y+D2y','x') %求解微分方程运行结果为:y = C1+C2*x+C3*exp(x)+C4*exp(x)*x 2、 用ode23、ode45等求解常微分方程(组)的初值问题的数值解(近似解)问题八：求解微分方程初值问题的数值解，求解范围为区间0, 0.5。解：MATLAB程序为：fun=inline('-2*y+2*x2+2*x','x'

14、,'y'); %定义函数x,y=ode23(fun,0,0.5,1); %在区间0,0.5求出微分方程的数值解x' %输出y(0)=1时x 的值y' %输出y(0)=1时y 的值plot(x,y,'o-') ; %画出函数的图像 运行结果为：x'ans = 0.0000 0.0400 0.0900 0.1400 0.1900 0.24000.2900 0.3400 0.3900 0.4400 0.4900 0.5000 y'ans = 1.0000 0.9247 0.8434 0.7754 0.7199 0.6764 0.6440

15、 0.6222 0.6105 0.6084 0.6154 0.6179图形结果如图 2所示。图2 函数图象问题九：用 Euler 折线法求解微分方程初值问题的数值解(步长h取0.4)，求解范围为区间0,2解：本问题的差分方程为相应的Matlab 程序为： Clear; %清屏 f=sym('y+2*x/y2'); %定义变量 a=0; %求解区间左端点 b=2; %求解区间右端点 h=0.4; %步长 n=(b-a)/h+1; %在求解区间内取n个值 x=0; %x的初值为0y=1; %y的初值为1szj=x,y; %定义szj函数for i=1:n-1 %利用Euler方法求

16、解 确定i的取值范围y=y+h*subs(f,'x','y',x,y); %y的表达式x=x+h; %x的表达式szj=szj;x,y; %调用szj函数end %结束Szj %输出sziplot(szj(:,1),szj(:,2) %画出图像运行结果为：szj = 0 1.0000 0.4000 1.4000 0.8000 2.1233 1.2000 3.1145 1.6000 4.4593 2.0000 6.3074图形结果如图3所示。 图3 函数图像问题十：求解微分方程先求解析解，再求数值解，并进行比较。由clear; %清屏s=dsolve('D

17、y=-y+t+1','y(0)=1','t') ; %求微分方程的解析解simplify(s) %以最简形式显示s运行结果：ans = t+exp(-t)可得解析解为。下面再求其数值解，先编写M文件fun8.mM函数fun8.mfunction f=fun8(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;  %保留已经画好的图形，如果下面再画图，两

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

19、e45('f',a b,1); %调用龙格库塔求解函数求解数值解；y(1)=1;N=100;h=(b-a)/N; % 利用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('数值解与真解图'); %标注图形名称legend('RK4','Euler&#

20、39;,'真解'); %用红色*在折线上标出真解xlabel('x');ylabel('y'); %标出x,y轴运行结果如图5所示 。图5 数值解与真解图问题十二：设有初值问题：试求其数值解，并与精确解相比较（精确解为）。解：（1）建立函数文件funt.m function y=funt(t,y) %建立函数文件 y=(y2-t-2)/4/(t+1); %y的表达式(2)求解微分方程 t0=0;tf=10; %确定t的范围y0=2; %t=0时y的值t,y=ode23('funt',t0,tf,y0); %求数值解y1=sqrt(

21、t+1)+1; %求精确解plot(t,y,'b.',t,y1,'r-'); %绘出图形并比较 运行结果如图6所示。图6 微分方程精确解与数值解的比较问题十三：已知一个二阶线性系统的微分方程为：其中a=2，绘制系统的时间响应曲线和相平面图。解：令则得到系统的状态方程;建立一个函数文件sys.m: function xdot=sys(t,x) %建立函数文件xdot=-2*x(2);x(1); % xdot的表达式取t0=0,tf=20,求微分方程;t0=0;tf=20; %确定t的值t,x=ode45('sys',t0,tf,1,0); %求数值

22、解t,x %输出结果subplot(1,2,1);plot(t,x(:,2); %以子图形式绘出解的曲线 subplot(1,2,2);plot(x(:,2),x(:,1); %以子图形式绘出相平面曲线 axis equal运行结果为：ans = 0 1.0000 0 0.0001 1.0000 0.0001 0.0001 1.0000 0.0001 0.0002 1.0000 0.0002 0.0002 1.0000 0.0002 0.0005 1.0000 0.0005 0.0007 1.0000 0.0007 0.0010 1.0000 0.00100.0012 1.0000 0.001

23、2 19.1332 -0.3498 0.661819.2670 -0.5196 0.603619.4007 -0.6708 0.523819.5344 -0.7980 0.425319.6681 -0.8968 0.311619.7511 -0.9422 0.235219.8340 -0.9747 0.155619.9170 -0.9937 0.073820.0000 -0.9991 -0.0090方程的时间相应及相平面曲线如图7所示。图7 时间相应及相平面曲线四、结论利用MATLAB语言求解微分方程的分析我们不难得出以下结论：1、 利用MATLAB的dsolve函数求微分方程的解析解非常方便快捷。2、 龙格-库塔法求常微分方程的数值解准确方便。3、 在求解常微分方程的过程中，函数文件的使用可是程序简单许多，用起来方便快捷。4、 MATLAB有很强大的数学计算和绘图功能，把MATLAB的的数学计算功能与绘图功能充分的应用，可以方便快捷的解决很多数学方面的问题。 5、在应用MATLAB的过程中，不但可以普及计算机的利用，充分地利用资源，更可以极大的提高学习者的兴趣，为其他学科的学习提供了较好的计算工具，并且对未来更进一的学习打下了良好的基础。五、课程体会1、通过运用MATLAB可以很轻松的计算出一些方程（组）的解，绘制图形，以及更多