5.3 Iterative Lösungsverfahren für lineare Gleichungssysteme

5.3.3 Relaxationsverfahren: das SOR-Verfahren

Das vorangehende Beispiel zeigt, dass Jacobi- sowie Gauß-Seidel-Verfahren sehr langsamkonvergieren. Für die Modellmatrix A ∈ R

n×n steigt die Anzahl der notwendigen Iterati-onsschritte (zum Erreichen einer vorgegebenen Genauigkeit) quadratisch O(n2). Obwohljeder einzelne Schritt sehr einfach ist und äußerst effizient in O(n) Operationen durch-geführt werden kann, sind diese Verfahren den direkten nicht überlegen. Oft, etwa beiTridiagonalsystemen sind direkte Löser mit einem Aufwand von O(n) sogar unschlagbarschneller.

Das SOR-Verfahren ist eine Weiterentwicklung der Gauß-Seidel-Iteration durch die Ein-führung eines Relaxationsparameters ω > 0. Das i-te Element berechnet sich laut Satz 5.24als:

xk,GSi =

1

aii

bi −∑

j<i

aijxk,GSj −

∑

j>i

aijxk−1j

, i = 1, . . . , n.

Zur Bestimmung der SOR-Lösung verwenden wir nicht unmittelbar diese Approximationxk,GS

i , sondern führen einen Relaxationsparameter ω > 0 ein und definieren

xk,SORi = ωxk,GS

i + (1 − ω)xk−1i , i = 1, . . . , n

als einen gewichteten Mittelwert zwischen Gauß-Seidel Iteration und alter Approximati-on. Dieser Relaxationsparameter ω kann nun verwendet werden, um die Konvergenzeigen-schaften der Iteration wesentlich zu beeinflussen. Im Fall ω = 1 ergibt sich gerade dasGauß-Seidel-Verfahren. Im Fall ω < 1 spricht man von Unterrelaxation, im Fall ω > 1von Überrelaxation. Das SOR-Verfahren steht für Successive Over Relaxation, verwendetalso Relaxationsparameter ω > 1. Successive (also schrittweise) bedeutet, dass die Relaxa-tion für jeden einzelnen Index angewendet wird. Die Vorstellung zunächst die kompletteGauß-Seidel Approximation xk,GS zu berechnen und xk,SOR = ωxk,GS + (1 − ω)xk−1 zubestimmen ist falsch! Stattdessen definieren wir in Indexschreibweise:

xk,SORi = ω

1

aii

bi −∑

j<i

aijxk,SORj −

∑

j>i

aijxk−1j

+ (1 − ω)xk−1i , i = 1, . . . , n.

In Vektorschreibweise gilt:

xk,SOR = ωD−1(b − Lxk,SOR − Rxk−1) + (1 − ω)xk−1.

Trennen der Terme nach xk,SOR sowie xk−1 ergibt

(D + ωL)xk,SOR = ωb + [(1 − ω)D − ωR]xk−1,

also die Iteration

xk,SOR = Hωxk−1 + ω[D + ωL]−1b, Hω := [D + ωL]−1[(1 − ω)D − ωR].

Dieses Verfahren, mit der Iterationsmatrix Hω passt wieder in das Schema der allgemeinenFixpunktiterationen und gemäß Satz 5.21 hängt die Konvergenz des Verfahrens an ρω :=spr(Hω) < 1.

211

5 Numerische Iterationsverfahren

Die Schwierigkeit bei der Realisierung des SOR-Verfahrens ist die Bestimmung von gutenRelaxationsparametern, so dass die Matrix Hω einen möglichst kleinen Spektralradiusbesitzt. Es gilt die erste Abschätzung:

Satz 5.31 (Relaxationsparameter des SOR-Verfahrens). Es sei A ∈ Rn×n mit regulärem

Diagonalteil D ∈ Rn×n Dann gilt:

spr(Hω) ≥ |ω − 1|, ω ∈ R.

Für spr(Hω) < 1 muss gelten ω ∈ (0, 2).

Beweis: Wir nutzen die Matrix-Darstellung der Iteration:

Hω = [D + ωL]−1[(1 − ω)D − ωR] = (I + w D−1L︸ ︷︷ ︸

=:L̃

)−1 D−1D︸ ︷︷ ︸

=I

[(1 − ω)I − ω D−1R︸ ︷︷ ︸

=:R̃

].

Die Matrizen L̃ sowie R̃ sind echte Dreiecksmatrizen mit Nullen auf der Diagonale. D.h., esgilt det(I+ωL̃) = 1 sowie det((1−ω)I−ωR̃) = (1−ω)n, also Nun gilt für die Determinantevon Hω

det(Hω) = 1 · (1 − ω)n.

Für die Eigenwerte λi von Hω gilt folglich

n∏

i=1

λi = det(Hω) = (1 − ω)n ⇒ spr(Hω) = max1≤i≤n

|λi| ≥

(n∏

i=1

|λi|

) 1

n

= |1 − ω|.

Die letzte Abschätzung nutzt, dass das geometrische Mittel von n Zahlen kleiner ist, alsdas Maximum. �

Dieser Satz liefert eine erste Abschätzung für die Wahl des Relaxationsparameters, hilftjedoch noch nicht beim Bestimmen eines Optimums. Für die wichtige Klasse von positivdefiniten Matrizen erhalten wir ein sehr starkes Konvergenzresultat:

Satz 5.32 (SOR-Verfahren für positiv definite Matrizen). Es sei A ∈ Rn×n eine symme-

trisch positiv definite Matrix. Dann gilt:

spr(Hω) < 1 für 0 < ω < 2.

SOR-Verfahren und auch Gauß-Seidel-Verfahren sind konvergent.

Beweis: Siehe [9]. �

Für die oben angegebene Modellmatrix ist die Konvergenz von Jacobi- sowie Gauß-Seidel-Iteration auch theoretisch abgesichert. Für diese Matrizen (und allgemein für die Klasseder konsistent geordneten Matrizen, siehe [9]) kann für die Jacobi- J und Gauß-Seidel-Iteration H1 der folgende Zusammenhang gezeigt werden:

spr(J)2 = spr(H1). (5.15)

212

5.3 Iterative Lösungsverfahren für lineare Gleichungssysteme

0

0.2

0.4

0.6

0.8

1

0 0.5 1 1.5 2

sp

r

omega

Abbildung 5.1: Bestimmung des optimalen Relaxationsparameters ωopt.

Ein Schritt der Gauß-Seidel-Iteration führt zu der gleichen Fehlerreduktion wie zwei Schrit-te der Jacobi-Iteration. Dieses Resultat findet sich in Beispiel 5.30 exakt wieder. Wei-ter kann für diese Matrizen ein Zusammenhang zwischen Eigenwerten der Matrix Hω

sowie den Eigenwerten der Jacobi-Matrix J hergeleitet werden. Angenommen, es giltρJ := spr(J) < 1. Dann gilt für den Spektralradius der SOR-Matrix:

spr(Hω) =

ω − 1 ω ≤ ωopt,14(ρJω +

√

ρ2Jω2 − 4(ω − 1))2 ω ≥ ωopt

Ist der Spektralradius der Matrix J bekannt, so kann der optimale Parameter ωopt kannals Schnittpunkt dieser beiden Funktionen gefunden werden, siehe Abbildung 5.1. Es gilt:

ωopt =2(1 −

√

1 − ρ2J)

ρ2J

. (5.16)

Beispiel 5.33 (Modellmatrix mit SOR-Verfahren). Wir betrachten wieder die vereinfachteModellmatrix aus Beispiel 5.30. Für die Jacobi-Matrix J = −D−1(L + R) gilt:

J =

0 12

12

0 12

. . .. . .

. . .12

0 12

12

0

.

Zunächst bestimmen wir die Eigenwerte λi und Eigenvektoren wi für i = 1, . . . , n dieserMatrix. Hierzu machen machen wir den Ansatz:

wi = (wik)k=1,...,n, wi

k = sin(

πik

n + 1

)

.

Dann gilt mit dem Additionstheoremen sin(x ± y) = sin(x) cos(y) ± cos(x) sin(y):

(Jwi)k =1

2wi

k−1 +1

2wi

k+1 =1

2

(

sin(

πi(k − 1)

n + 1

)

+ sin(

πi(k + 1)

n + 1

))

=1

2sin(

πik

n + 1

)(

cos(

πi

n + 1

)

+ cos(

πi

n + 1

))

= wik cos

(πi

n + 1

)

.

213

5 Numerische Iterationsverfahren

Man beachte, dass diese Gleichung wegen wi0 = wi

n+1 = 0 auch für die erste und letzteZeile, d.h. für i = 1 sowie i = n gültig ist. Es gilt λi = cos(πi/(n+1)) und der betragsmäßiggrößte Eigenwert von J wird für i = 1 sowie i = n angenommen. Hier gilt mit derReihenentwicklung des Kosinus:

λmax = λ1 = cos(

π

n + 1

)

= 1 −π2

2(n + 1)2+ O

(1

(n + 1)4

)

.

Der größte Eigenwert geht mit n → ∞ quadratisch gegen 1. Hieraus bestimmen wirmit (5.16) für einige Schrittweiten aus Beispiel 5.30 die optimalen Relaxationsparame-ter:

n λmax(J) ωopt

320 0.9999521084 1.980616162640 0.9999879897 1.9902456641280 0.9999969927 1.995107064

Schließlich führen wir für diese Parameter das SOR-Verfahren mit optimalem Relaxati-onsparameter durch und fassen die Ergebnisse in folgender Tabelle zusammen:

Matrixgröße Jacobi Gauß-Seidel SORSchritte Zeit (sec) Schritte Zeit (sec) Schritte Zeit (sec)

320 147 775 0.92 73 888 0.43 486 ≪ 1640 588 794 7.35 294 398 3.55 1 034 0.02

1 280 2 149 551 58 1 074 776 29 1937 0.052 560 4 127 0.225 120 zu aufwendig 8 251 0.90

10 240 16 500 3.56

Die Anzahl der notwendigen Schritte steigt beim SOR-Verfahren nur linear in der Problem-größe. Dies ist im Gegensatz zum quadratischen Anstieg beim Jacobi- sowie beim Gauß-Seidel-Verfahren ein wesentlicher Fortschritt. Da der Aufwand eines Schrittes des SOR-Verfahrens mit dem von Jacobi- und Gauß-Seidel vergleichbar ist für das SOR-Verfahrenzu einem Gesamtaufwand von nur O(n2) Operationen. Dieses positive Resultat gilt jedochnur dann, wenn der optimale SOR-Parameter bekannt ist.

5.3.4 Praktische Aspekte

Wir fassen zunächst die bisher vorgestellten Verfahren zusammen:

Beispiel 5.34 (Einfache Iterationsverfahren). Es gilt in allgemeiner Darstellung

xk+1 = xk + C−1(b − Axk) = (I − C−1A)︸ ︷︷ ︸

=B

xk + C−1b.

214

5.3 Iterative Lösungsverfahren für lineare Gleichungssysteme

Ausgehend von der natürlichen Aufspaltung A = L + D + R sowie mit einem Relaxati-onsparameter ω sind Richardson-, Jacobi-, Gauß-Seidel- sowie SOR-Verfahren gegeben alsals:

• Gedämpftes Richardson Verfahren:

C−1 = ωI, B = I − ωA,

mit der Index-Schreibweise:

xki = ωbi + xk−1

i − ωn∑

j=1

aijxk−1j , i = 1, . . . , n.

• Jacobi-Verfahren:C−1 = D−1, B = −D−1(L + R),

mit der Index-Schreibweise:

xki =

1

aii

bi −n∑

j=1,j 6=i

aijxk−1j

, i = 1, . . . , n.

• Gauß-Seidel-Verfahren

C−1 = [D + L]−1, B = −(D + L)−1R

mit der Index-Schreibweise:

xki =

1

aii

bi −∑

j<i

aijxkj −

∑

j>i

aijxk−1j

, i = 1, . . . , n.

• SOR-Verfahren (englisch. Successive Over-Relaxation):

C = [D + ωL]−1, B = [D + ωL]−1[(1 − ω)D − ωR], ω = ωopt ∈ (0, 2),

mit der Index-Schreibweise:

xki = ω

1

aii

bi −∑

j<i

aijxkj −

∑

j>i

aijxk−1j

+ (1 − ω)xk−1i , i = 1, . . . , n.

Zur einfachen Durchführung der Verfahren eignet sich stets die Index-Schreibweise. DieMatrix-Form dient insbesondere der einfachen Charakterisierung sowie zum Herleiten vonKonvergenzaussagen.

Als iterative Verfahren werden die Gleichungssysteme nur im (praktisch irrelevanten) Falln → ∞) wirklich gelöst. Üblicherweise muss die Iteration nach einer bestimmten Anzahlvon Schritten abgebrochen werden. Als Kriterium für ein Abbrechnen kann zunächst die

215

5 Numerische Iterationsverfahren

asymptotische Konvergenzaussage aus Satz 5.21 herangezogen werden. Mit ρ := spr(B)gilt im Grenzfall:

‖xk − x‖ ≈ ρk‖x0 − x‖.

Und bei vorgegebener Toleranz TOL kann die notwendige Zahl an Iterationsschrittenabgeschätzt werden:

‖xk − x‖ < TOL ⇒ k =log

(T OL

‖x0−x‖

)

log(ρ).

Die Toleranz TOL gibt hier an, um welchen Faktor der Anfängliche Fehler ‖x0 − x‖reduziert wird. Dieses Vorgehen ist in der praktischen Anwendung wenig hilfreich, da derSpektralradius ρ der Iterationsmatrix B im Allgemeinen nicht bekannt ist.

Ein alternatives allgemeines Kriterium liefert die Abschätzung aus Satz 4.44 für den Defektdk := b − Axk:

‖xk − x‖

‖x‖≤ cond(A)

‖b − Axk‖

‖b‖.

Hier entsteht jedoch ein ähnliches Problem: die Konditionszahl der Matrix A ist im All-gemeinen nicht bekannt, so kann auch keine quantitativ korrekte Abschätzung hergeleitetwerden. Dieser einfache Zusammenhang zwischen Defekt und Fehler kann jedoch genutztwerden um eine relative Toleranz zu erreichen:

Bemerkung 5.35 (Relative Toleranz). Bei der Durchführung von iterativen Lösungsver-fahren werden als Abbruchskriterium oft relative Toleranzen eingesetzt. Die Iteration wirdgestoppt, falls gilt:

‖xk − x‖ ≤ TOL ‖x0 − x‖.

Als praktisch durchführbares Kriterium werden die unbekannten Fehler durch die Defekteersetzt:

‖b − Axk‖ ≤ TOL ‖b − Ax0‖.

216