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

M文件和Simulink求解连续微分系统实例分析Matlab 自带的一个 S 函数源代码：D:MATLAB7toolboxsimulinkblockssfuntmpl 例1.常微分方程Lorenze混沌系统：*3 -人 X2 ex 一 , x其中 a =10, r =28,b =8/ 3。(1) m文件实现：文件名为 examl.mfun cti on examl x0=0;0;1e-3;t,x=ode45(lore nzfu n,0,100,x0); figure(1)plot(t,x)figure(2)plot3(x(:,1),x(:,2),x(:,3)% fun cti on dx=lore nzfun (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);-10 20 103010 2030405060708090100-20 0(2) (l) Simulink 模块实现:(见 lorenzblok )A T ，dinteqratar>2integratorGain?Prozuct'其中三个积分模块的初始值设置与examl相同，仿真时长为100s。精度设置：Simulation-Configuration Parameter Relative toleranee, 1e-3 改为 1e-5 试试 不 作此修改的结果比拟。运行后双击示波器scope后可看到：或在matlab命令窗口输入画图命令:>> plot(tout,yout) >>plot3(yout(:,2),yout(:,3),yout(:,1)5040302010-10-20'i '- j i JI Fi J2iif50403020100.40 '20's.201020-0 f 10-40-2060657075808590951o0ClockTo Workspace1(ll)Simulink 向量模块实现：(见 lorenzevector)FenScopeTo WorkspaceFcn1Fcn2ClockTo Workspace1画图语句：>>plot(t,x)>>plot3(x(:,1),x(:,2),x(:,3)(3)Simulink 中 S 函数的实现：(见 lorenzsfun 和 lorenzsystem.m)ScopeorenzsystemS-Function To Workspace3例2.常时滞微分方程：yi =yi(t 1), y2 = yi(t-1) * y2 (t Q2) y3 = y2 ，(1) m文件需调用dde23来求解：(见exam2.m ) fun cti on 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)');%fun cti on 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);ClockClock1To Workspace180160140120100806040200.51.522.5time t3.54.5200 <Example 2(2)Simulink 中 S 函数来实现：(见 exam2sfun 和 exam2mfun.m )将时滞变量作为S注：用Simulink中S函数求解时滞微分方程的核心思想在于:函数的外部输入Transporte 悟y T 倒 I TransportDelay 1;|ClockTo WorkspaceTo Workspacel画图语句为：>>plot(t,y)例3.用Simulink解决一个变时滞的例子。(见 logisticvariable 和 logisticconstant )单时滞量的Logistic 一阶时滞微分方程x'(t 卜 rx(tb1<(tK-t)x'(t)取r =1和K =2。对于(3)和(4)而言，假设分别取？ (t)二exP( t) 1.5和.=1.5, exp( 1)+1易知(t)?-当t充分大时。Simulink模块框图分别如下:Scopef(u)FcnVariableTransport Delay3xTo Workspace1Scopel_l<rIntegratorGainlTranlpor t DelayClockl取相同的初始条件，运行后分别双击3xio vvorkspacel-0.GainConstantProductTo WorkspaceScope可得或画图命令：>>plot(t,x)思考：a.改变时滞量的值，看有什么变化；b.试试看用S函数替代上面的框图会更简单