系统仿真实验报告

1、系统仿真技术实验报告 信息科学与工程学院 目录实验一 MATLAB中矩阵与多项式的基本运算1实验二 MATLAB绘图命令2实验三 MATLAB程序设计3实验四 MATLAB的符号计算与SIMULINK的使用5实验五 MATLAB在控制系统分析中的应用8实验六 连续系统数字仿真的基本算法15实验一对不同命令操作训练1.eye(3)ans = 1 0 0 0 1 0 0 0 12.ones(2)ans = 1 1 1 1ones(3,2)ans = 1 1 1 1 1 1 3.zeros(3,2)ans = 0 0 0 0 0 04.rand(3,2)ans = 0.8147 0.9134 0.9

2、058 0.6324 0.1270 0.09755.diag(3) ans = 36. A=ones(3,3);B=rand(3,3); ans=Bans = 0.0357 0.6787 0.3922 0.8491 0.7577 0.6555 0.9340 0.7431 0.17127. A/Bans = 0.1195 1.5417 -0.3355 0.1195 1.5417 -0.33550.1195 1.5417 -0.3355 AB警告: 矩阵为奇异工作精度。 ans = NaN NaN NaN NaN NaN NaN Inf Inf -Inf inv(A)*B、B*inv(A)同上8.

3、roots(1)ans = Empty matrix: 0-by-1A*Bans = 1.8188 2.1796 1.2189 1.8188 2.1796 1.2189 1.8188 2.1796 1.2189>> B*Aans = 1.1067 1.1067 1.1067 2.2623 2.2623 2.2623 1.8483 1.8483 1.84839.who您的变量为:A B ans th x whos Name Size Bytes Class Attributes A 3x3 72 double B 3x3 72 double ans 3x3 72 double th

4、1x21 168 double x 1x21 336 double complex 10.size(A)ans = 3 3 length(A)ans = 3实验二1、二维直角坐标 t=-pi:0.3:pi;y=1./(1+exp(-t);subplot(221),plot(t,y,'b');xlabel('t');title('plot(t,y)');subplot(222),stem(t,y,'g');xlabel('t');title('plot(t,y)');subplot(223),semi

5、logy(t,y,'m');xlabel('t');title('plot(t,y)');subplot(224),stairs(t,y,'k');xlabel('t');title('plot(t,y)')实验三 function f=SY3()while 1 b=input('输入一个自然数:'); if b<2 fprintf('%d既不是质数也不是合数',b); elseif isprime(b)=1 fprintf('%d是一个素数',

6、b); else fprintf('%d是一个合数，可以因式分解为:',b) fprintf('n') for n=2:sqrt(b) if rem(b,n)=0 num1=n; num2=b/n; fprintf('%d = %d×%d ',b,num1,num2) fprintf('n') end end end fprintf('n') y=input('继续判断？继续按1，退出按0:'); if y=0 break; endend 命令窗口输入 SY3输入一个自然数:7878是一个

7、合数，可以因式分解为:78 = 2×39 78 = 3×26 78 = 6×13 继续判断？继续按1，退出按0:实验四1求矩阵对应的行列式和特征根a=sym('a11 a12;a21 a22');da=det(a)ea=eig(a)da = a11*a22 - a12*a21 ea = a11/2 + a22/2 - (a112 - 2*a11*a22 + a222 + 4*a12*a21)(1/2)/2 a11/2 + a22/2 + (a112 - 2*a11*a22 + a222 + 4*a12*a21)(1/2)/2 2 求方程的解（包括精

8、确解和一定精度的解）r1=solve('x2-x-1')rv=vpa(r1)rv4=vpa(r1,4)rv20=vpa(r1,20)r1 = 5(1/2)/2 + 1/2 1/2 - 5(1/2)/2 rv = 1.6180339887498948482045868343656 -0.61803398874989484820458683436564 rv4 = 1.618 -0.618 rv20 = 1.6180339887498948482 -0.61803398874989484823 a=sym('a');b=sym('b');c=sym(

9、'c');d=sym('d'); %定义4个符号变量w=10;x=5;y=-8;z=11; %定义4个数值变量A=a,b;c,d %建立符号矩阵AB=w,x;y,z %建立数值矩阵Bdet(A) %计算符号矩阵A的行列式det(B) %计算数值矩阵B的行列式A = a, b c, d B = 10 5 -8 11 ans = a*d - b*c ans = 150 4 syms x y;s=(-7*x2-8*y2)*(-x2+3*y2);expand(s) %对s展开collect(s,x) %对s按变量x合并同类项(无同类项)factor(ans) % 对an

10、s分解因式5 对方程 AX=b求解A=34,8,4;3,34,3;3,6,8;b=4;6;2;X=linsolve(A,b) %调用linsolve函数求解X = 138/2045 66/409 212/2045 Ab %用另一种方法求解ans = 138/2045 66/409 212/2045 6syms a b t x y z;f=sqrt(1+exp(x);diff(f) %未指定求导变量和阶数，按缺省规则处理f=x*cos(x);diff(f,x,2) %求f对x的二阶导数diff(f,x,3) %求f对x的三阶导数f1=a*cos(t);f2=b*sin(t);diff(f2)/d

11、iff(f1) %按参数方程求导公式求y对x的导数 ans = exp(x)/(2*(exp(x) + 1)(1/2) ans = - 2*sin(x) - x*cos(x) ans = x*sin(x) - 3*cos(x) ans = -(b*cos(t)/(a*sin(t) 结果： 实验五1、求下面系统的单位阶跃响应%程序如下：num=5 ; den=1 , 1 , 4 ;(分子4改为5)step(num , den)y , x , t=step(num , den) ;tp=spline(y , t , max(y) %计算峰值时间max(y) %计算峰值 tp = 1.6579>

12、;> max(y) %计算峰值ans =1.80402. 求如下系统的单位阶跃响应%程序如下：a=0,1;-6,-5;b=0;1;c=1,0;d=0;y,x=step(a,b,c,d);plot(y) 3. 求下面系统的单位脉冲响应：%程序如下：num=4 ; den=1 , 1 ,4 ;impulse(num,den) 4. 已知二阶系统的状态方程为：求系统的零输入响应和脉冲响应。 %程序如下：a=0 , 1 ; -10 , -2 ; b=0 ; 1 ; c=1 , 0 ; d=0 ; x0=1 ,0; subplot(1 , 2 , 1) ; initial(a , b , c ,d

13、,x0) subplot(1 , 2 , 2) ; impulse(a , b , c , d) 6. 有一二阶系统，求出周期为4秒的方波的输出响应 %程序如下： num=2 5 1;den=1 2 3;t=(0:.1:10);period=4;u=(rem(t,period)>=period./2);%看rem函数功能lsim(num,den,u,t);7. 已知开环系统传递函数，绘制系统的根轨迹，并分析其稳定性 %程序如下： num=1 2;den1=1 4 3;den=conv(den1,den1);figure(1)rlocus(num,den)k,p= rlocfind(num

14、,den) figure(2)k=55;num1=k*1 2;den=1 4 3;den1=conv(den,den);num,den=cloop(num1,den1,-1);impulse(num,den)title('impulse response (k=55) ')figure(3)k=56;num1=k*1 2;den=1 4 3;den1=conv(den,den);num,den=cloop(num1,den1,-1);impulse(num,den)title('impulse response(k=56)') 8. 作如下系统的bode图 %程

15、序如下： n=1 , 1 ; d=1 , 4 , 11 , 7 ; bode(n , d) 9. 系统传函如下 求有理传函的频率响应，然后在同一张图上绘出以四阶伯德近似表示的系统频率响应 %程序如下： num=1;den=conv(1 2,conv(1 2,1 2); w=logspace(-1,2); t=0.5;m1,p1=bode(num,den,2);p1=p1-t*w'*180/pi;n2,d2=pade(t,4);numt=conv(n2,num);dent=(conv(den,d2);m2,p2=bode(numt,dent,w);subplot(2,1,1);semil

16、ogx(w,20*log10(m1),w,20*log10(m2),'g-');grid on ; title('bode plot');xlabel('frequency');ylabel('gain');subplot(2,1,2);semilogx(w,p1,w,p2,'g-');grid on;xlabel('frequency');ylabel('phase'); 11. 二阶系统为：令wn=1,分别作出=2 , 1 , 0.707 , 0.5时的nyquist曲线。 %程

17、序如下：n=1 ; d1=1 , 4 , 1 ; d2=1 , 2 , 1 ; d3=1 , 1.414 , 1; d4=1,1,1; nyquist(n,d1) ; hold on nyquist(n,d2) ; nyquist(n,d3) ; nyquist(n,d4) ; 12. 已知系统的开环传递函数为 绘制系统的Nyqusit图，并讨论系统的稳定性 %程序如下： G=tf(1000,conv(1,3,2,1,5);nyquist(G);axis('square') 13. 分别由w的自动变量和人工变量作下列系统的nyquist曲线： %程序如下： n=1 ; d=1

18、, 1 ,0 ;nyquist(n ,d) ; %自动变量n=1 ; d=1 , 1 ,0 ; w=0.5 : 0.1 : 3 ;nyquist(n , d , w) ; %人工变量 14. 一多环系统，其结构图如下，使用Nyquist频率曲线判断系统的稳定性。 %程序如下： k1=16.7/0.0125;z1=0;p1=-1.25 -4 -16;num1,den1=zp2tf(z1,p1,k1);num,den=cloop(num1,den1);z,p,k=tf2zp(num,den);p figure(1)nyquist(num,den)figure(2)num2,den2=cloop(n

19、um,den);impulse(num2,den2); 15. 已知系统为：作该系统的nichols曲线。 %程序如下： n=1 ; d=1 , 1 , 0 ; ngrid(new) ; nichols(n , d) ; 16. 已知系统的开环传递函数为：当k=2时，分别作nichols曲线和波特图。 %程序如下：num=1;den=conv(conv(1 0,1 1),0.5 1);subplot(1,2,1);nichols(num,den);grid; % nichols曲线subplot(1,2,2);g=tf(num,den);bode(feedback(g,1,-1);grid;

20、%波特图17. 系统的开环传递函数为： 分别确定k=2和k=10时闭环系统的稳定性。%程序如下：d1=1 , 3 , 2 , 0 ; n1=2 ;nc1 , dc1=cloop(n1 , d1 ,-1) ;roots(dc1)d2=d1 ; n2=10 ;nc2 , dc2=cloop(n2 , d2,-1) ; roots(dc2)ans = -2.5214 + 0.0000i -0.2393 + 0.8579i -0.2393 - 0.8579ians = -3.3089 + 0.0000i 0.1545 + 1.7316i 0.1545 - 1.7316i18. 系统的状态方程为： 试确

21、定系统的稳定性。 %程序如下： a=-4,-3,0 ; 1,0,0 ; 0,1,0 ; b=1;0;0 ; c=0,1,2 ; d=0 ;eig(a) %求特征根rank(ctrb(a,b) ans = 0 -1 -3ans = 3实验六1. 取h=0.2，试分别用欧拉法、RK2法和RK4法求解微分方程的数值解，并比较计算精度。 注：解析解： %程序如下cleart(1)=0; % 置自变量初值y(1)=1; y_euler(1)=1; y_rk2(1)=1; y_rk4(1)=1; % 置解析解和数值解的初值h=0.2; % 步长% 求解析解for k=1:5 t(k+1)=t(k)+h;

22、y(k+1)=sqrt(1+2*t(k+1);end% 利用欧拉法求解for k=1:5 y_euler(k+1)=y_euler(k)+h*(y_euler(k)-2*t(k)/y_euler(k);end% 利用RK2法求解for k=1:5 k1=y_rk2(k)-2*t(k)/y_rk2(k); k2=(y_rk2(k)+h*k1)-2*(t(k)+h)/(y_rk2(k)+h*k1); y_rk2(k+1)=y_rk2(k)+h*(k1+k2)/2;end% 利用RK4法求解for k=1:5 k1=y_rk4(k)-2*t(k)/y_rk4(k); k2=(y_rk4(k)+h*k

23、1/2)-2*(t(k)+h/2)/(y_rk4(k)+h*k1/2); k3=(y_rk4(k)+h*k2/2)-2*(t(k)+h/2)/(y_rk4(k)+h*k2/2); k4=(y_rk4(k)+h*k3)-2*(t(k)+h)/(y_rk4(k)+h*k3); y_rk4(k+1)=y_rk4(k)+h*(k1+2*k2+2*k3+k4)/6;end% 输出结果disp(' 时间 解析解 欧拉法 RK2法 RK4法')yt=t', y', y_euler', y_rk2', y_rk4'disp(yt)程序运行结果如下： 时间

24、 解析解 欧拉法 RK2法 RK4法 0 1.0000 1.0000 1.0000 1.0000 0.2000 1.1832 1.2000 1.1867 1.1832 0.4000 1.3416 1.3733 1.3483 1.3417 0.6000 1.4832 1.5315 1.4937 1.4833 0.8000 1.6125 1.6811 1.6279 1.6125 1.0000 1.7321 1.8269 1.7542 1.73212. 考虑如下二阶系统： 在上的数字仿真解（已知：，），并将不同步长下的仿真结果与解析解进行精度比较。 说明：已知该微分方程的解析解分别为： 采用RK4法

25、进行计算，选择状态变量： 则有如下状态空间模型及初值条件 采用RK4法进行计算。程序如下：clearh=input('请输入步长h='); % 输入步长M=round(10/h); % 置总计算步数t(1)=0; % 置自变量初值y_0(1)=100; y_05(1)=100; % 置解析解的初始值(y_0和y_05分别对应于为R=0和R=0.5)x1(1)=100; x2(1)=0; % 置状态向量初值y_rk4_0(1)=x1(1); y_rk4_05(1)=x1(1); % 置数值解的初值 % 求解析解for k=1:M t(k+1)=t(k)+h; y_0(k+1)=1

26、00*cos(t(k+1); y_05(k+1)=100*exp(-t(k+1)/2).*cos(sqrt(3)/2*t(k+1)+100*sqrt(3)/3*exp(-t(k+1)/2).*sin(sqrt(3)/2*t(k+1);end% 利用RK4法求解% R=0for k=1:M k11=x2(k); k12=-x1(k); k21=x2(k)+h*k12/2; k22=-(x1(k)+h*k11/2); k31=x2(k)+h*k22/2; k32=-(x1(k)+h*k21/2); k41=x2(k)+h*k32; k42=-(x1(k)+h*k31); x1(k+1)=x1(k)+h*(k11+2*k21+2*k31+k41)/6; x2(k+1)=x2(k)+h*(k12+2*k22+2*k32+k42)/6; y_rk4_0(k+1)=x1(k+1);end% R=0.5for k=1:M k11=x2(k); k12=-x1(k)-x2(k); k21=x2(k)+h*k12/2; k22=-(x1(k)+h*k11/2)-(x2(k)+h*k12/2); k31=x2(k)+h*