常微分方程组的MATLAB求解方法

1、一、常微分方程组(ODEs)简介 11. 简谐振动12 电路 Vander Pol方程 13. 生物种群的Volterra-Lotka方程 24. 蝴蝶效应Lorenz方程2二、MATLAB数值求解0DEs的方法 31. 多变量常微分方程组的求解42. 高阶常微分方程如何表示？ 43. 相图和极限环怎么绘制？ 4个人在学习自动控制原理、现代控制理论、非线性动力学等课程时，经常遇到求解常微分方程组的问 题。很多人知道MATLAB是简便易行的一个工具，但是不会调用它自带的ode求解器，往往还在自己 编写单步欧拉法的程序，不仅求解精度差，而且程序不规范，还浪费了大量时间。以下我就工程中常 见的一些非

2、线性系统，利用MATLAB自带的求解器，说明一下如何求解ODE方程组、以及如何绘制 相轨迹和极限环的问题。供相尖专业工科大学生参考和借鉴。一、常微分方程组(ODEs)简介以下列出了一些较为著名的非线性动力学系统的数学表达式，大都是由常微分方程组表达的。这种形式 在工程中应用非常广泛，如力学中的非线性振动、航天领域的弹道计算、控制工程中的非线性系统等， 由于自然界的大多数现象都表现出非线性，因此对于该种动力系统的研究以及微分方程的求解也具有重大 的意义。以下列出一些工程应用中常见的一些由ODE方程组所描述的动力系统。1. 简谐振动该式是一个2阶非线性常微分方程。2. 电路Vander Pol方程

3、2Fig 1VanderPol系统时域响应3. 生物种群的Volterra丄otka方程Fig 2. Volterra丄otka方程时域响应（左）Fig 3非线性动力学方程的极限环（右）左图的捕食者猎物随时间变化的曲线表现出强烈的非线性，而状态变量x、y的变化却呈现出一个规则的鹅卵石状。4. 蝴蝶效应Lorenz方程01 0203040 SO 60604020x020-40Fig 4.洛伦兹方程系统3维状态的时域响应(左亠s小“亠亠、一Fig 5系统的3维空间相图：蝴蝶效应(右)世微分方程组的魅力就在于，它可以表现出使人意想不到的规律性与自然之美。按照数值积分的方法对3个状态 变量xyz进行求

4、解，由于该动力系统较为敏感，因此需要以尽可能小的积分步长来求 解。求解所得的相图在3维呈现出一个蝴蝶状。MATLAB数值求解ODEs的方法对于单个的常微分方程，给定初值和求解区间，同学们一定会求解。如以下形式的一阶非 线性微分方程，典型的方法有：单步Euler法、4阶Runge-Kutta积分法等。大多数 的数值分析或计算方法教材上均有讲述。dyf ( x，y)dxy (xo) yo在MATLAB中，编写被积函数f(x,y)的脚本，通过调用命令ode45(f,yO,xO,xf)可以得到较好的数值解。然而，难倒大多数初学者的情况却是：多变量微分方程，状态如何表示，函数的输入和输出该怎么用？单变量

5、二阶或高阶微分方程，只能求解一阶ODE的函数命令如何用在高阶系统上？。要求绘制非线性动力学系统的相轨迹图，该如何绘制？以下我就这两个问题进行解释1-多变量常微分方程组的求解多变量常微分方程组的求解，在MATLAB中仍然可 以按照函数ode模板进行求解，其形式与一阶单变量是一样的。由于MATLAB的ode 函数只能接受形参表形式为(t, x)的被积函数，所以我们要将状态变量通过多维矢量x传 递进去。如14所述的生物种群2变量微分方程组Valter Lotka，将方程中的状态变量X与y分别表示为一个二维向量的2维，从而按照类似方法进行求解。function dx=ValterraLotka(t,x

6、) A=2/3;B=4/3;C=1;D=1; dx=x(1)*(A-B*x(2); dx(2)=x(2)*(C*x(1)-D); dx=dx(:);end2. 高阶常微分方程如何表示？对于第2种形式的ODE，般来说我们可以引入一些附加变量而将其化为第1种形式，从而通过第1种形式的求解而得到高阶微分的结果。如Vander Pol方程是一个 二阶常微分方程，我们可以将其化为：在MATLAB中，同样可以编写多变量的被积函数VanderPol(t,x)，需要注意的是这里 的状态变量x是以上我们降阶之后的，即为矢量形式。同样函数输出的也是2变量的矢量形式，即.function dx=VanderPol(

7、t，x)%经过测试可得mu影响极限环的形状；初值不影响最终形状mu=0.8;dx(1)=x(2);dx(2)=mu*(1 -x(1)*x(1 )*x(2)-x;dx=dx(:);end3. 相图和极限环怎么绘制？对于该二阶单变量微分方程，控制工程中有极限环的概念。 它的定义是：在时间t趋于正无穷时，相轨迹都逐渐靠近某一条封闭的曲线，这种被称为 稳定的极限环。绘制极限环或者相图，首先得求解ODE方程组，得到状态的时域响应。 在MATLAB中，极限环的绘制方法很简单，因为以上降阶的过程中有：所以对微分方程 组的两个状态变量x1和x2是函数输出的结果，对其直接作plot即可得到相平面上状 态变量的运动轨迹，即极限环。t,y=ode45(VanderPol,0,15,-0.5,0); figure(1 ),plot(t,y(:,1 ),t,y(:,2),xlabel(f ),ylabel(刈; figure(2),plot(y(:,1),y(:,2),axis equal, xlabel(x),ylabel(dx/