信息科学与工程学院系统仿真实验报告

1、中南大学系统仿真实验报告 指导老师 胡杨 实验者 学 号 专业班级 实验日期 2014.6.4 学 院 信息科学与工程学院 目录实验一 matlab中矩阵与多项式的基本运算.3实验二 matlab绘图命令.7实验三 matlab程序设计.9实验四 matlab的符号计算与simulink的使用.13 实验五 matlab在控制系统分析中的应用.17 实验六 连续系统数字仿真的基本算法.30实验一 matlab中矩阵与多项式的基本运算一、实验任务1了解matlab命令窗口和程序文件的调用。2熟悉如下matlab的基本运算： 矩阵的产生、数据的输入、相关元素的显示； 矩阵的加法、乘法、左除、右除；

2、 特殊矩阵：单位矩阵、“1”矩阵、“0”矩阵、对角阵、随机矩阵的产生和运算； 多项式的运算：多项式求根、多项式之间的乘除。二、基本命令训练1eye(m)m=3;eye(m)ans = 1 0 0 0 1 0 0 0 12ones(n)、ones(m,n)n=1;m=2;ones(n)ones(m,n)ans = 1ans = 1 13zeros(m,n)m=1,n=2;zeros(m,n)m = 1ans = 0 04rand(m,n)m=1;n=2;rand(m,n)ans = 0.8147 0.90585diag(v)v=1 2 3;diag(v)ans = 1 0 0 0 2 0 0 0

3、 36ab 、a/b、 inv(a)*b 、b*inv(a)a=1 2;3 4;b=5 6;7 8;a=ab b=a/bc=inv(a)*b d=b*inv(a)a = -3 -4 4 5b = 3.0000 -2.0000 2.0000 -1.0000c = -3.0000 -4.0000 4.0000 5.0000d = -1.0000 2.0000 -2.0000 3.00007roots(p)syms x;a=3*x.3+2*x+1;p=3,0,2,1;roots(p)ans = 0.2012 + 0.8877i 0.2012 - 0.8877i -0.4023 8polya=1 2;

4、3 4;poly(a)ans = 1.0000 -5.0000 -2.00009conv 、deconva=1 2;b=5 6;a=conv(a,b)b=deconv(a,b)a = 5 16 12b = 0.200010a*b 与 a.*b的区别a=1 2;b=5 6;a=a*b a=1 2;b=5 6;b=a.*ba = 17b = 5 1211who与whos的使用a=1 2;3 4;whowhosyour variables are:a name size bytes class attributes a 2x2 32 double12disp、size(a)、length(a)的使用

5、a=a b c d e f;disp(a)a=1 2 3 4;b=size(a)c=length(a)a b c d e fb = 1 4c = 4三、实验要求根据实验内容和相关命令进行实验，自拟输入元素，将上述各命令的输入和输出结果写成实验报告一（全部实验完成后交实验报告）。实验二 matlab绘图命令一、实验任务熟悉matlab基本绘图命令，掌握如下绘图方法：1坐标系的选择、图形的绘制；2图形注解（题目、标号、说明、分格线）的加入；3图形线型、符号、颜色的选取。二、基本命令训练1plot 2loglog 3semilogx 4semilogy5polar 6title 7xlabel 8y

6、label9text 10grid 11bar 12stairs13contour三、实验举例1t=0:pi/360:2*pi*22/3; x=93*cos(t)+36*cos(t*4.15); y=93*sin(t)+36*sin(t*4.15); plot(y,x),grid; %绘制二维坐标网格图 2 t=0:0.05:100; x=t;y=2*t;z=sin(2*t); plot3(x,y,z,r-.) %绘制三维坐标图3 t=0:pi/20:2*pi; y=sin(x); stairs(x,y) %绘制阶梯图4 th=pi/200:pi/200:2*pi; r=cos(2*th);

7、polar(th,r),grid %在网格里画极坐标图5th=0:pi/10:2*pi; x=exp(j*th); %x为复数 plot(real(x),imag(x),r*); %以实部为横轴，虚部为纵轴画图 grid;四、实验要求 在两种或两种以上坐标系绘制35个图形，要有颜色、图形种类、注解的不同实验结果写成实验报告二（全部实验完成后交实验报告）。实验三 matlab程序设计一、实验任务1熟悉matlab程序设计的方法和思路；2掌握循环、分支语句的编写，学会使用look for、help命令。二、程序举例1计算11000之内的斐波那契亚数列 f=1,1;i=1;while f(i)+f(

8、i+1)1000 f(i+2)=f(i)+f(i+1); i=i+1;end f,i f = columns 1 through 14 1 1 2 3 5 8 13 21 34 55 89 144 233 377 columns 15 through 16 610 987i = 152 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/63 m=3; n=4;for i=1:m for j=1:n a(i,j)=1/(i+

9、j-1); endendaa = 1 1/2 1/3 1/4 1/2 1/3 1/4 1/5 1/3 1/4 1/5 1/6 请比较程序2与程序3的区别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/647 5去掉多项式或数列开头的零项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;p p = columns 1 through 5 1 3 0

10、 2 0 columns 6 through 7 0 96 建立matlab的函数文件，程序代码如下，以文件名ex2_4.m存盘function f=ffibno(n)%ffibno 计算斐波那契亚数列的函数文件%n可取任意自然数%程序如下f=1,1;i=1;while f(i)+f(i+1)editex2_4(200)ans = columns 1 through 5 1 1 2 3 5 columns 6 through 10 8 13 21 34 55 columns 11 through 12 89 144 lookfor ffibnoex2_4 - ffibno 计算斐波那契亚数列的

11、函数文件ex2_4 - ffibno 计算斐波那契亚数列的函数文件 help ex2_4ffibno 计算斐波那契亚数列的函数文件n可取任意自然数程序如下输入完毕后在matlab的命令窗口输入ex2_4(200)，得到运行结果。在matlab的命令窗口输入lookfor ffibno，得到结果：ex2_4.m: %ffibno 计算斐波那契亚数列的函数文件在matlab的命令窗口输入help ex2_4，得到结果： ffibno 计算斐波那契亚数列的函数文件 n可取任意自然数 程序如下三、程序设计题 用一个matlab语言编写一个程序：输入一个自然数，判断它是否是素数，如果是，输出“it is

12、 one prime”，如果不是，输出“it is not one prime.”。要求通过调用子函数实现。最好能具有如下功能：设计较好的人机对话界面，程序中含有提示性的输入输出语句。能实现循环操作，由操作者输入相关命令来控制是否继续进行素数的判断。如果操作者希望停止这种判断，则可以退出程序。如果所输入的自然数是一个合数，除了给出其不是素数的结论外，还应给出至少一种其因数分解形式。例：输入 6， 因为6不是素数。则程序中除了有“it is not one prime”的结论外，还应有：“6=2*3”的说明。close all;c=1;c=input(是否进行素数运算 1为是 0为否: );wh

13、ile c=1a=input(请输入一个自然数： );if factor(a)=a disp(it is one prime)else disp(it is not one prime); b=factor(a); fprintf(%3d =,a) for j=1:(length(b)-1) fprintf(%3d *,b(j) end fprintf(%3d n,b(length(b)endc=input(是否进行素数运算 1为是 0为否: );end是否进行素数运算 1为是 0为否: 1请输入一个自然数： 20it is not one prime 20 = 2 * 2 * 5 是否进行素

14、数运算 1为是 0为否: 1请输入一个自然数： 17it is one prime是否进行素数运算 1为是 0为否: 四、实验要求1参照已知程序，改动程序中的参数和输入量，验证输出结果。2使用lookfor、help命令，验证输出结果3实验结果写成实验报告三（全部实验完成后交实验报告）。 实验四 matlab的符号计算与simulink的使用一、实验任务1掌握matlab符号计算的特点和常用基本命令；2掌握simulink的使用。二、程序举例1求矩阵对应的行列式和特征根a=sym(a11 a12;a21 a22);da=det(a)ea=eig(a)da =a11*a22-a12*a21ea

15、= 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)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/2 rv = 1.6180339887498948482045868343656 -.618033988749894848204586834

16、36560 rv4 = 1.618 -.6180 rv20 = 1.6180339887498948482 -.618033988749894848203 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 %建立符号矩阵a b=w,x;y,z %建立数值矩阵b det(a) %计算符号矩阵a的行列式 det(b) %计算数值矩阵b的行列式a = a, b c, db = 10 5 -8 11ans =a*d-b*cans = 1504 syms x y; s=(-7*x2

17、-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*y4ans =(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.1037 ans = 0.0675 0.1614 0.

18、10376 对方程组求解 a11*x1+a12*x2+a13*x3=b1 a21*x1+a22*x2+a23*x3=b2 a31*x1+a32*x2+a33*x3=b3 syms 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 %用左除运算求解 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-

19、a21*a12*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-a21*a12*a33+a32*a21*a13-a22*a31*a13+a31*a12*a23) (a32*a21*b1-a11*a32*b2+a11*a22*b3-a22*a31*b1-a21*a12*b3+a31*a12*b2)/(a11*a22*a33-a11*a32*a23-a21*a12*a33+a32

20、*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) an

21、s = -b*cos(t)/a/sin(t)三、simulink的使用1在命令窗口中输入：simulink（回车）得到如下simulink模块： 2双击打开各模块，选择合适子模块构造控制系统，例如：3双击各子模块可修改其参数，选择simulation菜单下的start命令运行仿真，在示波器（scope）中观察结果。四、实验要求1符号计算部分：参照所给示例，自拟输入量，验证求矩阵的特征根、行列式及方程求解命令2simulink部分：逐一熟悉各模块的使用，对下面的随动系统模型进行仿真，实验报告中包含：连好的系统模型及用scope观测的结果g1(s)g2(s)r(s)c(s) 其中：r(s)为阶跃输

22、入，c(s)为输出 3实验结果写成实验报告四备注：全部实验完成后，在实验报告的封面注明：课程、班级、姓名、学号等信息。按时交实验报告。 实验五 matlab在控制系统分析中的应用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

23、,b,c,d);plot(y);grid;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);grid;subplot(1 , 2 , 2) ;impulse(a , b , c , d);grid; 5：系统传递函数为：输入正弦信号时，观察输出信号的相位差。

24、%程序如下：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: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,

25、den1);figure(1)rlocus(num,den);grid;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);grid;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);grid;tit

26、le(impulse response(k=56)select a point in the graphics windowselected_point = 2.0142 - 8.1739ik = 766.4105p = -11.2238 2.6126 + 7.8623i 2.6126 - 7.8623i -2.0013 8. 作如下系统的bode图：n=1 , 1 ; d=1 , 4 , 11 , 7 ;bode(n , d);grid;9. 系统传函如下求有理传函的频率响应，然后在同一张图上绘出以四阶伯德近似表示的系统频率响应%程序如下： num=1;den=conv(1 2,conv(1

27、 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

28、-);grid on;xlabel(frequency);ylabel(phase);10. 已知系统模型为求它的幅值裕度和相角裕度%程序如下： n=3.5; d=1 2 3 2; gm,pm,wcg,wcp=margin(n,d) gm = 1.1433 pm = 7.1688 wcg = 1.7323 wcp = 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)

29、;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. 一多环系统

30、，其结构图如下，使用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.056215. 已知系统为： ,作该系统的nichols

31、曲线。%程序如下：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; %波特图17. 系统的开环传递函数为：分别确定k=2和k=10时闭环系统的稳定性。%程序如下：d1=1 , 3 ,

32、2 , 0 ; n1=2 ;nc1 , dc1=cloop(n1 , d1 ,-1) ;a=roots(dc1)if (a0) disp(the system is stable)else disp(the system is not stable)endd2=d1 ; n2=10 ;nc2 , dc2=cloop(n2 , d2,-1);b=roots(dc2)if (b0) disp(the system is stable)else disp(the system is not stable)enda = -2.5214 -0.2393 + 0.8579i -0.2393 - 0.857

33、9ithe system is stableb = -3.3089 0.1545 + 1.7316i 0.1545 - 1.7316ithe system is not stable18. 系统的状态方程为：试确定系统的稳定性。%程序如下：a=-4,-3,0 ; 1,0,0 ; 0,1,0 ; b=1;0;0 ; c=0,1,2 ; d=0 ;eig(a) %求特征根nc=rank(ctrb(a,b)n=length(a);if n=ncdisp(system is completely state controllable)elsedisp(system is not completely

34、state controllable)endans = 0 -1 -3nc = 3system is completely state controllable四、实验要求1参照已知程序，改动程序中的参数和输入量，验证输出结果。2掌握相关命令的用法3实验结果写成实验报告五（全部实验完成后交实验报告）。实验六 连续系统数字仿真的基本算法一、实验任务1理解欧拉法和龙格-库塔法的基本思想；2理解数值积分算法的计算精度、速度、稳定性与步长的关系；二、程序举例1. 取h=0.2，试分别用欧拉法、rk2法和rk4法求解微分方程的数值解，并比较计算精度。 注：解析解：%程序如下：t(1)=0; % 置自变量

35、初值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; 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

36、); 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*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.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. 考虑如下二阶系统