% % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
%                                                         %
%        SS 2007 - Numerik 2 Programmieraufgaben          %
%                                                         %
%                Aufgabe 8 - Maehly-Alg.                  %   
%                                                         % 
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % %

% Hinweis: hier werden nur die ersten beiden Eigenwerte der Matrix A 
% berechnet und dann die Aufgabe 5c ausgenutzt. 


clear all

A = [ 16 1 0 0 0 ; 1 12 1 0 0 ; 0 1 8 1 0 ; 0 0 1 4 1; 0 0 0 1 0];
[n,m] = size(A); 

if n ~= m 
    error('Matrix ist nicht quadratisch!')
end

disp(' ')
disp('Berechnung der Eigenwerte der Matrix')
disp(A)
disp(' ')
disp('mit Hilfe des Maehly Algorithmus. ')
disp(' ')
disp(' ')


% -------------------------------
% einige Initialisierungen vorab: 
% -------------------------------

% Nutze Bezeichnungen aus der Vorlesung: 
delta = diag(A);            % Hauptdiagonale
gamma(2:n) = diag(A,-1);    % untere Nebendiagonale

% Angestrebte Genauigkeit des Newton-Verfahrens:
epsilon = 10^-8;

% Vektor, der spaeter die Eigenwerte beinhalten soll: 
z = zeros(n,1);

% Berechnung eines Startwerts lambda^(0), als Schaetzung fuer den ersten 
% Eigenwert:
lambda=norm(A,inf);

% Anzahl der benoetigten Iterationen: 
it = 0; 

% ------------------------------------
% Berechnung der fehlenden Eigenwerte: 
% ------------------------------------

for j=1:2     % Schleife ueber die ersten beiden Eigenwerte: 
    
    err=1;      % zu Beginn muss der Fehler err groesser als epsilon sein
    
    % Der erste Eigenwert wird mit dem normalen Newtonverfahren berechnet: 
    if j==1
         while err>epsilon
              lambda_old=lambda;
              
              % Berechnung der Werte p(lambda), p'(lambda)
              [p,p_prime]=charPoly(lambda,delta,gamma,n);
              lambda = lambda -p/(p_prime);
              err = abs(lambda-lambda_old);  
              it = it+1; 
         end
         z(j)=lambda;
    else
        lambda = lambda-1;
        while err>epsilon

            lambda_old=lambda;
            % Berechnung der Werte p(lambda) (=p), p'(lambda) (=p_prime):
            [p,p_prime]=charPoly(lambda,delta,gamma,n);

            hp=0;
            for k=1:j-1
                hp = hp + p/(lambda-z(k));       % siehe Formel (11) der Datei maehly.pdf
            end
            
            lambda = lambda -p/(p_prime-hp);     % siehe Formel (11) der Datei maehly.pdf
            err = abs(lambda-lambda_old);        % Fehlerberechnung
            it = it+1; 

        end
        z(j)=lambda;
    end
end

% Es gilt delta(i)+delta(n+1+i)=2c=16, daraus folgt c=8. Die Eigenwerte
% liegen somit symmetrisch um c, da n=5 ungerade ist, ist c auch ein
% Eigenwert von A. Ist lambda ein Eigenwert von A, so auch 2c - lambda.
% Die verbleibenden Eigenwerte koennen daher berechnet werden: 

z(3) = 8;
z(4) = 2*8-z(2);
z(5) = 2*8-z(1);


% sortiere die Eigenwerte in absteigender Reihenfolge
% (fuer den Vergleich mit den exakten Eigenwerten).
% z_neu ist also der Vektor, der die Eigenwerte enthaelt. 
reihenfolge = 'ascend';
z_neu=sort(z,1,reihenfolge);

% ------------------------------------------
% Berechnung der zugehoerigen Eigenvektoren: 
% ------------------------------------------
EV=zeros(n,n);          % Matrix, die die Eigenvektoren spaltenweise 
                        % enthalten soll.
for i=1:n
    q=zeros(n,1);       % Eigenvektor zum Eigenwert z_neu(i)
    q(1)=1;             % siehe Formel (5) der Datei maehly.pdf
    prod = gamma(2);
    for k=1:n-1
        % Berechnung des Eigenvekors q zum Eigenwert z_neu(i)
        % (siehe Formel (5) der Datei maehly.pdf):
        prod = prod*gamma(k+1);
        q(k+1)=(-1)^k*charPoly(z_neu(i),delta,gamma,k)/prod; 
    end
    % Speichern des normierten Eigenvektors (Normierungs nur zum Vergleich
    % mit den exakten Eigenvektoren, welche mit Matlab berechnet werden
    % koennen):
    EV(:,i)=q/norm(q);
end


% Berechnung der exakten Eigenvektoren v und Eigenwerte l der Matrix A:
% ---------------------------------------------------------------------
[v,l] = eig(A);


disp('Die exakten Eigenwerte der Matrix A sind:')
disp(diag(l))

disp('Die approiximierten Eigenwerte der Matrix A sind:')
disp(z_neu)


disp('Die exakten Eigenvektoren der Matrix A sind (spaltenweise):')
disp(v)

disp('Die approximierten Eigenvektoren der Matrix A sind (spaltenweise):')
disp(EV)

disp(['Dabei wurden ',num2str(it),' Iterationen benoetigt, um die Genauigkeit'])

disp(['epsilon = ',num2str(epsilon),' zu erreichen.'])
