系统仿真实验报告Word文件下载.docx
《系统仿真实验报告Word文件下载.docx》由会员分享，可在线阅读，更多相关《系统仿真实验报告Word文件下载.docx（30页珍藏版）》请在冰豆网上搜索。
B3x372double
ans3x372double
th1x21168double
x1x21336doublecomplex
10.size（A）
33
length（A）
3
实验二
1、二维直角坐标
t=-pi:
0.3:
pi;
y=1./（1+exp（-t））;
subplot（221）,plot（t,y,'
b'
）;
xlabel（'
t'
title（'
plot（t,y）'
subplot（222）,stem（t,y,'
g'
subplot（223）,semilogy（t,y,'
m'
subplot（224）,stairs（t,y,'
k'
）
实验三
functionf=SY3（）
while1
b=input（'
输入一个自然数:
'
ifb<
2
fprintf（'
%d既不是质数也不是合数'
b）;
elseifisprime（b）==1
%d是一个素数'
else
%d是一个合数，可以因式分解为:
b）
\n'
forn=2:
sqrt（b）
ifrem（b,n）==0
num1=n;
num2=b/n;
%d=%d×
%d'
b,num1,num2）
end
y=input（'
继续判断？
继续按1，退出按0:
ify==0
break;
end
命令窗口输入SY3
78
78是一个合数，可以因式分解为:
78=2×
39
78=3×
26
78=6×
13
实验四
1．求矩阵对应的行列式和特征根
a=sym（'
[a11a12;
a21a22]'
da=det（a）
ea=eig（a）
da=
a11*a22-a12*a21
ea=
a11/2+a22/2-（a11^2-2*a11*a22+a22^2+4*a12*a21）^（1/2）/2
a11/2+a22/2+（a11^2-2*a11*a22+a22^2+4*a12*a21）^（1/2）/2
2．求方程的解（包括精确解和一定精度的解）
r1=solve（'
x^2-x-1'
rv=vpa（r1）
rv4=vpa（r1,4）
rv20=vpa（r1,20）
r1=
5^（1/2）/2+1/2
1/2-5^（1/2）/2
rv=
1.6180339887498948482045868343656
-0.61803398874989484820458683436564
rv4=
1.618
-0.618
rv20=
1.6180339887498948482
-0.6180339887498948482
3．a=sym（'
a'
b=sym（'
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
a*d-b*c
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分解因式
5．对方程AX=b求解
A=[34,8,4;
3,34,3;
3,6,8];
b=[4;
6;
2];
X=linsolve（A,b）%调用linsolve函数求解
X=
138/2045
66/409
212/2045
A\b%用另一种方法求解
6．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=
exp（x）/（2*（exp（x）+1）^（1/2））
-2*sin（x）-x*cos（x）
x*sin（x）-3*cos（x）
-（b*cos（t））/（a*sin（t））
结果：
实验五
1、求下面系统的单位阶跃响应
%程序如下：
num=[5];
den=[1,1,4];
（分子4改为5）
step（num,den）
[y,x,t]=step（num,den）;
tp=spline（y,t,max（y））%计算峰值时间
max（y）%计算峰值
tp=
1.6579
max（y）%计算峰值
1.8040
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）
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）;
impulseresponse（k=55）'
figure（3）
k=56;
impulseresponse（k=56）'
8.作如下系统的bode图
n=[1,1];
d=[1,4,11,7];
bode（n,d）
9.系统传函如下
求有理传函的频率响应，然后在同一张图上绘出以四阶伯德近似表示的系统频率响应
num=[1];
den=conv（[12],conv（[12],[12]））;
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--'
gridon;
title（'
bodeplot'
frequency'
ylabel（'
gain'
subplot（2,1,2）;
semilogx（w,p1,w,p2,'
gridon;
phase'
11.二阶系统为：
令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）;
holdon
nyquist（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）;
%自动变量
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
nyquist（num,den）
figure
（2）
[num2,den2]=cloop（num,den）;
impulse（num2,den2）;
15.已知系统为：
作该系统的nichols曲线。
d=[1,1,0];
ngrid（‘new’）;
nichols（n,d）;
16.已知系统的开环传递函数为：
当k=2时，分别作nichols曲线和波特图。
num=1;
den=conv（conv（[10],[11]）,[0.51]）;
subplot（1,2,1）;
nichols（num,den）;
grid;
%nichols曲线
subplot（1,2,2）;
g=tf（num,den）;
bode（feedback（g,1,-1））;
%波特图
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）
-2.5214+0.0000i
-0.2393+0.8579i
-0.2393-0.8579i
-3.3089+0.0000i
0.1545+1.7316i
0.1545-1.7316i
18.系统的状态方程为：
试确定系统的稳定性。
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））
0
-1
-3
实验六
1.取h=0.2，试分别用欧拉法、RK2法和RK4法求解微分方程的数值解，并比较计算精度。
注：
解析解：
%程序如下
clear
t
（1）=0;
%置自变量初值
y
（1）=1;
y_euler
（1）=1;
y_rk2
（1）=1;
y_rk4
（1）=1;
%置解析解和数值解的初值
h=0.2;
%步长
%求解析解
fork=1:
5
t（k+1）=t（k）+h;
y（k+1）=sqrt（1+2*t（k+1））;
%利用欧拉法求解
y_euler（k+1）=y_euler（k）+h*（y_euler（k）-2*t（k）/y_euler（k））;
%利用RK2法求解
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;
%利用RK4法求解
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;
%输出结果
disp（'
时间解析解欧拉法RK2法RK4法'
yt=[t'
y'
y_euler'
y_rk2'
y_rk4'
];
disp（yt）
程序运行结果如下：
时间解析解欧拉法RK2法RK4法
01.00001.00001.00001.0000
0.20001.18321.20001.18671.1832
0.40001.34161.37331.34831.3417
0.60001.48321.53151.49371.4833
0.80001.61251.68111.62791.6125
1.00001.73211.82691.75421.7321
2.考虑如下二阶系统：
在
上的数字仿真解（已知：
，
），
并将不同步长下的仿真结果与解析解进行精度比较。
说明：
已知该微分方程的解析解分别为：
采用RK4法进行计算，选择状态变量：
则有如下状态空间模型及初值条件
采用RK4法进行计算。
程序如下：
h=input（'
请输入步长h='
%输入步长
M=round（10/h）;
%置总计算步数
%置自变量初值
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）;
%置数值解的初值
M
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+1））+100*sqrt（3）/3*exp（-t（k+1）/2）.*sin（sqrt（3）/2*t（k+1））;
%R=0
k11=x2（k）;
k12=-x1（k）;
k21=x2（k）+h*k12/2;
k22=-（x1（k）+h*k11/2）;
k31=x2（k）+h*k22/2;
k32=-（x1（k）+h*k21/2）;
k41=x2（k）+h*k32;
k42=-（x1（k）+h*k31）;
x1（k+1）=x1（k）+h*（k11+2*k21+2*k31+k41）/6;
x2（k+1）=x2（k）+h*（k12+2*k22+2*k32+k42）/6;
y_rk4_0（k+1）=x1（k+1）;
%R=0.5
k12=-x1（k）-x2（k）;
k22=-（x1（k）+h*k11/2）-（x2（k）+h*k12/2）;
k32=-（x1（k）+h*k21/2）-（x2（k）+h*k22/2）;
k42=-（x1（k）+h*k31）-（x2（k）+h*k32）;
y_rk4_05（k+1）=x1（k+1）;
%求出误差最大值
err_0=max（abs（y_0-y_rk4_0））;
err_05=max（abs（y_05-y_rk4_05））;
最大误差（R=0）最大误差（R=0.5）'
err_max=[err_0,err_05];
disp（err_max）
步长h取0.0001到0.5之间变化，并且将数值解直接与解析解比较，求出最大误差数据列入下表中。
步长h
0.0001
0.0005
0.001
0.005
0.01
0.05
0.1
0.5
R=0
最大误差
5.4330×
10-10
1.6969×
1.0574×
4.1107×
10-9
6.6029×
10-8
4.1439×
10-5
6.6602×
10-4
4.2988×
10-1
R=0.5
2.7649×
10-11
6.8123×
10-12
5.3753×
4.0902×
6.5425×
4.1365×
10-6
6.7152×
4.5976×
10-2
从上表中可以看出，当步长h=0.001时，总误差最小；
当步长h小于0.001时，由于舍入误差变大而使总误差增加；
当步长h大于0.001时，则由于截断误差的增加也使得总误差加大。
另外，当系统的解变化激烈时（如R=0），误差对步长的变化较为敏感；
当系统的解变化平稳时，步长的变化对误差的影响就要缓和得多。
数值积分算法确定以后，在选择步长时，需要综合考虑。
升级会员 