matlab仿真课后习题

第一章习题
3.请指出以下的变量名（函数名、M文件名）中，哪些是合法的Abc 2004x lil-1 wu_2004 a&b _xyz 解：合法的变量名有：Abc wu_2004
4．指令窗操作
（1）求[12+2×（7-4）]÷32的运算结果
解：>> [12+2*(7-4)]/3^2
ans =
2
（2）输入矩阵A=[1，2，3；4，5，6；7，8，9]，观察输出。

解：>> A=[1,2,3;4,5,6;7,8,9]
A =
1 2 3
4 5 6
7 8 9
（3）输入以下指令，观察运算结果；
clear;x=-8::8;
y=x';
X=ones(size(y))*x;
Y=y*ones(size(x));
R=sqrt(X.^2+Y.^2)+eps;
Z=sin(R)./R;
mesh(X,Y,Z);
colormap(hot)
xlabel('x'),ylabel('y'),zlabel('z')
解：
7．指令行编辑
（1）依次键入以下字符并运行：y1=2*sin*pi)/(1+sqrt(5))
解：>>y1=2*sin*pi)/(1+sqrt(5))
y1 =
(2)通过反复按键盘的箭头键，实现指令回调和编辑，进行新的计算；y2=2*cos*pi)/(1+sqrt(5))
解：>>y2=2*cos*pi)/(1+sqrt(5))
y2 =
11.编写题4中（3）的M脚本文件，并运行之。

解：
第二章习题
1.在指令窗中键入x=1::2和y=2::1,观察所生成的数组。

解:>> x=1::2 x =
>> y=2::1 y =
Empty matrix: 1-by-0
2．要求在[0，2π]上产生50个等距采样数据的一维数组，试用两种不同的指令实现。

解： y1=0:2*pi/49:2*pi y2=linspace(0,2*pi,50)
3.计算e -2t sint ，其中t 为[0，2π]上生成的10个等距采样的数组。

解：>> t=linspace(0,2*pi,10); x=exp(-2*t).*sin(t) x =
4.已知A=⎥⎦⎤⎢⎣⎡4321 , B=⎥
⎦⎤
⎢⎣⎡8765,计算矩阵A 、B 乘积和点乘. 解：>> A=[1,2;3,4];
B=[5,6;7,8]; x=A*B x =
19 22 43 50 >> x=A.*B x =
5 12
21 32
5.已知A=⎥⎦⎤⎢⎣⎡05314320，B=⎥
⎦⎤
⎢⎣⎡05314320，计算A&B, A|B, ~A, A==B, A>B. 解：>> A=[0,2,3,4;1,3,5,0];
B=[1,0,5,3;1,5,0,5]; a1=A&B a2=A|B
a3=~A
a4=(A==B)
a5=(A>B)
a1 =
0 0 1 1
1 1 0 0
a2 =
1 1 1 1
1 1 1 1
a3 =
1 0 0 0
0 0 0 1
a4 =
0 0 0 0
1 0 0 0
a5 =
0 1 0 1
0 0 1 0
7.将题5中的A阵用串转换函数转换为串B，再size指令查看A、B的结构，有何不同
解：>> A=[0,2,3,4;1,3,5,0]
B=num2str(A)
size(A)
size(B)
A =
0 2 3 4
1 3 5 0
B =
0 2 3 4
1 3 5 0
ans =
2 4
ans =
2 10
第三章习题
1.已知系统的响应函数为)sin(1
1)(θββ
ε+-
=-t e t y t ,其中
⎪⎪⎭
⎫
⎝⎛-=-=εεθεβ22
1arctan ,1 ,要求用不同线型或颜色,在同一张图上绘制ε取值分别为、、、时,系统在t ∈[0,18] 区间内的响应曲线,并要求用ε=和 ε=对他们相应的两条曲线进行文字标志。

解： clc
close all clear all t=0::18; xi=[,,,]';
sxi=sqrt(1-xi.^2); sita=atan(sxi./xi);
y=1-exp(-xi*t).*sin(sxi*t+sita*ones(1,901))./(sxi*ones(1,901))
plot(t,y(1), 'r-', t,y(2), ' b*', t,y(3), ' g+', t,y(4), ' k.')
text,,'\xi =') text,,'\xi=')
2.用plot3、mesh 、surf 指令绘制
()()2
22
2111
y x y x z +++
+-=
三维图(x,y 范围自定)。

解：
clc;close all ;clear all ; x=-5::5;y=-5::5; [X,Y]=meshgrid(x,y); a=sqrt((1-X).^2+Y.^2); b=sqrt((1+X).^2+Y.^2); Z=1./(a+b); a1=sqrt((1-x).^2+y.^2); b1=sqrt((1+x).^2+y.^2); z=1./(a1+b1);
subplot(1,3,1),plot3(x,y,z),xlabel('x'),ylabel('y'),zlabel('z');box on ; subplot(1,3,2),surf(X,Y,Z),xlabel('x'),ylabel('y'),zlabel('z');box
on ; subplot(1,3,3),mesh(X,Y,Z),xlabel('x'),ylabel('y'),zlabel('z');box on ;
3.对向量t 进行以下运算可以构成三个坐标的值向量:x=sin(t),y=cos(t),z=t.利用指令plot3,并选用绿色的实线绘制相应的三维曲线. 解：
t=(0::2)*pi;
x=sin(t);
y=cos(t);
z=t;
plot3(x,y,z,'b-');box on
1.请分别用for 和while 循环语句计算K=∑=63
02i i 的程序，再写出一种避免循环的
计算程序。

(提示：可考虑利用MATLAB 的sum(X,n)函数，实现沿数组X 的第n 维求和。

) 解：
1）K=0; for i=0:63; K=K+2^i; end K
K =+019 2）i=0;K=0; while i<=63; K=K+2^i; i=i+1; end; K
K =+019 3）i=0; X=0:63; for i=0:63;
X(i+1)=2^i; end
sum(X,2) ans =+019
1.将下列系统的传递函数模型用MATLAB 语言表达出来。

(1))170046842541254289()
1700109329135()(23452341+++++++++=s s s s s s s s s s G 解：
num=[1,35,291,1093,1700];
den=[1,289,254,2541,4684,1700]; sys=tf(num,den) (2))
15).(5).(1()
3(15)(2++++=
s s s s s G 解：
z=-3;
p=[-1,-5,-15]; k=15;
sys=zpk(z,p,k)
(3))252).(1).(1()
23.()2.(.100)(23223+++-++++=s s s s s s s s s s G 解：
z=[0,-2,-2]; p=[-1,1]; k=100;
sys1=zpk(z,p,k); num=[1,3,2]; den=[1,2,5,2]; sys2=tf(num,den); sys=series(sys1,sys2)
4.求题3中的系统模型的等效传递函数模型和零极点模型。

解：
A=[3,2,1;0,4,6;0,-3,-5]; B=[1,2,3]' ; C=[1,2,5]; D=0;
sys=ss(A,B,C,D); systf=tf(sys) syszpk=zpk(sys)
Transfer function: 20 s^2 - 83 s + 138 --------------------- s^3 - 2 s^2 - 5 s + 6 Zero/pole/gain: 20 (s^2 - + ----------------------- (s-3) (s-1) (s+2)
5.已知系统的动力学方程如下，试用MATLAB 语言写出它们的传递函数。

(1))(2)()(500)(50)(15)(.
.....)3(t r t r t y t y t y t y +=+++ 解：
num=[1,2,0];
den=[1,15,50,500]; sys=tf(num,den) Transfer function: s^2 + 2 s -------------------------
s^3 + 15 s^2 + 50 s + 500
(2) )(4)(4)(6)(3)(.
..
t r dt t y t y t y t y =+++⎰ 解：
num=[4,0]; den=[1,3,6,4]; sys=tf(num,den) Transfer function: 4 s --------------------- s^3 + 3 s^2 + 6 s + 4
6.试用MATLAB 语言表示图5-13所示系统。

当分别以y =x 2和f 为系统输出、输
入时的传递函数模型和状态空间模型(图中k =7N/m,c 1=, c 2=,m 1=3.5kg, m 2=5.6kg)。

f
解：)(t
k=7;
c1=;
c2=;
m1=;
m2=;
num=[m1,c1,k];
den=[m1*m2,c1*m1+c2*m1+c1*m2,c1*c2+m2*k,c1*k+c2*k,0];
sys=tf(num,den)
Transfer function:
s^2 + s + 7
--------------------------------------
s^4 + s^3 + s^2 + s
7.试用MATLAB语言分别表示图5-14所示系统质量m1,m2的位移x1,x2对输入f 的传递函数X2(s)/F(s)和X1(s)/F(s),其中m1=12kg, m2=38kg,k=1000N/m, c=。

解：
m1=12;
m2=38;
k=1000;
c=;
num=[c,k];
den=[m1*m2,m1*c+m2*c,m1*k+m2*k,0,0];
sys1=tf(num,den)
num=[m1,c,k];
den=[m1*m2,m1*c+m2*c,m1*k+m2*k,0,0];
sys2=tf(num,den)
Transfer function:
s + 1000
---------------------------
456 s^4 + 5 s^3 + 50000 s^2
Transfer function:
12 s^2 + s + 1000
---------------------------
456 s^4 + 5 s^3 + 50000 s^2
补充题
求图示传递函数
sys1=tf([1,2],[1,3,4]);
sys2=tf([1,4,5] ,[1,6,7,8]);
sys3=tf([1,0],[1,2]);
sys4=tf([1],[1,3]);
sys5=parallel(sys3,sys4);
sys=feedback(sys1*sys2*sys5,1,-1)
结果
s^5 + 10 s^4 + 39 s^3 + 74 s^2 + 66 s + 20
-----------------------------------------------------------------
s^7 + 14 s^6 + 81 s^5 + 262 s^4 + 530 s^3 + 684 s^2 + 538 s + 212
第六章习题2.将例6-2中的微分方程改写为以下形式：
1
)0(
,0
)0(0
.) 1.(
..
2
..
=
==
+
-
-
y
y y
y
y
yμ
求μ分别为1、2时，在时间区间t=[0,20]微分方程的解。

解：
M函数文件
function dx=wffc(t,x,flag,ps)
dx=zeros(2,1);
dx(1)=x(2);
dx(2)=ps*(1-x(1)^2)*x(2)-x(1);
调用程序
clc;close all;clear all;
tspan=[0,20];
x0=[0,1];
ps=1;
[T1,X1]=ode45('wffc',tspan,x0,odeset,ps);
ps=2;
[T2,X2]=ode45('wffc',tspan,x0,odeset,ps);
plot(T1,X1(:,1),'r',T2,X2(:,1),'b-.')
X1(:,1)
X2(:,1)
3.对图6-18所示反馈系统进行单位阶跃响应和方波响应(方波周期为30s)仿真。

要求：
(1)利用MATLAB模型连接函数求出系统闭环传递函数。

(2)利用step 函数求单位阶跃响应。

(3)利用gensig 函数产生方波信号，利用lsim 函数求方波响应。

解：
clc;close all ;clear all ; % (1)
sys1=tf([1,],[1,]);
sys2=ZPK([],[0,-2,-10],20); sys3=series(sys1,sys2); sys4=feedback(sys3,1,-1); % (2)
subplot(1,2,1) step(sys4); % (3)
[u,t]=gensig('square',30,60); subplot(1,2,2) lsim(sys4,'r',u,t)
20 (s+
-------------------------------------------- (s+ (s+ (s^2 + +
4.已知系统传递函数01
.12.01
)(2++=
s s s G ;
(1)绘制系统阶跃响应曲线。

(2)绘出离散化系统阶跃响应曲线，采样周期T s =。

解：
clc;close all ;clear all ; % (1)
sys=tf([1],[1,,]); subplot(1,2,1)
step(sys) % (2)
sys=tf([1],[1,,]); sys1=c2d(sys,,'zoh');
[num,den]=tfdata(sys1,'v'); subplot(1,2,2) dstep(num,den)
附加题
1、已知二阶微分方程0342=-+-y y y y y ，其初始条件为0)0(=y ，1)0(=y ，求在时间范围t=[0 5]内该微分方程的解。

M 函数为：
function dy=vdp(t,y)
dy=zeros(2,1); dy(1)= y(2); dy(2)= 4*y(2)-(y(1)^2)*y(2)+3*y(1); 调用函数为：
[T,Y]=ode45('vdp',[0 5],[0,1]); plot(T,Y(:,1),'r-',T,Y(:,2),'b:')
2、已知系统模型为
7
22)(3
+++=s s s s G ，计算系统在周期10s 的方波信号作用下5
个周期内的时间响应，并在同一图形窗口中绘制输入信号和时间响应曲线。

sys=tf([1,2],[1,0,2,7]);
[u,t]=gensig('square',10,50); %产生方波信号数据
lsim(sys,'r',u,t) , hold on %产生方波响应并绘曲线 plot(t,u,'-.') %在同一坐标系绘方波波形 hold off
第七章习题
1.绘制下列各单位反馈系统开环传递函数的Bode 图和Nyquist 图，并根据其稳定裕度判断系统的稳定性。

(1))
31).(21).(1(10
)(s s s s G k +++=
解：
clc;clear all ;close all ; % (1)
Gk=zpk([],[0,,-1/3],5/3); subplot(1,2,1) margin(Gk) grid on
subplot(1,2,2) nyquist(Gk)
由上图的稳定裕度知系统临界稳定。

(2))
101).(1.(10
)(s s s s G k ++=
解：
clc;clear all ;close all ; % (2)
Gk=zpk([],[0,-1,],1); subplot(1,2,1) margin(Gk) grid on
subplot(1,2,2) nyquist(Gk)
由上图的稳定裕度知系统不稳定。

(3))
2.01).(1.01.(10
)(2s s s s G k ++=
解：
clc;clear all ;close all ; % (3)
Gk=zpk([],[0,0,-10,-5],500); subplot(1,2,1) margin(Gk) grid on
subplot(1,2,2) nyquist(Gk)
由上图的稳定裕度知系统不稳定。

(4) )
101).(1.01.(2
)(2s s s s G k ++=
解：
clc;clear all;close all; % (4)
Gk=zpk([],[0,0,-10,],2); subplot(1,2,1) margin(Gk) grid on
subplot(1,2,2) nyquist(Gk)
由上图的稳定裕度知系统不稳定。

2.设单位反馈系统的开环传递函数为
)
12.()(22++=
n n
k w s w s s K
s G ξ,其中无阻尼固有频率w n =90rad/s ，阻尼比ξ=,试确定使系统稳定的K 的范围。

解： 方法1
g=tf(1,[1/90^2 90 1 0]);%系统开环模型
w=logspace(0,3,1000); %生成频率向量 bode(g,w)
[mag,phase,w]=bode(g,w); %产生幅值（非分贝）和相位向量
mag1=reshape(mag,1000,1); %重构幅值向量（1000*1）
phase1=reshape(phase,1000,1);%重构相频向量（1000*1）
wc=interp1(phase1,w,-180) %插值求-180度所对应的频率——wc
gk=interp1(w,mag1,wc) %插值求wc所对应的增益
gkk=1/gk %该增益的倒数即为可增加的最大增益
wc =
gk =
gkk =
方法2
wc=0;wg=;k=1;
while wc<wg
sys=tf(k,[1/(90*90),2*90,1,0]);
[gm,pn,wg,wc]=margin(sys);
k=k+;
end
ans =
方法3
xi=;omega=90;w=90;
sys1=tf(1,[1,0]);
sys2=tf(1,[1/w^2,2*xi/w,1]);
sys=series(sys1,sys2);
[Gm,Pm,Wcg,Wcp]=margin(sys);
k=Gm
k =
36
3.设系统结构如图7-22所示，试用LTI Viewer分析系统的稳定性,并求出系统的稳定裕度及单位阶跃响应峰值。

clc;close all;clear all;
G11=;
G12=zpk([0],[],1);
G1=G11-G12;
G2=tf(1,[1 2 0]);
Gk=G1*G2;
Gb=feedback(Gk,1,-1);
[Gm,Pm,Wcg,Wcp]=margin(Gb)
step(Gb)
[y,t]=step(Gb);
[yp,k]=max(y)
yp
Gm =
Pm =
yp =
4.设闭环离散系统结构如图7-23所示，其中G(s)=10/(s.(s+1))，H(s)=1,绘制T=、1s时离散系统开环传递函数的Bode图和Nyquist图,以及系统的单位阶跃响应曲线。

解：
clc;close all;clear all;
ts=,ts1=1;
Gk=zpk([],[0,-1],10);
Gz1=c2d(Gk,ts,'zoh');
Gz2=c2d(Gk,ts1,'zoh');
[num1,den1,ts]=tfdata(Gz1,'v');
[num2,den2,ts1]=tfdata(Gz2,'v');
figure(1)
subplot(1,3,1)
dbode(num1,den1,ts);
grid
subplot(1,3,2)
dnyquist(num1,den1,ts);
subplot(1,3,3)
dstep(num1,den1)
figure(2)
subplot(1,3,1)
dbode(num2,den2,ts1);
grid
subplot(1,3,2)
dnyquist(num2,den2,ts1);
subplot(1,3,3)
dstep(num2,den2)
第九章习题
3.构建图9-63所示的仿真模型。

图中的PID模块为图9-39所示的积分可分离式PID子系统，取kp=5,kd=,ki=5,分别取delta为、时比较系统的单位阶跃响应性能。

解：
% delta=
delta=;
kp=5;
kd=;
ki=5;
% delta=
delta=;
kp=5;
kd=;
ki=5;
第九章补充习题
2、建立下图所示系统的动力学方程，并绘制用于模拟单位阶跃输入作用下系统响应的simulink模型，系统输出y接至示波器，并给出仿真结果，仿真时间为10s。

解：
对如图所示的系统进行受力分析有：
y c
t
y
m=
+
+
t
t
)(
)(
)(
ky
f
)(t
Simulink仿真模型为：
仿真结果略
3、建立下图所示微分方程所对应的simulink模型(u为输入、y为输出)，u为单位斜坡输入作用下系统响应的输出y接至示波器，并给出仿真结果，仿真时
间为10s。

6
)(
+
+
)(2=
3
y
u y3
t
y
t
解：
仿真结果略。