1、M 文件和 Simulink 求解连续微分系统实例分析Matlab 自带的一个 S 函数源代码： D:MATLAB7toolboxsimulinkblockssfuntmpl例 1. 常微分方程（Lorenze 混沌系统）：12323,xbxac（1）其中 。0,28,/3arb（1） m 文件实现：文件名为 exam1.mfunction exam1x0=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)
2、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);0 10 20 30 40 50 60 70 80 90 100-30-20-10010203040500 10 2030 40 50-20-1001020-30-20-100102030(2) (I)Simulink 模块实现：（见 lorenzblok）x1-8/3x1x2x328x2x1*x2x2x23Out22x21x1ScopeProduct1Product1sIntegrat
3、or21sIntegrator11sIntegrator28Gain4-8/3Gain310Gain2-10Gain1其中三个积分模块的初始值设置与 exam1 相同，仿真时长为 100s。精度设置：Simulation-Configuration ParametersRelative tolerance, 1e-3 改为 1e-5（试试不作此修改的结果比较） 。运行后双击示波器 scope 后可看到：或在 matlab 命令窗口输入画图命令： plot(tout,yout) plot3(yout(:,2),yout(:,3),yout(:,1)60 65 70 75 80 85 90 95
4、100-30-20-1001020304050-20 -10 010 20-40-200204001020304050(II)Simulink 向量模块实现：(见 lorenzevector)tTo Workspace1xTo WorkspaceScope1sIntegrator-8*u(1)/3+u(2)*u(3)Fcn2-10*u(2)+10*u(3)Fcn1-u(1)*u(2)+28*u(2)-u(3)FcnClock画图语句：plot(t,x) plot3(x(:,1),x(:,2),x(:,3)(3)Simulink 中 S 函数的实现：（见 lorenzsfun 和 lorenzs
5、ystem.m）tTo Workspace1xTo WorkspaceScopelorenzsystemS-FunctionClock例 2. 常时滞微分方程:1223(),(0.),ytytt（2）（1） m 文件需调用 dde23 来求解：（见 exam2.m）function exam2sol = 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 =
6、 Z(:,2);v = zeros(3,1);v(1) = ylag1(1);v(2) = ylag1(1) + ylag2(2);v(3) = y(2);0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5020406080100120140160180200 Example 2time ty(t)(2)Simulink 中 S 函数来实现：（见 exam2sfun 和 exam2mfun.m）注：用 Simulink 中 S 函数求解时滞微分方程的核心思想在于：将时滞变量作为S 函数的外部输入。TransportDelay1TransportDelaytTo Workspace1
7、yTo WorkspaceScope2Scope1exam2mfunS-FunctionClock画图语句为： plot(t,y)例 3. 用 Simulink 解决一个变时滞的例子。 （见 logisticvariable 和logisticconstant）单时滞量的 Logistic 一阶时滞微分方程：()()1xtxtrK（3）()()xtxtr（4）取 和 。对于（3）和（4）而言，若分别取 和 ，1r2Kep()1.5)xtt.易知 当 充分大时。Simulink 模块框图分别如下：()tt3xVariableTransport DelayxTo Workspace1tTo WorkspaceScopeProduct1sIntegrator1Gain1-0.5Gainf(u)Fcn1 ConstantClock1Clock3xTransportDelayxTo Workspace1tTo WorkspaceScopeProduct1sIntegrator1Gain1-0.5Gain1 ConstantClock1取相同的初始条件，运行后分别双击 Scope 可得或画图命令：plot(t,x)思考：a. 改变时滞量的值，看有什么变化;b. 试试看用 S 函数替代上面的框图会更简单。
