Matlab - powers of a matrix

%
% Graphical representation of A^K
%
% input: A ... 2x2 matrix 
%
% output: 
%
%   fig. 1 - before transformation:
%              unit circle and unit square (blue)
%              eigenvectors (red)
%              a few other vectors (different colors)
%   fig. 2 - after K repeated transformations 
%              - blue circle -- transformed unit circle
%              - transformed vectors (the same colors)
%              - black dashed circle -- spectral radius of A
%              - blue dashed circle -- the unit circle

close all
clear all
%                                                               
 A = [1 0.1; 0.3 1]; h=4; % real eigenv., ro > 1
% A = [ 0.9 0.6; 0 0.6]; h=1.5;   % real eigenv., ro < 1, norm > 1
% A = 0.3*[ 3 -1; -1 2]; h=2;  % symmetric transformation, ro > 1
% A = 0.8*[1 -1; 1 2]; h=10;   % imag. eigenv., ro > 1
% A = 0.8*[ 1 1; 0 1]; h=2;   % shear, ro < 1, norm > 1
% A = [ 1 1; 0 1]; h=8;   % shear, ro = 1
% A = 0.8*[1 0.1; 0.2 1]; h=1.5; % real eigenv., ro < 1, norm < 1

K = 12; % number of iterations
uc = 1; % plot of transf. unit circle (0 = no, 1 = yes)

% check if eigenvalues are real
if ( (A(1,1)-A(2,2))^2 + 4*A(2,1)*A(1,2) ) < 0
  ie = 0; % imaginary eigenvectors - will not be displayed 
  disp('eigenvectors are imaginary')
else
  ie = 1; % auxiliary value
end

% X - unit circle
t = linspace(0,2*pi,40);
X = [sin(t);cos(t)];
hold off
plot(X(1,:),X(2,:),'linewidth',2)
axis equal;
hold on;

% V - eigenvectors, E - eigenvalues of A
[V,E] = eig(A);

if (ie > 0)
  % plot the eigenvectors
  UP = [ [0;0] V(:,1),  [0;0] V(:,2) ];
  plot(UP(1,:),UP(2,:),'r','linewidth',2)
end

% a, b, c - different unit vectors given by 2 points
t = 1.5:2:5.5; % angles of the vectors
U = [cos(t);sin(t)]; % endpoints of the vectors
a = [ [0;0] U(:,1)];
b = [ [0;0] U(:,2)];
c = [ [0;0] U(:,3)];
plot(a(1,:), a(2,:),'linewidth',2);
plot(b(1,:), b(2,:), 'k','linewidth',2);
plot(c(1,:), c(2,:), 'g','linewidth',2);

U=X; % unit circle

figure(2);

for k=1:K

  % unit circle
  plot(U(1,:),U(2,:), '--','linewidth',2)
  axis([-h h -h h] ,'equal');
  hold on

  % spectral radius
  s = max(abs(diag(E)));
  plot(s*U(1,:),s*U(2,:), '--k','linewidth',2) % ':'

  % transformation of the unit circle
  if uc == 1
    X = A*X;
    plot(X(1,:),X(2,:),'linewidth',2)
  end

  if (ie > 0)
    % transformation of the eigenvectors
    UP = A*UP;
    plot(UP(1,:),UP(2,:),'r','linewidth',2)
  end

  % transformation of the other vectors
  a = A*a;  b = A*b;  c = A*c; 
  plot(a(1,:), a(2,:),'linewidth',2);
  plot(b(1,:), b(2,:), 'k','linewidth',2);
  plot(c(1,:), c(2,:), 'g','linewidth',2);
  hold off
  pause
end