%《MATLAB在電子信息課程中的應用》第六章例6.21程序q621
% 求狀態轉移矩陣及齊次狀態方程解
% 電子工業出版社出版  陳懷琛 吳大正 高西全合著 2001年10月

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,	
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('系統輸出狀態變量')