%Look at convergence of Power Method for nonsymmetric matrix A with real
%eigenvalues
%NOTE: FOR FUN, MODIFY THIS CODE TO DO A SYMMETRIC MATRIX!
%Power method with random initial guess, eigenvector initial guess
%Aitken's del^2 speedup with random initial guess 
%Plot and compare with theoretical convergence rates 
%(lambda1/lambda2)^k for plain power method
%(lambda1/lambda2)^k for Aitken del^2 speedup

clear all
close all
n = 10;
A = diag(1:n);   %Make a matrix with eigenvalues 1:n
S  = rand(n,n);
A = S\(A*S);   %do a random similarity transform (eigenvalues stay the same!) 
[V,D] = eig(A); %Find true eigenstuff
W = (sortrows([diag(D)';V]'))';  %organize e-vals and e-vecs in ascending order
V = W(2:end,:);  %colunms are eigenvectors
d = W(1,:);     %entries are eivenvalues   
maxiter = 100;
%%%%%%%%%%%%%%%%%
%Power Method

v0 = rand(n,1);  %random initial guess
v0 = v0/norm(v0);
for k = 1:maxiter
    w = A*v0;
    v = w/norm(w);
    lambda1(k) = v'*A*v;
    v0 = v;
end

%%%%%%%%%%%%%%%%%%%%%%
%Power Method again
v0 = V(:,4);  %Eigenvector (from one of the smaller eigenvalues) as initial guess
v0 = v0/norm(v0);

for k = 1:maxiter
    w = A*v0;
    v = w/norm(w);
    lambda2(k) = v'*A*v;
    v0 = v;
end
%%%%%%%%%%%%%%%%%%%%%
%Aitken's Del^2 Acceleration
v0 = rand(n,1);
v0 = v0/norm(v0);
lambda3(1) = lambda1(1);  %initialize algorithm with values from power iterations
lambda3(2) = lambda2(2);
for k = 3:maxiter+2
    w = A*v0;
    v = w/norm(w);
    lambda3(k) = v'*A*v;
    del_lambda = lambda3(k-1)-lambda3(k-2);
    del2_lambda = lambda3(k)-2*lambda3(k-1)+lambda3(k-2);
    lambdahat(k-2) = lambda3(k-2) - (del_lambda)^2/del2_lambda;
    v0 = v;
end

%%%%%%%%%%%%%%%%%
%Plot everything
semilogy(1:maxiter, (d(end-1)/d(end)).^([1:maxiter]),'k--', 'LineWidth',2)
hold on
semilogy(1:maxiter, (d(end-1)/d(end)).^(2*[1:maxiter]),'m--', 'LineWidth',2)
semilogy(1:maxiter, abs(lambda1-d(end)),'g')
semilogy(1:maxiter, abs(lambda2-d(end)),'r')
semilogy(1:maxiter, abs(lambdahat-d(end)),'b')
hleg = legend('(\lambda_2/\lambda_1)^k','(\lambda_2/\lambda_1)^{2k}','power method, random v_0','power method, aother e-vector as v_0','Aiktens \Delta^2, random v_0')
set(hleg, 'Location','SouthWest')
xlabel('k = iteration #')
ylabel('|\lambda_k-\lambda_{act}|')

