%%
clear
clc

%%
% We construct a symmetric matrix with a known EVD.
n = 16;
lambda = (1:n)';  
D = diag(lambda);
[V, ~] = qr(randn(n));   % get a random orthogonal V
A = V*D*V';

%%
% The Rayleigh quotient of an eigenvector is its eigenvalue.
format long
R = @(x) (x'*A*x)/(x'*x);
R(V(:, 7))

%% 
% The Rayleigh quotient's value is much closer to an eigenvalue than its
% input is to an eigenvector. In this experiment, each additional digit of
% accuracy in the eigenvector estimate gives two more digits to the
% eigenvalue estimate.
seed = 2;
rng(seed);
np = 6;
delta = (0.1).^(1:np)';
quotient = 0*delta;
for k = 1:np
    e = randn(n,1); 
    e = e/norm(e);
    x = V(:,7) + delta(k)*e;
    quotient(k) = R(x);
end

loglog(delta,abs(quotient - D(7,7)), '.-')
hold on
loglog(delta, delta.^2, '--')
grid on
xlabel('\delta')
ylabel('error of RQ')
title('Accuracy of RQ vs accuracy of eigenvector')
legend('Nemerical experiment', '\delta^2', 'location', 'best')