系统建模与仿真习题6及答案.docx
《系统建模与仿真习题6及答案.docx》由会员分享，可在线阅读，更多相关《系统建模与仿真习题6及答案.docx（12页珍藏版）》请在冰豆网上搜索。
系统建模与仿真习题6及答案
系统建模与仿真习题六及答案
1.一个离散的系统：
分别利用dimpulse（）函数、filter（）函数、impz（）求解系统的脉冲响应。
解：
clc;clear;
b=[1201];
a=[1-0.50.25]
u=[1,zeros（1,9）]
num=[1201];
den=[1-0.50.250];%补零
subplot（311）
[y,x]=dimpulse（num,den,10）;
stem（y）
title（'dimpulse求解'）
subplot（312）
y=filter（b,a,u）;
stem（y）
title（'filter求解'）
subplot（313）
y=impz（b,a,10）;
stem（y）
title（'impz求解'）
2.假设信号：
分别以采样频率为2000Hz、1500Hz、1200Hz三种情况对该信号进行采样。
然后用理想的低通滤波器恢复该连续信号，观察恢复波形与实际波形的差异，分别计算出最大恢复误差。
解：
clear;clc;
fork=1:
3
ifk==1
Fs=2000;
elseifk==2
Fs=1500;
else
Fs=1200;
end
t1=0:
1/Fs:
0.002;%采样的时间间隔
n=0:
0.002*Fs;%采样的点数
x=2*sin（1000*pi*t1）;
T=1/Fs;
dt=T/10;%每个采样周期取多个样点，保证g（t）的连续性
t=0:
dt:
0.002;
M=ones（length（n）,1）*t-n'*T*ones（1,length（t））;
xa=x*sinc（Fs*M）;
x0=2*sin（1000*pi*t）;%实际函数
subplot（3,1,k）;
stem（t,xa,'.'）;
holdon
plot（t,x0,'r','linewidth',1.5）
wucha（k）=max（abs（xa-x0））;
end
wucha
wucha=
0.39930.28471.5271
3.系统传递函数为：
假设系统输入为
，请用卷积、lsim函数方法求其零状态响应的结果。
解：
clc;clear;
b=[21];a=[321];
t=0:
0.001:
10
y1=impulse（b,a,t）;%求脉冲响应
u=1+sin（2*t）;
Y1=0.001*conv（y1,u）;%利用卷积求响应输出
Y2=lsim（b,a,u,t）%直接求零状态响应
subplot（211）
plot（t,Y1（1:
length（t）））
title（'卷积计算的结果'）
grid
subplot（212）
plot（t,Y2）
title（'直接求零状态响应的结果'）
grid
4.已知微分方程
（1）用Euler方法求解常微分方程初值问题，并将数值解和该问题的解析解（
）比较。
（2）用二阶Runge-Kutta方法求解常微分方程初值问题，并将数值解和该问题的解析解（
）比较。
（3）利用四阶Runge-Kutta方法编程仿真；
（4）利用ode45函数求解并仿真
解：
（1）
clc;clear;
y
（1）=1;
h=0.1
x=0:
h:
20;
forn=1:
length（x）-1
xn=x（n）;
yn=y（n）;
y（n+1）=yn+h*sin（xn）;
end
x0=0:
h:
20;
y0=-cos（x0）+2;
plot（x0,y0,'bo',x,y,'r*'）
legend（'解析解','数值解'）
（2）
clc;clear;
y
（1）=1;
h=0.1;
t=0:
h:
20;
forn=1:
length（t）-1
%k1=sin（t（n））;%该语句不需要
k2=sin（t（n）+h/2）;
y（n+1）=y（n）+h*k2;
end
y1=2-cos（t）;%精确解
subplot（211）
plot（t,y）
subplot（212）
plot（t,y1）
（3）
clear;clc;
t0=0;tN=20;y0=1;h=0.1;
t=t0:
h:
tN;
N=length（t）;
fori=1:
N-1
t1=t0+h;
K1=Runge（t0,y0）;
K2=Runge（t0+h/2,y0+h*K1/2）;
K3=Runge（t0+h/2,y0+h*K2/2）;
K4=Runge（t0+h,y0+h*K3）;
y=y0+（h/6）*（K1+2*K2+2*K3+K4）;
t0=t1;
y0=y;
yy1（i）=y;
end
plot（t,[1,yy1]）;
结果：
（4）
子程序：
Runge.m
functiondy=Runge（t,y）
dy=sin（t）;
主程序：
clc;clear;
y0=1;
[t,y]=ode45（'Runge',[0,20],y0）;
plot（t,y）
结果：
5.小型火箭初始质量为900千克，其中包括600千克燃料。
火箭竖直向上发射时燃料以15千克/秒的速率燃烧掉，由此产生30000牛顿的恒定推力。
当燃料用尽时引擎关闭。
设火箭上升的整个过程中，空气阻力与速度平方成正比，比例系数为0.4（千克/米）。
重力加速度取9.8米/秒2.
（1）建立火箭升空过程的数学模型（微分方程）；
（2）分别利用四阶Runge-Kutta编程方法，ode45函数方法求解引擎关闭瞬间火箭的高度、速度及火箭到达最高点的高度。
解：
（1）
建立模型：
——
时刻的火箭高度；
=30000（牛顿）——火箭推力，当时间
>40秒时，T=0；
——火箭飞行过程中的质量，
>40秒时，
千克
=900（千克）——火箭初始质量；
=600（千克）——燃料质量；
=15（千克/秒）——燃料的燃烧速率；
=0.4（千克/米）——空气阻力系数；
=9.8（米/秒
）——重力加速度
由牛顿第二定律，可得到火箭飞行过程的方程：
这是一个初值问题，初始条件为
设
，则问题化为求下列微分方程组的初值问题：
关闭引擎时
（秒），所求的是此时火箭的高度，速度，以及火箭到最高点的高度。
（2）
四阶Runge-Kutta编程方法：
子程序：
functiony=huojian（t,x）
globalm0kgT
g=9.8;m0=900;T=30000;k=0.4;
m=m0-15*t;
ift>40
T=0;
m=300;
end
y=[x
（2）;-（k/m）*x
（2）^2+T/m-g];
主程序：
clear;clc;
globalm0kgT
y0=[0;0];h=0.01;
t=0:
h:
55;
t0=0;
N=length（t）;
fori=1:
N
t1=t0+h;
K1=huojian（t0,y0）;
K2=huojian（t0+h/2,y0+h*K1/2）;
K3=huojian（t0+h/2,y0+h*K2/2）;
K4=huojian（t0+h,y0+h*K3）;
y=y0+（h/6）*（K1+2*K2+2*K3+K4）;
yy1（i）=y
（1）;
yy2（i）=y
（2）;
t0=t1;
y0=y;
end
subplot（211）
plot（t,yy1）
subplot（212）
plot（t,yy2）
a=[t',yy1',yy2'];
x40=a（40/h+1,2）%燃料用尽时的高度
v40=a（40/h+1,3）%燃料用尽时的速度
xmax=max（yy1）
结果：
x40=
8.3254e+003
v40=
257.8302
xmax=
9.1906e+003
Ode45方法：
子程序：
huojian.m
functiony=huojian（t,x）
k=0.4;g=9.8;m0=900;T=30000;m=m0-15*t;
ift>40
T=0;
m=300;
end
y=[x
（2）;-（k/m）*x
（2）^2+T/m-g];
主程序：
feixing.m
clear;clc;
k=0.4;g=9.8;m0=900;T=30000;
x0=[0,0];
h=0.01;
ts=0:
h:
55;
[t,x]=ode45（'huojian',ts,x0）;
a=[t,x];
x40=a（40/h+1,2）%燃料用尽时的高度
v40=a（40/h+1,3）%燃料用尽时的速度
xmax=max（x（:
1））%火箭到达最高点的高度
subplot（2,1,1）,plot（t,x（:
1））,title（'altitude'）
subplot（2,1,2）,plot（t,x（:
2））,title（'speed'）
结果：
x40=
8.3234e+003
v40=
260.2437
xmax=
9.2151e+003
升级会员 