% % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
%                                                         %
%        SS 2007 - Numerik 2 Programmieraufgaben          %
%                                                         %
%                Aufgabe 4 - Potenzmethode                %   
%                                                         %                                               
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
clear all

c1=8;
c2=3;
c3=11;
A = [c1+c2 -c2 0; -c2 c2+c3 -c3; 0 -c3 c3];
x = [ 1; 1; 1];     % Startvektor
err=1;              % Fehler zwischen zwei Iterationen, zu Beginn "gross" gesetzt
epsilon = 10^-6;    % Genauigkeit


disp('Fuehre die Potenzmethode durch, um den betragsgroessten')
disp('Eigenwert lambda der Matrix A = ')
disp(' ')
disp(A)
disp(' ')
disp('zu berechnen. Waehle den Startvektor x^0=')
disp(x)
disp(' ')


% einige Initialisierungen: 
err=1;
iteration = 0;

[r,index]=max(abs(x));  % findet die im abs-Betrag groesste Komponente von x (r)
                   % und gibt auch die Zeile aus (index)

while err>epsilon
    
    x_old = x;
    r_old = r; 
    index_old = index; 
    if iteration >=1
        lam_old=lam;
    end
    
    x = A*x_old; 

    [r,index]=max(abs(x));
    
    if r==0
        error('Alle Komponenten von x sind 0!') ;
    else
        lam = x(index_old)/x_old(index_old);
        x=x/x(index);
    end
    
    
    iteration = iteration + 1;
    
    if iteration >= 2 & lam~=0
        err = abs(lam-lam_old);
    else 
        err = 1; 
    end
    
end


disp(['Der groesste Eigenwert von A wird approximiert durch lambda = ', num2str(lam),'.'])
disp(' ')
disp(['Der zugehoerige normierte Eigenvektor ist '])
disp(x/norm(x))
disp(' ')
disp(['Es wurden ', num2str(iteration),' Iterationen benoetigt, um eine Genauigkeit'])
disp('von 10^-6 zu erreichen. ')

% Test: 

[v_exact,lam_exact] = eig(A);
disp(' ')
disp(' ')
disp(['Der exakte groesste Eigenwert von A ist lambda = ', num2str(lam_exact(3,3)),'.'])
disp(['Der exakte zugehoerige normierte Eigenvektor ist '])
disp(v_exact(:,3))
disp(' ')
