%% Dokumentaion des prizipiellen Ablaufs eines Lanczos Verfahrens
% wenig geeignet fuer "echte Anwendungen"
% Besser z. Bsp. MATLAB Funktion eigs direkt auf Hamiltonmatrix anwenden

% W. Schweizer 
% Vers. 3.2 05/2016 (Vers. 1 03/2000)
%
% Simulation physikalischer Systeme - Computational Physics mit MATLAB 
% 2016

clear, close all, clc

%% Berechnung der Testmatrix
[E,V,Emat] = anharmoscmatrix(100,0.1,0.1,1,3);
%% Startvektor
nv = length(Emat);
psi = zeros(nv,1);
psi(1,1) = 1; 
%% psi-2
zwi = Emat*psi(:,1);
psi(:,2) = (zwi -zwi'*psi(:,1)*psi(:,1))/sqrt(zwi'*zwi - (psi(:,1)'*zwi)^2);
HLan=psi'*Emat*psi;
%%
for n=3:60
    zwi = Emat*psi(:,n-1);
    psi(:,n) = zwi - (psi(:,n-2)'*zwi)*psi(:,n-2) - (psi(:,n-1)'*zwi)*psi(:,n-1);
    psin = sqrt(zwi'*zwi - (psi(:,n-2)'*zwi)^2 - (psi(:,n-1)'*zwi)^2);
    psi(:,n) = psi(:,n)/psin;
    HLan(n,n) = psi(:,n)'*Emat*psi(:,n);
    HLan(n-1,n) = psi(:,n-1)'*Emat*psi(:,n);
    HLan(n,n-1) = HLan(n-1,n);
end
eigs(HLan,2,'sm')
%% Visualisierung
for nv = 3:length(HLan)
    eiwert = eigs(HLan(1:nv,1:nv),3,'sm');
    eiw1(nv) = eiwert(3);
    eiw2(nv) = eiwert(2);
    eiw3(nv) = eiwert(1);
end
xn = 3:length(HLan);
figure
plot(xn,eiw1(3:end),'k'), hold on
plot(xn,eiw2(3:end),'k')
plot(xn,eiw3(3:end),'k')  
%% Orthogonalitaet - Visualisierung
Otest = NaN*ones(60,60);
for ni=1:59
    for nj=ni+1:60
        Otest(ni,nj)=abs(psi(:,ni)'*psi(:,nj));
    end
end
figure,semilogy(Otest,'.k'),shg
%% Orthogonalit"at Tridiagonalelemente
Otesttri = NaN*ones(60,60);
for ni=1:59
    for nj=ni+1
        Otesttri(ni,nj)=abs(psi(:,ni)'*psi(:,nj));
        Otesttri(nj,ni)=Otesttri(ni,nj);
    end
end   