2022年中国矿业大学《控制系统计算机仿真》实验试题及仿真程序及结果 .pdf

《2022年中国矿业大学《控制系统计算机仿真》实验试题及仿真程序及结果 .pdf》由会员分享，可在线阅读，更多相关《2022年中国矿业大学《控制系统计算机仿真》实验试题及仿真程序及结果 .pdf（34页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、实验一1、（教材 P37 页例2.1 ）已知一阶微分方程：试编写程序，用欧拉法求，即的值；并画出图形。解：程序如下：clcclearTf=input( 仿真时间 Tf=); % 输入仿真时间h=input( 计算步长 h= ); %输入仿真步长x0=1/3; t=0; x=x0; % 输入初始值for i=1:Tf/h K1=-30*x0; x1=x0+h*K1; x=x;x1; % 变量以向量形式保存 t=t;t(i)+h; % 对应时刻以向量形式保存 x0=x1;endt,x % 以数据形式输出plot(t,x) % 以曲线形式输出所绘图形如下：00.20.40.60.811.21.41.

2、6-0.2-0.100.10.20.30.40.5X: 1.5Y: 3.104e-0102、2. （教材 P79习题 2.5）已知系统的状态方程和输出方程为：式中 u(t)=1(t)。初始条件为：x1(0)=x2(0)=x3(0)=0。取h=0.05，试用RK4法求t=0.5时的y 值。解：程序如下：05.0,3/1)0(,30hyyy)5.1(ty30y321321321001104000121019018xxxyuxxxxxx精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 1 页，共 34 页clcclearu=1;x0=0; t=0;y=0;

3、x0=0;0;0;A=-8 1 0;-19 0 1;-12 0 0;B=0;4;10;C=1 0 0;Tf=input( 仿真时间 Tf=);h=input( 计算步长 h= );x=x0;for i=1:Tf/h K1=A*x+B*u; K2=A*(x+h*K1/2)+B*u; K3=A*(x+h*K2/2)+B*u; K4=A*(x+h*K3)+B*u; x=x+h*(K1+2*K2+2*K3+K4)/6; y=y;C*x; t=t;t(i)+h;endyplot(t,y) 所绘图形如下：0510152025303500.10.20.30.40.50.60.70.80.900.05 0.1

4、0.15 0.20.25 0.30.35 0.40.45 0.500.050.10.150.20.25X: 0.5Y: 0.2201精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 2 页，共 34 页3、已知系统的传递函数，试求其状态方程及输出方程。（教材 P26页 习题1.7 ）解：程序如下：den=1 3;num=1 3 3 1;A,B,C,D=tf2ss(den,num) 运行结果如下：A = -3 -3 -1 1 0 0 0 1 0 B = 1 0 0 C = 0 1 3 D = 0 4、已知线性定常系统的状态方程及输出方程，试求系统的传递

5、函数。（教材 P26页 习题1.8 ）解：程序如下：A=-1 1 0;0 -1 0;0 0 -2;B=0;0;1;C=1 0 1;den,num=ss2tf(A,B,C,D) 运行结果如下：1333)(23sssssG321321321 101 1002-00010011xxxyuxxxxxx精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 3 页，共 34 页den = 0 1 2 1 num = 1 4 5 2 5、 （教材P48页 例2.4 ）如图所示的直流电机拖动系统，试在SIMULINK 环境下建立系统仿真模型：1）采用欧拉法与四阶龙格库塔

6、法进行仿真，并确定其步长的稳定区间，取不同步长比较仿真结果；2）采用RK4研究外环 PI控制器参数的改变 (0.17 0.5 1 1.5)，观其阶跃响应对系统动态参数的影响。（1）、欧拉法clcclearTf=10; % 输入仿真时间h=input( 计算步长 h= ); %输入仿真步长P=1 0.01 0.1 0; 0 0.85 1 0.17; 1 0.01 1 0; 0 0.051 1 0.15; 1 0.0067 70 0; 1 0.15 0.21 0; 0 1 130 0; 1 0.01 0.1 0; 1 0.01 0.0044 0; % 输入各环节参数W0=zeros(9,1);W0

7、(1,1)=1; % 输入外部链接矩阵W=zeros(9,9);W(2,1)=1;W(2,9)=-1;W(3,2)=1;W(4,3)=1;W(4,8)=-1;W(5,4)=1;W(6,5)=1;W(6,7)=-0.212;W(7,6)=1;W(8,6)=1;W(9,7)=1; % 输入系统连接矩阵精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 4 页，共 34 页A=diag(P(:, 1);B=diag(P(:, 2);C=diag(P(:, 3);D=diag(P(:, 4); % 生成系数矩阵Q=B-D*W;Q1=inv(Q);R=C*W-A

8、; V=C*W0;A1=Q1*R; B1=Q1*V; % 生成闭环系数矩阵C1=0 0 0 0 0 0 1 0 0;u=1; %阶跃输入幅值x0=zeros(9,1); t=0;y=0; % 设置初值x=x0;for i=1:Tf/h K1=A1*x+B1*u; x=x+h*K1; y=y;C1*x; t=t;t(i)+h;endplot(t,y)计算步长 h=0.011 计算步长 h=0.012 0123456789100510152025303540012345678910-2-1.5-1-0.500.511.522.5x 1050四阶龙格库塔法：clcclearTf=input( 仿真时

9、间 Tf=); % 输入仿真时间h=input( 计算步长 h= ); %输入仿真步长P=1 0.01 0.1 0; 0 0.85 1 0.17; 1 0.01 1 0; 0 0.051 1 0.15; 1 0.0067 70 0; 1 0.15 0.21 0; 0 1 130 0;精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 5 页，共 34 页 1 0.01 0.1 0; 1 0.01 0.0044 0; % 输入各环节参数W0=zeros(9,1);W0(1,1)=1; % 输入外部链接矩阵W=zeros(9,9);W(2,1)=1;W(2

10、,9)=-1;W(3,2)=1;W(4,3)=1;W(4,8)=-1;W(5,4)=1;W(6,5)=1;W(6,7)=-0.212;W(7,6)=1;W(8,6)=1;W(9,7)=1; % 输入系统连接矩阵A=diag(P(:, 1);B=diag(P(:, 2);C=diag(P(:, 3);D=diag(P(:, 4); % 生成系数矩阵Q=B-D*W;Q1=inv(Q);R=C*W-A; V=C*W0;A1=Q1*R; B1=Q1*V; % 生成闭环系数矩阵C1=0 0 0 0 0 0 1 0 0;u=1; %阶跃输入幅值x0=zeros(9,1); t=0;y=0; % 设置初值x

11、=x0;for i=1:Tf/h K1=A1*x+B1*u; K2=A1*(x+h*K1/2)+B1*u; K3=A1*(x+h*K2/2)+B1*u; K4=A1*(x+h*K3)+B1*u; x=x+h*(K1+2*K2+2*K3+K4)/6; y=y;C1*x; t=t;t(i)+h;endplot(t,y)仿真时间 Tf=20 仿真时间Tf=5 仿真时间Tf=100 计算步长 h=0.015 计算步长 h=0.015 计算步长h=0.016 02468101214161820051015202530354000.5 11.5 22.5 33.5 44.5 505101520253035

12、40020406080100120-18-16-14-12-10-8-6-4-202x 10302精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 6 页，共 34 页（2）、clcclearTf=input( 仿真时间 Tf=); % 输入仿真时间h=input( 计算步长 h= ); %输入仿真步长P=1 0.01 0.1 0; 0.17 0.5 1 1.5; 1 0.01 1 0; 0 0.051 1 0.15; 1 0.0067 70 0; 1 0.15 0.21 0; 0 1 130 0; 1 0.01 0.1 0; 1 0.01 0.00

13、44 0; % 输入各环节参数W0=zeros(9,1);W0(1,1)=1; % 输入外部链接矩阵W=zeros(9,9);W(2,1)=1;W(2,9)=-1;W(3,2)=1;W(4,3)=1;W(4,8)=-1;W(5,4)=1;W(6,5)=1;W(6,7)=-0.212;W(7,6)=1;W(8,6)=1;W(9,7)=1; % 输入系统连接矩阵A=diag(P(:, 1);B=diag(P(:, 2);C=diag(P(:, 3);D=diag(P(:, 4); % 生成系数矩阵Q=B-D*W;Q1=inv(Q);R=C*W-A; V=C*W0;A1=Q1*R; B1=Q1*V;

14、 % 生成闭环系数矩阵C1=0 0 0 0 0 0 1 0 0;u=1; %阶跃输入幅值x0=zeros(9,1); t=0;y=0; % 设置初值x=x0;for i=1:Tf/h K1=A1*x+B1*u; K2=A1*(x+h*K1/2)+B1*u; K3=A1*(x+h*K2/2)+B1*u; K4=A1*(x+h*K3)+B1*u; x=x+h*(K1+2*K2+2*K3+K4)/6; y=y;C1*x; t=t;t(i)+h;endplot(t,y)精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 7 页，共 34 页仿真时间 Tf=20

15、 计算步长h=0.015 02468101214161820051015202530仿真时间 Tf=5 计算步长h=0.015 00.511.522.533.544.55051015202530精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 8 页，共 34 页SIMULINK仿真0.00440.01s+10.10.01s+1130s0.210.15s+1700.00167s+10.15s+10.051s10.01s+10.17s+10.85s0.10.01s+1y1To WorkspaceSaturation10.212精选学习资料 - - -

16、- - - - - - 名师归纳总结 - - - - - - -第 9 页，共 34 页0.00440.01s+10.10.01s+1130s0.210.15s+1700.00167s+10.15s+10.051s10.01s+10.17s+10.85s0.10.01s+1y1To Workspace0.212精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 10 页，共 34 页实验二1、已知某四阶非线性系统结构框图如图1 所示。当系统输入幅值为10 的阶跃信号时，试仿真分析系统输出响应：1）在 SIMULINK环境下建立仿真模型，输出仿真结果；2

17、）采用 MATLAB 语言编程实现仿真程序，与1)进行结果比较；3）在第一环节前加饱和非线性环节，分析限幅C1=6、2 时系统的输出响应，并与没限幅前比较其超调、峰值时间及调节时间。解： （1）SIMULINK图形如下1ss+0.785s+0.1610.1s+110.5s+11ss+0.785s+0.1610.1s+110.5s+110.1s+110.5s+11ss+0.785s+0.16Saturation1Saturation10.1s+110.5s+11ss+0.785s+0.16精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 11 页，共

18、34 页（2）子程序： satura.mfunction uc=satura(ur,s1)%ur=1;s1=5;if (ur=s1) uc=s1;elseif (ur=s1)if (ur0) uc=ur-s1;else uc=ur+s1;endelse uc=0;End 主程序：clcclearTf=20; % 输入仿真时间h=0.01; %输入仿真步长z=input( 非线性标志向量z= ); %输入非线性标志向量c1=input( 饱和非线性环节限幅c1=); % 输入饱和非线性环节参数a1=6;% 输入死区非线性环节参数精选学习资料 - - - - - - - - - 名师归纳总结 -

19、- - - - - -第 12 页，共 34 页P=0.16 1 0.785 1; 0 1 1 0; 1 0.5 1 0; 1 0.1 1 0; % 输入各环节参数A=P(:,1);B=P(:,2);C=P(:,3);D=P(:,4);W0=zeros(4,1);W0(1,1)=1; % 输入外部链接矩阵W=zeros(4,4);W(1,4)=-1;W(2,1)=1;W(3,2)=1;W(4,3)=1;for i=1:4if(A(i)=0); F(i)=1; G(i)=h*C(i)/B(i); H(i)=0.5*h*h*C(i)/B(i);if (D(i)=0) C1(i)=1;D1(i)=0

20、; else C1(i)=1;D1(i)=D(i)/B(i); % 求积分、比例积分环节离散系数endelse F(i)=exp(-h*A(i)/B(i);if (D(i)=0) G(i)=(1- F(i)*C(i)/A(i); H(i)=h*C(i)/A(i)-G(i)*B(i)/A(i); C1(i)=1;D1(i)=0; else G(i)=(1- F(i)*D(i)/A(i); H(i)=h*D(i)/A(i)-G(i)*B(i)/A(i); C1(i)=C(i)/D(i)-A(i)/B(i); D1(i)=D(i)/B(i); % 求惯性、比例惯性环节离散系数endendend%F1

21、=diag(F(:, 1);G1=diag(G(:, 1);H1=diag(H(:, 1);%C2=diag(C1(:, 1);D2=diag(D1(:, 1);% 求各环节输入Y=zeros(4,1);X=Y;y=zeros(1,4);r=10;Uk=zeros(4,1);t=0:h:Tf; N=length(t); % 求采样点个数for k=1:N-1 Ub=Uk; %保存前一次输入 Uk=W*Y+W0*r; % 求当前个环节输入v for i=1:4;if(z(i)=0)if(z(i)=1) Uk(i,1)=satura(Uk(i,1),c1);end精选学习资料 - - - - -

22、- - - - 名师归纳总结 - - - - - - -第 13 页，共 34 页if(z(i)=2) Uk(i,1)=dead(Uk(i,1),a1);endendend Udot=(Uk-Ub)/h; % 求当前各环节输入导数 Uf=2*Uk-Ub; %求下一拍输入%求各环节状态变量及输出 X=F.*X+G.*Uk+H.*Udot; Y=C1.*X+D1.*Uf;for i=1:4;if(z(i)=0)if(z(i)=3) Y(i,1)=satura(Y(i,1),c1);endif(z(i)=4) Y(i,1)=dead(Y(i,1),a1);endendend y=y;Y;endy4=

23、y(:,4);plot(t,y4)非线性标志向量z=0 0 0 0 饱和非线性环节限幅c1=6 024681012141618200246810121416X: 2.74Y: 15.17精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 14 页，共 34 页（3）非线性标志向量z=1 0 0 0 饱和非线性环节限幅c1=6 02468101214161820051015X: 3.12Y: 14.7非线性标志向量z=1 0 0 0 饱和非线性环节限幅c1=2 0246810121416182002468101214X: 4.91Y: 13.182、 某

24、控制系统如图所示，选择增益K 的值，使系统阶跃响应的超调整量小于20%，且调节时间小于5 s 。 （教材 P80 习题 2.11）K(s+1)(s+1)(s+8)(s+20)Zero-Pole1s +3.2s +3.56s32Transfer FcnStepScope精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 15 页，共 34 页解： SIMULINK图形如下：1550(s+1)(s+1)(s+8)(s+20)Zero-Pole1s +3.2s +3.56s32Transfer FcnStepScope670(s+1)(s+1)(s +8)(

25、s+20)Zero-Pole1s +3.2s +3.56s32Transfer FcnStepScope精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 16 页，共 34 页3、 某非线性系统如图所示， 试求时系统的动态特性。 (P52页例题 2.7) 解： SIMULINK图形如下：2s+1s +2s +s32Transfer FcnStepScopeSaturation1Gain)( 12)(ttr2s +1s +2s +s32Transfer FcnStepScopeSaturation5Gain2s +1s +2s +s32Transfer

26、 FcnStepScopeMATLABFunctionMATLAB Fcn5Gain精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 17 页，共 34 页2s +1s +2s +s32T ransfer Fcn12s +1s +2s +s32T ransfer FcnStep1StepScopeSaturation5Gain15GainuyfcnEmbeddedMATLAB Function其中：function y = fcn(u)%#emlif u=1 y=1;elseif u=s1)if (ur0) uc=ur-s1;else uc=ur+s

27、1;endelse uc=0;end6. clcclear精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 27 页，共 34 页Tf=input( 仿真时间 Tf=); % 输入仿真时间h=input( 采样周期 h= ); %输入仿真步长z=input( 非线性标志向量z= ); %输入非线性标志向量c1=input( 饱和非线性环节限幅c1=); % 输入饱和非线性环节参数a1=input( 死区非线性环节参数a1= ); % 输入死区非线性环节参数P=0.16 1 0.785 1; 0 1 1 0; 1 0.5 1 0; 1 0.1 1 0;

28、 % 输入各环节参数A=P(:,1);B=P(:,2);C=P(:,3);D=P(:,4);W0=zeros(4,1);W0(1,1)=1; % 输入外部链接矩阵W=zeros(4,4);W(1,4)=-1;W(2,1)=1;W(3,2)=1;W(4,3)=1;for i=1:4if(A(i)=0); F(i)=1; G(i)=h*C(i)/B(i); H(i)=0.5*h*h*C(i)/B(i);if (D(i)=0) C1(i)=1;D1(i)=0; else C1(i)=1;D1(i)=D(i)/B(i); % 求积分、比例积分环节离散系数endelse F(i)=exp(-h*A(i)

29、/B(i);if (D(i)=0) G(i)=(1- F(i)*C(i)/A(i); H(i)=h*C(i)/A(i)-G(i)*B(i)/A(i); C1(i)=1;D1(i)=0; else G(i)=(1- F(i)*D(i)/A(i); H(i)=h*D(i)/A(i)-G(i)*B(i)/A(i); C1(i)=C(i)/D(i)-A(i)/B(i); D1(i)=D(i)/B(i); % 求惯性、比例惯性环节离散系数endendend%F1=diag(F(:, 1);G1=diag(G(:, 1);H1=diag(H(:, 1);%C2=diag(C1(:, 1);D2=diag(

30、D1(:, 1);% 求各环节输入Y=zeros(4,1);X=Y;y=zeros(1,4);r=10;Uk=zeros(4,1);t=0:h:Tf; N=length(t); % 求采样点个数for k=1:N-1 Ub=Uk; %保存前一次输入 Uk=W*Y+W0*r; % 求当前个环节输入v 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 28 页，共 34 页for i=1:4;if(z(i)=0)if(z(i)=1) Uk(i,1)=satura(Uk(i,1),c1);endif(z(i)=2) Uk(i,1)=dead(Uk(i,1)

31、,a1);endendend Udot=(Uk-Ub)/h; % 求当前各环节输入导数 Uf=2*Uk-Ub; %求下一拍输入%求各环节状态变量及输出 X=F.*X+G.*Uk+H.*Udot; Y=C1.*X+D1.*Uf;for i=1:4;if(z(i)=0)if(z(i)=3) Y(i,1)=satura(Y(i,1),c1);endif(z(i)=4) Y(i,1)=dead(Y(i,1),a1);endendend y=y;Y;endy4=y(:,4);plot(t,y4)7. clcclearnum=input(num=); %输入分子向量den=input(den=); % 输

32、入分母向量T=input(T=); % 输入采样周期a,b,c,d=tf2ss(num,den); % 转化为状态空间表达式f,g,c1,d1=c2dm(a,b,c,d,T,zoh) %求离散状态方程系数8. function uc,urb=hysteresis(urb,ur,ucb,s1,h1)%urb 前一次输入， ur 当前输入， ucb 前一次输出，s1 滞环环宽精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 29 页，共 34 页if (ururb) % 当输入增大if (urs1) %当输入大于环宽 uc=h1; % 输出为 h1els

33、e uc=ucb; % 输入没有大于环宽，保持前一次输出endelseif (uurb) %当输入减小if (ur=s1) uc=s1;elseif (ur=-1*s1) uc=-1*s1;else uc=ur;end11.clcclearTf=input( 仿真时间 Tf=); % 输入仿真时间h=input( 计算步长 h= ); %输入仿真步长a=input(ASR 参数 a= ); % 输入转速调节器参数P=1 0.01 0.1 0; 0 0.85 1 a;精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 31 页，共 34 页 1 0.01

34、 1 0; 0 0.051 1 0.15; 12.clcclearu=1;x0=0; t=0;y=0;x0=0;0;0;A=0 1 0;0 0 1;-22.06 -27 -10;B=0;0;1;C=40.6 0 0;Tf=input( 仿真时间 Tf=);h=input( 计算步长 h= );x=x0;for i=1:Tf/h K1=A*x+B*u; K2=A*(x+h*K1/2)+B*u; K3=A*(x+h*K2/2)+B*u; K4=A*(x+h*K3)+B*u; x=x+h*(K1+2*K2+2*K3+K4)/6; y=y;C*x; t=t;t(i)+h;endplot(t,y) 1

35、0.0067 70 0; 1 0.15 0.21 0; 0 1 130 0; 1 0.01 0.1 0; 1 0.01 0.0044 0; % 输入各环节参数W0=zeros(9,1);W0(1,1)=1; % 输入外部链接矩阵W=zeros(9,9);W(2,1)=1;W(2,9)=-1;W(3,2)=1;W(4,3)=1;W(4,8)=-1;W(5,4)=1;W(6,5)=1;W(6,7)=-0.212;W(7,6)=1;W(8,6)=1;W(9,7)=1; % 输入系统连接矩阵A=diag(P(:, 1);B=diag(P(:, 2);C=diag(P(:, 3);D=diag(P(:,

36、 4); % 生成系数矩阵Q=B-D*W;Q1=inv(Q);R=C*W-A; V=C*W0;A1=Q1*R; B1=Q1*V; % 生成闭环系数矩阵C1=0 0 0 0 0 0 1 0 0;u=1; %阶跃输入幅值x0=zeros(9,1); t=0;y=0; % 设置初值x=x0;for i=1:Tf/h K1=A1*x+B1*u;精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 32 页，共 34 页 K2=A1*(x+h*K1/2)+B1*u; K3=A1*(x+h*K2/2)+B1*u; K4=A1*(x+h*K3)+B1*u; x=x+h

37、*(K1+2*K2+2*K3+K4)/6; y=y;C1*x; t=t;t(i)+h;endplot(t,y) 13.clcclearu=1;x0=0 0 0;C=40.6 0 0;tspan=0,10;t,x=ode45(funexam2,tspan,x0);y=x*C;plot(t,y);14.clcclearTf=input( 仿真时间 Tf=); % 输入仿真时间h=input( 计算步长 h= ); %输入仿真步长x0=1/3; t=0; x=x0; % 输入初始值for i=1:Tf/h K1=-30*x0; x1=x0+h*K1; x=x;x1; % 变量以向量形式保存 t=t;

38、t(i)+h; % 对应时刻以向量形式保存 x0=x1;endt,x % 以数据形式输出plot(t,x) % 以曲线形式输出15.clearx0=620,10,70; % 置状态变量初值tspan=0,30; % 置仿真时间区间t,x=ode45(myf,tspan,x0); % 调用 ode45 求仿真解plot(t,x(:,1),t,x(:,2),-,t,x(:,3),:); % 用不同的线型绘制仿真结果曲线xlabel(time(天) t0=0, tf=30); % 对t-x轴进行标注ylabel(x(人):x1(0)=620,x2(0)=10;x3(0)=70);精选学习资料 - -

39、 - - - - - - - 名师归纳总结 - - - - - - -第 33 页，共 34 页legend(x1, x2, x3);grid; 16. function dx=myf(t,x)dx = zeros(3,1);dx(1)=-0.001*x(1)*x(2); % 第一个微分方程dx(2)=0.001*x(1)*x(2)-0.072*x(2); % 第二个微分方程dx(3)=0.072*x(2);17.clear allnum=0.005,0.005,0; % 脉冲传递函数分子多项式按z 的降幂系数% 排列的行向量den=1,-2.4,1.863,-0.453; % 脉冲传递函数分

40、母多项式按z 的降幂系数% 排列的行向量yk,x,n=dstep(num,den); % yk 为存放输出离散序列的数组,n 为dstep% 函数自动设定的采样点数T=0.1; % 已知系统采样周期为0.1sfor k=1:n plot(k*T,yk(k),*k); % k 为采样序列号,k*T为第 k 次采样对应的时刻 hold onend xlabel( 时间 (s); grid hold off18. function dy = vdp1000(t,y)dy = zeros(2,1); % a column vectordy(1) = y(2);dy(2) = 1000*(1 - y(1)2)*y(2) - y(1);精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 34 页，共 34 页