中南大学系统仿真实验报告.docx
《中南大学系统仿真实验报告.docx》由会员分享,可在线阅读,更多相关《中南大学系统仿真实验报告.docx(32页珍藏版)》请在冰豆网上搜索。
![中南大学系统仿真实验报告.docx](https://file1.bdocx.com/fileroot1/2023-1/23/d68ce6df-04eb-4518-b9e6-08e142b31ad9/d68ce6df-04eb-4518-b9e6-08e142b31ad91.gif)
中南大学系统仿真实验报告
实验一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=
>>rand(2,3)
ans=
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=
>>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=
ans=
ans=
ans=
ans=
7、>>p=[1,-6,-72,-27],roots(p)
p=
1-6-72-27
ans=
>>p=[2,3,6],roots(p)
p=
236
ans=
+
-
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:
:
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:
:
2;
8、y=-2:
:
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:
:
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:
:
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=
.560
rv4=
.6180
rv20=
.6484820
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=
ans=
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的使用
其中:
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=
ans=
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:
:
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=
-
k=
p=
+
-
8.作如下系统的bode图
>>n=[1,1];d=[1,4,11,7];
bode(n,d),gridon
9.系统传函如下
求有理传函的频率响应,然后在同一张图上绘出以四阶伯德近似表示的系统频率响应
>>num=[1];den=conv([12],conv([12],[12]));
w=logspace(-1,2);t=;
[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)