%% Einführung in die numerische Mathematik

%% Blatt 5

%% Bsp_6

clear; clc;

%%nicht relevant
%%beliebige SPD Testmatrix A  wird erstellt (ähnlich wie bei Blatt 4, Bsp 6.)
M=2;
N=2;
N1=N+1;
N1N1=N1*N1;
h=1/N;
tau=1/M;
r=2;
d=3;
D=diag([4*ones(N1N1-1,1);3])-diag(ones(N1N1-1,1),1)-diag(ones(N1N1-1,1),-1)-diag(ones(N1N1-N1,1),N1)-diag(ones(N1N1-N1,1),-N1);
for k=N1:N1:N1N1-N1  
    D(k:k+1,k)=[3,0];
    D(k,k+1)=0;
    F(k+1,k)=0;    
end
B=(d/(h^2))*D+r*eye(N1N1);
A=eye(N1N1)+tau*B;

n=N1N1;

%% Bsp 6

%%a)-----------------------------------------------------------------------


%% Transformation der SPD-Matrix auf tridiagonale Form

%% A SPD Testmatrix
%% n  Anzahl der Zeilen bzw. Spalten der Matrix A
%% ad, my, sigma Hilfsgrößen zur Housholder-Transformation, u Vektor zur
%% Householder Transformation P=I-u*u'/sigma
%% Q Matrix zur Speicherung des Produktes der orthogonalen
%% Transformationsmatrizen
%% H tridiagonale Form der Matrix A

H=A;
Q=eye(n);
    
for j=1:n-2
    
ad=H(j+1:1:n,j);

my=norm(ad);

e=[1;zeros(n-(j+1),1)];
 
u=ad+sign(ad(1))*my*e;

sigma=(1/2)*(norm(u)^2);
 
u=[zeros(j,1);u];

P=eye(n)-(u*u')/sigma;

Q=Q*P;
 
 %%Q=P*Q; Speicherung der Transformationsmatrizen
 H=P*H*P;

end

%%Ausgabe zur Überprüfung 

Tridiagonal_A=H;

norm(A-Q*H*Q')

%%b)-----------------------------------------------------------------------

%% Bestimmung der Eigenwerte einer tridiagonalen Matrix

%% A SPD Testmatrix, H tridiagonale Form der Matrix A
%% n  Anzahl der Zeilen bzw. Spalten der Matrix A
%% kmax maximale Anzahl der Iterationen zur Absicherung
%% a, b, x, y, g, c, s, d, e Analog zum Pseudo-Code der VO

for i=1:n
    
    a(i)=H(i,i);
    
end

for i=1:n-1
   
    b(i)=H(i,i+1);
end
    
kmax=100000;            %% maximale Anzahl der Iterationen(zur Absicherung)

    for k=1:kmax
        
       x(1)=a(1);
       y(1)=b(1);
       
       for j=1:n-1
           
           g(j)=sqrt(x(j)^2+b(j)^2);
           c(j)=x(j)/g(j);
           s(j)=b(j)/g(j);
           d(j)=c(j)*y(j)+s(j)*a(j+1);
           x(j+1)=-s(j)*y(j)+c(j)*a(j+1);
           
           if j <n-1
               e(j)=s(j)*b(j+1);
               y(j+1)=c(j)*b(j+1);       
           end
       end
       g(n)=x(n);
       
       
       a(1)=s(1)*d(1)+c(1)*g(1);
       b(1)=s(1)*g(2);
       
       for j=2:n-1
           a(j)=s(j)*d(j)+c(j)*c(j-1)*g(j);
           b(j)=s(j)*g(j+1);
           
       end
       
       a(n)=c(n-1)*g(n);
       
       if norm(b)<10e-23
           break
       end
    end
       
%% Ausgabe der Eigenwerte

Eigenwerte_A=a;

norm(sort(eig(H))-sort(a)')
