q621.m

MATLAB及在电子信息课程中的应用,电子工业出版社出版 陈怀琛 吴大正 高西全合著

%《MATLAB及在电子信息课程中的应用》第六章例6.21程序q621
% 求状态转移矩阵及齐次状态方程解
% 电子工业出版社出版  陈怀琛 吴大正 高西全合著 2001年10月初版，2003年7月第二版
%
A=[-2,1;-17,-4];		        % 输入状态方程系数矩阵
x0=[3;4];				        % 输入初始条件
t=0:.02:3; Nt=length(t);		% 设定自变量数组并确定其长度
F=zeros(2,2,Nt);x=zeros(2,Nt);	% 状态转移矩阵F及状态变量初始化
for k=1:Nt							% 对时间循环
   F(:,:,k)=expm(A*t(k));		% 计算各时刻的状态转移矩阵F
end
z=reshape(F,[4,Nt]);			% 把F变为二维矩阵以便绘图
% 第一张图是系统状态转移矩阵,plot语句只接受二维变量
% z(1,:)=F(1,1,:),z(2,:)=F(2,1,:),z(3,:)=F(1,2,:),z(2,:)=F(2,2,:)
subplot(2,1,1),plot(t,z(1,:),'-.',t,z(2,:),':',t,z(3,:),'-',t,z(4,:),'--'), grid,	
set(gcf,'color','w')            % 设置图形背景色为白色
legend('F(1,1)','F(2,1)','F(1,2)','F(2,2)')
title('系统的状态转移矩阵')
for k=1:Nt						% 对时间循环，求各点状态变量
   % 矩阵乘法只能用于二维，因此对每一时刻的F,用squeeze函数缩去长度为1的第三维
   x(:,k)=squeeze(F(:,:,k))*x0;
end
% 第二张图是在给定初始条件下的输出
subplot(2,1,2),plot(t,x(1,:),'-',t,x(2,:),'--'),grid
legend('x(1,:)','x(2,:)')
title('系统输出状态变量')