% Matlab code powerplot.m % For "Applied Numerical Linear Algebra", Algorithm 4.1 % Written by James Demmel, Oct 23, 1995 % Modified, Jun 5, 1997 % % Plot convergence of the power method, Algorithm 4.1 % function [x,vv,verror,eerror] = powerplot(a,v,m) % % Inputs: % a = matrix % v = starting vector % m = number of steps % % Outputs: % Plot of actual and predicted error versus iteration % x = result of 30+m steps of power method % vv = matrix of all intermediate results of power method % verror = difference between x and each column of vv % eerror = difference between x'*a*x and each vv'*a*vv % function [x,vv,verror,eerror] = powerplot(a,v,m) vv = v/norm(v); for i = 1:m vv(:,i+1)=a*vv(:,i); vv(:,i+1)= vv(:,i+1)/norm(vv(:,i+1)); end x=vv(:,m+1); for i = 1:2*m+30, x = a*x; x = x/norm(x); end vv=vv*diag(sign(vv(1,:))); x=x*sign(x(1)); % compute error = angle between vectors % (note that it does not converge to macheps!) % verror = acos(x'*vv); % compute error = sine of angle between vectors % (note that it does not converge to macheps!) % cosine = x'*vv; % verror = sqrt((1- cosine).*(1+cosine)) ; % compute error = norm of difference of vectors verror = (max(abs(vv-x*ones(1,m+1))))'; figure(1) hold off, clf % plot error in approximate eigenvectors subplot(211) semilogy(verror,'-b') title('Eigenvector error, blue = actual, red = predicted') xlabel('Iteration number') grid hold on e=sort(abs(eig(a))); [n,n]=size(a); rat=e(n-1)/e(n); rat=(rat*ones(1,m+1)).^(0:m); semilogy(rat,'-r') % % plot error in approximate eigenvalues subplot(212) eerror = abs( ( diag(vv'*a*vv) - x'*a*x ) / (x'*a*x) ); semilogy(eerror,'-b') title('Eigenvalue error, blue = actual, red = predicted') xlabel('Iteration number') grid hold on semilogy(rat,'-r')