% Ein Skript zum Testen verschiedener iterativer Solver
% fuer die tridiagonale  (-1 2 -1)-Matrix.
% 
% Die Experimente zeigen:
% - Anzahl der Iterationen beim Jacobi-Verfahren waechst ~ 1/h^2
% - Gauss-Seidel konvergiert etwa doppelt so schnell wie Jacobi
% - Anzahl der Iterationen beim SOR-Verfahren mit 
%   optimal gewaehltem omega waechst ~ 1/h 
% - SOR funktioniert nur sinnvoll fuer 1<omega<2, wobei der 
%   optimale omega-Wert mit zunehmendem n immer knapper bei 2
%   liegt.
% - geringfuegige Unterschreitung des omega_opt ist weniger kritisch
%   als Ueberschreitung.
%
% Diese Ergebnisse sind qualitativ auch fuer andere Gleichungssysteme mit schwach 
% diagonaldominanter Matrix typisch.


n = 200;
A = zeros(n);
for i=1:n
    A(i,i) = 2;
end
for i=2:n
    A(i-1,i) = -1;
    A(i,i-1) = -1;
end
    
b = ones(n,1);

% Das ist die direkte Loesung
% mit eingebautem Matlab-Befehl
%x = A\b;
%plot(x)

%Iterative Lösung: 
x = b;
xnew = x;
itmax = 100000;
epsilon = 1.e-5;

% theoretisch optimales omega fuer diese Beispiel-Matrix :
%omega_opt = 2/(1 + sqrt(1-cos(pi/(n+1))^2))
% fuer grosse n ist das annaehernd
omega_opt = 2/(1 + pi/(n+1))  
% Experimente zeigen trotzdem: geringste It-Anzahl bei leicht groesseren
% Werten....

omega = 1.975

for it = 1:itmax
    xnew(1) = (b(1)+ xnew(2))/2;
    for i=2:n-1
        %xnew(i) = (b(i)+ x(i-1)+ x(i+1))/2; % Jacobi
        %xnew(i) = (b(i)+ xnew(i-1)+ xnew(i+1))/2; % Gauss-Seidel
        tmp = (b(i)+ xnew(i-1)+ xnew(i+1))/2; % Gauss-Seidel 
        xnew(i) = omega * tmp + (1-omega)* x(i);
    end
    xnew(n) = (b(n)+ xnew(n-1))/2;
    diff = norm(x-xnew,'inf');
    x = xnew;
    if(diff<epsilon) 
        it
        break; 
    end
    
end
    