《MATLAB及应用》上机作业
学院名称： 机械工程学院 专业班级： 测控1201 学生姓名： 学生学号：
201 年 4 月
1
《MATLAB及应用》上机作业要求及规范
一、作业提交方式：word文档打印后提交。 二、作业要求：
1．封面：按要求填写学院、班级、姓名、学号，不要改变封面原有字体及大小。 2．内容：只需解答过程（结果为图形输出的可加上图形输出结果），不需原题目；为便于批阅，题与题之间应空出一行；每题答案只需直接将调试正确后的M文件内容复制到word中（不要更改字体及大小），如下所示：
%作业1_1 clc
A=[1 2 3 4;2 3 5 7;1 3 5 7;3 2 3 9;1 8 9 4];
B=[1+4*i 4 3 6 7 7;2 3 3 5 5 4+2*i;2 6+7*i 5 3 4 2;1 8 9 5 4 3]; C=A*B D=C(4:5,4:6)
三、大作业评分标准：
1．提交的打印文档是否符合要求； 2．作业题的解答过程是否完整和正确；
3．答辩过程中阐述是否清楚，问题是否回答正确；
4．作业应独立完成，严禁直接拷贝别人的电子文档，发现雷同者都以无成绩论处。
2
作业1
1、用MATLAB可以识别的格式输入下面两个矩阵
?1??2A??1??3?1?234?4?1?4i?357??23357?，B???26?7i?239???18??894?7??3554?2i? ?5342?9543??367再求出它们的乘积矩阵C，并将C矩阵的右下角2?3子矩阵赋给D矩阵。赋值完成后，调用相应的命令查看MATLAB工作空间的占有情况。 解：
A=[1 2 3 4;2 3 5 7;1 3 5 7 ;3 2 3 9 ;1 8 9 4;]
B=[1+4i 4 3 6 7 7;2 3 3 5 5 4+2i;2 6+7i 5 3 4 2;1 8 9 5 4 3;] B=[1+4i 4 3 6 7 7;2 3 3 5 5 4+2i;2 6+7i 5 3 4 2;1 8 9 5 4 3;] C=A*B D=C(4:5,4:6);
?162313???511108?，求A，A?1，A3，2A?A?1，3A?1?A'，并求矩阵A的2、设矩阵??97612???414152??特征值和特征向量。 解：
A=[16 2 3 13;5 11 10 8;9 7 6 12;4 14 15 2;] det(A) inv(A)
3
A.^3 2*A+inv(A) 3*A-A' [V,D]=eig(A) abs (A)
3、解下列矩阵方程:
?010??100??1?43??????? 100X001?20?1
???????001??010??1?20???????解：A=[0 1 0;1 0 0;0 0 1;]; B=[1 0 0;0 0 1;0 1 0;]; C=[1 -4 3; 2 0 -1;1 -2 0;]; X=inv(A)*C*inv(B)
4、求多项式f(x)?(x?5x?35x?13)(x?45x?3)(31x?3)当x?3时的值，并求
322f(x)的导数。
解：a1=[1 5 35 13]
a2=[0 1 45 3] a3=[0 0 31 3] p1=conv(a1,a2) p1=conv(p1,a3) polyval(p1,3) polyder(p1)
4
5、求多项式f?x?5x?14x?10x?3的根和导数。 解:
432p1=[1 -5 -14 -10 -3] roots(p1) polyder(p1)
10(s4?4s2?5s2?6s?7)6、对于有理多项式 f(s)? (1)计算该多项式相除的结果;(2)
(s?1)(s?2)(s?3)将该多项式展开为部分分式的形式;(3)计算
df。 ds解(1)A=[1 0 4 5 6 7]; B=A.*10; C=[1 1]; D=[1 2]; E=[1 3];
F=conv(conv(C,D),E); G=deconv(B,F); (2)[r,p,k]=residue(B,F) (3)[p,q]=ployder(B,F)
7、在某次传感器输入输出特性实验中测得输入输出的一组数据如下表所示:
x（输入） 1 1。3 2 1。8 3 2。2 4 2。9 5 3。5 y（输出） 5
已知输入x和输出y可以近似成线性关系，即y?kx?b，求系数k和b，并求当输入x?8时输出y的值。 解x=[1 2 3 4 5]; y=[1.3 1.8 2.2 2.9 3.5]; p=polyfit(x,y,1) a=polyval(p,8)
8、根据人口理论的马尔萨斯模型可知，人口总数N(t)可以采用指数函数N(t)= e 口数据进行拟合。据统计，六十年代世界人口数据如下(单位：亿)
a + b t对人
t y 1960 1961 1962 1963 1964 1965 1966 1967 1968 3。3918 3。4213 3。4503 3。4698 3。4763 3。4920 3。5133 3。5322 3。5505 试求马尔萨斯模型中的a，b值，并画出拟合曲线图，同时预测一下2010年的人口数值。 解year=[1960:1:1968];
n=[3.3918 3.4213 3.4503 3.4698 3.4763 3.4920 3.5133 3.5322 3.5505]; y=log(n);
p=polyfit(year,y,1); plot(year,y,'g*'); a=p(2) b=p(1)
y2010=exp(polyval(p,2010))
6
9、某实验测得强度随时间变化的一组数据:
t y 0 0 0。5 0。4794 1 0。8415 1。5 0。9915 2 0。9093 2。5 0。5985 3 0。1411 1)利用二次曲线拟合求出t?2.25秒处强度指标y1。
2)利用样条曲线插值求出t?2.25秒处强度指标y1。 解：
t=[0 0.5 1 1.5 2 2.5 3]
y=[0 0.4794 0.8415 0.9915 0.9093 0.5985 0.1411] p1=polyfit(t,y,2) y1=polyval(p1,2.25)
y2=interp1(t,y,2.25,'spline')
10、试用MATLAB求解下面的问题:
sin2x?sin2a(1) 求极限lim；
x?ax?a7
(2) 不定积分(x2?e2x?sinx?ex)dx；
?(3) 对（2）的不定积分结果进行微分，看是否能还原原函数；
(4) 对xsin(cosx)cosx函数做20项的Taylor幂级数展开；
22(5)求级数S(n)=1?2?3?......?n的和。 解：（1）f=((sin(x))^2-(sin(a))^2)/(x-a); limit(f,x,a)
（2）f=[x^2*e^(2*x)+sinx*e^x]; int(f)
(3) x=sym('x');
f=1/3*x^3+2*x+log(x^2+1)+3*atan(x); diff(f)
(4) f=x^2*sin(cos(x^2))*cos(x); taylor(f,x,21)
(5) S=symsum(n^2,n,1,inf) S=symsum(n^2,n,1,100)
11、求解下列方程组的解
2222?sin(x?y)?exy?0(1)? (2) 2?x?cosy?0?x2?xy?y?0 ?2?x?4x?3?0解（1）[x y]=solve('sin(x+y)-exp(x)*y=0','x^2-cos(y)=0','x,y') （2）[x y]=solve('x^2+x*y+y=0','x^2-4*x+3=0','x,y')
8
12、求微分方程(组)的解
(1)y(3)?y''?x, y(1)?，8 y'(1)?7， y''(2)?4
(2)(x)?x?1， x(0)?0
'22dy?2x?ydt(3)
dx?4x?2ydt解：(1) y=dsolve('D3y+Dy-x','y(1)=8','Dy(1)=7','Dy(2)=4','x') (2) x=dsolve('(Dx)^2+x^2-1','x(0)=0')
(3) [x,y]=dsolve('Dy-2*x+y','Dx-4*x+2*y','t')
作业2
1、一球从100米高度自由落下，每次落地后反跳回原高度的一半，再落下。求它在第10次落地时，共经过多少米?第10次反弹有多高? 解：y=0.5*x; sum=0; for i=2:10, sum=sum+x+y; x=0.5*x; y=0.5*x; end sum y
9
2、用MATLAB的M函数文件定义如下分段函数：
?5,x?10?1?y?f(x)??x, x?10
?2x??10?5,??解：function y=f(x)
if x>10 y=5; else
if x>=-10&x<=10 y=1/2*x; else y=-5; end
3、分别用for和while循环编写程序，求出
K??2i?1?2?22?23?i?163?262?263
并考虑一种避免循环的简洁方法来进行求和，并比较各种算法的运行时间。 解：x=2; K=0;
for i=1:63, K=K+x^i; end K
4、应用MATLAB语言及二分法编写求解一元方程f(x)?x?14x?59x?70?0在区间[3，6]的实数解的算法，要求绝对误差不超过0.001。 解：A=[1 -14 59 -70]; a=3; b=6; c=0.001;
while c<0.5*(b-a), x1=(a+b)/2; f1=polyval(A,x1); fa=polyval(A,a); fb=polyval(A,b); if f1==0; x=x1;
10
32else if f1*fa>0; a=x1;
else f1*fb>0; b=x1; end end end x=x1
5、二阶系统的单位阶跃响应为y?1?11??2e??tsin(1??2t?acos?)，在同一平面绘
制?分别为0，0。3，0。5，0。707的单位阶跃响应曲线。要求：
(1) 四条曲线的颜色分别为蓝、绿、红、黄，线型分别为“——”、“……”、“oooooo”、“++++++”；(2) 添加横坐标轴和纵坐标轴名分别为“时间t”和“响应y”，并在平面图上添加标题“二阶系统曲线”和网格；
(3) 在右上角添加图例（即用对应的字符串区分图形上的线），并分别在对应的响应曲线的第一个峰值处标示“zeta＝0”、“zeta＝0。3”、“zeta＝0。5”、“zeta＝0。707”。
解：t=[0:0.1:10]; kos=0;
y=1-1/sqrt(1-kos^2)*exp(-kos*t).*sin(sqrt(1-kos^2)*t+2*cos(kos)); plot(t,y,'b-') hold on; kos=0.3;
y=1-1/sqrt(1-kos^2)*exp(-kos*t).*sin(sqrt(1-kos^2)*t+2*cos(kos)); plot(t,y,'g.') hold on; kos=0.5;
y=1-1/sqrt(1-kos^2)*exp(-kos*t).*sin(sqrt(1-kos^2)*t+2*cos(kos)); plot(t,y,'ro') hold on; kos=0.707;
y=1-1/sqrt(1-kos^2)*exp(-kos*t).*sin(sqrt(1-kos^2)*t+2*cos(kos)); plot(t,y,'y+') hold on;
xlabel('时间t'); ylabel('响应y'); grid on;
title('二阶系统曲线');
11
legend('kos=0','kos=0.3','kos=0.5','kos=0.707') gtext('kos=0') gtext('kos=0.3') gtext('kos=0.5') gtext('kos=0.707')
6、绘制如下图所示的图形，把图形窗口分割为2行2列，窗口左上角画一正弦曲线
y?sin(2?t),t?[0,2]；窗口
右上角画3条单边指数曲线
y?e?t,y?e?2t,y?e?3t,t?[0,2]；窗口左下角画一矩形脉冲信号，脉冲宽度为1，高为2，开始时间为1；窗口右下角画一单位圆。
解：t=0:0.01:2;
subplot(2,2,1);plot(t,sin(2*pi*t));title('plot(x,2*PI*t)');grid on; subplot(2,2,2);
plot(t,exp(-t));hold on;
plot(t,exp(-2*t),'g');hold on; plot(t,exp(-3*t),'r');hold on;
title('exp(-t),exp(-2*t),exp(-3*t)');grid on; x=[0 1 1 2 2 3 4]; y=[0 0 2 2 0 0 0]; subplot(2,2,3); plot(x,y);
axis([0 4 -0.5 3]); title('pulse signal');
12
syms x y
subplot(2,2,4); ezplot(x^2+y^2-1); axis([-1.5 1.5 -1.2 1.2]); title('circle'); xlabel('') ylabel('')
7、已知函数z?f(x,y)?1(1?x)?y22?1(1?x)?y22，试分别应用三维曲线图绘制命令plot3、三维网线图绘制命令mesh、三维曲面图绘制命令surf在同一窗口中绘制出3个子
图。
解：x=-2:0.1:2; [X,Y]=meshgrid(x);
Z=1./sqrt((1-X).^2+Y.^2)+1./sqrt((1+X).^2+Y.^2); subplot(311); plot3(X,Y,Z); grid on;
axis([-1 1 -1 1 -1 9]); subplot(312); mesh(X,Y,Z);
axis([-1 1 -1 1 -1 9]); subplot(313); surf(X,Y,Z);
axis([-1 1 -1 1 -1 9]);
8、求解下面两个方程构成的联立方程组在区间x?[?2,6],y?[?3,5]内的解，并用绘图的方法绘出两曲线在同一坐标上的图，以验证求得的解的正确性。
13
(x?2)2?3(y?1)2?3， y?6(x?22)
解：clc
[x,y]=solve('(x-2)^2+3*(y-1)^2-3=0','y-6*(x-2)^2=0') ezplot('(x-2)^2+3*(y-1)^2-3');hold on ezplot('y-6*(x-2)^2=0');grid on axis([0 4 -1 4]);
作业3
4(s?2)(s2?6s?6)21、已知系统的传递函数为G(s)?，试在MATLAB中建立其传332s(s?1)(s?3s?2s?5)递函数模型，并将其转化为零极增益模型和部分分式模型。 解%[n,d]=numden(G);num=sym2poly(n);den=sym2poly(d); syms s ;
G=4*(s+2)*(s^2+6*s+6)^2/(s*(s+1)^3*(s^3+3*s^2+2*s+5)); [n,d]=numden(G); num=sym2poly(n); den=sym2poly(d); [z,p,k]=tf2zp(num,den)
14
[r1,p1,k1]=residue(num,den)
2、一反馈控制系统的前向通道传递函数和反馈通道传递函数分别为
G(s)?3(s?6)5s?1和，试在MATLAB中实现两个子系统的负反馈H(s)?2(s?1)(s?3s?5)15s?1联结，给出该控制系统的传递函数和以零极点形式表示的数学模型，并绘出单位阶跃响应曲线图。
解：G1=tf(3*[1 6],conv([1,1],[1,3,5])); H=tf([5 1],[15 1]) Gtf=feedback(G1,H) Gzpk=zpk(Gtf) Gzpk=zpk(Gtf) step(Gtf); hold on; step(Gzpk,'r');
?n23、二阶系统的传递函数为H(s)?2，设其固有频率?n?10，在阻尼
s?2??ns??n2??[0.1,0.3,0.7,1]时，分别画出其单位阶跃响应和单位脉冲响应曲线。
15
解：wn=10; syms ks
for ks=[0.1 0.3 0.7 1] num=[1];
den=[1 20*ks 100]; figure(1); step(num,den); hold on; figure(2); impulse(num,den); hold on; end figure(1); grid on; figure(2); grid on;
16
4、已知某控制系统的开还传递函数为G(s)?环频率特性曲线，即系统的Bode图。 解：K=1.5; num=[K];
den=conv([1 0],conv([1 1],[1 2])); bode(num,den)
5、二阶系统的传递
K，K?1.5，试绘制系统的开
s(s?1)(s?2)?n2函数为H(s)?2，设其固有频率?n?10，在阻尼??[0.1,0.3,0.7,1]时，
s?2??ns??n2分别画出其Bode图，研究阻尼系数对二阶系统频率响应的影响。 解wn=10;
for ks=[0.1 0.3 0.7 1] num=[1];
den=[1 20*ks 100];
17
bode(num,den); hold on; end grid on;
3s3?16s2?41s?286、系统模型如下所示：G(s)?6，判断5432s?14s?110s?528s?1494s?2117s?112系统的稳定性，以及系统是否为最小相位系统。 解：num=[3 16 41 28];
den=[1 14 110 528 1494 2117 112];
[z,p,k]=tf2zp(num,den)%求系统的零极点，并显示 ii=find(real(p))>0 n1=length(ii); jj=find(real(z)>0) n2=length(jj); if(n1>0)
disp('the system is Unstable') else
disp('the system is stable') end
%判断系统是否为最小相位系统 if(n2>0)
disp('the system is a nonminimal phase one') else
disp('the syetem is a minimal phase one')
18
end
%绘制零极点图 pzmap(p,z)
作业4
1、由周期信号的分解与合成可知方波分解为多次正弦波之和。如图所示的周期方波，其傅
立叶级数为
f(t)?1[sint?sin3t??34?1sin(2k?1)t?2k?1]
用MATLAB演示谐波合成情况。 解：t=-2*pi:0.1:2*pi;
19
sum=0;
for a=1:2:300 a1=1./a;
sum=sum+4*a1*sin(a*t)/pi; plot(t,sum) grid on pause(0.25) end
0 pi 2pi 2、考虑简单的线性微分方程
dydx?2yx?4x?0，且方程的初值为y(1)?2，解析解和在区间[1,2]上的数值解，并比较二者得出的曲线。
解：syms x1 y1 x y;
y1=dsolve('Dy1+2*y1/x-4*x','y1(1)=2','x') f=inline('4*x-2*y/x','x','y') [x,y]=ode45(f,[1,2],2) subplot(121); ezplot(y1);
axis([1 2 0 4]) subplot(122); plot(x,y,'r')
axis([1 2 0 4])
20
试求该方程的 d2ydy?8y?u，求其冲击响应。3、设二阶连续系统，其特性可用常微分方程表示为2?2dtdt若输入为u?3t?cos(0.1t)，求其零状态响应。
解：
num=[1];den=[1 2 8];impulse(num,den)symssg1XGtG=1/(s^2+2*s+8);X=laplace(3*t+cos(0.1*t)) ;g1=G*X;[n,d]=numden(g1);num=sym2poly(n);den=sym2poly(d);figure(2);impulse(num,den)axis([0 5 0 2]);
4、一阶低通电路的频率响应 下图所示是一阶RC低通电路，若以Uc为响应，求频率响应函数画出其幅频响应(幅频特性)H(j?)和相频的响应(相频特性)?(?)。
21
解：G(jw)=1/1+jwrcden=[2 1];num=[1];bode(num,den);gridon;
5、n级放大器，每级的传递函数均为性。
?0
，求阶跃响应，画出n不同时的波形和频率特s??0
解：
w0=input('??ê?è?w0μ??μ:')num=[w0];den=[1 w0];G0=tf(num,den);forn=1:6G=G0^n; figure(1); step(G); gridon; holdon; figure(2); bode(G); gridon; holdon;end
22
23
6、下图试典型的二阶动态电路，其零输入响应有过阻尼，临界阻尼和欠阻尼三种情况。如
iL(0)?0，已知L?0.5H，初始值uc(0)?1V，求t?0时的uc(t)C?0.02F，R?12.5?，
和iL(t)的零输入响应，并画出波形。
解: clear,formatcompactL=0.5;
C=0.02;R=12.5; uc0=1; iL0=0; dt=0.01;t=0:dt:3; num=[uc0,R/L*uc0+iL0/C]; den=[1,R/L,1/L/C]; [r,p,k]=residue(num,den);
ucn=r(1)*exp(p(1)*t)+r(2)*exp(p(2)*t); iLn=C*diff(ucn)/dt;
figure(1),plot(t,ucn),holdon;gridon;
figure(2),plot(t(2:end),iLn),holdon;gridon;
24
25
作业5
?n21、二阶系统的传递函数为H(s)?2，设其固有频率?n?10，在阻尼2s?2??ns??n??[0.1,0.3,0.7,1]，输入信号为单位阶跃信号时，分别按以下要求利用Simulink绘制系统
的响应曲线：（1）在同一幅图的同一窗口；（2）在同一幅图的不同窗口。 解：
26
...2、考虑简单的线性微分方程y(4)?3y(3)?3y?4y?5y?e?3t?e?5tsin(4t?)，
3?27
且方程的初值为y(0)?1，y(0)?y(0)?仿真模型，并绘制出仿真结果曲线。 解：
...1(3)，y(0)?0.2，试用Simulink搭建起系统的2
28
3、考虑上面的模型，假设给定的微分方程变化成时变线性微分方程
...y(4)?3ty(3)?3ty?4y?5y?e?3t?e?5tsin(4t?)
3...2?而方程的初值仍为y(0)?1，y(0)?y(0)?的仿真模型，并绘制出仿真结果曲线。
1(3)，y(0)?0.2，试用Simulink搭建起系统2解：
29
dx1
?x2dt
4、已知一个系统的微分方程表示为： ，其中，状态变量初始条件
dx2
?u?5x1dtx1(0)?0,x2(0)?0,输入u为阶跃函数，要求：（1）利用Simulink对系统建立仿真模型，
并绘制时域响应曲线；（2）将x1和x2输出到工作空间，并绘制其曲线。
解：
syms t figure(1)
30
ezplot(t,x1); axis([0 10 0 0.5]); figure(2) ezplot(t,x2); axis([0 10 -0.5 0.5]);
s2?5s?25、假设已知线性系统模型为G(s)?，试输入该系统模型，并求出系统在4(s?4)?4s?4阶跃输入和斜坡输入下的解析解，并和仿真曲线相比较，检验得出的结果。
解：
31
syms s t;
Y=(s^2+5*s+2)/((s+4)^4+4*s+4); [n,d]=numden(Y); num=sym2poly(n) den=sym2poly(d) figure (1) impulse(num,den); figure (2) step(num,den); figure (3)
den1=[1 16 96 260 260 0]; step(num,den1);%斜坡响应
32
6、分别使用解微分方程方法、控制工具箱和Simulink求解传递函数为
G(s)?10的系统的阶跃响应。
s4?8s3?36s2?40s?1033
解：
num=[10];
den=[1 8 36 40 10]; step(num,den);
34
num=[10];
den=[1 8 36 40 10]; step(num,den);
35
7、(选做)已知一级倒立摆的数学模型为
..?f/m?l?2sin??gsin?cos?y??M/m?sin2?? ?..2????fcos??(M?m)gsin?/m?l?sin?cos??l(M/m?sin2?)?其中?为摆体与垂直方向的夹角(单位为rad)，y为小车的位移(单位为m)，f为电机对小车的作用力(单位为N)，M和m分别为小车和摆体的质量(单位为kg)，l为摆长的一半(单位为m)，g为重力加速度(9.81m/s)。试建立起倒立摆的Simulink模型。若取m?0.21kg，
2M?0.455kg，l?0.61/2m，并取f为系统的输入信号，试在平衡点y???0处对系统
进行线性化，并比较原系统和线性化系统的阶跃响应曲线。
36