%%
clear
clc
rng(1)

%%
% We set up a 5x5 matrix with prescribed eigenvalues, then apply
% the power iteration.
lambda = [1 -0.75 0.6 -0.4 0];
A = triu(ones(5),1) + diag(lambda)   % triangular matrix, eigs on diagonal. 
                                     % triu(X, k) are the elements on and 
                                     % above the K-th diagonal of X. K = 0
                                     % is the main diagonal, K > 0 is above
                                     % the main diagonal and K < 0 is below 
                                     % the main diagonal.

%%
% We run the power iteration 60 times. The best estimate of the dominant
% eigenvalue is the last entry of gamma.
np = 60;
[gamma, x] = poweriter(A, np);
eigval = gamma(np)

%%
% We check linear convergence using a log-linear plot of the error. We use
% our best estimate in order to compute the error at each step.
err = eigval - gamma;
semilogy(1:np, abs(err), '.-')
title('Convergence of power iteration')
xlabel('k')
ylabel('|\lambda_1 - \gamma_k|')
grid on

%%
% The trend is clearly a straight line asymptotically. We can get a refined
% estimate of the error reduction in each step by using the exact
% eigenvalues.
hold on
theor = abs(lambda(2)/lambda(1))
semilogy(1:np, theor.^(1:np), '--')