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

《MATLAB语言》课程论文MATLAB在求解常微分方程中的应用姓名：袁学婷学号：12010245278专业：通信工程班级：通信一班指导老师：汤全武学院：物理电气信息学院完成日期：2012年12月10日MATLAB在求解常微分方程中的应用（袁学婷120102452782010通信班）【摘要】MATLAB成为许多学科的解题工具，将MATLAB融入其它课程的学习中，可以大大提高运算效率和准确性。随着计算机的普及和国民整体素质的提高，科学计算将会更加的普及。MATLAB在矩阵及数值计算、多项式和线形代数、符号数学的基本方法等方面都有较好的应用。本文概括地介绍了MATLAB在求解常微分方程中的应用。【关键词】MATLAB求解常微分方程解析解数值解一、问题的提出自20世纪80年代以来，出现了多种科学计算语言，亦称数学软件，比较流行的有MATLAB、Mathematica、Maple等。因为他们具有功能强、效率高、简单易学等特点，在在许多领域等到广泛应用。MATLAB便是一种影响大、流行广的科学计算语言。MATLAB的语法规则简单，更加贴近人的思维方式。MATLAB是英文MATrixLABoratory(矩阵实验室)的缩写。自1984年由美国MathWorks公司推向市场以来，得到了广泛的应用和发展。在欧美各高等院校MATLAB已经成为线性代数、自动控制理论、数字信号处理、时间序列分析、动态系统仿真、图像处理等诸多课程的基本教学工具，成为大学生、硕士生以及博士生必须掌握的基本技能。在设计研究单位和工业部门，MATLAB已被广泛的应用于研究和解决各种具体的工程问题。近年来，MATLAB在我国也开始流行，应用MATLAB的单位和个人急剧增加。可以预见，MATLAB将在我国科学研究和工程应用中发挥越来越大的作用。现在简单的介绍一下MATLAB在解微分方程中的应用。众所周知，只有对一些典型的常微分方程，才能求出它们的一般解表达式并用初始条件确定表达式中的任意常数。所以能够求解的微分方程是十分有限的，特别是高阶方程（组）．这就要求我们必须研究微分方程（组）的解法，既要研究微分方程（组）的解析解法（精确解），更要研究微分方程（组）的数值解法（近似解）．对微分方程（组）的解析解法(精确解)，Matlab有专门的函数可以用，本文将作一定的介绍．二、常微分方程的概述1、微分方程的概念未知的函数以及它的某些阶的导数连同自变量都由一已知方程联系在一起的方程称为微分方程。如果未知函数是一元函数，称为常微分方程。常微分方程的一般形式为(1)2、常微分方程的解析解有些微分方程可直接通过积分求解.例如,一解常系数常微分方程1ydtdy可化为dtydy1,两边积分可得通解为1tcey.其中c为任意常数.有些常微分方程可用一些技巧,如分离变量法,积分因子法,常数变异法,降阶法等可化为可积分的方程而求得解析解.线性常微分方程的解满足叠加原理,从而他们的求解可归结为求一个特解和相应齐次微分方程的通解.一阶变系数线性微分方程总可用这一思路求得显式解。高阶线性常系数微分方程可用特征根法求得相应齐次微分方程的基本解，再用常数变异法求特解。0),,,',,()(nyyyytF一阶常微分方程与高阶微分方程可以互化，已给一个n阶方程),,,',()1()(nnyyytfy(2)设)1(21,,',nnyyyyyy，可将上式化为一阶方程组),,,,(''''2113221nnnnyyytfyyyyyyy(3)反过来，在许多情况下，一阶微分方程组也可化为高阶方程。所以一阶微分方程组与高阶常微分方程的理论与方法在许多方面是相通的，一阶常系数线性微分方程组也可用特征根法求解。3、微分方程的数值解法除常系数线性微分方程可用特征根法求解，少数特殊方程可用初等积分法求解外，大部分微分方程无限世界，应用中主要依靠数值解法。考虑一阶常微分方程初值问题000)()),(,()('ytyttttytftyf(4)其中)'.,,,(,)',,,(,)',,,(0201002121mmmyyyyffffyyyy所谓数值解法，就是寻求)(ty在一系列离散节点fntttt10上的近似值nkyk,,1,0,称kkktth1为步长，通常取为常量h。最简单的数值解法是Euler法。Euler法的思路极其简单：在节点出用差商近似代替导数htytytykkk)()()('1(5)这样导出计算公式（称为Euler格式）,2,1,0),,(1kythfyykkkk(6)他能求解各种形式的微分方程。Euler法也称折线法。4、解微分方程的MATLAB命令MATLAB中主要用dsolve求符号解析解,0de45,0de23,ode15s求数值解。在MATLAB中，由函数dsolve()解决常微分方程（组）的求解问题，其具体格式如下：r=dsolve('eq1,eq2,...','cond1,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,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求微分方程的解析解问题一：求解微分方程的通解2222xyxdxdy．解：方程求解的Matlab程序为y=dsolve('Dy-(x^2+y^2)/x^2/2','x');%求出的微分方程的解运行结果为：y=x*(-log(x)+2+C1)/(-log(x)+C1)问题二：求xexydxdyx22的通解。解：方程求解的Matlab程序为y=dsolve('Dy*x^2+2*x*y-exp(x)','x')%求出微分方程的解运行结果为：y=1/x^2*exp(x)+1/x^2*C1问题三：求微分方程0'xeyxy在初始条件ey2)1(下的特解，并画出解函数的图形．解：方程求解的Matlab程序为：symsxy;%定义x,y为符号变量y=dsolve('x*Dy+y-exp(x)=0','y(1)=2*exp(1)','x');%求出的微分方程在初始条件ey2)1(下的特解ezplot(y);%作出解函数的图象微分方程的特解为：y=1/x*exp(x)+1/x*exp(1)即xeeyx，解函解的图形如图1所示。图1解函数的图问题四：求yxdtdyyxdtdx224的通解。解：方程求解的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）43'yy（2）1)0(',0)0(,)2sin(''yyyxy（3）1)0(',1)0(',','gffgggff方程（1）求解的MATLAB程序为：clear;%清屏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;%清屏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求常微分方程的初值问题2(1)''2'xyxy，01xy，0'3xy的解。解：MATLAB程序为：y=dsolve('D2y*(1+x^2)-2*x*Dy=0','y(0)=1,Dy(0)=3','x')%求解微分方程运行结果：y=1+3*x+x^3问题七：利用MATLAB求常微分方程(4)2'''''0yyy的解。解：MATLAB程序为：y=dsolve('D4y-2*D3y+D2y','x')%求解微分方程运行结果为:y=C1+C2*x+C3*exp(x)+C4*exp(x)*x2、用ode23、ode45等求解常微分方程(组)的初值问题的数值解(近似解)．问题八：求解微分方程初值问题1)0(2222yxxydxdy的数值解，求解范围为区间[0,0.5]。解：MATLAB程序为：fun=inline('-2*y+2*x^2+2*x','x','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.00000.04000.09000.14000.19000.24000.29000.34000.39000.44000.49000.5000y'ans=1.00000.92470.84340.77540.71990.67640.64400.62220.61050.60840.61540.6179图形结果如图2所示。图2函数图象问题九：用Euler折线法求解微分方程初值问题1)0(,22yyxydxdy的数值解(步长h取0.4)，求解范围为区间[0,2]．解：本问题的差分方程为1,,2,1,0).2),(),(,,4.0,1,021100nkyxyyxfyxhfyyhxxhyxkkkkkk（其中：相应的Matlab程序为：Clear;%清屏f=sym('y+2*x/y^2');%定义变量a=0;%求解区间左端点b=2;%求解区间右端点h=0.4;%步长n=(b-a)/h+1;%在求解区间内取n个值x=0;%x的初值为0y=1;%y的初值为1szj=[x,y];%定义szj函数fori=1:n-1%利用Euler方法求解确定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=01.00000.40001.40000.80002.12331.20003.11451.60004.45932.00006.3074图形结果如图3所示。图3函数图像问题十：求解微分方程,1)0(,1'ytyy先求解析解，再求数值解，