%% Convergence of Rayleigh Quotient iteration

A=diag([4:0.2:10],0); 
[m,n]=size(A);
T=rand(n);
A=T*A*inv(T); % Similarity transform, expect quadratic convergence
A=A*A'; % Normal matrix, expect cubic convergence

lmax=max(eig(A))

%% Solve using Rayleigh Quotient iteration
tol=1e-8;
%t=cputime;[l,v,k,error] = ray_it(A,lmax+0.2,tol);tend=cputime-t;
t=cputime;[l,v,k,error] = ray_it(A,norm(A,inf),tol);tend=cputime-t;
disp('Rayleigh Quotient iteration');
disp(['Eigenvalue = ' sprintf('%0.5e',l) ' , number of iterations = ' sprintf('%d',k) ' , cpu time = ', sprintf('%0.5e',tend)]);
disp(' ');

t=cputime;[l,v,k,errpit] = pow_it_err(A,tol);tend=cputime-t;
disp('Power iteration');
disp(['Eigenvalue = ' sprintf('%0.5e',l) ' , number of iterations = ' sprintf('%d',k) ' , cpu time = ', sprintf('%0.5e',tend)]);
disp(' ');


loglog(error(1:end-1),error(2:end),'ro-'); hold on;grid on;
loglog(errpit(1:end-1),errpit(2:end),'r*-');
loglog(error(1:end-1),error(1:end-1).^3,'b-');
loglog(error(1:end-1),error(1:end-1).^2,'g-');
loglog(error(1:end-1),error(1:end-1).^1,'k-');
legend('Rayleigh iteration error','Power itation error','cubic','quadratic','linear');
xlabel('err_k');ylabel('err_{k+1}');