中南大学系统仿真实验报告.docx
《中南大学系统仿真实验报告.docx》由会员分享，可在线阅读，更多相关《中南大学系统仿真实验报告.docx（28页珍藏版）》请在冰点文库上搜索。
中南大学系统仿真实验报告
实验一MATLAB中矩阵与多项式的基本运算
实验任务
1．了解MATLAB命令窗口和程序文件的调用。
2．熟悉如下MATLAB的基本运算：
①矩阵的产生、数据的输入、相关元素的显示；
②矩阵的加法、乘法、左除、右除；
③特殊矩阵：
单位矩阵、“1”矩阵、“0”矩阵、对角阵、随机矩阵的产生和运算；
④多项式的运算：
多项式求根、多项式之间的乘除。
基本命令训练
1、>>eye
（2）
ans=
10
01
>>eye（4）
ans=
1000
0100
0010
0001
2、>>ones
（2）
ans=
11
11
>>ones（4）
ans=
1111
1111
1111
1111
>>ones（2,2）
ans=
11
11
>>ones（2,3）
ans=
111
111
>>ones（4,3）
ans=
111
111
111
111
3、>>zeros
（2）
ans=
00
00
>>zeros（4）
ans=
0000
0000
0000
0000
>>zeros（2,2）
ans=
00
00
>>zeros（2,3）
ans=
000
000
>>zeros（3,2）
ans=
00
00
10
4、随机阵>>rand（2,3）
ans=
0.27850.95750.1576
0.54690.96490.9706
>>rand（2,3）
ans=
0.95720.80030.4218
0.48540.14190.9157
5、>>diag（5）
ans=
5
>>diag（5,5）
ans=
000005
000000
000000
000000
000000
000000
>>diag（2,3）
ans=
0002
0000
0000
0000
6、
（inv（A）为求A的逆矩阵）>>B=[531;238;111],inv（B）
B=
531
238
111
ans=
0.62500.2500-2.6250
-0.7500-0.50004.7500
0.12500.2500-1.1250
>>A=[23;44],B=[53;38],inv（A）,inv（B）;A\B,A/B,inv（A）*B,B*inv（A）
A=
23
44
B=
53
38
ans=
-1.00000.7500
1.0000-0.5000
ans=
-2.75003.0000
3.5000-1.0000
ans=
0.22580.2903
0.64520.2581
ans=
-2.75003.0000
3.5000-1.0000
ans=
-2.00002.2500
5.0000-1.7500
7、>>p=[1,-6,-72,-27],roots（p）
p=
1-6-72-27
ans=
12.1229
-5.7345
-0.3884
>>p=[2,3,6],roots（p）
p=
236
ans=
-0.7500+1.5612i
-0.7500-1.5612i
8、（A为n*n的方阵）>>A=[010;-440;-212],poly（A）,B=sym（A）,poly（B）
A=
010
-440
-212
ans=
1-612-8
B=
[0,1,0]
[-4,4,0]
[-2,1,2]
ans=
x^3-6*x^2+12*x-8
9,、（conv是多项式相乘，deconv是多项式相除）
>>u=[1246],v=[500-67],conv（u,v）
u=
1246
v=
500-67
ans=
5102024-5-10-842
>>v=[1246],u=[500-67],deconv（u,v）
v=
1246
u=
500-67
ans=
5-10
10、（点乘是数组的运算，没有点的乘是矩阵运算）
>>a=[25;34],b=[31;47],a.*b,a*b
a=
25
34
b=
31
47
ans=
65
1228
ans=
2637
2531
>>a=[23];b=[47];
a.*b=[821];
a*b%错误
a*b'=29;
11、（who可以看到你用过的一些变量，whos是把该变量及所存储的大小等信息都显示出来了）>>who
Yourvariablesare:
ABaansbpuv
>>whos
NameSizeBytesClassAttributes
A2x232double
B2x232double
a1x216double
ans1x216double
b1x216double
p1x324double
u1x540double
v1x432double
12、>>A=[253;654],disp（A）,size（A）,length（A）
A=
253
654
253
654
ans=
23
ans=
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．掌握循环、分支语句的编写，学会使用lookfor、help命令。
程序举例
1、>>f=[1,1];
i=1;
whilef（i）+f（i+1）<1000
f（i+2）=f（i）+f（i+1）;
i=i+1;
end
f,i
f=
Columns1through14
1123581321345589144233377
Columns15through16
610987
i=
15
2、>>m=3;
n=4;
fori=1:
m
forj=1:
n
a（i,j）=1/（i+j-1）;
end
end
formatrat
a
a=
11/21/31/4
1/21/31/41/5
1/31/41/51/6
（分数格式形式。
用有理数逼近显示数据）
>>m=5;
n=4;
fori=1:
m
forj=1:
n
a（i,j）=1/（i+j-1）;
end
end
formatrat
a
a=
11/21/31/4
1/21/31/41/5
1/31/41/51/6
1/41/51/61/7
1/51/61/71/8
3、程序中没有formatrat命令时，如果上次运行结果没有清除，输出的结果就是上次运行的结果！
但是运用clear命令清楚之前的运行结果之后就会正常运行。
4、>>x=input（'请输入x的值:
'）;
ifx==10
y=cos（x+1）+sqrt（x*x+1）;
else
y=x*sqrt（x+sqrt（x））;
end
y
请输入x的值:
2
y=
2391/647
x=input（'请输入x的值:
'）;
ifx==10
y=fprintf（'不在定义域内，请重新输入：
'）;return
else
y=1/（x-10）;
end
y
请输入x的值:
2
y=
-1/8
5、>>p=[0001302009];
fori=1:
length（p）,ifp
（1）==0,p=p（2:
length（p））;
end;
end;
p
p=
Columns1through5
13020
Columns6through7
09
>>p=[0001302009];p（p==0）=[];p
p=
1329
6、>>e2（500）
ans=
1123581321345589144233377
>>lookforffibno
e2-ffibno计算斐波那契亚数列的函数文件
>>helpe2
ffibno计算斐波那契亚数列的函数文件
n可取任意自然数
程序如下
（用法：
lookfor关键词
在所有M文件中找“关键词”，比如：
lookformax（即寻找关键词“max”）
其实就和我们平时用CTRL+F来查找“关键词”是一样的
而help是显示matlab内置的帮助信息
用法：
help命令，比如helpinv，作用就是调用inv这个命令的帮助）
程序设计题
用一个MATLAB语言编写一个程序：
输入一个自然数，判断它是否是素数，如果是，输出“Itisoneprime”，如果不是，输出“Itisnotoneprime.”。
要求通过调用子函数实现。
最好能具有如下功能：
①设计较好的人机对话界面，程序中含有提示性的输入输出语句。
②能实现循环操作，由操作者输入相关命令来控制是否继续进行素数的判断。
如果操作者希望停止这种判断，则可以退出程序。
③如果所输入的自然数是一个合数，除了给出其不是素数的结论外，还应给出至少一种其因数分解形式。
例：
输入6，因为6不是素数。
则程序中除了有“Itisnotoneprime”的结论外，还应有：
“6=2*3”的说明。
functionsushu
while1
x=input（'请输入一个自然数'）;
ifx<2
disp（'既不是质数又不是合数'）;
Elseifisprime（x）==1
disp（'这是一个素数'）;
Elsedisp（'这是一个合数，可以因式分解为：
'）
forn=2:
sqrt（x）
ifrem（x,n）==0
num3=x;
num1=n;
num2=x/n;
disp（[num2str（num3）,'=',num2str（num1）,'x',num2str（num2）]）
End
End
End
end
y=input（'是否继续判断？
继续请按1，按任意键退出:
'）;
ify~=1
break
end
end
实验四MATLAB的符号计算与SIMULINK的使用
实验任务
1．掌握MATLAB符号计算的特点和常用基本命令；
2．掌握SIMULINK的使用。
程序举例
1．求矩阵对应的行列式和特征根
>>a=sym（'[a11a12;a21a22]'）;
da=det（a）
ea=eig（a）
da=
a11*a22-a12*a21
ea=
1/2*a11+1/2*a22+1/2*（a11^2-2*a11*a22+a22^2+4*a12*a21）^（1/2）
1/2*a11+1/2*a22-1/2*（a11^2-2*a11*a22+a22^2+4*a12*a21）^（1/2）
a=sym（'[23;15]'）;
da=det（a）
ea=eig（a）
da=
7
ea=
7/2+1/2*21^（1/2）
7/2-1/2*21^（1/2）
2.求方程的解（包括精确解和一定精度的解）
>>r1=solve（'x^2+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.8343656
rv4=
.6180
-1.618
rv20=
.4820
-1.482
3．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,d]
B=
105
-811
ans=
a*d-b*c
ans=
150
4.>>symsxy;s=（-7*x^2-8*y^2）*（-x^2+3*y^2）;
expand（s）%对s展开
collect（s,x）%对s按变量x合并同类项（无同类项）
factor（ans）%对ans分解因式
ans=
7*x^4-13*x^2*y^2-24*y^4
ans=
7*x^4-13*x^2*y^2-24*y^4
ans=
（8*y^2+7*x^2）*（x^2-3*y^2）
5.对方程AX=b求解>>A=[34,8,4;3,34,3;3,6,8];
b=[4;6;2];
X=linsolve（A,b）%调用linsolve函数求解
A\b%用另一种方法求解
X=
0.0675
0.1614
0.1037
ans=
0.0675
0.1614
0.1037
6．对方程组求解
a11*x1+a12*x2+a13*x3=b1
a21*x1+a22*x2+a23*x3=b2
a31*x1+a32*x2+a33*x3=b3
>>symsa11a12a13a21a22a23a31a32a33b1b2b3;
A=[a11,a12,a13;a21,a22,a23;a31,a32,a33];
b=[b1;b2;b3];
XX=A\b%用左除运算求解
（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）
7．symsabtxyz;
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.step2.impulse3.initial4.lsim5.rlocfind
6.bode7.margin8.nyquist9.Nichols10.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.6062
ans=
1.4441
2.求如下系统的单位阶跃响应
>>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）;holdon
plot（t,u,'r'）
lsim（num,den,u,t）
6.有一二阶系统，求出周期为4秒的方波的输出响应
>>num=[251];
den=[123];
t=（0:
.1:
10）;
period=4;
u=（rem（t,period）>=period./2）;
%看rem函数功能
lsim（num,den,u,t）;
7.已知开环系统传递函数，绘制系统的根轨迹，并分析其稳定性
>>num=[12];
den1=[143];
den=conv（den1,den1）;
figure
（1）
rlocus（num,den）
[k,p]=rlocfind（num,den）
figure
（2）
k=55;
num1=k*[12];
den=[143];
den1=conv（den,den）;
[num,den]=cloop（num1,den1,-1）;
impulse（num,den）
title（'impulseresponse（k=55）'）
figure（3）
k=56;
num1=k*[12];
den=[143];
den1=conv（den,den）;
[num,den]=cloop（num1,den1,-1）;
impulse（num,den）
title（'impulseresponse（k=56）'）
Selectapointinthegraphicswindow
selected_point=
-2.5924-0.0248i
k=
0.7133
p=
-3.4160
-2.5918
升级会员 