function y = powers(A,v,n)
%POWERS  Execute the power method on a 2 by 2 matrix A
%
%       y = powers(A,v,n) executes the power method on A (2 by 2) 
%       starting from the vector v.  The number of iterations is n.
%	For examples that use power see expower

% Written by Ismail Khalil Bustany                            07/20/93
%
% Function :  powers(A,v,n)
%
% Purpose  :  Computes a candidate for the largest eigenvalue (in 
%             magnitude) of a 2 by 2 matrix. The sequence 
%             of vectors {v(i)} is plotted. The input vector v is 
%	      v(1) and then  v(i+1) = Av(i).
%
% Arguments:  A  -- 2 by 2 matrix
%             v  -- 2 by 1 vector
%             n  -- number of iterations 
%
% Returns  :  A candidate for the largest eigenvalue (in magnitude).
%
% Comments :  1. Contrast the speed of convergence with the inverse power 
%                method, where A inverse is used instead of A.
%             2. Explain the pictures obtained with the following:
%
%               a. Try A = [1.0000000e+00   0.0000000e+00
%                           5.0000000e-03   1.0500000e+00]
%                      v = [1 0]';
%
%               b. Try A = [1.0000000e+00   0.0000000e+00
%                           5.0000000e-03  -1.0500000e+00]
%                      v = [1 0]';
%
%               c. Try [0 1; 1 0], v =  [1 1] 
%             
% Credits  :  Idea contributed by Carlos Kirjner Neto (UC Berkeley). 
%             Edited by Tom Bryan 24 August 1993. 
%-----------------------------------------------------------------------
axis('square');
scale = 1.15;
axis([-scale scale -scale scale]);
clf;
xlabel('');
ylabel('');
history(:,1) = v;
v = v/sqrt(v'*v);
R1 = [cos(5*pi/6) -sin(5*pi/6);  sin(5*pi/6) cos(5*pi/6)];
R2 = [cos(5*pi/6)  sin(5*pi/6); -sin(5*pi/6) cos(5*pi/6)];
tip1 = .1*R1*v+v;
tip2 = .1*R2*v+v;
plot([0 0], [-scale scale], '--', [-scale scale],[0 0],'--',...
   [0 v(1)], [0 v(2)], '-', [v(1) tip1(1)], [v(2) tip1(2)], '-',...
                              [v(1) tip2(1)], [v(2) tip2(2)], '-');
xlabel('v(1)')
ylabel('v(2)')
title('Sequence of vectors generated by power method')
hold on 
fprintf('%4.0f.  <v,Av>/<v,v> = %e\n',0,v'*A*v/(v'*v));

text(v(1),v(2),'0')

for i = 1:n 
   v = A*v;
   v = v/norm(v);
   history(:,i) = v;
   tip1 = .1*R1*v+v;
   tip2 = .1*R2*v+v;
   fprintf('%4.0f.  <v,Av>/<v,v> = %e\n',i,v'*A*v/(v'*v));
   plot([0 0], [-scale scale], '--', [-scale scale],[0 0],'--',...
      [0 v(1)], [0 v(2)], '-', [v(1) tip1(1)], [v(2) tip1(2)], '-',...
                                  [v(1) tip2(1)], [v(2) tip2(2)], '-');   
   text(v(1),v(2),num2str(i))   
end %for 
% plot(history(1,:),history(2,:),'+');
y = v'*A*v/(v'*v);
fprintf('Candidate for largest eigenvalue (in magnitude) = %e\n',y);
hold off