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