% 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')