例1.：多子图演示：

t=(pi*(0:1000)/1000)’ ; 

y1=sin(t); 

y2=sin(10*t);

y12=sin(t).*sin(10*t);

subplot(2,2,1),plot(t,y1);axis([0,pi,-1,1]) 

subplot(2,2,2),plot(t,y2); axis([0,pi,-1,1]) 

subplot(‘position’,[0.2, 0.05, 0.6, 0.45 ])  

plot(t, y12,’b-’, t , [y1, -y1], ‘r :’) ; axis([0, pi ,-1, 1])

例2:：已知某系统运动方程及初始条件：y1’=y2y3, y2’=-y1y3,y3’=-2y1y2,y1(0)=0,y2(0)=0.5,

y3(0)=-0.5,求时间区间t=[0,20]微分方程的解。

解：（1）建立m函数文件：

function dy=right(t,y)  

dy=zeros(3,1);

dy(1)=y(2)*y(3);

dy(2)=-y(1)*y(3);

dy(3)=-2*y(1)*y(2);

（2 ）编写调用函数right（）的M文件，并执行：

[T,y]=ode45(‘right’,[0,20],[0,0.5,-0.5]);

plot(T,y(:,1),’r’,T,y(:,2),’b*’,T,y(:,3),’k-.’)  

 legend(‘y1’,’y2’,’y3

例3：已知二阶微分方程：y’’-(1-y^2)y’+y=0,y(0)=0,y(0)’=1,求时间t=[0,20]微分方程解。

解：（1）将微分方程表示为一阶微分方程组：

y1=y,

y1’=y2,

y2’=(1-y12)y2-y1；

（2）建立描述系统微分方程的m-函数文件

function dy=vdp(t,y)  

dy=zeros(2,1);

dy(1)=y(2);

dy(2)=(1-y(1)^2)*y(2)-y(1)    

（3)编写matlab主程序，并执行。

[T,Y]=ode45(‘vdp’,[0 20],[0,1]);

plot(T,Y(:,1),’r-’,T,Y(:,2),’b:’)  

 legend(‘y1’,’y2’)

例4：已知系统模型G(s)=（3s+100）/(s^3+10s^2+40s+100),计算系统在周期为5s的方波信号作用下的响应。

解：sys=tf([3,100],[1,10,40,100]);

[u,t]=gensig(‘square’,5,10);

hold on  

plot(t,u,’-.’) 

hold off   

text(1.3,0.8,‘输入\rightarrow’)  

text(5.4,0.8,‘\leftarrow输出’)

例5：已知系统传递函数为G(s)=5(s+1)/(s^3+4s^2+6s+5),求最大超调量Mp，调整时间tp

解：num=[5 5];den=[1 4 6 5];

sys=tf(num,den) 

finalvalue=polyval(num,0)/polyval(den,0)  

 [y,t]=step(sys);

[yp,k]=max(y)  

tp=t(k)  

Mp=100*(yp-finalvalue)/finalvalue 

 len=length(t);

while(y(len)>0.98*finalvalue)&(y(len)<1.02*finalvalue)  

len=len-1;

end  

ts=t(len)  

step(sys)

例6：已知单位负反馈系统的开环传递函数为Gk(s)=10/[s(s+1)(s+5)],求其单位斜坡输入时，系统的稳态误差。

解：Gk=zpk([],[0 -1 -5],10);

Xi=zpk([],[0 0],1);

sys=1/(1+Gk);

Es=sys*Xi  

ess=dcgain(tf([1 0],[1])*Es)

t=[0:0.05:10];

xi=t;

y=lsim(sys*Gk,xi,t);

plot(t,xi,’r-.’,t,y,t,xi-y’,’k:’)

 legend(‘输入’,’输出’,’误差’,0) 

xlabel(‘t(s)’), ylabel(‘幅值、误差’)

例7：已知某一伺服机构的开环传递函数为G(s)=1/[s(0.5s+1)(0.1s+1)],设计PID控制器，使得系统的加速度误差系数Ka>=10，幅值穿越频率wc>=4rad/s,相位裕度>=50度。

解.Kv=10;

G=zpk([],[0 -2 -10],20);

Ki=Kv/dcgain(G*tf([1 0],1))  

 wc=4.1;

[num,den]=tfdata(G,’v’);

numc=polyval(num,j*wc);

denc=polyval(den,j*wc);

Gjwc=numc/denc;theta=(-180+50.1)*pi/180;

Ejwc=cos(theta)+j*sin(theta);

sum=Ejwc/Gjwc+j*Ki/wc;

Kp=real(sum)   

Kd=imag(sum)/wc  

Gc=tf([Kd Kp Ki],[1 0]); 

figure(1) bode(Gc),grid  

figure(2)  bode(G,G*Gc),grid   

 [kg,ph,wp,wc]=margin(G*Gc)

例8：已知单位反馈系统开环传递函数。系统输入信号为三角波(T=4s,幅值为1)。试求系统的输出响应，并将输入输出信号对比显示。G(s)=(s+2)/(s^2+10s+1)

解：num = [1,2];den = [1 10 1];

[numc ,denc]=cloop(num,den,-1);

v1 = [0:0.1:1];

v2 = [0.9:-0.1:-1];

v3 = [-0.9:0.1:0];

t = [0:0.1:4];

u = [v1,v2,v3];

[y,x] = lsim(numc,denc,u,t);

plot(t,y,t,u),grid

例9：求1+2+3+…+100＝？

解：for循环：

sum = 0;

for i = 1:100    

sum = sum + i; 

 end   

sum    

双重循环：

m = 3;

n = 3;

for i = 1:m   

  for j = 1:n 

   H(i,j) = 1/(i+j); 

  end   

end

while循环：

sum = 0;

i = 1;

while(i <= 100)  

  sum = sum +i;  

  i = i + 1;

end  

sum

例10：求1+2+3+…+N>10000的最小N值。

解：(1)while语句实现：

sum = 0;

i = 0;

while(sum <= 10000) 

 i = i + 1;  

sum = sum +i;

end  

[sum,i]  

(2)用for语句实现：

sum = 0;

for i = 1:10000 

   sum = sum + i;

   if sum >10000,

     break; 

   end  

end   

[sum,i]

例11：已知某一伺服机构的开环传递函数为G(s)=1/(s(0.5s+1)(0.1s+1)), 设计PID控制器，使得系统的加速度误差系数Ka>=10，幅值穿越频率wc>=4rad/s，相位裕度>=50^o 。

解：Ka = 10;

G = zpk([], [0 -2 -10], 20);

Ki = Ka/dcgain(G*tf([1 0], 1))       

wc = 4.1;

[num, den]= tfdata(G, ‘v’);

numc = ployval(num, j*wc);

denc = ployval(den, j*wc);

Gjwc = numc/denc;

theta = (-180 + 50.1)*pi/180;

Ejwc = cos(theta) + j*sin(theta); 

sum = Ejwc/Gjwc + j*Ki/wc;   

Kp = real(sum)     

Kd = imag(sum)/wc      

Gc = tf([Kd Kp Ki], [1 0]);

figure(1);bode(Gc), grid  

figure(2);bode(G, G*Gc), grid         

[kg, ph, wp, wc] = margin(G*Gc)

例12：试用分别用函数调用、子函数的方式画出所给两个函数的图形：y₁=10/(1+x²),y2=5+4sinx .

解：函数1：fun1.m：        

function y = fun1(x)   

 y = 10./(1+x.^2);

函数2：fun2.m：        

function y = fun2(x)   

 y = 5 +4.*sin(x);

调用程序：

t = 1:2:50;  

f1 = fun1(t); 

 f2 = fun2(t);  

plot(t,f1,t,f2)

例13：G(s)=(s²+10s+20)/(s³+15s²+10s+2)，G_c(s)=(10s+1)/s，H(s=1)。编写MATLAB程序，绘制系统的单位阶跃响应曲线。 

解：num = [1 10 20]; 

den = [1 15 10 2]; 

G = tf(num,den); 

Gc = tf([10 1],[1 0]);

GB = feedback(G*Gc,1,-1);

step(GB)

例14：已知系统开环传递函数：G(s)=20/(s⁴+8s³+36s²+40s)。求该系统在单位反馈下的（1）单位阶跃响应曲线（2）单位脉冲响应曲线

解：G = tf([20],[1 8 36 40 0]);

GB = feedback(G,1,-1);

figure(1),step(GB) 

figure(2),impulse(GB)

例15：已知Cummins电喷柴油机外特性数据，试绘制其扭矩特性曲线。

ne(rpm)	800	1000	1200	1400	1500	1700	1800	1900	2000	2200
Ne(kW)	57	88	113	136	139	144	145	145	145	144

解：ne=[800,1000,1200,1400,1500,1700,1800,1900,2000,2200];

Ne=[57,88,113,136 ,139,144,145,145,145,144];

 I=Ne./ne; 

plot(ne,I)

例16：已知系统开环传递函数：G(s)=5(s+1)(s+2)/((s+0.5)(s+3)(s+1+2j)(s+1-2j)),编写matlab程序，在三个图形窗口中分别绘制系统以下曲线：（1）单位反馈下的阶跃响应曲线；（2）开环Nyquist曲线；（3）开环Bode图

解：G=zpk([-1,-2],[-0.5,-3,-1-2j,-1+2j],5);

 sys=feedback(G,1,-1)

figure(1)  step(sys)   

figure(2)  nyquist(G)   

igure(3)   bode(G)

例17：n=10,用for循环求n！。

解：y=1;

for n=1:1:10  

   y=y*n  

end  

y

例18：在[0,50]区间产生101个等距采样数据，一维数组，两种方法实现。

解：a=0:0.5:50     

    a=linspace(0,50,101)