clear
rng(1)
format long

%%
% We set up a 5x5 triangular matrix with prescribed eigenvalues on
% its diagonal.
lambda = [1 -0.75 0.6 -0.4 0];
A = triu(ones(5),1) + diag(lambda); 

%%
% We run inverse iteration with the shift s=0.7 and take the final
% estimate as our ``exact'' answer to observe the convergence. 
np = 30;
[gamma, x] = inviter(A, 0.7, np);
eigval = gamma(end)

%%
% As expected, the eigenvalue that was found is the one closest to 0.7.
% The convergence is again linear (i.e exponential)
err = eigval - gamma;
semilogy(abs(err),'.-')
title('Convergence of inverse iteration')
xlabel('k')
ylabel('|\lambda_1 - \gamma_k|')
grid on

%%
% In the numbering of this example, the eigenvalue closest to s=0.7 is
% lambda_3 and the next-closest is lambda_1.
theoretical_rate = abs((lambda(3)-0.7) / (lambda(1)-0.7))
hold on
semilogy(1:np, theoretical_rate.^(1:np), '--')
legend('numerical experiment', 'theory', 'location', 'best')