微分方程在MATLAB中的实现
作者：吴建宏
时间：2011.09.20
*********************************************
作者简介：
吴建宏，男，毕业于哈尔滨工业大学（威海），现就读于同济大学，攻读汽车电子方向，有两年的MATLAB实践经验，个人喜欢建模，编程和电控方面的。真诚愿意和各位志同道合的朋友一起探讨交流，一起搭建广阔的知识平台。
*********************************************
PART ONE 微分方程简介
一、微分方程基本概念
微分方程：未知函数以及它们的某些阶的导数连同自变量都由一已知方程联系在一起的方程称为微分方 程。其中，未知函数的最高阶数称为微分方程的阶。
偏微分方程：如果未知函数是多元函数，那么这个微分方程就是偏微分方程。常微分方程：如果未知函数是一元函数，那么这个微分方程就是常微分方程。
高等数学里讲过，微分方程的解法有分离变量法，常数变易法，特征值法……而在解决工程实际问题时候，利用MATLAB 求解微分方程的方法主要有符号解法和数值解法。其中，符号解法主要可以求解可以用特征值法求解的常系数微分方程和少数特殊的微分方程，有很大的局限性，而且运行时间较长，优点就是可以得到精确的数学表达式。大部分实际的微分方程不得不通过数值解法来求解，优点就是运行时间快，可以解决复杂的线性或者非线性微分方程。缺点就是解是近似值。 PART TWO 微分方程的符号解法（dsolve ）
微分方程的符号解法主要是函数dsolve ，这个函数用起来很简单，最重要的是要知道神马时候能够用dsolve ，神马时候不能用，当运行出错的时候，该怎么处理。接下来进入正文。
一、语法
r =dsolve('eq1,eq2,...','cond1,cond2,...','v')
r =dsolve('eq1','eq2',...,'cond1','cond2',...,'v')
maxstep二、语法详解
1.eq1,eq2……用来代表看上去能够用高等数学介绍过的手段求解出来的微分方程，如果实在不太熟悉或者忘了，没关系，后面会提供出现错误时候的处理措施；v 代表自变量，默认值为时间t 。边界条件/初始值用cond1,cond2,...来表示。
2.在eq1,eq2……中，用D 来代表变量的微分，例如•y 用Dy 表示，而•
•y 用y D 2来表示。
3.初始值的边界条件用b a y =)(或者b a Dy =)(来表示。如果初始值的个数小于变量的个数，输出的结果中就会出现待定系数。 4.Dsolve 能够接受的最大输入初始条件是12个
5.如果是一个等式一个输出，输出的结果是非线性的符号表达式；如果是多个等式和多个输出，将会按照左边输出参数[y1,y2……]中的顺序输出相应的符号表达式；如果是多等式，单输出，返回一个结构体数组。
出现错误时候怎么处理？如果dsolve 不能到有限解，它会试图去方程的隐式解，当到隐式解得时候，就会返回警告标志；如果既不到有限解，也不到隐式解，就会返回empty 的错误警告"warning Warning:explicit solution could not be found"。在这两种情况下都只能通过使用数值解法来解决。
三、实例求解
1.求解方程b
ay dt dy +=/>>dsolve('Dy=a*y+b')
ans =
-b/a+exp(a*t)*C1
2.求解方程y x dx y d −=)2sin(/22，0)0(=y ，0
/0==x dx
dy >>dsolve('D2y=sin(2*x)-y','y(0)=0','Dy(0)=1','x')
ans =
5/3*sin(x)-1/3*sin(2*x)
3.求解方程g f dt df +=/，f g dt dg −=/，1/0==t dt df ，1
/0==t dt dg >>y=dsolve('Df=f+g','Dg=g-f','Df(0)=1','Dg(0)=1')
y =
f:[1x1sym]
g:[1x1sym]
>>y.f
ans =
exp(t)*sin(t)>>y.g
ans =
exp(t)*cos(t)
PART TWO 微分方程的数值解
由于微分方程的类型很多，所以就会出现很多类型的解法，我们不仅要掌握这些解决方程的函数，更要掌握这些函数应用的条件。函数odexx 主要用来解决ODES 一类的微分方程组，这类微分方程组主要是给定初值（自变量=0时的初值）的微分方程或者微分方程组，其中ode45可以解决大多数ODES 问题，ode15s 可以解决DAE （微
分方程+代数方程），ode15i 可以解决IDE （隐式微分方程）。而bvp4c 则可以解决含边界条件的微分方程或者微分方程组。
I.含边界条件的微分方程（组）解法bvp4c
一、bvp4c 解决问题的类型
函数bvp4c 用来解决含有边界条件的微分方程（组），其中bvp ：boundary value problems
二、语法
sol =bvp4c(odefun,bcfun,solinit)
sol =bvp4c(odefun,bcfun,solinit,options)
sol =bvp4c(odefun,bcfun,solinit,options,)
三、语法详解
1.odefun 描述的是微分方程（组）),(y x f ，其标准形式：
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 是与)(a y 和)(b y 对应的列向量，输出res 是列向量。这儿要注意的是啊，a ，b 分别是微分方程的积分区间的左边界和右边界。即积分区间是[a b]。
（2）多边界条件的bcfun
函数bcfun 可以解决多边界条件的微分方程（组）。设b a a a a a n =<<<<=⋯210是积分区间[a b]的边界点。此时，函数标准形式为：
res =bcfun(yleft,yright)
其中yleft(:,k)代表],[1k k a a −左边界对应的函数值，而yright(:,k)代表],[1k k a a −右边界对应的函数值。例如，给出式子的区间210<<，且有：
4)0(=y ，5.4)1(=y ，5)1(=y ，5
.5)2(=y 则bcfun 就可以写成：
function res =bc(yleft,yright)
res =[yleft(1)-4
yright(1)-4.5
yleft(2)-5
yright(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'-相对误差，默认值是（31−e ），可以进行设定；'AbsTol'-绝对误差，默认值是（61−e ）。
5.p1,p2是已知的参数，传到odefun 和bcfun 函数中。
6.sol =bvp4c(odefun,bcfun,solinit)求解形如),(/y x f dx dy =的普通微分方程，给方程在区间[a b]，满足单边界条件：
))(),((=b y a y bc 函数bvp4c 也可以解决形如),,(/p y x f dx dy =，0)),(),((=p b y a y bc 的未知参数p 。
7.利用函数deval 和输出的sol 可以估计出区间内任意一点x 对应的函数值。sxint =deval(sol,xint)
四、范例分析
例1：接下面的微分方程：
0/22=+y dt y d ，0)0(=y ，2
)4(−=y 解：先化成标准形式：
⎩⎨⎧−==)
1()2()2()1(y dy y dy 下来开始编程：
function dydx=myfun1(x,y)
dydx=[y(2)
-abs(y(1))];
function res=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...)
豫公网安备41160202000603
站长QQ:729038198
关于我们
投诉建议