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

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

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

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

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

5、常数变异法，降阶法等可化为可积分的方程而求得解析解.线性常微分方程的解满足叠加原理，从而他们的求解可归结为求一个特解和相应齐次微分方程的通解.一阶变系数线性微分方程总可用这一思路求得显式解。 高阶线性常系数微分方程可用特征根法求得相应齐次微分方程的基本解，再用常数变异法求特解。一阶常微分方程与高阶微分方程可以互化，已给一个n n阶方程设yi=y,y2=y,yn=y(y(n)=f(t,y,y,y(n),可将上式化为一阶方程组y y1=y2y y2=y3反过来，在许多情况下,yn=ynyn=f(t,y1,y2,yn)一阶微分方程组也可化为高阶方程。所以一阶微分方程组与高阶常微分方程的理论与方法在许

6、多方面是相通的， 一阶常系数线性微分方程组也可用特征根法求解。3、微分方程的数值解法除常系数线性微分方程可用特征根法求解，少数特殊方程可用初等积分法求解外，大部分微分方程无限世界，应用中主要依靠数值解法。考虑一阶常微分方程初值问题,y(t)=f(t,y(t)t0ttf,.,.y(ty(t0)= =y y0其中y y= =（y y1,y y2,，y ym）,f f= =（3 3。,，f fm）,y y。= =（y y10,y y2。,，y ym0）.所谓数值解法，就是寻求y在一系列离散节点t0t1tnwtf上的近似值丫卜*=01称限=。由-八为步长，通常取为常量h0最简单的数值解法是Euler法

7、。Euler法的思路极其简单：在节点出用差商近似代替导数y(tk1)-y(tk)y(tk)h这样导出计算公式（称为Euler格式）yki=ykhf(tk,yk),k=0,1,2,(6)他能求解各种形式的微分方程。Euler法也称折线法。4、解微分方程的MATLAB命令MATLAB中主要用dsolve求符号解析解，0de45,0de23,ode15s求数值解在MATLA中，由函数dsolve（）解决常微分方程（组）的求解问题，其具体格式如下:r=dsolve（eq1,eq2,，cond1,cond2,，v）eq1,eq2,.为微分方程或微分方程组，cond1,cond2,.,是初始条件或边界条件

8、，v是独立变量，默认的独立变量是t。函数dsolve用来解符号常微分方程、方程组，如果没有初始条件，则求出通解，如果有初始条件，则求出特解。基于龙格-库塔法，MATLA魄供了求常微分方程的数值解的函数，一般调用格式为：t,y=ode23(filename,tspan,y0)t,y=ode45(filename,tspan,y0)其中filename是定义f(t,y)的函数文件名，该函数文件必须返回一个列向量。,tspan形式为t0,tf,表示求解区间，y0是初始状态列向量。t和y分别给出时间限量和相应的状态向量。ezplot(x,y,tmin,tmax)：符号函数的作图命令.x,y为关于参数t

9、的符号函数，tmin,tmax为t的取值范围.inline。:建立一个内联函数.格式：inline(expr,var1,var2,-),注意括号里的表达式要加引号.三、实例应用1、直接用Matlab求微分方程的解析解.22问题一：求解微分方程的通解.dx2x解：方程求解的Matlab程序为y=dsolve(Dy-(xA2+yA2)/xA2/2,x);%求出的微分方程的解运行结果为：y=x*(-log(x)+2+C1)/(-log(x)+C1)问题二：求 x2电 2xy=ex的通解。dx解：方程求解的Matlab程序为y=dsolve(Dy*xA2+2*x*y-exp(x),x)%求出微分方程的

10、解运行结果为：y=1/xA2*exp(x)+1/xA2*C1问题三：求微分方程xy+y-ex=0在初始条件 y(1)=2e 下的特解，并画出解函数的图形.解：方程求解的Matlab程序为：symsxy;%定义x,y为符号变量y=dsolve(x*Dy+y-exp(x)=0,y(1)=2*exp(1),x);%求出的微分方程在初始条件y(1)=2e下的特解ezplot(y);%作出解函数的图象微分方程的特解为：y=1/x*exp(x)+1/x*exp(1)x即 y=e+e，解函解的图形如图1所示 x生=4x-2y问题四：求dj的通解。-2x-y.dt解：方程求解的Matlab程序为：x,y=ds

11、olve(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)y=3y4(2)y=sin(2x)-y,y(0)=0,y(0)=1f=fg,g=g-f,f(0)=1,g(0)=1方程(1)求解的MATLA勰序为：%求解方程组图1解函数的图问题六：利用MATLAB求常微分方程的初值问题(1(1+ +x2x2)y”=2xy,y yx x4=1,4=1,y yx4=x4=3 3的解。解：MATLA毓序为：y=dsolve(D2y*(1+xA2)-2*x*Dy=0,y(0)=1,Dy(0

12、)=3,x,)%求解微分方程运行结果：y=1+3*x+xA3问题七：利用MATLA家常微分方程y-2y“+y=0的解。解：MATLA毓序为：y=dsolve(D4y-2*D3y+D2y,x)%求解微分方程运行结果为：y=C1+C2*x+C3*exp(x)+C4*exp(x)*x2、用ode23、ode45等求解常微分方程(组)的初值问题的数值解(近似解).dy-2V2x22x问题八：求解微分方程初值问题 jdxy 的数值解，求解范围为区间0,0.5o.y(0)=1clear;%s=dsolve(Dy=3*y+4);%运行结果为s=-4/a+exp(3*t)*C1方程(2)求解的MATLA醒序为

13、：clear;%s=dsolve(D2y-sin(2*x)+y,y(0)=0,Dy(0)=1,x)%simplify(s)%运行结果为s=5/3*sin(x)-1/3*sin(2*x)ans=-1/3*(-5+2*cos(x)*sin(x)方程(3)求解的MATLA勰序为：clear;%s=dsolve(Df-f-g,Dg-g+f,f(0)=1,g(0)=1);%simplify(s.f);%ssimplify(s.g);运行结果为：ans=exp(t)*cos(t)+exp(t)*sin(t)ans=-exp(t)*sin(t)+exp(t)*cos(t)清屏求出微分方程的解清屏求出微分方程

14、的解以最简形式显示s清屏求出微分方程的解是一个结构0.00000.04000.09000.14000.19000.24000.29000.34000.39000.44000.49000.5000yans=1.00000.92470.84340.77540.71990.67640.64400.62220.61050.60840.61540.6179图形结果如图2所示。口小金。fe 入5立q/嚏n=_o图2函数图象问题九：用Euler折线法求解微分方程初值问题电学,dxyy(0)=1的数值解(步长h取0.4),求解范围为区间0,2.解：本问题的差分方程为解：MATLAB程序为：fun=inline

15、(-2*y+2*xA2+2*x,x,y);%x,y=ode23(fun,0,0,5,1);%x;%y;%plot(x,y,o-);%运行结果为：xans=定义函数在区间0,0.5求出微分方程的数值解输出y(0)=1时x的值输出y(0)=1时y的值画出函数的图像00050101O202503Q.3504046D501.00000.40001.40000.80002.12331.20003.11451.60004.45932.00006.3074图形结果如图3所示X0=0,y0=1,h=0.4,jX=Xk+h,yk+=yk+hf(Xk,yk)I相应的Matlab程序为：Clear;%f=sym(y

16、+2*x/yA2);%a=0;%b=2;%h=0.4;%n=(b-a)/h+1;%x=0;%xy=1;%yszj=x,y;%fori=1:n-1%y=y+h*subs(f,x,y,x,y);%yx=x+h;%xszj=szj;x,y;%end%Szj%plot(szj(:,1),szj(:,2)%运行结果为：szj=k=0,1,2,n-1(其中：f(x,y)=y2x).y清屏定义变量求解区间左端点求解区间右端点步长在求解区间内取n个值的初值为0的初值为1定义szj函数利用Euler方法求解确定i的取值范围的表达式的表达式调用szj函数结束输出szi画出图像问题十：求解微分方程y-yt1,y(0

17、)=1,先求解析解，再求数值解，并进行比较。由clear;%s=dsolve(Dy=-y+t+1,y(0)=1,t);%simplify(s)%运行结果：ans=t+exp(-t)可得解析解为y=t+e 上。 下面再求其数值解， 先编写M文件fun8.m编写M文件fun8.m的表达式清屏确定求解区间解析解画解析解的图形保留已经画好的图形，如果下面再画图，两个图形和并在一起调用龙格库塔求解函数求解数值解画数值解图形图3函数图像清屏求微分方程的解析解以最简形式显示sM函数fun8.mfunctionf=fun8(t,y);%f=-y+t+1;%f再用命令clear;%close;t=0:0.1:1

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

19、十二：设有初值问题：1_y2-t-2y=4(t1)，0，t三1y(0)=2试求其数值解，并与精确解相比较(精确解为解：(1)建立函数文件funt.mfunctiony=funt(t,y)y=(yA2-t-2)/4/(t+1);(2)求解微分方程t0=0；tf=10；y0=2;t,y=ode23(funt,t0,tf,y0);y1=sqrt(t+1)+1;plot(t,y,b.,t,y1,r-);运行结果如图6所示。y(t)=vt+1+1)。%建立函数文件%y的表达式%确定t的范围%t=0时y的值%求数值解%求精确解%绘出图形并比较图6微分方程精确解与数值解的比较X2=Xix1=-ax2x2(0

20、)=0,x1(0)=1建立一个函数文件sys.m:functionxdot=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);axisequal%建立函数文件%xdot的表达式%确定t的值%求数值解%输出结果%以子图形式绘出解的曲线%以子图形式绘出相平面曲线问题十三: 已知一个二阶线性系统的微分方程为:出2ax=0,a0 x(0)=0,x=1其中a=2,绘

21、制系统的时间响应曲线和相平面图。解：令x2=x,x1=x,则得到系统的状态方程;运行结果为：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.31

22、1619.7511-0.94220.235219.8340-0.97470.155619.9170-0.99370.073820.0000-0.9991-0.0090方程的时间相应及相平面曲线如图7所示。Qfigure1I圆FileEditVies*nsertTochWind网Hdp口/)Im 眼口0Noten的itoglbijrbutton郭 E 折1iq&lindddWX觉P.诉M口乂图7时间相应及相平面曲线四、结论利用MATLA第言求解微分方程的分析我们不难得出以下结论：1、利用MATLAB勺dsolve函数求微分方程的解析解非常方便快捷。2、龙格-库塔法求常微分方程的数值解准确方便。3、在求解常微分方程的过程中，函数文件的使用可是程序简单许多，用起来方便快捷。4、MATLABT很强大的数学计算和绘图功能，把MATLAB勺的数学计算功能与绘图功能充分的应用，可以方便快捷的解决很多数学方面的问题。5、在应用MATLAB勺过程中，不但可以普及计算机的利用，充分地利用资源，更可以极大的提高学习者的兴趣，为其他学科的学习提供了较好的计算工具，并且对未来更进一的学习打下了良好的基础。五、课程体会1、通过运用MATLABT以很轻松的计算出一些方程（组）的解，绘制图形，以及更多的如对高等数学上的应用等。只要在经过合理的编排后运