clear all, close all

format compact
format long

A = laplacian2d(8);
[n, ~] = size(A);

tol = 1.e-14;
rng(1);
x0 = rand(n, 1);

maxiter = 400;

[lambda, x, conv1] = mypower(A, x0, tol, maxiter);

abs(lambda)
eigen = max(abs(eig(A)))
abserr = abs( eigen - abs(lambda) )

niter1 = length(conv1);
nps1 = 1:niter1;

lams = eigs(A, 2);
sigma = 2*log(abs(lams(2)/lams(1))); % theoretical convergence rate

semilogy(nps1, conv1, 'r.-');
grid on
hold on
semilogy(nps1, exp(sigma*nps1), '--');
legend('Power iterations', 'Theory')
xlabel('Iteration number');
ylabel('|Delta|');
title('Convergence of power iteration process');
