中南大学系统仿真 实验报告.doc

实验一 MATLAB中矩阵与多项式的基本运算实验任务1了解MATLAB命令窗口和程序文件的调用。2熟悉如下MATLAB的基本运算： 矩阵的产生、数据的输入、相关元素的显示； 矩阵的加法、乘法、左除、右除； 特殊矩阵：单位矩阵、“1”矩阵、“0”矩阵、对角阵、随机矩阵的产生和运算； 多项式的运算：多项式求根、多项式之间的乘除。基本命令训练1、 eye(2)ans = 1 0 0 1 eye(4)ans = 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 12、 ones(2)ans = 1 1 1 1 ones(4)ans = 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ones(2,2)ans = 1 1 1 1 ones(2,3)ans = 1 1 1 1 1 1 ones(4,3)ans = 1 1 1 1 1 1 1 1 1 1 1 13、 zeros(2)ans = 0 0 0 0 zeros(4)ans = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 zeros(2,2)ans = 0 0 0 0 zeros(2,3)ans = 0 0 0 0 0 0 zeros(3,2)ans = 0 0 0 00 04、随机阵 rand(2,3)ans = 0.2785 0.9575 0.15760.5469 0.9649 0.9706 rand(2,3)ans =0.9572 0.8003 0.4218 0.4854 0.1419 0.91575、 diag(5)ans = 5 diag(5,5)ans = 0 0 0 0 0 5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 diag(2,3)ans = 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 06、（inv（A）为求A的逆矩阵） B=5 3 1;2 3 8;1 1 1,inv(B)B = 5 3 1 2 3 8 1 1 1ans = 0.6250 0.2500 -2.6250 -0.7500 -0.5000 4.7500 0.1250 0.2500 -1.1250 A=2 3;4 4,B=5 3;3 8,inv(A),inv(B);AB,A/B,inv(A)*B,B*inv(A)A = 2 3 4 4B = 5 3 3 8ans = -1.0000 0.7500 1.0000 -0.5000ans = -2.7500 3.0000 3.5000 -1.0000ans = 0.2258 0.2903 0.6452 0.2581ans = -2.7500 3.0000 3.5000 -1.0000ans = -2.0000 2.25005.0000 -1.75007、 p =1,-6,-72,-27, roots(p)p = 1 -6 -72 -27ans = 12.1229 -5.7345 -0.3884 p=2,3,6,roots(p)p = 2 3 6ans = -0.7500 + 1.5612i -0.7500 - 1.5612i8、（A为n*n的方阵） A=0 1 0;-4 4 0;-2 1 2,poly(A),B=sym(A),poly(B)A = 0 1 0 -4 4 0 -2 1 2ans = 1 -6 12 -8 B = 0, 1, 0 -4, 4, 0 -2, 1, 2 ans = x3-6*x2+12*x-89,、（conv是多项式相乘，deconv是多项式相除） u=1 2 4 6 ,v=5 0 0 -6 7,conv(u,v)u = 1 2 4 6v = 5 0 0 -6 7ans = 5 10 20 24 -5 -10 -8 42 v=1 2 4 6 ,u=5 0 0 -6 7,deconv(u,v)v = 1 2 4 6u = 5 0 0 -6 7ans = 5 -1010、（点乘是数组的运算，没有点的乘是矩阵运算） a = 2 5;3 4, b =3 1;4 7,a.*b,a*ba = 2 5 3 4b = 3 1 4 7ans = 6 5 12 28ans = 26 3725 31 a = 2 3; b = 4 7;a.*b = 8 21;a*b %错误a*b = 29;11、（who 可以看到你用过的一些变量，whos是把该变量及所存储的大小等信息都显示出来了） whoYour variables are:A B a ans b p u v whos Name Size Bytes Class Attributes A 2x2 32 double B 2x2 32 double a 1x2 16 double ans 1x2 16 double b 1x2 16 double p 1x3 24 double u 1x5 40 double v 1x4 32 double 12、 A=2 5 3;6 5 4,disp(A),size(A),length(A)A = 2 5 3 6 5 4 2 5 3 6 5 4ans = 2 3ans = 3实验二 MATLAB绘图命令实验任务 熟悉MATLAB基本绘图命令，掌握如下绘图方法： 1坐标系的选择、图形的绘制； 2图形注解（题目、标号、说明、分格线）的加入； 3图形线型、符号、颜色的选取。基本命令训练1、t=0:pi/360:2*pi;x=cos(t)+ cos(t*4);y=sin(t)+ sin(t*4);xlabel(x轴);ylabel(y轴);plot(y,x),grid; 2、t=0:0.1:100;x=3*t;y=4*t;z=sin(2*t);plot3(x,y,z,g:)3、x = linspace(-2*pi,2*pi,40);y=sin(x);stairs(x,y) 4、t=0:pi/360:2*pi;x=cos(t)+ cos(t*4) + sin(t*4);y=sin(t)+ sin(t*4);plot(y,x,r:);xlabel(x轴);ylabel(y轴); 5、th=0:pi/1000:2*pi;r=cos(2*th);polar(th,r);title(四叶草图)6、th=0:pi/20:2*pi;x=exp(j*th);plot(real(x),imag(x),r-.) ;grid; text(0,0,中心) ; 7、 x=-2:0.01:2;8、 y=-2:0.01:2;9、 X,Y = meshgrid(x,y);Z = Y.*exp(-X.2-Y.2);C,h = contour(X,Y,Z);set(h,ShowText,on,TextStep,get(h,LevelStep)*2)8、x = 0:0.2:10;y = 2*x+3;subplot(411);plot(x,y); grid;title(y的原函数);subplot(412) ;semilogy(x,y); grid;title(对y取对数);subplot(413) ;semilogx(x,y); grid;title(对x取对数);subplot(414) ;loglog(x,y);grid;title(对xy均取对数); 9、x = -3:0.3:3;bar(x,exp(-x.*x),g) 实验三 MATLAB程序设计实验任务 1熟悉MATLAB程序设计的方法和思路； 2掌握循环、分支语句的编写，学会使用look for、help命令。程序举例1、f=1,1;i=1;while f(i)+f(i+1)m=3;n=4;for i=1:m for j=1:n a(i,j)=1/(i+j-1); endendformat rataa = 1 1/2 1/3 1/4 1/2 1/3 1/4 1/5 1/3 1/4 1/5 1/6 (分数格式形式。用有理数逼近显示数据)m=5;n=4;for i=1:m for j=1:n a(i,j)=1/(i+j-1); endend format rataa = 1 1/2 1/3 1/4 1/2 1/3 1/4 1/5 1/3 1/4 1/5 1/6 1/4 1/5 1/6 1/7 1/5 1/6 1/7 1/8 3、程序中没有format rat命令时，如果上次运行结果没有清除，输出的结果就是上次运行的结果！但是运用clear命令清楚之前的运行结果之后就会正常运行。4、x=input(请输入x的值:); if x=10 y=cos(x+1)+sqrt(x*x+1); else y=x*sqrt(x+sqrt(x); end y请输入x的值:2y = 2391/647x=input(请输入x的值:); if x=10 y=fprintf(不在定义域内，请重新输入：);return else y=1/(x-10); end y请输入x的值:2y = -1/85、p=0 0 0 1 3 0 2 0 0 9;for i=1:length(p),if p(1)=0,p=p(2:length(p); end;end;pp = Columns 1 through 5 1 3 0 2 0 Columns 6 through 7 0 9 p=0 0 0 1 3 0 2 0 0 9;p(p=0)=;pp = 1 3 2 96、 e2(500)ans = 1 1 2 3 5 8 13 21 34 55 89 144 233 377 lookfor ffibnoe2 - ffibno 计算斐波那契亚数列的函数文件 help e2 ffibno 计算斐波那契亚数列的函数文件 n可取任意自然数 程序如下（用法： lookfor 关键词在所有M文件中找“关键词”，比如：lookfor max（即寻找关键词“max”）其实就和我们平时用CTRL+F来查找“关键词”是一样的而help是显示matlab内置的帮助信息 用法：help 命令，比如 help inv ，作用就是调用inv这个命令的帮助）程序设计题用一个MATLAB语言编写一个程序：输入一个自然数，判断它是否是素数，如果是，输出“It is one prime”，如果不是，输出“It is not one prime.”。要求通过调用子函数实现。最好能具有如下功能：设计较好的人机对话界面，程序中含有提示性的输入输出语句。能实现循环操作，由操作者输入相关命令来控制是否继续进行素数的判断。如果操作者希望停止这种判断，则可以退出程序。如果所输入的自然数是一个合数，除了给出其不是素数的结论外，还应给出至少一种其因数分解形式。例：输入 6， 因为6不是素数。则程序中除了有“It is not one prime”的结论外，还应有：“6=2*3”的说明。function sushuwhile 1 x=input( 请输入一个自然数);if x a=sym(a11 a12;a21 a22);da=det(a)ea=eig(a)da = a11*a22-a12*a21 ea = 1/2*a11+1/2*a22+1/2*(a112-2*a11*a22+a222+4*a12*a21)(1/2) 1/2*a11+1/2*a22-1/2*(a112-2*a11*a22+a222+4*a12*a21)(1/2)a=sym(2 3;1 5);da=det(a)ea=eig(a)da =7ea = 7/2+1/2*21(1/2) 7/2-1/2*21(1/2)2. 求方程的解（包括精确解和一定精度的解） r1=solve(x2+x-1)rv=vpa(r1)rv4=vpa(r1,4)rv20=vpa(r1,20) r1 = 1/2*5(1/2)-1/2 -1/2*5(1/2)-1/2rv = .61803398874989484820458683436560 -1.6180339887498948482045868343656 rv4 = .6180 -1.618rv20 = .61803398874989484820 -1.61803398874989484823 a=sym(a);b=sym(b);c=sym(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, dB = 10 5 -8 11 ans = a*d-b*cans = 1504. syms x y;s=(-7*x2-8*y2)*(-x2+3*y2);expand(s) %对s展开collect(s,x) %对s按变量x合并同类项(无同类项)factor(ans) % 对ans分解因式ans =7*x4-13*x2*y2-24*y4ans = 7*x4-13*x2*y2-24*y4 ans = (8*y2+7*x2)*(x2-3*y2)5. 对方程 AX=b求解 A=34,8,4;3,34,3;3,6,8;b=4;6;2;X=linsolve(A,b) %调用linsolve函数求解Ab %用另一种方法求解X = 0.0675 0.1614 0.1037ans = 0.0675 0.16140.10376 对方程组求解a11*x1+a12*x2+a13*x3=b1a21*x1+a22*x2+a23*x3=b2a31*x1+a32*x2+a33*x3=b3syms a11 a12 a13 a21 a22 a23 a31 a32 a33 b1 b2 b3;A=a11,a12,a13;a21,a22,a23;a31,a32,a33;b=b1;b2;b3;XX=Ab %用左除运算求解(X=linsolve(A,b) %调用linsolve函数求的解)XX = (a12*a23*b3-a12*b2*a33+a13*a32*b2-a13*a22*b3+b1*a22*a33-b1*a32*a23)/(a11*a22*a33-a11*a32*a23-a12*a21*a33+a32*a21*a13-a22*a31*a13+a31*a12*a23) -(a11*a23*b3-a11*b2*a33-a21*a13*b3-a23*a31*b1+b2*a31*a13+a21*b1*a33)/(a11*a22*a33-a11*a32*a23-a12*a21*a33+a32*a21*a13-a22*a31*a13+a31*a12*a23) (a32*a21*b1-a11*a32*b2+a11*a22*b3-a22*a31*b1-a12*a21*b3+a31*a12*b2)/(a11*a22*a33-a11*a32*a23-a12*a21*a33+a32*a21*a13-a22*a31*a13+a31*a12*a23)7syms 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)/diff(f1) %按参数方程求导公式求y对x的导数ans = 1/2/(1+exp(x)(1/2)*exp(x) ans = -2*sin(x)-x*cos(x) ans = -3*cos(x)+x*sin(x) ans = -b*cos(t)/a/sin(t)三、SIMULINK的使用G1(s)G2(s)R(s)C(s)其中：R(s)为阶跃输入，C(s)为输出 仿真图：波形图：实验五 MATLAB在控制系统分析中的应用实验任务1掌握MATLAB在控制系统时间响应分析中的应用；2掌握MATLAB在系统根轨迹分析中的应用； 3. 掌握MATLAB控制系统频率分析中的应用； 4. 掌握MATLAB在控制系统稳定性分析中的应用基本命令 1. step 2. impulse 3. initial 4. lsim 5. rlocfind 6. bode 7. margin 8. nyquist 9. Nichols 10. cloop程序举例1. 求下面系统的单位阶跃响应 num=4 ; den=1 , 1 , 4 ;step(num , den)y , x , t=step(num , den) ;tp=spline(y , t , max(y) %计算峰值时间max(y) %计算峰值tp = 1.6062ans =1.44412. 求如下系统的单位阶跃响应 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,x0)subplot(1 , 2 , 2) ; impulse(a , b , c , d)5：系统传递函数为：输入正弦信号时，观察输出 信号的相位差。 num=1 ; den=1 , 1 ;t=0 : 0.01 : 10 ;u=sin(2*t) ; hold onplot(t,u, r)lsim(num,den,u,t)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,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)Select a point in the graphics windowselected_point = -2.5924 - 0.0248ik = 0.7133p = -3.4160 -2.5918 -0.9961 + 0.4306i -0.9961 - 0.4306i8. 作如下系统的bode图 n=1 , 1 ; d=1 , 4 , 11 , 7 ; bode(n , d),grid on9. 系统传函如下 求有理传函的频率响应，然后在同一张图上绘出以四阶伯德近似表示的系统频率响应 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);semilogx(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);10. 已知系统模型为 求它的幅值裕度和相角裕度 n=3.5; d=1 2 3 2; Gm,Pm,Wcg,Wcp=margin(n,d)Gm = 1.1433Pm = 7.1688Wcg = 1.7323Wcp =1.654111. 二阶系统为：令wn=1,分别作出=2 , 1 , 0.707 , 0.5时的nyquist曲线。 n=1 ; d1=1 , 4 , 1 ; d2=1 , 2 , 1 ; d3=1 , 1.414 , 1; d4=1,1,1;nyquist(n,d1) ;hold onnyquist(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 , 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(num,den);impulse(num2,den2);p = -10.5969 +36.2148i -10.5969 -36.2148i -0.0562 15. 已知系统为：作该系统的nichols曲线。 n=1 ; d=1 , 1 , 0 ; 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; %波特图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.2393 + 0.8579i -0.2393 - 0.8579ians = -3.3089 0.1545 + 1.7316i 0.1545 - 1.7316i18. 系统的状态方程为： 试确定系统的稳定性。 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理解欧拉法和龙格-库塔法的基本思想； 2理解数值积分算法的计算精度、速度、稳定性与步长的关系；程序举例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.001; % 步长修改为0.001for k=1:5 t(k+1)=t(k)+h; y(k+1)=sqrt(1+2*t(k+1);endfor k=1:5 y_euler(k+1)=y_euler(k)+h*(y_euler(k)-2*t(k)/y_euler(k);endfor 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;endfor k=1:5 k1=y_rk4(k)-2*t(k)/y_rk4(k); k2=(y_rk4(k)+h*k1/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) 时间 解析解 欧拉法 RK2法 RK4法 0 1.0000 1.0000 1.0000 1.0000 0.0010 1.0010 1.0010 1.0010 1.0010 0.0020 1.0020 1.0020 1.0020 1.0020 0.0030 1.0030 1.0030 1.0030 1.0030 0.0040 1.0040 1.0040 1.0040 1.0040 0.0050 1.0050 1.0050 1.0050 1.0050 2. 考虑如下二阶系统： 在上的数字仿真解（已知：，），并将不同步长下的仿真结果与解析解进行精度比较。 说明：已知该微分方程的解析解分别为： 采用RK4法进行计算，选择状态变量： 则有如下状态空间模型及初值条件 采用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)=100*cos(t(k+1); y_05(k+1)=100*exp(-t(k+1)/2).*cos(sqrt(3)/2*t(k