function [x] = demo_jacobi_sor(matrix,omega_jac,omega_sor) 

% demo_jacobi_sor ist zur Demonstration der Arbeitsweise
% des Jacobi-Verfahrens und des SOR-Verfahrens in der
% Vorlesung
% Aufruf demo_jacobi_sor(matrix,omega_jac,omega_sor)
% matrix == 0 : Zufallsmatrix
% sonst, diagonal dominante Matrix
% omega_jac: Dämpfung im Jacobi-Verfahren
% omega_sor: Relaxation im SOR-Verfahren

n = 20;
b = rand(n,1);
A=rand(n,n);
if matrix ~= 0
   A = A + (n-1)*eye(20,20);
end
A
pause()
disp('Lösung ');
x_gauss = A\b
pause()

% jacobi l"ost das System
%   A x = b
% mittels der ged"ampften Jacobi-Iteration
% (falls das Verfahren konvergiert)
% Der Aufruf ist [x] = jacobi(A,b,omega)
% A - (n x n) Matrix
% b - n-dimensionaler Vektor rechte Seite
% omega aus (0,1] - D"ampfungsparameter 
% als Startwert wird der Nullvektor genommen

% kontrolliere, ob alle Dimensionen passen
[n,m] = size(A);
if n ~= m
  error('A ist keine quadratische Matrix')
end
if n ~= size(b)
  error('Die Dimension von A und b passen nicht')
end
%
omega = omega_jac;

% maximale Anzahl der Iterationen
maxit = 100;

% Startvektor
x_alt = zeros(n,1);
x_neu = x_alt;

% Abbruchkriterium
eps = 1e-8;

% die Iteration
for i= 0:maxit
  % Residuum
  r = b - A*x_alt;
  % berechne neue Iterierte
  for j = 1:n
    x_neu(j) = x_alt(j) + omega*r(j)/A(j,j);
  end

  norm_aenderung = norm(x_neu-x_alt);
  s = sprintf('Jacobi-Iteration %d  Aenderung der Iterierten %g',...
               i, norm_aenderung);
  disp(s);
  % falls die Aenderung der Iterierten klein genug, dann breche ab
  if (norm_aenderung < eps)
    break;
  end

  % speichere Iterierte f"ur n"achtste Iteration um
  x_alt = x_neu;

end

% vergleiche iterative L"osung mit direkter L"osung
fehler = x_gauss - x_neu;
s = sprintf('Fehler zur L"osung mit Gauss-Verfahren %g',...
             norm(fehler));
disp(s);
pause()

% nun SOR-Verfahren
omega = omega_sor;

% sor l"ost das System
%   A x = b
% mittels der ged"ampften sor-Iteration
% (falls das Verfahren konvergiert)
% Der Aufruf ist [x] = sor(A,b,omega)
% A - (n x n) Matrix
% b - n-dimensionaler Vektor rechte Seite
% omega aus (0,2) - D"ampfungsparameter 
% als Startwert wird der Nullvektor genommen

% maximale Anzahl der Iterationen
maxit = 100;

% Startvektor
x_alt = zeros(n,1);

% Abbruchkriterium
eps = 1e-8;

% die Iteration
for i= 0:maxit
  % berechne neue Iterierte
  x_neu = x_alt;
  for j = 1:n
    s = b(j) - A(j,1:n) * x_neu;
    x_neu(j) = x_neu(j) + omega*s/A(j,j);
  end
  norm_aenderung = norm(x_neu-x_alt);
  s = sprintf('SOR(%g) Iteration %d  Aenderung der Iterierten %g',...
               omega, i, norm_aenderung);
  disp(s);
  % falls die Aenderung der Iterierten klein genug, dann breche ab
  if (norm_aenderung < eps)
    break;
  end

  % speichere Iterierte f"ur n"achtste Iteration um
  x_alt = x_neu;
  
end

% vergleiche iterative L"osung mit direkter L"osung
fehler = x_gauss - x_neu;
s = sprintf('Fehler zur L"osung mit Gauss-Verfahren %g',...
             norm(fehler));
disp(s);

  









