Das Ausgangsprobleme Ax = b, A ∈ ℝnxn regulär, b ∈ ℝn wird mit einer regulären Matrix B ∈ ℝnxn als Bx+(A-B)x=b geschrieben und das, aufgrund der Regularität von B, zum Fixpunktproblem xk+1=(I-B-1A)xk+B-1b umgeformt. Man kann zeigen, dass das Verfahren gegen den Fixpunkt x* genau dann konvergiert, wenn ρ(I-B-1A)<1, mit ρ(C) bezeichnet den Spektralradius von C, bzw. wenn ∥I-B-1A∥<1 für eine beliebige induzierte Matrixnorm. Die Matrix I-B-1A wird meist mit M bezeichnet und bildet die Iterationsmatrix des Verfahrens. Es gibt verschiedene Splitting-Verfahren, die ich hier kurz vorstellen werde. Sie unterscheiden sich in der Iterationsmatrix, können aber alle auf die allgemeine Form von M zurückgeführt werden.
Dazu wird häufig die Matrix A in folgender Weise zerlegt: A=D-L-R, wobei D die Diagonalmatrix, L die (negative) untere Dreiecksmatrix und R die (negative) obere Dreiecksmatrix von A darstellt. Hier die Iterationsmatrizen zu den gängigen Verfahren:

Jacobi: MJ := D-1(L+R)
Gauß-Seidel (GS): MGS := (D-L)-1R
SOR: MSOR := (D-ωL)-1((1-ω)D+ωR) wobei ω ∈ (0,2)
SSOR: MSSOR := (D-ωR)-1((1-ω)D+ωL)(D-ωL)-1((1-ω)D+ωR)

Um nicht mit der Matrix-Vektor-Notation arbeiten zu müssen, wird die Iteration meist als xk+1 = Φ(xk) aufgefasst. Für die einzelnen Verfahren sehen dann die Iterationsvorschriften wie folgt aus (komponentenweise für i=1,...,n):

Jacobi:

xik+1 =

1 aii

(

bi -

n ∑ j=1,j≠i

aijxjk

)

GS:

xik+1 =

1 aii

(

bi -

i-1 ∑ j=1

aijxjk+1 -

n ∑ j=i+1

aijxjk

)

SOR:

xik+1 = (1-ω)xik +

ω aii

(

bi -

i-1 ∑ j=1

aijxjk+1 -

n ∑ j=i+1

aijxjk

)

Die Iterationsvorschrift für das SSOR-Verfahrens (symmetrisches SOR-Verfahren) ist noch komplexer als die des SOR-Verfahrens (successive-overreaxation). Deshalb, und weil in der Anwendung, wenn überhaupt, nur die Iterationsmatrix eine Rolle spielt, wurde die Vorschrift hier weggelassen. Interessant ist, dass im GS- und SOR-Verfahren auch auf Komonenten des (k+1)ten Schritts zurückgegriffen wird und zwar genau auf die, die schon berechnet sind. Dadurch werden bessere Nährungswerte zur Berechnung der nächsten Iterierten herangezogen, was die Konvergenz verbessert.
Die Bestimmung eines optimalen ω im (S)SOR-Verfahren ist auch recht kompliziert. Es sollte, wie schon geschrieben, aus dem Intervall (0,2) kommen. Das optimale ω für SOR wird wie folgt berechnet:

ω* :=

2 1+√(1-ρJ2)

∈ (1,2)

Um es auszurechnen benötigt man den Spektralradius der Iterierten des Jacobi-Verfahrens (ρJ). Das bedeutet man muss den betragsmäßig größten Eigenwerte dieser Matrix bestimmen. Das ist oftmals recht aufwendig.

Hier nun die implementierten Verfahren:

• function [x] = jacobi(A,b,x,maxfehler)
• % Jacobi Verfahren zur iterativen Lösung eines linearen GLS Ax=b
• %
• % Als Eingabe wird die Matrix und der Vektor b erwartet
• % Als 3. Eingabewert kann man einen Startvektor x0 angeben,
• % wird dieser nicht angegeben wird der Nullvektor (0,0,...,0) angenommen
• % Als 4. Parameter kann eine Fehlerschranke angegeben werden,
• % wird diese nicht angegeben wird der Wert 10^(-10) angenommen
• % Die Funktion liefert eine Näherung für den Vektor x
•  
• n=lx
•  
• n=length(b);
•  
• if(nargin<3)
• x=zeros(n,1);
• else
• [n1 m1] = size(x);
• if(n1==1) x=x'; end
• end
• if(nargin<4)
• maxfehler=10^(-10);
• end
•  
• [n1 m1] = size(b);
• if(n1==1) b=b'; end
•  
• xneu = zeros(n,1);
• m=0;
• while true
• m=m+1;
• for k=1:n
• komp=0;
• for i=1:n
• if(i~=k)
• komp = komp+A(k,i)*x(i);
• end
• end
• xneu(k) = (b(k)-komp)/A(k,k);
• end
• x=xneu;
• fehler = sqrt(sum(abs(A*x-b).^2));
• if(fehler
• break;
• end
• end
• IterationsSchritte = m
• Fehler = fehler

• function [x] = gs(A,b,x,maxfehler)
• % Gauß-Seidel Verfahren zur iterativen Lösung eines linearen GLS Ax=b
• % A
• % Als Eingabe wird die Matrix und der Vektor b erwartet
• % Als 3. Eingabewert kann man einen Startvektor x0 angeben,
• % wird dieser nicht angegeben wird der Nullvektor (0,0,...,0) angenommen
• % Als 4. Parameter kann eine Fehlerschranke angegeben werden,
• % wird diese nicht angegeben wird der Wert 10^(-10) angenommen
• % Die Funktion liefert eine Näherung für den Vektor x
•  
• n=lx
•  
• n=length(b);
•  
• if(nargin<3)
• x=zeros(n,1);
• else
• [n1 m1] = size(x);
• if(n1==1); x=x'; end
• end
• if(nargin<4)
• maxfehler=10^(-10);
• end
•  
• [n1 m1] = size(b);
• if(n1==1); b=b'; end
•  
• m=0;
• while true
• m=m+1;
• for k=1:n
• komp=0;
• for i=1:n
• if(i~=k)
• komp = komp+A(k,i)*x(i);
• end
• end
• x(k) = (b(k)-komp)/A(k,k);
• end
• fehler = sqrt(sum(abs(A*x-b).^2));
• if(fehler
• break;
• end
• end
• IterationsSchritte = m
• Fehler = fehler

• function [x] = sor(A,b,w,x,maxfehler)
• % SOR-Verfahren zur iterativen Lösung eines linearen GLS Ax=b
• %
• % Als Eingabe wird die Matrix und der Vektor b erwartet
• % Zusätzlich kann ein Relaxionsparameter w angegeben werden, sonst wird w=1 angenommenn
• % Als 4. Eingabewert kann man einen Startvektor x0 angeben,
• % wird dieser nicht angegeben wird der Nullvektor (0,0,...,0) angenommen
• % Als 5. Parameter kann eine Fehlerschranke angegeben werden,
• % wird diese nicht angegeben wird der Wert 10^(-10) angenommen
• % Die Funktion liefert eine Näherung für den Vektor x
•  
• if(x
•  
• if(nargin<3)
• w=1
• end
• if(nargin<4)
• n=length(b);
• x=zeros(n,1);
• else
• [n m] = size(x);
• if(n==1)
• x=x';
• n=m;
• end
• end
• if(nargin<5)
• maxfehler=10^(-10);
• end
•  
• [n1 m1] = size(b);
• if(n1==1)
• b=b';
• end
•  
• m=0;
• while true
• m=m+1;
• for k=1:n
• x(k) = (1 - w)*x(k) + (w/A(k,k))*(b(k) - (A(k,:)*x-A(k,k)*x(k)));
• end
• fehler = sqrt(sum(abs(A*x-b).^2));
• if(fehler
• break;
• end
• end
•  
• IterationsSchritte = m
• Fehler = fehler