Sie sind hier: Numerik > Lösung linearer Gleichungssysteme > Splitting-Verfahren zur Lösung ...

Splitting-Verfahren zur Lösung von LGS

Mittwoch 27. Juni 2007 von
Simon Praetorius
Splitting-Verfahren zu Lösung linearer Gleichungssysteme sind eine Klasse iterativer Verfahren, die in den letzten Jahren stark an Bedeutung verloren haben, da z.B. die Konvergenzeigenschaften mit anderen Verfahren viel besser sind, trotzdem dienen sie noch als Vorkonditionierer oder z.B. in einem modifizierten Newtonverfahren zur Bestimmung der Lösung einer nichtlinearen Gleichung.
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:
matlab Code
  • function [x] = jacobi(A,b,x,maxfehler)
  • % Jacobi Verfahren zur iterativen Lsung 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 fr 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
matlab Code
  • function [x] = gs(A,b,x,maxfehler)
  • % Gau-Seidel Verfahren zur iterativen Lsung 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 fr 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
matlab Code
  • function [x] = sor(A,b,w,x,maxfehler)
  • % SOR-Verfahren zur iterativen Lsung 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 fr 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
Besucher: 14772 | Permalink | Kategorie: Numerik
Tags: , , ,

Kommentar hinzufügen

Dieses Feld bitten nicht ausfüllen: