作者：吴建宏时间：2011.09.20
*********************************************作者简介：
吴建宏，男，毕业于哈尔滨工业大学（威海），现就读于同济大学，攻读汽车电子方向，有两年的MATLAB实践经验，个人喜欢建模，编程和电控方面的。真诚愿意和各位志同道合的朋友一起探讨交流，一起搭建广阔的知识平台。
*********************************************
PARTONE微分方程简介
一、微分方程基本概念
微分方程：未知函数以及它们的某些阶的导数连同自变量都由一已知方程联系在一起的方程称为微分方程。其中，未知函数的最高阶数称为微分方程的阶。
偏微分方程：如果未知函数是多元函数，那么这个微分方程就是偏微分方程。常微分方程：如果未知函数是一元函数，那么这个微分方程就是常微分方程。
二、微分方程的解法
高等数学里讲过，微分方程的解法有分离变量法，常数变易法，特征值法……而在解决工程实际问题时候，利用MATLAB求解微分方程的方法主要有符号解法和数值解法。其中，符号解法主要可以求解可以用特征值法求解的常系数微分方程和少数特殊的微分方程，有很大的局限性，而且运行时间较长，优点就是可以得到精确的数学表达式。大部分实际的微分方程不得不通过数值解法来求解，优点就是运行时间快，可以解决复杂的线性或者非线性微分方程。缺点就是解是近似值。
PARTTWO微分方程的符号解法（dsolve）
微分方程的符号解法主要是函数dsolve，这个函数用起来很简单，最重要的是要知道神马时候能够用dsolve，神马时候不能用，当运行出错的时候，该怎么处理。接下来进入正文。
一、语法
r=dsolve('eq1,eq2,...','cond1,cond2,...','v')r=dsolve('eq1','eq2',...,'cond1','cond2',...,'v')
二、语法详解
1.eq1,eq2……用来代表看上去能够用高等数学介绍过的手段求解出来的微分方程，如果实在不太熟悉或者忘了，没关系，后面会提供出现错误时候的处理措施；v代表自变量，默认值为时间t。边界条件/初始值用cond1,cond2,...来表示。
2.在eq1,eq2……中，用D来代表变量的微分，例如y用Dy表示，而y用D2y来表示。
3.初始值的边界条件用y(a)=b或者Dy(a)=b来表示。如果初始值的个数小于变量的个数，输出的结果中就会出现待定系数。
4.Dsolve能够接受的最大输入初始条件是12个
5.如果是一个等式一个输出，输出的结果是非线性的符号表达式；如果是多个等式和多个输出，将会按照左边输出参数[y1,y2……]中的顺序输出相应的符号表达式；如果是多等式，单输出，返回一个结构体数组。
出现错误时候怎么处理？如果dsolve不能找到有限解，它会试图去找方程的隐式解，当找到隐式解得时候，就会返回警告标志；如果既找不到有限解，也找不到隐式解，就会返回empty的错误警告\"warningWarning:explicitsolutioncouldnotbefound\"。在这两种情况下都只能通过使用数值解法来解决。
•
••
三、实例求解
1.求解方程dy/dt=ay+b>>dsolve('Dy=a*y+b')ans=
-b/a+exp(a*t)*C1
2.求解方程d2y/dx2=sin(2x)−y，y(0)=0，dy/dx>>dsolve('D2y=sin(2*x)-y','y(0)=0','Dy(0)=1','x')ans=
5/3*sin(x)-1/3*sin(2*x)
3.求解方程df/dt=f+g，dg/dt=g−f，df/dtt=0=1，dg/dtt=0=1>>y=dsolve('Df=f+g','Dg=g-f','Df(0)=1','Dg(0)=1')y=
f:[1x1sym]g:[1x1sym]>>y.fans=exp(t)*sin(t)
x=0
=0
>>y.gans=exp(t)*cos(t)
PARTTWO微分方程的数值解
由于微分方程的类型很多，所以就会出现很多类型的解法，我们不仅要掌握这些解决方程的函数，更要掌握这些函数应用的条件。函数odexx主要用来解决ODES一类的微分方程组，这类微分方程组主要是给定初值（自变量=0时的初值）的微分方程或者微分方程组，其中ode45可以解决大多数ODES问题，ode15s可以解决DAE（微
分方程+代数方程），ode15i可以解决IDE（隐式微分方程）。而bvp4c则可以解决含边界条件的微分方程或者微分方程组。
I.含边界条件的微分方程（组）解法bvp4c一、bvp4c解决问题的类型
函数bvp4c用来解决含有边界条件的微分方程（组），其中bvp：boundaryvalueproblems
二、语法
sol=bvp4c(odefun,bcfun,solinit)
sol=bvp4c(odefun,bcfun,solinit,options)
sol=bvp4c(odefun,bcfun,solinit,options,p1,p2...)
三、语法详解
1.odefun描述的是微分方程（组）f(x,y)，其标准形式：
dydx=odefun(x,y)
dydx=odefun(x,y,p1,p2,...)dydx=odefun(x,y,parameters)
dydx=odefun(x,y,parameters,p1,p2,...)
其中，参数x是与相应的x对应的标量，参数y是与y对应的列向量，parameters是个未知参数的向量，p1,p2,...是已知参数，输出dydx是一个列向量。2.bcfun用来表示边界条件。分为单边界条件和多边界条件。
（1）单边界条件的bcfun具有如下的标准形式：res=bcfun(ya,yb)
res=bcfun(ya,yb,p1,p2,...)res=bcfun(ya,yb,parameters)
res=bcfun(ya,yb,parameters,p1,p2,...)
其中ya和yb是与y(a)和y(b)对应的列向量，输出res是列向量。这儿要注意的是啊，a，b分别是微分方程的积分区间的左边界和右边界。即积分区间是[ab]。（2）多边界条件的bcfun
函数bcfun可以解决多边界条件的微分方程（组）。设a=a0 y(0)=4，y(1)=4.5，y(1)=5，y(2)=5.5 则bcfun就可以写成： functionres=bc(yleft,yright) res=[yleft(1)-4 yright(1)-4.5yleft(2)-5yright(2)-5.5]; 3.Solinit是包含初始猜测值的结构体，可以通过solinit=bvpinit(x,yinit,params)来产生solinit。其中x是初始的网格节点，y是输出的猜测值，如solinit.y(:,i)是节点solinit.x(i)处的猜测值。 4.Options是积分参数，可以通过bxpset设定，语法：options=bvpset('name1',value1,'name2',value2,...)options=bvpset(oldopts'name1',value1,...)options=bvpset(oldopts,newopts) 就是将'name1'等特性设定为value1等值。如果没有任何设定，则取默认值。常见的有'RelTol'-相对误差，默认值是（1e−3），可以进行设定；'AbsTol'-绝对误差，默认值是（1e−6）。 5.p1,p2是已知的参数，传到odefun和bcfun函数中。 6.sol=bvp4c(odefun,bcfun,solinit)求解形如dy/dx=f(x,y)的普通微分方程,给方程在区间[ab]，满足单边界条件： bc(y(a),y(b))=0 函数bvp4c也可以解决形如dy/dx=f(x,y,p)，bc(y(a),y(b),p)=0的未知参数p。7.利用函数deval和输出的sol可以估计出区间内任意一点x对应的函数值。sxint=deval(sol,xint) 四、范例分析 例1：接下面的微分方程： d2y/dt2+y=0，y(0)=0，y(4)=−2 解：先化成标准形式：⎧dy(1)=y(2)⎨ ⎩dy(2)=−y(1) 下来开始编程： functiondydx=myfun1(x,y)dydx=[y(2) -abs(y(1))]; functionres=myfun(ya,yb)res=[ya(1) yb(1)+2]; >>solinit=bvpinit(linspace(0,4,5),[10]);>>sol=bvp4c(@myfun1,@myfun,solinit);>>x=linspace(0,4);y=deval(sol,x); plot(x,y(1,:)); 运行结果图见下： solinit=bvpinit(linspace(0,4,5),[-10]); sol=bvp4c(@myfun1,@myfun,solinit);x=linspace(0,4);y=deval(sol,x);plot(x,y(1,:));运行结果图： 从上例中可以看出，对于多解得微分方程，设置合适的猜测值是必要的，来求出自己想要的正确结果。 II.含初值条件的微分方程（组）解法odexx 一、odexx解决问题的类型 函数odexx主要解决含初始条件的微分方程（组），常见的odexx有ode45，ode23，ode113，ode15s，ode23s，ode23t等，它们的用法大同小异，不同的是它们解决的问题不同。接下来会讲讲它们的共同用法，别且举几个特殊的odexx。 二、语法 [t,Y]=solver(odefun,tspan,y0) [t,Y]=solver(odefun,tspan,y0,options) [t,Y,TE,YE,IE]=solver(odefun,tspan,y0,options)sol=solver(odefun,[t0tf],y0...) 其中solver就是函数odexx。 三、语法详解 1.odefun用来表示微分方程（组）的函数，所有的odexx都可以用来解决形如 dy/dt=f(t,y)或者包含质量矩阵的微分方程（组）。定义的一般形式： functiondx=odefun(t,x) 2.tspan是积分区间向量[t0tf]。Odexx会在tspan(1)处赋初值，然后在tspan(1)到 tspan(end)区间上进行积分求解。要想的到指定时间点得函数值，主需要使 tspan=[t0,t1,⋯,tf]。对于形如[t0tf]的tspan，odexx会返回每个积分步长处的函数值， 对于超过两个元素的tspan，odexx会给出指定时间点处的函数值。Tspan中元素必须是按递增或者递减的顺序排列的。 3.语法中的y0是初始向量。 4.语法中options可以改变默认的积分特性。通过options=odeset('name1',value1,'name2',value2,...)来设定相应的特征。 在opyions中，'RelTol'和'AbsTol'与前边讲述的相同，这儿就不讲了；'Jacobian'对微分方程的可靠性和有效性有很重要的作用，所以可以对它进行设置，其值分别可以是constant和matrix。 Odexx也可以解决形如M(t,y)dy/dt=f(t,y)的微分方程，这时候就需要对质量矩阵进行设置，如果质量矩阵是常数，则应该使用'Mass'这个特征；如果质量矩阵与状态变量y无关，只与t有关，则将'MStateDependence'设置为'none';如果质量矩阵受y的影响微弱，将'MStateDependence'设置为'weak'，否则设置为'strong'；如果M是非奇异的质量矩阵，则是DAES问题，可以使用ode15s，ode23t来解决。 5.下面列出了odexx的输出参数：t:时间点控制向量 y：解数组，y的每列与t相对应。 6.[t,Y]=solver(odefun,tspan,y0)在积分区间tspan=[t0tf]上，对含有初值向量y0，形如dy/dt=f(t,y)的微分方程进行求解的。 7.各种算法的优缺点见下： 特性 RelTol, AbsTol MaxStep, InitialStepJacobianMass MStateDependence MvPatternMassSingularode45ode23ode113ode15sode23t可以设定可以设定可以设定可以设定可以设定可以设定可以设定可以设定可以设定可以设定 可以设定可以设定 可以设定可以设定可以设定可以设定可以设定可以设定可以设定可以设定可以设定 可以设定可以设定可以设定可以设定 ode23tb可以设定可以设定可以设定可以设定可以设定可以设定 ode23s可以设定可以设定可以设定可以设定 从上边的特性设置来看，凡是不能设定的就是不能够参与解决相应类型的微分方程。问题类型准确度 非刚性中等 非刚性低 非刚性低到高 刚性刚性低到中等低 刚性低 刚性低 四、几种典型的odexx 1.ode45，在大多数情况下是作为第一选择的，下面看个例子：解微分方程：d2x/dt2−1000(1−x2)dx/dt−x=0，x(0)=0，dx/dt首先需要化成odes的标准形式：⎧x1=x⎨ ⎩x2=dx/dt⎧dx1/dt=x2⎨2 ⎩dx2/dt=x1+1000(1−x1)x2 t−0 =1 接下来开始编写函数：functiondx=myfun(t,x)dx=[x(2) x(1)+1000*(1-x(1)^2)*x(2)];>>[t,y]=ode45(@myfun,[010],[01]);>>plot(t,y(:,1))运行结果： 2.ode15s可以用来解决ode45解决不了的一些问题，包括微分代数方程DAE。微分代数方程其实就是在微分方程的基础上加入了代数方程。期一般形式可写成M(t,y)dy/dt=f(t,y)。接下来举个例子。 求解方程 ⎧dx1/dt=−0.2x1+4x2x3+0.5x1x2⎪2 ⎨dx2/dt=−2x2−6x2x3+3x1x2⎪x+x+x=4⎩123 的数值解，其中 x1(0)=0.8,x2(0)=x3(0)=0.4 首先把它转化为标准形式： ⎡100⎤⎡dx1/dt⎤⎡−0.2x1+4x2x3+0.5x1x2⎤⎢010⎥⎢dx/dt⎥=⎢−2x2−6xx+3xx⎥ 22312⎥⎢⎥⎢2⎥⎢ ⎢⎥x1+x2+x3−4⎣000⎥⎦⎢⎣dx3/dt⎥⎦⎢⎣⎦ 编写函数程序>> myfun=@(t,x)[-0.2*x(1)+4*x(2)*x(3)+0.5*x(1)*x(2);-2*x(2)^2-6*x(2)*x(3)+3*x(1)*x(2);x(1)+x(2)+x(3)-4];>>M=[100;010;000];>>options=odeset('mass',M); >>[t,x]=ode15s(myfun,[010],[0.8,0.4,0.4],options);>>plot(t,x) >>legend('x1(t)','x2(t)','x3(t)')运行结果见下： 3.当遇到隐式微分方程的时候，怎么办呢？所谓的隐式微分方程就是没法化成上述标准型的微分方程，我推荐先考虑使用solve函数解出二阶的显示解。这样便可套用标准型了。举个例子： 2222⎧⎪dx/dt+4xdy/dt=7dy/dt解⎨2微分方程222 ⎪⎩dx/dt+2(dx/dt)dy/dt+xdy/dt−y=8 我们可以首先使用solve函数解出d2x/dt2和d2y/dt2的符号表达式，见下：>>[d2x,d2y]=solve('d2x+4*x*dy-7*d2y','d2x+2*dx*d2y+x*dy-y-8','d2x','d2y')d2x= -(7*x*dy+8*x*dy*dx-7*y-56)/(7+2*dx) d2y= (3*x*dy+y+8)/(7+2*dx) 然后便可利用ode45进行求解了。当然solve也不是万能的，有时候根本求不出二阶的显式来的，这时候就需要fsolve来大展身手了，举个例子： 2222222⎧⎪dx/dtsin(dy/dt)+(dy/dt)=−4xy+6x(dx/dt)dy/dt求解⎨222222 ⎪⎩x(dx/dt)(dy/dt)+cos(dy/dt)=3ydx/dt初始值为： x0=[1;0;0;1] 很明显对于这种方程来说，用solve是无能为力了。这时候只能启用数值解法fsolve， 解法见下： 首先写出状态变量表达式：⎧dx1/dt=x2⎧x1=x⎪⎪x=dx/dt2 ⎪2⎪(dx2/dt)(sinx3)+(dx4/dt)=−4x1x3+6x1(dx2/dt)x4⎨⎨x=y⎪3⎪x3=x(4)⎪⎪x1(dx2/dt)(dx4/dt)+cos(dx4/dt)=3x3x2⎩x4=dy/dt⎩编写程序见下： functiondx=odefun(t,x) fun=@dxy[dxy(1)*sin(x(3))+dxy(2)^2+4*x(1)*x(3)-6*x(1)*dxy(1)*x(4);x(1)*dxy(1)dxy(2)+cos(dxy(2))-3*x(3)*x(2)]; options=optimset('display','off');y=fsolve(fun,x([1,3]),options);dx=[x(2);y(1);x(4);y(2)];主函数程序： >>[t,y]=ode45(@odefun,[03],[1;0;0;1]);>>plot(t,y) >>legend('x','x''','y','y''')运行结果见下： 因篇幅问题不能全部显示，请点此查看更多更全内容