Matlab - powers of a matrix

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

close all
clear all

  % choose one of the matrices (h is graph window size):
                                                 
 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.8*[1 0.1; 0.2 1]; 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

K = 10; % number of iterations
uc = 1; % plot of transformed unit circle (0 = no, 1 = yes)

  % norms of the matrix
normc = num2str(norm(A,1));
normr = num2str(norm(A,inf));
norms = num2str(norm(A,2));
normfro = num2str(norm(A,'fro'));
disp(['norms - col, row, spe, fro: ' normc ', ' normr ', ' norms ', ' normfro]);

  % 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);
realV = isreal(V);
if realV
    % plot the eigenvectors
  UP = [ [0;0] V(:,1),  [0;0] V(:,2) ];
  plot(UP(1,:),UP(2,:),'r','linewidth',2)
else
  disp('eigenvectors are imaginary')
end

  % rho / spectral radius of A
rho = max(abs(diag(E)));
disp(['spect. radius: ', num2str(rho)])

% 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
    % transformation of the eigenvectors
  if realV
    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