matlab程序大全

频率特性类题目

1 一个系统的开环传递函数为 ，试绘制其当K=5、30时系统的开环频率特性Nyquist 图，并判断系统的稳定性。

w=linspace(0.5,5,1000)*pi; sys1=zpk([],[0 -10 -2],100) sys2=zpk([],[0 -10 -2],600) figure(1) nyquist(sys1,w)

title('system nyquist charts with k=5') figure(2)

nyquist(sys2,w)

title('system nyquist charts with k=30')

由图可知K=5时，开环Nyquist 曲线没有包围（-1，j0）点，所以系统稳定。

K=30时，开环Nyquist 曲线包围（-1，j0）点，所以系统不稳定。

2系统开环传递函数为 ，建立其零极点增益模型，然后分别绘制当K=5、K=30时系统的开环频率特性Bode 图，并判断系统的稳定性。

sys1=zpk([],[0 -10 -2],100) sys2=zpk([],[0 -10 -2],600)

figure(1)

[Gm1,Pm1,Wcg1,Wcp1]=margin(sys1) bode(sys1)

title('system bode charts with k=5'),grid figure(2)

[Gm2,Pm2,Wcg2,Wcp2]=margin(sys2) bode(sys2)

title('system bode charts with k=30'),grid 因为K=5时，Wcg>Wcp,所以系统稳定。

K=10时，Wcg

3 系统开环传递函数为 ，计算K=5和K=30时系统

的幅值裕度与相位裕度。

sys1=zpk([],[0 -10 -2],100)

sys2=zpk([],[0 -10 -2],600)

[Gm1,Pm1,Wcg1,Wcp1]=margin(sys1) [Gm2,Pm2,Wcg2,Wcp2]=margin(sys2)

4 已知某系统的闭环传递函数()s Φ如下，试用roots （）命令来判断系统的稳定性。

)

5.01)(1.01()(s s s k

s G ++=)5.01)(1.01()(s s s k

s G ++=)

5.01)(1.01()(s s s k

s G ++=

2

5

4

3

2

325

()24576

s s s s s s s s ++Φ=

+++++

num=[0 0 0 3 2 5]; den=[1 2 4 5 7 6]; D=num+den; roots(D)

5 某单位负反馈系统的开环控制系统的传递函数为

2

k (0.80.64)()(0.05)(5)(40)

K s s G s s s s s ++=

+++

（1）绘制系统的根轨迹；

num=[1 0.8 0.64]; x=[0 -0.05 -5 -40]; den=poly(x); rlocus(num,den)

（2）当10K =时，绘制系统的Bode 图，判断系统的稳定性，并且求出幅值裕度和相角裕度。

num=[1 0.8 0.64]; x=[0 -0.05 -5 -40]; den=poly(x); sys=tf(num,den) bode(sys),grid

[Gm,Pm,Wcg,Wcp]=margin(sys)

6 已知系统的状态空间模型如下：

???=1

1x

??

?-31x +

??

????01u

[1=y ]1x

（1）绘制系统的Bode 图和nyquist 图； （2）求系统的幅值裕度和相位裕度；

A=[1,-1;1,3];B=[1,0]'; C=[1,1];D=0;

w=linspace(0.5,5,1000)*pi; sys=ss(A,B,C,D)

figure(1),bode(sys),grid figure(2),nyquist(sys,w) [Gm,Pm,Wcg,Wcp]=margin(sys)

7 已知单位负反馈系统的开环传递函数为

)

1(12

++s s s ，试绘制系统的单位阶跃响

应、开环Bode 图和Nyquist 曲线，并求系统的幅值裕度和相位裕度。

sys=zpk(-1,[0 0 -1],1) sys1=feedback(sys,1,-1) figure(3) step(sys1) figure(1)

bode(sys),grid figure(2)

nyquist(sys),grid

[Gm,Pm,Wcg,Wcp]=margin(sys)

编写程序题

1 编写两个m 文件，分别使用for 和while 循环语句计算200

31

k k =∑。

y=0; for k=1:200 y=y+k^3; end y

y=0;k=1; while k<=200 y=y+k^3; k=k+1; end y

2 一个三位整数各位数字的立方和等于该数本身则称该数为水仙花数，编写M 文件，输出全部水仙花数。

for m=100:999 m1=fix(m/100);

m2=fix((m-100*m1)/10); m3=fix(m-100*m1-10*m2); if m==m1^3+m2^3+m3^3 disp(m) end end

3 编写M 文件，实现以下功能：从键盘输入若干个数，当输入0时结束输入，求这些数的平均值和它们之和。

4 编写M 脚本文件，使得运行该文件后能求出[100，200]之间第一个能被21整除的整数。

for i=100:200 if mod(i,21)==0 break end end i

5 先编写一个函数文件，文件名为tran ，使得该函数文件实现直角坐标(x,y)与极坐标(ρ,θ)之间的转换。然后编写脚本文件，文件名为main1，在脚本文件中实现如下功能：要求用户从键盘输入直角坐标X 和Y ，调用函数文件实现直角坐标到极坐标的转换，输出极坐标的值。

tran.m ：

function [rho,theta]=tran(x,y) rho=sqrt(x*x+y*y); theta=atan(y/x); main1.m ：

x=input('Please input x=:') y=input('Please input y=:') [rho,theta]=tran(x,y); rho theta

6 猜数游戏。首先由计算机产生[1,100]之间的随机整数，然后由用户猜测所产生的随机数。根据用户猜测的情况给出不同提示，如猜测的数大于产生的数，则显示“High ”，小于则显示“Low ”，等于则显示“You won ”，同时退出游戏。用户最多可以猜7次。

x=fix(100*rand(1)); i=1 while i<=7 y=input('') if x==y

disp('Y ou won') break; elseif x>y

disp('Low') else disp('High') end i=i+1; end

7 编写M 脚本文件，分别使用for 和 while 循环语句计算∑==10

1

i i i sum 的值,并

且显示出sum 的结果。

sum=0;

for i=1:10

sum=sum+i^i;

end

sum

sum=0;

while i<=10

sum=sum+i^i;

i=i+1;

end

sum

8 用筛选法求某自然数范围内的全部素数。素数是大于1，且除了1和它本身以外，不能被其他任何整数所整除的整数。用筛选法求素数的基本思想是：要找出2～m之间的全部素数，首先在2～m中划去2的倍数(不包括2)，然后划去3的倍数(不包括3)，由于4已被划去，再找5的倍数 (不包括5)，…，直到再划去不超过的数的倍数，剩下的数都是素数。

function y=sushu(x)

y=2:x;

n=2:x;

for j=2:x

m=j*n;

y=setdiff(y,m);

end

9 已知∑

=-

=

n

k k

y

1

1

2

1

，当n=100时，求y的值。

y=0;n=100;

for k=1:n

y=y+(1/(2*k-1));

k=k+1;

end

y

10 利用for循环计算1!+2!+3!+….. +20!的值。

sum=0;

for i=1:20

y=1;

for n=1:i

y=y*n;

end

sum=sum+y;

end

sum

11用while循环来实现求1+2+……+100的值。

y=0;n=1;

while n<=100 y=y+n; n=n+1; end y

12 编一函数计算以下分段函数值。

2

32(0)21

(01)321(1)x

x x x x x x

≤≥ function y=f(x); if x<0

y=x; elseif x>=0&x<1 y=2*x^2+1; else y=3*x^3+2*x^2+1; end

13 我国新税法规定：个体工商户的生产、经营所得和对企事业单位的承包经营、承租经营所得应缴纳的个人所得税为

全年收入中应缴纳所得税部分 税率/% (1) 不超过5000元的 5 (2) 超过5000元至10000元的部分 10 (3) 超过10000元至30000元的部分 20 (4) 超过30000元至50000元的部分 30 (5) 超过50000元的部分 35 试编程加以计算。

x=input('x=') if x<=5000 y=0.05*x;

elseif x>5000&x<=10000 y=0.1*x;

elseif x>10000&x<=30000 y=0.2*x;

elseif x>30000&x<=50000 y=0.3*x; else y=0.35*x; end y

14 矩阵乘法运算要求两矩阵的维数相容，否则会出错。设两个矩阵为A=[1,2; 5,6]; B=[7,9;20,11]；先求两矩阵的乘积，若出错，则自动转去求两矩阵的点乘。

A=[1,2;5,6];B=[7,9;20,11]; try

C=A*B; catch C=A.*B;

C

15 编制M脚本文件，要求从键盘逐个输入数值（input），然后判断输入的数是大于0还是小于0，并输出提示（使用disp命令）是正数（positive one）还是负数（negative one），同时记录输入的正数、负数的个数。当输入0时，终止此M文件的运行（使用break命令）；当输入第10个数字时，显示记录的正、负数个数并终止程序。

i=0;j=0;

for n=1:10

x=input('The num is:')

if x>0

i=i+1;

disp('positive one')

elseif x<0

j=j+1;

disp('negative one')

else break

end

end

i,j

16 一个三位整数各位数字的立方和等于该数本身则称该数为水仙花数，编写M 文件，输出全部水仙花数。

17 编写M文件，实现以下功能：从键盘输入若干个数，当输入0时结束输入，求这些数的平均值和它们之和。

sum=0;cnt=0;

x=input('x=')

while x~=0

sum=sum+x;

cnt=cnt+1;

x=input('x=')

end

if cnt>0

average=sum/cnt;

end

sum,average

18 一数组的元素满足规则： a(1)=a(2)=１, a(3)=a(1)+a(2),

a(4)=a(2)+a(3),…依次类推，求该数组中第一个大于10000的元素。

a(1)=1;a(2)=1;i=1;

while a(i)<=10000

a(i+2)=a(i+1)+a(i);

i=i+1;

end

i

19 某商场对顾客所购买的商品实行打折销售，标准如下(商品价格用price 来表示)：

price<200 没有折扣 200≤price<500 3%折扣 500≤price<1000 5%折扣 1000≤price<2500 8%折扣 2500≤price<5000 10%折扣

5000≤price 14%折扣

请编写程序，使得当从计算机输入所售商品的打折前价格后，能输出其实际销售价格。

x=input('x=') if x<200 y=x;

elseif x>=200&x<500 y=0.97*x;

elseif x>=500&x<1000 y=0.95*x;

elseif x>=1000&x<2500 y=0.92*x;

elseif x>=2500&x<5000 y=0.9*x; else y=0.86*x; end y

20 编写M 脚本文件，分别使用for 和 while 循环语句计算∑==10

1

i i i sum 的值,

并且显示出sum 的结果。

21 已知 ∑=-=

n

k k y 1

121

，当n=100时，求y 的值。

22编写函数文件使其实现如下功能：当在命令窗口输入数字后，根据输入参数的个数决定输出的结果；当没有参数输入时，显示“无输入参数”；当有一个参数时，则以该参数为边绘制正方形；当有两个参数时，则以两个参数为长和宽绘制矩形。

function y=H(varargin) % varargin 表示可变的参数 if nargin==0

disp('无输入参数') elseif nargin==1 r=varargin{1};

rectangle('position',[0 0 r r]) %表示从（0，0）点开始画宽为R,高为R的矩形else nargin==2

r1=varargin{1};

r2=varargin{2};

rectangle('position',[0 0 r1 r2])

end

23 编制M脚本文件，要求从键盘逐个输入数值（input），然后判断输入的数是大于0还是小于0，并输出提示（使用disp命令）是正数（positive one）还是负数（negative one），同时记录输入的正数、负数的个数。当输入0时，终止此M文件的运行（使用break命令）；当输入第10个数字时，显示记录的正、负数个数并终止程序。

i=0;j=0;

for x=1:10

n=input('')

if n<0

i=i+1;

disp('negative one')

elseif n>0

j=j+1;

disp('position one')

else break

end

end

i

j

24 若一个数等于它的各个真因子之和，则称该数为完数，如6=1+2+3，所以6是完数。求[1,500]之间的全部完数。

for m=1:500

s=0;

for k=1:m/2

if rem(m,k)==0

s=s+k;

end

end

if m==s

disp(m)

end

end

25 今有多项式P1(x)=x^4-2x+1,P2(x)=x^2+4x-0.5，要求先求得P(x)=P1(x)+P2(x)，然后计算xi=0.2*i各点上的P(xi)（i=0,1,2,…,5）值。

A=[1 0 0 -2 1];B=[0 0 1 4 -0.5];

C=A+B

i=0:5;

xi=0.2*i;

y=polyval(C,xi)

26 试编一个m 程序，将一维数组x 中的N 个数按颠倒的次序重新存储。要求运行程序后，在命令窗口提示“请用[]输入一个一维数组；”，输入完毕后若回车，就能在输出窗口显示刚才所输入的数组的颠倒次序。（如N=5，原来x 为：

x=[ 1 3 5 7 9 ]

而经过颠倒处理后x 中数据的次序应该为： x=[ 9 7 5 3 1 ]）

x=input('请用[]输入一个一维数组:x=') x1=fliplr(x)

模型建立及转换题

1 某系统的传递函数为

2

3

21.323()0.5 1.21

s s G s s s s ++=

+++

使用MATLAB 求出状态空间表达式和零极点模型。

num=[1.3 2 3];den=[1 0.5 1.2 1]; sys=tf(num,den) nsys=zpk(sys)

2 求出以下系统的传递函数

10

1012000

031u -????

????

=-+????

????-????X

X ， []1

1

0y =X

A=[-1,0,1;1,-2,0;0,0,-3];B=[0,0,1]';C=[1 1 0];D=1;

sys=ss(A,B,C,D) nsys=tf(sys)

3 已知某系统的传递函数为：2

264220122412)(2

3

4

2

3+++++++=

s s s s s s s s G ，试用MATLAB 语言

求出该系统的传递函数模型、状态空间模型和零极点增益模型。

num=[12 24 12 20];den=[2 4 6 2 2]; sys=tf(num,den) sys1=zpk(sys) sys2=ss(sys)

4 已知双闭环调速系统电流环内的前向通道3个模块传递函数分别为：

试求串联连接的等效传递函数及其等效状态空间模型。

num1=[0.0128 1];den1=[0.04 0]; num2=[30];den2=[0.00167,1]; num3=[2.5];den3=[0.0128 1];

s

s s G 04.010128.0)(1+=

1

00167.030)(2+=

s s G 1

0128.05.2)(3+=

s s G

sys1=tf(num1,den1) sys2=tf(num2,den2) sys3=tf(num3,den3) sys=sys1*sys2*sys3 nsys=ss(sys)

5 已知两个环节的传递函数分别为：

试求两环节并联连接等效传递函数的num 与den 向量及等效的状态空间模型。 num1=[1];den1=[1 2];

num2=[2 1];den2=[1 1 2]; sys1=tf(num1,den1) sys2=tf(num2,den2) sys=sys1+sys2

[num,den]=tfdata(sys,'v') nsys=ss(sys)

6 将传递函数模型2

3)(2

++=

s s s s G 转换为状态空间模型和零极点模型。

num=[1 0];den=[1 3 2]; sys=tf(num,den) sys1=zpk(sys) sys2=ss(sys)

7 系统方程为u y y y

y 616116=+++ ,求出其状态空间表达式，写出A 、B 、C 、D 阵，建立其MATLAB 状态空间模型，然后绘制其单位阶跃响应曲线。

A=[0 1 0;0 0 1;-16 -11 -6];B=[0 0 6]'; C=[1 0 0];D=0; sys=ss(A,B,C,D) step(sys)

8 设有一质量弹簧阻尼系统，如下图所示。F(t)为输入的外力，y(t)为质量块的

位移，弹簧的弹性系数为k ，阻尼器的阻尼系数为f 。系统的动态方程组为：

)()()()(t F t ky t y f t y

m =+?+ 求出其状态空间表达式，写出A 、B 、C 、D 阵，建立其MATLAB 状态空间模型，然后绘制其单位阶跃响应曲线（只要写出MATLAB 语句即可）。

A=[0 1;-k/m -f/m];B=[0 F(t)/m]; C=[1 0];D=0;

2

1)(1+=

s s G 2

12)(2

2+++=

s s s s G

sys=ss[A,B,C,D] step(sys)

9 对下图中所示系统进行建模、模型转换及仿真：

⑴ 用MATLAB 语言表示图示系统的模型，并给出传递函数、零极点和状态空间三种模型的具体表达。

num1=[1 0.5];den1=[1 0.1]; z=[];p=[0 -2 -10];k=[20]; sys1=tf(num1,den1) sys2=zpk(z,p,k) sys3=sys1*sys2 sys=feedback(sys3,1,-1) nsys1=tf(sys) nsys2=ss(sys)

⑵ 用MATLAB 语言画出系统单位阶跃响应的波形。

10 对图中所示的质量—弹簧—阻尼机械系统，系统的动力学方程为：

ku

ky y c y m =++?

?

?,其中，m=5kg ，k=2N/m ，c=0.1N/m.s 1-，建立其MATLAB 状态空

间模型，写出A 、B 、C 、D 阵，然后绘制其单位阶跃响应曲线。

A=[0,1;-2/5,-0.1/5]; B=[0,2/5]'; C=[1,0]; D=0;

sys=ss(A,B,C,D) step(sys)

y

k

11求下面系统在阶跃信号为0.11(t)时系统的响应。

20

4036820

2

3

4

++++s s s s

并求系统性能指标：稳态值、上升时间、调节时间、超调量 num=[20];den=[1 8 36 40 20]; sys=tf(num,den)

final=polyval(num,0)/polyval(den,0) [y,t]=step(sys) [yp,k]=max(y) tp=t(k)

chaotiao=100*(yp-final)/final n=1;

while y(n)

risetime=t(n) i=length(t);

while (y(i)>0.98*final&y(i)<1.02*final) i=i-1; end ts=t(i)

step(0.1*sys)

12 已知系统传递函数01

.12.01

)(2

++=

s s G s

，绘制系统阶跃响应曲线；绘出离

散化系统阶跃响应曲线，采样周期s T s 3.0=。

num=[1];den=[1 0.2 1.01]; sys1=tf(num,den) sys2=tf(num,den,0.3) figure(1) step(sys1)

figure(2) step(sys2)

数值计算类题目

1 求解以下线性代数方程：

123102211313

1

21x x x ??????

??????=??????????????????

A=[1,0,2;1,1,3;3,1,2]; B=[2,1,1]'; X=inv(A)*B

2 已知矩阵

013=121542??

??

??????

A ，218=414332??

????????

B 试分别求出A 阵和B 阵的秩、转置、行列式、逆矩阵以及特征值。

A=[0 1 3;1 2 1;5 4 2];B=[2 1 8;4 1 4;3 3 2];

A1=rank(A) B1=rank(B) A2=A.' B2=B.' A3=det(A) B3=det(B) A4=inv(A) B4=inv(B) A5=eig(A) B5=eig(B)

3 已知矩阵

013=121542??

??

??????

A ，218=

414332??

????????

B 在Command 窗口中分别求出2

C =A 、

D 矩阵为A 中每个元素平方组成的矩阵、

E

矩阵为A 阵乘以B 阵、F 矩阵为A 阵和B 阵数组乘积（即：对应元素分别相

乘的积构成的矩阵）。

A=[0 1 3;1 2 1;5 4 2];B=[2 1 8;4 1 4;3 3 2]; C=A^2 D=A.^2 E=A*B F=A.*B

4 求复数矩阵13

573

62

32

43

i i i i i i +-+??=?

?++

-??C 的转置C1与共轭转置C2。

A=[1 5 7;6 3 4];B=[3 -1 3;2 2 -3];

C=A+B*i C1=C.' C2=C'

5 假设a 、b 、c 和d 的定义如下：

???=2

1

a ?

??10； ???-=0

1b ???12； ???

???=23c ； 5=d

分别写出下列表达式的运算结果。

(1) a*b (2) a*c (3) a.*d (4) a*d

ans = -1 2 -2 5 ans = 3 8 ans = 5 0 10 5 ans = 5 0

6 求解由下面四个方程组成的方程组。

5

22778638

232432143214214321=+-+=++-=++=++-x x x x x x x x x x x x x x x

A=[2 -3 1 2;1 3 0 1;1 -1 1 8;7 1 -2 2]; B=[8 6 7 5]'; C=inv(A)*B

7 求解由下面四个方程组成的方程组。

51

2279823317232432143214214321=+-+=++-=++=++-x x x x x x x x x x x x x x x

A=[2 -3 1 2;1 3 0 1;1 -1 1 8;7 1 -2 2]; B=[17 23 9 51]'; C=inv(A)*B

8 已知某班的5名学生的三门课成绩列表如下：

学生序号 1 2 3 4 5 高等数学 78 89 64 73 68 外语 83 77 80 78 70 MATLAB 语言 82 91 78 82 68

试写出有关命令，先分别找出三门课的最高分及其学生序号；然后找出三门课总分的最高分及其学生序号。

A=[78 89 64 73 68;83 77 80 78 70;82 91 78 82 68]; [Y ,L]=max(A,[],2) B=sum(A,1); [Y1,L1]=max(B)

SIMULINK 建模分析题

1 某一单位负反馈控制系统，其开环传递函数为

1()(1)

G s s s =

+

给它输入阶跃量为2的阶跃信号，试使用Simulink 构造其闭环仿真模型，并且绘制其响应曲线。然后将已经建立的闭环控制系统封装成一个子系统。

2 典型的直流电动机双闭环调速系统的结构图如图所示，编写MATLAB 程序（或利用SIMULINK 工具），画出系统的单位阶跃响应曲线。

3 双环调速的电流环系统的方框图模型为：

建立该系统的SIMULINK 模型，然后在单位阶跃输入信号作用下对系统进行仿真。根据仿真结果计算超调量、上升时间、调节时间、稳态误差。

4 利用sisotool 工具完成下面的题目： 系统开环传递函数为：

绘制当K=5、30时系统的Bode 图，并判断系统的稳定性；计算K=5、30时系统的幅值和相位裕度；绘制K=5、30时系统的Nyquist 图；绘制K=5、30时系统的Nichols 图。

5 已知二阶系统微分方程：.

0)1(2=+--y y y y

1)0(,0)0(==y

y 求在时间区间t=[0,20]内的微分方程的解。

Vf.m

Function dy=vf(t,y) Dy=zeros(2,1); Dy(1)=y(2);

Dy(2)=(1-y(1)^2)*y(2)-y(1); Main.m

[T,Y]=ode45(‘vf ’,[0 20],[0,1]);%调用离散ode45产生离散点时间向量和解向量 Plot(T,Y(:,1),’r-’,T,Y(:,2),’b:’) Legend(‘y1’,’y2’)

)

5.01)(1.01()(s s s K

s G ++=

6具有死区非线性特性的控制系统方框图如题下图所示。试绘制当输入信号

t

At

(t

A

t

=时系统响应曲线。

?

+

r?

)

(1

)

(1

)

(t

A

=和)(1

)

r?

t

图形绘制题

1绘制z=sin(x)*cos(y)的三维网格和三维曲面图，x,y变化范围均为[0，2π]。 x=0:pi/100:2*pi;

y=0:pi/100:2*pi;

[X,Y]=meshgrid(x,y)

Z=sin(X).*cos(Y);

subplot(1,2,1),surf(X,Y,Z); %曲面图

subplot(1,2,2),mesh(X,Y,Z); %网线图

2用MATLAB语言绘制［－2π，2π］区间内sinθ、cosθ曲线，数据点间隔

0.1弧度。要求分上下两个子图分别绘制sinθ、cosθ，并且sinθ用红色、实

线，cosθ用兰色、虚线，对图形标题及横纵坐标轴进行标注（图形标题分别为sinθ、cosθ，横纵坐标轴分别为x,y）。

t=-2*pi:0.1:2*pi;

y1=sin(t);

y2=cos(t);

subplot(2,1,1),plot(t,y1,'r')

title('sin(t)'),xlabel('x'),ylabel('y')

subplot(2,1,2),plot(t,y2,'b:')

title('cos(t)'),xlabel('x'),ylabel('y')

3 在[0 2π]范围内绘制二维曲线图y=sin(x)*cos(5x)。

x=0:pi/50:2*pi;

y=sin(x).*cos(5*x);

plot(x,y)

4 在[0 2π]范围内绘制二维曲线图。

y=|1000sin(4x)|+1 x=0:pi/50:2*pi;

y=abs(1000*sin(4*x))+1; plot(x,y)