%   A matlab script to illustrate conditioning of eigenvalue computation.
%   Written for the course Scientific Computing, homework 4, see 
%   http://www.math.nyu.edu/faculty/goodman/teaching/SciComp2002/


clear               % So leftover matrices don't get in the way.
startn = 5;         % The smallest matrix, the smallest n.
stride = 5;         % The difference between n values.  I put this in 
                    % in case your computer is too slow.  Use stride=1
                    % unless you can't stand the wait.
nMax   = 80;        % The largest matrix.
run    = 1;         % Keep track of which run, for indexing arrays of results.

for n=startn:stride:nMax  % Loop over the runs, go through n values by stride

              % Clear matrices between runs so the dimensions will be right.
   clear A    % The matrix whose eigenvalues we seek.
   clear B    % The same matrix conjugated by the DFT matrix
   clear AR   % The eigenvectors of A.
   clear BR   % The eigenvectors of B, would be related in exact arithmetic.
   clear Q    % The DFT matrix.
   clear lA   % The eigenvalues of A.
   clear lB   % The eigenvalues of B, would be equal in exact arithmetic.
   clear D
   z = exp(2*i*pi/n);  % The primative root of unity for the DFT matrix
   Q=zeros(n);         % Get memory for an array of the correct size.
   nf = 1/sqrt(n);     % The normalization constant.
   for j=1:n           % Construct Q.
      w = z^(j-1);     % The root of unity for row j.
      for k=1:n
         Q(j,k) = nf*w^(k-1);
       end
    end
   A = zeros(n);       % Construct the matrix A.
   for j=1:n
      A(j,j) = .25;    % .25 on the diagonal
    end
   for j=1:n-1
      A(j,j+1) = .25;  % .25 above the diagonal
      A(j+1,j) = .5;   % .5 below the diagonal.
    end
   [AR,D] = eig(A);    % Compute the eigenvalues and eigenvectors of A.
   for j=1:n
     lA(j) = D(j,j);   % Copy the eigenvalues to a one index array.
   end
   B=Q'*A*Q;           % Get B by conjugating by Q.
   [BR,D] = eig(B);
   for j=1:n
     lB(j) = D(j,j);
   end
   me   = max( abs( imag (lB) ) );   % The eigenvalues of A are real.  Imaginary
                                     % parts are a symptom of roundoff.
   mplot(run) = me;                  % Store in an array for plotting.
   c   = cond(AR);                   % Compute the condition number of AR.
   cplot(run) = c;
   nplot(run) = n;
   run = run + 1;
   sprintf('n = %5d,  cond = %12.4e,   maxImag = %12.4e',n,c,me)
end
figure(1)
grid
semilogy(nplot,cplot)
title('something here')
grid
title('something different here')
figure(2)
semilogy(nplot,mplot)
grid
title('something different here')