  function Demo1
% This script looks at the convergence of the power method for the
% matrix M given that the start vector has centroid zero.
  
  clc
  close all
  n  = 7;

% Construct M...
  I = eye(n,n);
  S = [I(:,n) I(:,1:n-1)];
  M = I + S;
  
% and its schur decomp ordering eigenvalues from big to little: 
  [Q,T] = schur(M);
  format short
  [d,idx] = sort(-diag(T));
  Q = Q(:,idx)
  T = T(idx,idx)
% Notice how T is block diagonal. The complex conjugate eigenvalues 
% come in pairs. The first column of Q is a multiple of ones(n,1) and
% is the dominant eigenvector. Note that if x has centroid zero then
% x'*ones(n,1) = sum(x) = 0. Thus, the dominant eigenvector is "out of the
% ppicture if we start the power method with a zero centroid vector.
%
% Note that span{Q(:,2),Q(:,3)}, span{Q(:,4),Q(:,5)}, span{Q(:,6),Q(:,7)}
% are invariant associated with
%
%     lambda2 = 1+c1 + i*s1   c1 = cos(2pi/n), s1 = sin(2pi/n)
%     lambda3 = 1+c1 - i*s1

%     lambda4 = 1+c2 + i*s2   c2 = cos(4pi/n), s2 = sin(4pi/n)
%     lambda5 = 1+c2 - i*s2
% and
%     lambda6 = 1+c3 + i*s3   c3 = cos(6pi/n), s3 = sin(6pi/n)
%     lambda7 = 1+c3 - i*s3

% Let's look at the power method iterates in terms of the basis
% defined by the columns of Q

% Starting vector...
  x = randn(n,1); x = x-sum(x)/n; x = x/norm(x);
 
  alfaMatrix = [];
  for k=1:20
%    Solve Q*alfa = x to express x in terms of Q-columns
     alfa = Q'*x;
     alfaMatrix = [alfaMatrix;alfa'];
     x = M*x;
     x = x/norm(x);
  end
  subplot(2,1,1)
  plot(alfaMatrix(:,2:3))
  title('Coefficients of Q(:,2) and Q(:,3)')
  subplot(2,1,2)
  plot(alfaMatrix(:,[1 4:7]))
  title('Coefficients of Q(:,1), Q(:,4), Q(:,5), Q(:,6), and Q(:,7)')
  shg
  pause(1)
  
% The rate of convergence should resemble |lambda4/lambda2|^k which can be
% shown to equal ( (1+cos(4pi/n))/(1+cos(2pi/n)) )^k/2. Let's look at the
% critical ratio tau  = sqrt( (1+cos(4pi/n))/(1+cos(2pi/n))  ) as a function of n:

  figure
  n=5:50; z = sqrt( (1+cos(4*pi./n))./(1+cos(2*pi./n))); plot(n,z);
  ylabel('Convergence Factor tau');xlabel('n')
  title('| lambda4 |/| lambda2 | as a function of n')
  shg 
  pause(1)
  
% Noting that the convergence factor gets closer and closer to one as n
% increases, we ask, how big must k be so that tau^k < .001? This is a good
% tolerance value given screen granularity and human perception
 
  kmin = [];
  for nVal=5:50
      tau = sqrt( (1+cos(4*pi/nVal))/(1+cos(2*pi/nVal)));
      k=1;
      while tau^k>.001
          k=k+1;
      end
      kmin = [kmin k];
  end
  figure
  plot(n,kmin)
  ylabel('kmin'); xlabel('n'); title('Number Power Iterationss to Get O(.001) Error')
  shg
  pause(1)

% The power method iterates converge to span{Q(:,2),Q(:,3)} =
% span{cVec,sVec} where 
%               cVec = cos((0:n-1)*(2pi/n))' 
%               sVec = sin((0:n-1)*(2pi/n))'


% Let's look at a random power method example, reporting how the minimum
% value of norm( [cVec svec]*a - x) goes to zero where x is the kth iterate vector.
% This quantity measures how far x is to the invariant subspace
% span{cVec,sVec}.

  nVal = 30; itMax = 1000;
  x = randn(nVal,1); 
  x = x-(sum(x)/nVal);
  x = x/norm(x);
  [cVec sVec] = csVecs(nVal); 
  A = [cVec sVec];
  tau = sqrt( (1+cos(4*pi/nVal))/(1+cos(2*pi/nVal)));
  predErr = []; actualErr = [];
  for k=1:itMax
      x = Mprod(x); 
 % The following statement ensures that the iterates stay orthogonal
 % to M's dominant eigenvector. In exact arithmetic it is not necessary.
 % But for very large values of itMax you will notice that the dominant
 % eigenvector of all ones begins to creep in without this enforcement of
 % orthogonality.
      x = x-sum(x)/nVal;     
      x = x/norm(x);
      predErr = [predErr tau^k];
      actualErr = [actualErr norm(A*(A\x)-x)];
  end
  figure
  semilogy(1:itMax,predErr,1:itMax,actualErr,'r')
  legend('O(Predicted Error)','O(Actual Error)')
  title(sprintf('n = %1d',nVal))
  xlabel('Iterations')
  shg 
  
