% 
% Written by:
% -- 
% John L. Weatherwax                2009-02-24
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

close all; drawnow; 
clc; 
clear; 

F = [ 0.32857, -0.085714 ; 1.1429, 1.0714 ]; 
H = [ 10, 5 ]; 

% compute the eigenvalues of F (note all are less than 1 so the system is stable):
evalue = eig(F)

% the observability Grammian (note that this has one large eigenvalue and one very small one)
% this indicates that this system is not very observable since the ratio of small/large is very 
% close to zero.
% 
G = zeros(2); 
for j=1:2000, 
  G_add = ( H * F^j )' * ( H * F^j );
  if( mod(j,20)==0 ) fprintf('(%10d) norm G_add= %10.6f\n', j, norm(G_add) ); end; 
  G = G + G_add; 
end
G 
rank(G)
geig = eig(G)
min(geig)/max(geig)

% the observabilty matrix:
M=size(F,2);
obs_M = []; 
for j=0:(M-1)
  obs_M = [ obs_M; H * F^j ];
end
obs_M 
rank(obs_M)
obs_M_eig = eig(obs_M)
min(obs_M_eig)/max(obs_M_eig)