系统建模与仿真习题6及答案.docx

1、系统建模与仿真习题6及答案系统建模与仿真习题六及答案1.一个离散的系统：分别利用dimpulse（）函数、filter（）函数、impz（）求解系统的脉冲响应。解：clc;clear;b=1 2 0 1;a=1 -0.5 0.25u=1,zeros(1,9)num=1 2 0 1;den=1 -0.5 0.25 0; %补零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

2、);stem(y)title(impz求解)2. 假设信号：分别以采样频率为2000Hz、1500Hz、1200Hz三种情况对该信号进行采样。然后用理想的低通滤波器恢复该连续信号，观察恢复波形与实际波形的差异，分别计算出最大恢复误差。解：clear;clc;for k=1:3 if k=1 Fs=2000; elseif k=2 Fs=1500;else Fs=1200; endt1=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

3、.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,.);hold onplot(t,x0,r,linewidth,1.5)wucha(k)=max(abs(xa-x0);endwuchawucha = 0.3993 0.2847 1.52713.系统传递函数为：假设系统输入为，请用卷积、lsim函数方法求其零状态响应的结果。解：clc;clear; b=2 1;a=3 2 1;t=0:0.001:10y1=impuls

4、e(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(卷积计算的结果)gridsubplot(212)plot(t,Y2)title(直接求零状态响应的结果)grid 4. 已知微分方程（1）用Euler方法求解常微分方程初值问题，并将数值解和该问题的解析解（）比较。（2）用二阶Runge-Kutta方法求解常微分方程初值问题，并将数值解和该问题的解析解（）比较。（3）利用四阶Runge-Kutta方法编程

5、仿真；（4）利用ode45函数求解并仿真解：（1）clc;clear;y(1)=1;h=0.1x=0:h:20;for n=1:length(x)-1 xn=x(n); yn=y(n); y(n+1)=yn+h*sin(xn); endx0=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;for n=1:length(t)-1 % k1=sin(t(n);%该语句不需要 k2=sin(t(n)+h/2); y(n+1)=y(n)+h*k2;endy1=2-co

6、s(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);for i = 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 +

7、K4); t0=t1; y0=y; yy1(i)=y;end plot (t, 1,yy1); 结果：（4）子程序：Runge.mfunction dy=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) 建立火箭升

8、空过程的数学模型（微分方程）；(2) 分别利用四阶Runge-Kutta编程方法， ode45函数方法求解引擎关闭瞬间火箭的高度、速度及火箭到达最高点的高度。解：（1）建立模型：时刻的火箭高度； =30000（牛顿）火箭推力，当时间40秒时，T=0；火箭飞行过程中的质量， 40秒时，千克=900（千克）火箭初始质量； =600（千克）燃料质量；=15（千克/秒）燃料的燃烧速率； =0.4（千克/米）空气阻力系数；=9.8（米/秒）重力加速度由牛顿第二定律，可得到火箭飞行过程的方程：这是一个初值问题，初始条件为 设，则问题化为求下列微分方程组的初值问题：关闭引擎时（秒），所求的是此时火箭的高度，

9、速度，以及火箭到最高点的高度。（2）四阶Runge-Kutta编程方法：子程序：function y=huojian(t,x)global m0 k g Tg=9.8;m0=900;T=30000;k=0.4;m=m0-15*t;if t40 T=0; m=300;endy=x(2);-(k/m)*x(2)2+T/m-g;主程序：clear;clc;global m0 k g Ty0=0;0;h=0.01;t = 0: h :55;t0=0;N = length (t);for i = 1 : N t1 = t0 + h; K1 = huojian(t0, y0); K2 = huojian(

10、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; endsubplot(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

11、) 结果：x40 = 8.3254e+003v40 = 257.8302xmax = 9.1906e+003Ode45方法：子程序：huojian.mfunction y=huojian(t,x)k=0.4;g=9.8;m0=900;T=30000;m=m0-15*t;if t40 T=0; m=300;endy=x(2);-(k/m)*x(2)2+T/m-g;主程序：feixing.mclear;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+003v40 = 260.2437xmax = 9.2151e+003