M文件和Simulink求解连续微分系统实例分析

Matlab自带的一个S函数源代码： D:\\MATLAB7\oolbox\\simulink\\blocks\\sfuntmpl

例1.常微分方程（Lorenze混沌系统）：

                                                    （1）

其中。

（1）m文件实现：文件名为exam1.m

function exam1

x0=[0;0;1e-3];

[t,x]=ode45(@lorenzfun,[0,100],x0);

figure(1)

plot(t,x)

figure(2)

plot3(x(:,1),x(:,2),x(:,3))

%----------------------------------

function dx=lorenzfun(t,x)

a=10;c=28;b=8/3;

dx=zeros(3,1);

dx(1)=-b*x(1)+x(2)*x(3);

dx(2)=-a*x(2)+10*x(3);

dx(3)=-x(1)*x(2)+c*x(2)-x(3);

(2) (I)Simulink模块实现：（见lorenzblok）

其中三个积分模块的初始值设置与exam1相同，仿真时长为100s。精度设置：Simulation--Configuration Parameters—Relative tolerance, 1e-3改为1e-5（试试不作此修改的结果比较）。运行后双击示波器scope后可看到：

或在matlab命令窗口输入画图命令：

>> plot(tout,yout) >>plot3(yout(:,2),yout(:,3),yout(:,1))

(II)Simulink向量模块实现：(见lorenzevector)

画图语句：>>plot(t,x) >>plot3(x(:,1),x(:,2),x(:,3))

(3)Simulink中S函数的实现：（见lorenzsfun 和lorenzsystem.m）

例2.常时滞微分方程:

                                            （2）

（1）m文件需调用dde23来求解：（见exam2.m）

function exam2

sol = dde23('exam1f',[1, 0.2],ones(3,1),[0, 5]);

plot(sol.x,sol.y);

title('Example 2')

xlabel('time t');

ylabel('y(t)');

%-----------------------------------------------

function v = exam1f(t,y,Z)

ylag1 = Z(:,1);

ylag2 = Z(:,2);

v = zeros(3,1);

v(1) = ylag1(1);

v(2) = ylag1(1) + ylag2(2);

v(3) = y(2);

(2)Simulink中S函数来实现：（见exam2sfun和exam2mfun.m）

注：用Simulink中S函数求解时滞微分方程的核心思想在于：将时滞变量作为S函数的外部输入。

画图语句为： >>plot(t,y)

例3.用Simulink解决一个变时滞的例子。（见logisticvariable和logisticconstant）

单时滞量的Logistic一阶时滞微分方程：

                                           （3）

                                             （4）

取和。对于（3）和（4）而言，若分别取和，易知当充分大时。Simulink模块框图分别如下：

取相同的初始条件，运行后分别双击Scope可得

或画图命令：>>plot(t,x)

思考：a. 改变时滞量的值，看有什么变化;

      b. 试试看用S函数替代上面的框图会更简单。