Kleiner Einschub: Etwas Heuristik zu dieser Aussage Sei
(L x)(t) =f(t), (
(L x)(t) =x(m)(t) +∑m−1
k=0 ak(t)x(k)(t))
Rx= 0 (264)
eine Randwertaufgabe m-ter Ordnung und
h−mLhxh =fh (265)
eine Diskretisierung davon zur Schrittweiteh= 1/(n+ 1). Um die Argumenta-tion zu erleichtern, nehmen wir einmal für die MatrixLh ∈Rn×n an, dass ihre Diagonale einen konstanten Wert enthält, der auch noch h-unabhängig ist:
Diagonale(Lh) = αI.
Dann ist die Iterations-Matrix des Jacobi-Verfahrens zu (265) gegeben durch GJ =α−1(αI−Lh) =I−α−1hm(h−mLh).
Mit der Diskretisierung möchte man ja die Randwertaufgabe (264) approxi-mieren. Wir ziehen daraus den (zugegeben flinken) Schluß, dass die Eigenwerte λh,1, . . . , λh,nder Matrixh−mLh die erstenn Eigenwerteλ1, . . . , λnder zu (264) gehörigen Eigenwertaufgabe
(L x)(t) =λx(t),
Rx = 0
approximieren. Für die Eigenwerte solcher Eigenwertaufgaben der Ordnungm gilt normalerweise
λj =O(jm).
Damit schließen wir, dass tendenziell gelten sollte spec(GJ)≈{
1−α−1hm,1−α−1hm 2m, . . . ,1−α−1hm nm ≈1−α−1} . Dies ergibt (auch wenn wir noch versuchen α geeignet zu variieren (gedämpfte Jacobi-Iteration)) auf jeden Fall
ρ(GJ)−→h→0 1 wie hm −→0.
Weil sich die Eigenwerte und Eigenvektoren der Matrix (262) explizit angeben lassen, ist dies Problem ideal geeignet, um erste Ideen über das Verhalten von Iterationsverfahren zu erhalten. Man kann nämlich - wie oben gesehen - über die Spektralzerlegung der System-matrix selbst oft auch Spektralanalysen des Iterationsprozesses185 vornehmen.
Wir haben das obige System (262) aber aus zwei weiteren Gründen als Beispiel ausge-wählt. Denn es gehört gleichzeitig zu zwei wichtigen Klassen von Matrizen, für die Konver-genzaussagen zu den Standard-Iterationsmethoden bekannt sind. Einerseits ist die Matrix
„SPD“186 und andererseits ist sie „schwach diagonaldominant“ und „nicht zerfal-lend“ (bzw.anderslautend aber dasselbe bezeichnend: „irreduzibel“ ), vgl. hierzu Definitio-nen 5.5 und 5.7.
Die Konsequenzen dieser letzten Eigenschaften behandeln wir im nächsten Unterabschnitt.
Hier wird für die Systemmatrix weder Symmetrie noch positive Definitheit vorausgesetzt.
In den zuhörigen Konvergenzbeweisen für die beiden Splitting-Verfahren nach Jacobi und Gauss-Seidel, muss der Rechenfluss en detail analysiert werden.
Unter der Voraussetzung der positiven Definitheit vonAlassen sich für diese Verfahren als auch gleich für overrelaxierte oder gedämpfte Versionen elegantere Beweise führen. Dies geschieht im übernächsten Unterabschnitt.
5.1.1 Diagonaldominanz und schwache Diagonaldominanz
Wir wollen in diesem Abschnitt über Kriterien berichten, welche an der Systemmatrix erkennen lassen, ob die Standard-Verfahren für Gleichungssysteme wie (262) konvergieren.
Einerseits hatten wir als leicht kontrollierbare und für die Konvergenz der Jacobi-Iteration ausreichende Eigenschaft schon auf der Seite 52 die folgende Eigenschaft der starken Zeilen-Diagonaldominanz erkannt:
Definition 5.3 (Starke Zeilen Diagonaldominanz)
A= (aij)∈R(n,n) heißt „stark zeilen-diagonal-dominant“, wenn
|aii| >∑
k̸=i
|aik| ∀i= 1, . . . , n (266)
Auf Seite 52 hatten wir auch schon festgestellt, dass die starke Spaltendiagonaldominanz hinreichend für die Konvergenz der Jacobi-Iteration ist.
Definition 5.4 (Starke Spalten Diagonaldominanz)
A= (aij)∈R(n,n) heißt „stark spalten-diagonal-dominant“, wenn
|ajj|>
∑n i=1,i̸=j
|ai,j|, j = 1, . . . , n. (267)
Leider sind beide Eigenschaften bei der Matrix (262) nicht erfüllt. Hier gelten die Be-dingungen (266) und auch (267) nur noch in dem Sinne „schwach“, dass dort anstelle einen <-Zeichens fast überall das Gleichheitszeichen gilt. Nur in zwei Zeilen (bzw. Spal-ten) bleibt das<-Zeichen erhalten. Die einfach Bildung der ∞-Norm der Iterationsmatrix bringt uns hier nicht die beobachte Konvergenz sowohl vom Jacobi- als auch vom Gauss-Seidel-Verfahren.
Für Abhilfe sorgt hier das Zusammenwirken dieser „schwache Diagonaldominanz“ mit der sogenannten „Irreduzibilität“ der Systemmatrix:
Definition 5.5 (Schwache Zeilen-Diagonaldominanz)
A= (aij)∈R(n,n) heißt „schwach Zeilen-diagonal-dominant“, wenn
|aii| ≥∑
k̸=i
|aik| ∀i= 1, . . . , n (268)
aber auch
|app| >∑
k̸=p
|apk| für mindestens ein p∈ {1, . . . , n}, (269)
Anmerkung:Wenn die beiden Bedingungen (268, 269) gelten, so sagt man auch, dass die Matrix das „schwache Zeilensummenkriterium“ erfüllt.
Definition 5.6 (Schwache Spalten-Diagonaldominanz)
A = (aij) ∈ R(n,n) heißt „schwach Spalten-diagonal-dominant“, wenn AT schwach Zeilen-diagonal-dominant ist.
Definition 5.7 (Reduzibilität und Irreduzibilität)
Eine Matrix A ∈ R(n,n) heißt „reduzibel“ , wenn es eine Permutationsmatrix P gibt, mit der
PTAP =
(A˜11 A˜12 0 A˜22
)
wird, wobei A˜ii∈Rpi×pi, i= 1,2 sind mit p1 >0, p2 >0und p1+p2 =n.
Andernfalls ist A irreduzibel187. Anmerkungen 5.8
1. Anstelle der Bezeichnung „reduzibel“ verwendet man auch „zerfallend“ oder „zerleg-bar“.
2. Eine alternative äquivalente Definition der Reduzibilität ist die Folgende:
Eine Matrix A ∈ R(n,n) heißt reduzibel, wenn es nichtleere Teilmengen N1 und N2 der Indexmenge N :={1,2,3, . . . , n} gibt mit den Eigenschaften
– (i)N1∩N2 = ∅; – (ii) N1∪N2 = N;
– (iii) aij = 0 für alle i∈N1 und j ∈N2
3. Wir werden zeigen, dass schwache Diagonaldominanz zusammen mit der Irreduzibi-lität hinreichend für die Konvergenz sowohl der Jacobi- als auch der Gauss-Seidel-Iterationen sind.
Überzeugen Sie sich durch Konstruktuion von Beispielen, dass man auf keine der drei Eigenschaften (268), (269) und Irreduzibilität verzichten kann, ohne den Verlust der Konvergenz des Jacobi-Verfahrens zu verlieren. (Dies gilt auch für das Gauss-Seidel-Verfahren, ist aber schwieriger zu analysieren.)
4. Die Reduzibilität einer Matrix bedeutet, dass es eine Gruppe von Komponenten des Lösungsvektors gibt, die im System nicht von den komplementären Komponenten abhängen. (Nach Anwendung der obigen Permutation sind dies gerade die letzten p2 Komponenten.) Der folgende Satz sichert die Konvergenz der Jacobi-Iteration für schwach diagonal-dominante Matrizen, die gleichzeitig irreduzibel sind. Diese letzte Forderung ist für die Anwendungen nicht einschränkend. Ist nämlich eine Matrix reduzibel, so kann ein Teilsystem des Problems unabhängig vom Rest bearbeitet werden. Das ist sicher kein Nachteil. Nach Elimination der entsprechenden Variablen ist das Restsystem entweder irreduzibel, oder es kann weiter reduziert werden.
5. Es ist von Vorteil, die Irreduzibilität einer Matrix über einen zugehörigen gerichte-ten Adjazenz-Graphen zu interpretieren. Wie beim ungerichtegerichte-ten Adjazenzgraphen symmetrischer Matrizen188 wird der i–ten Zeile der Matrix A∈ Rn×n ein Knoten vi zugeordnet, i = 1, . . . , n. Im gerichteten Adjazenzgraphen GG(A) von A wird nun
der Knoten vi mit dem Knoten vj durch die gerichtete Kante (vi, vj) verbunden189, wennaij ̸= 0 ist190.
Eine Matrix ist nun genau dann irreduzibel, wenn es in ihrem gerichteten Adjazenz-Graphen zu je zwei Knoten vi und vj stets einen gerichteten verbindenden Weg191 gibt192.
6. Für die spätere Verwendung schließen wir, dass es bei irreduziben MatrizenM zu zwei nicht leeren und disjunkten KnotenmengeJ, K ⊂ {1, . . . , n} mit J∪K ={1, . . . , n} stets zwei Indizes j ∈ J und k ∈ K gibt mit mjk ̸= 0. Wählt man nämlich zwei beliebigen Indizesˆj ∈J und kˆ∈K, so gibt es einen gerichteten Weg inGG(M)von ˆj nach K. Dieser Weg muss irgendwannˆ J verlassen und in K eintreten. Hier findet man die Indizesj und k.
Beispiel 5.9 (Reduzibilität und Irreduzibilität)
∗ bedeutet im weiteren stets ein Nicht-Null-Element. Unter den Matrizen
A1 :=
∗ 0 ∗ 0 0 0 ∗ 0 ∗ 0
∗ 0 ∗ 0 ∗ 0 ∗ 0 ∗ 0 0 0 ∗ 0 ∗
, A2 :=
∗ 0 0 0 ∗
∗ ∗ 0 0 0 0 ∗ ∗ 0 0 0 0 ∗ ∗ 0 0 0 0 ∗ ∗
, A3 :=
∗ ∗ 0 0 0
∗ ∗ ∗ 0 0 0 ∗ ∗ ∗ 0 0 0 ∗ ∗ ∗ 0 0 0 ∗ ∗
istA1reduzibel, die anderen beiden sind irreduzibel. Ihre gerichteten GraphenGG(Ak), k = 1, ..,3 sind in der folgenden Abbildung dargestellt.
0 1 2 3 4 5 6
−1
−0.8
−0.6
−0.4
−0.2 0 0.2 0.4 0.6 0.8 1
1 3 5 2 4
GG(A1)
−1 −0.8 −0.6−0.4−0.2 0 0.2 0.4 0.6 0.8 1
−1
−0.8
−0.6
−0.4
−0.2 0 0.2 0.4 0.6 0.8 1
4 3
2
1
GG(A2) 5
0 1 2 3 4 5 6
−1
−0.8
−0.6
−0.4
−0.2 0 0.2 0.4 0.6 0.8 1
1 2 3 4 5
GG(A3)
Abbildung 96: Gerichtete Graphen
Wie man übrigens aus dem Graphen sofort abliest, zerfälltA1 in zwei unabhängige jeweils irreduzible Untersysteme.
Satz 5.10 (Eine zweite Konvergenz-Aussage für die Jacobi-Iteration )
Ist die System-Matrix A ∈ Cn×n schwach diagonal-dominant und irreduzibel so ist das Jacobi-Verfahren durchführbar und konvergent.
189Wichtig: Anders als beim ungerichteten Graphen kommt es hier bei der Kante(vi, vj)auf die Reihen-folge der Knoten an.
190Im unsymmetrischen Fall muß mit aij ̸= 0 nicht notwendig auch aji ̸= 0 sein. Dieser strukturellen Unsymmetrie wird im Graphen mit der Einführung einer Richtung der Kanten Rechnung getragen.
191Dabei erklären wir einen solchen gerichteten Weg nicht genauer, sondern gehen davon aus, dass diese Notation anschaulich klar ist.
192Die Forderung schließt ein, dass es sowohl einen Hinweg als auch einen Rückweg gibt; denn mitvi und vj sind auchvj undvi „ je zwei Knoten“.
Bevor wir auf die Konvergenz eingehen, wollen wir erst einmal zeigen, dass die Vorausset-zungen der letzen Konvergenzaussage implizieren, dass das Verfahren überhaupt ausführ-bar ist (dass also die Diagonalelemente ungleich Null sind). Dabei lernen wir zugleich auch noch, dass dann die Systemnmatrix zwingend regulär ist.
Satz 5.11 (Satz über schwach diagonaldominante, irreduzible Matrizen )
Eine quadratische schwach diagonaldominante und irreduzible Matrix M ist regulär und hat lauter nichtverschwindende Diagonalelemente.
Beweis: Wenn M nicht regulär ist, gibt es einen Vektor x ̸= 0 mit M x = 0. Das heißt, dass
miixi =−
∑n l=1,l̸=i
milxl, i= 1, . . . , n, und mit der Dreiecksungleichung folgt hieraus
|mii| |xi| ≤
∑n l=1,l̸=i
|mil| |xl|, i= 1, . . . , n. (270)
Wir definieren nun Indexmengen
J :={j | |xj|=∥x∥∞} und K :={k| |xk|<∥x∥∞}. Es ist trivialerweise J ̸=∅, und es ist ebenfalls K ̸=∅, weil sonst
|xi|=∥x∥∞ für alle i= 1, . . . , n wäre, was nach (270) bedeutete, dass
|mii| ≤
∑n l=1,l̸=i
|mil|, i= 1, . . . , n,
Dies widerspräche aber (269). Weil J und K beide nicht leer sind und M irreduzibel ist, gibt es j ∈J und k ∈K mit mjk ̸= 0. Damit ergibt sich
|mjj| ≤ ∑n
l=1,l̸=j|mjl|||xxjl|| =∑n
l=1,l̸=j|mjl|∥|xx∥l∞|
= ∑n
l∈J,l̸=j|mjl|∥∥xx∥∥∞∞ +∑n
l∈K,l̸=j,k|mjl|∥|xx∥l∞| +|mjk|∥|xx∥k∞|
< ∑n
l=1,l̸=j|mjl|, da ∥|xx∥k|
∞ <1 und mlk ̸= 0. Die gewonnene Ungleichung
|mjj|<
∑n l=1,l̸=j
|mjl|
widerspricht nun aber (268), weshalb M nicht singulär sein kann, was der erste Teil der Behauptung ist.
Dass die Diagonalelemente alle von Null verschieden sind, schließt man so. Wäre das
Dia-wie gerade gezeigt - nicht der Fall ist.
2 Einen eleganten Beweis für die Konvergenzaussagen findet man etwas weiter unten. Die-ser ist aber ein wenig indirekt und macht Gebrauch von der letzten Implikation „schwach diagonal-dominant und irreduzibel ⇒ regulär“. Es gibt andere, rechentechnisch aufwendi-gere Beweise, bei denen man aber die Wirkungsweise der Voraussetzungen besser versteht.
Wegen ihrer Länge werden wir hier die Grundidee dieser Beweise nur veranschaulichen.
Beweisidee zum letzten Konvergenzsatz:
Nach (96) gilt für die Fehlervektoren193 e[k] :=xk−x∗ die Rekursion e[k+1]=GJe[k],
wobei GJ := I −D−1A = (gij)ni,j=1 die Iterationsmatrix des Jacobi- oder Gesamtschritt-verfahrens bezeichne. Da die Betragssummen der Zeilen von GJ nach Voraussetzung der schwachen Diagonaldominanz alle kleiner oder gleich 1 sind, gilt für die i–te Komponente von e[k+1] sicherlich
|e[k+1]i |=
∑n j=1
gije[k]j ≤
∑n j=1
|gij| · ∥e[k]∥∞ ≤ ∥e[k]∥∞. Der Fehlervektor wird also in der Maximum-Norm sicher nicht größer.
Nach Voraussetzung gibt es mindestens eine Zeile vonGJ, deren Betragssumme echt kleiner als 1 ist. Sei dies die m–te Zeile:
∑n j=1
|gmj|=:W <1.
Dann wird die m–te Komponente des Fehlervektors e[k+1] sogar mindestens um den Wert (1−W)∥e[k]∥∞ kleiner als ∥e[k]∥∞. Diese Verkleinerung wirkt sich beim nächsten Iterati-onsschritt in allen Gleichungen aus, welche auf diem–te Komponente des Iterationsvektors zugreifen, deren zugeordnete Knoten im Adjazenzgraphen vonAalso direkt mitvm verbun-den sind. Hier bewirken sie garantierte Reduktion der entsprechenverbun-den Fehlerkomponenten unter∥e[k]∥∞. In den nächsten Iterationen breitet sich diese Reduktion nun entlang der We-ge im Adjazenzgraphen aus. Ist L die maximale Länge eines einfachen Verbindungsweges zwischen zwei Knoten in GG(A), so hat diese Reduktion nach spätestens LSchritten alle Komponenten erfaßt, und es tritt eine echte Fehlerverkleinerung∥e[k+L]∥∞<∥e[k]∥∞ein194. Wollen wir die Idee einprägsam formulieren, so können wir etwa festhalten:
In den „stark-diagonaldominanten Zeilen wird Fehler aus dem System direkt abgeführt“.
Die Fehler der anderen Zeilen „fließen über den zusammenhängenden Adjazenzgraphen in diese Fehlersenken ab.“
Es folgt nun der angekündigte indirekte Beweis:
193Die im Satz mitbewiesene Regularität von A (und damit die eindeutige Existenz der Lösung x∗) nehmen wir hier der Einfachheit halber einmal an.
194Hier muß man natürlich noch insofern Vorsicht walten lassen, als dass ja bekanntlich nicht jede monoton fallende nichtnegative Folge auch schon Nullfolge ist. Der Grundgedanke sollte aber so doch klar geworden sein.
Beweis von Satz 5.10:
Wir zielen darauf, für die Iterationsmatrix GJ =D−1(L+R) =D−1(D−A) der Jacobi-Iteration direkt ρ(GJ)<1nachzuweisen. Dazu werden wir zeigen, dass für alle λ∈C mit
|λ| ≥1 die Matrix
M(λ) :=GJ −λI
regulär ist. Damit wäre dann gezeigt, dass alle Eigenwerte von GJ im Inneren des komple-xen Einheitskreises liegen müssen.
Für den Nachweis der Regularität benutzen wir Satz 5.11, nach dem wir nur zu zeigen haben, dass M(λ) schwach diagonal-dominant und irreduzibel ist. Letzteres folg aber so-fort daraus, dass Irreduzibilität und Reduzibilität Eigenschaften sind, die nur von den Nichtdiagonalelementen einer Matrix bestimmt sind. Da die Nichtdiagonalelemente von D−1(D−A)für dieselben Indexpaare von Null verschieden sind wie die vonA, ist GJ mit A irreduzibel. Wenn wir die Nichtdiagonalelemente
mjk =−ajk ajj
über die Zeile aufsummieren, ergeben sich wegen der schwachen Diagonaldominanz vonA die Ungleichungen
∑n k=1,k̸=j
ajk ajj
≤1 für alle j = 1, . . . , n, mit dem <-Zeichen für mindestens einen Index aus {1, . . . , n}.
Da |λ| ≥ 1 ist und außerdem −λ überall auf der Diagonale von M(λ) steht, können wir die letzten Ungleichung wie folgt fortsetzen
∑n k=1,k̸=j
ajk ajj
≤1≤ |λ|=|mjj| für alle j = 1, . . . , n
und ebenfalls dem <-Zeichen für mindestens einen Index aus {1, . . . , n}. D.H. aber das M(λ) tatsächlich schwach diagonal-dominant ist, und wir sind fertig.
2 Korollar 5.12
Die Aussage des vorigen Konvergenzsatzes gilt natürlich auch wieder mit dem analogen schwachen Spaltensummenkriterium. Zum Nachweis geht man wieder genau so vor wie in Bemerkung 4 auf der Seite 52.
Die starke und die schwache Zeilen-Diagonal-Dominanz (mit Irreduzibilität) sind auch hin-reichend für die Konvergenz des Gauss-Seidel-Verfahrens. Wir zeigen dies in den nächsten beiden Sätzen.
Satz 5.13 (Konvergenz-Satz für das Gauss-Seidel-Verfahren) Erfüllt A das starke Zeilensummenkriterium, so gilt
∥GGS∥∞≤ ∥GJ∥∞<1,
wobei GJ und GGS die Iterationsmatrizen von Jacobi- bzw. Gauss-Seidel-Verfahren
be-Beweis:
Die Durchführbarkeit des Verfahrens haben wir hier gar nicht erst erwähnt, da die starke Diagonaldominanz nach Satz 5.11 sofort impliziert, dass alle Diagonalelemente vonA von Null verschieden sind.
Zum Beginn erinnern wir daran, dass bei starker Diagonaldominanz
∥GJ∥∞= max
1≤i≤n
∑n j=1,j̸=i
|aij| aii| <1 gilt.
Seiy∈Rn beliebig undz :=GGSy. Durch vollständige Induktion werden wir gleich zeigen, dass
|zi| ≤ ∥GJ∥∞∥y∥∞, für alle i= 1, . . . , n (271) gilt. Maximumbildung über i zeigt, dass dann
∥GGSy∥∞=∥z∥∞≤ ∥GJ∥∞∥y∥∞
also
∥GGS∥∞ ≤ ∥GJ∥∞ <1 ist.
Für i= 1 ist (271) klar, denn es ist
a11z1 =−
∑n k=2
a1kyk und mithin
|z1| ≤
∑n k=2
|a1k|
|a11||yk| ≤
∑n k=2
|a1k|
|a11|∥y∥∞ ≤ ∥GJ∥∞∥y∥∞.
Der Induktionsschluss ergibt sich, indem wir die i-te Gleichung von z =GGSy bzw. äqui-valent (D−L)z = Ry komponentenweise aufschreiben, wobei wir annehmen, dass (271) bisi−1 gezeigt ist. Aus
aiizi =−
i−1
∑
k=1
aikzk−
∑n k=i+1
aikyk ergibt sich dann
|zi| ≤ ∑i−1 k=1
|aik|
|aii| |zk|+∑n k=i+1
|aik|
|aii||yk|
≤ ∑i−1 k=1
|aik|
|aii| ∥| {z }GJ∥∞
<1
∥y∥∞+∑n k=i+1
|aik|
|aii|∥y∥∞
≤ ∑n
k=1,k̸=i |aik|
|aii|∥y∥∞≤ ∥GJ∥∞∥y∥∞ Damit ist alles gezeigt.
Anmerkungen 5.14
Man könnte geneigt sein, das Ergebnis des letzten Satzes dahingehend zu interpretie-ren, dass das Gauss-Seidel-Verfahren stets und bedingungslos schneller konvergiert als das Jacobi-Verfahren. Wir bemerken aber, dass dazu nach dem Satz 3.9 über die asymptotische Konvergenzrate von Seite 48 die Aussage ρ(GGS) ≤ ρ(GJ) <1 erforderlich wäre. Hierfür sind i.a. aber schärfere Voraussetzungen nötig. Wäre z.B. GJ elementweise größer oder
gleich Null, so erhielte man sie etwa aus dem Satz von Stein-Rosenberg195. Dieser sagt, dass unter dieser Voraussetzung stets genau eine der folgenden vier Aussagen gilt:
(1) ρ(GGS) =ρ(GJ) = 0;
(2) 0< ρ(GGS)< ρ(GJ)<1;
(3) ρ(GGS) =ρ(GJ) = 1;
(4) ρ(GGS)> ρ(GJ)>1.
Es ist naheliegend, nun auch noch nach der Konvergenz des Gauss-Seidel-Verfahrens un-ter den Voraussetzungen der schwachen Diagonaldominaz und Irreduzibilität zu fragen.
Tatsächlich läst sich aus dies mit einem Beweis ähnlich dem zum Satz 5.11 zeigen.
Satz 5.15 (GS-Konvergenz bei Diagonaldominanz und Irreduzibilität)
Ist die System-Matrix A ∈ Cn×n schwach diagonal-dominant und irreduzibel so ist das Gauss-Seidel-Verfahren durchführbar und konvergent.
Beweis:
Die Wohldefiniertheit des Verfahrens ist klar, weil die Diagonalelemente nach Satz 5.11 von Null verschieden sind. Wie beim Beweis von Satz 5.10 betrachten wir die Matrix
M(λ) = (D−L)−1R−λI für |λ| ≥1 und wollen wieder zeigen, dassM(λ) regulär ist.
Nun ist M(λ) = (D −L)−1R − λI genau dann regulär wenn M˜ := (D −L)M(λ) = R−λD+λL regulär ist. Offenbar istM˜ =R−λD+λL mitA=−R+D−Lirreduzibel.
Da nach
∑n k=1,k̸=j
|m˜jk|=|λ|
j−1
∑
k=1
|ajk|+
∑n k=j+1
|ajk| ≤λ|
∑n k=1,k̸=j
|ajk| ≤ |λ||ajj|=|m˜jj|für j = 1, . . . , n
mit strikter Ungleichung für mindestens ein j gilt, ist M˜ auch schwach diagonaldominent,
nach Satz 5.11 also regulär. 2
5.1.2 SPD-Systeme
Bei der Wahl der Splitting Matrix eines Splitting-Verfahrens besteht die Kunst darin, unter den leicht numerisch behandelbaren Matrizen B eine solche auszusuchen, die alle Spektralwerte von A über G = I −B−1A in das Innere des Einheitskreises treibt und möglichst dicht bei Null versammelt.
Dies wird umso einfacher machbar sein, je netter das Spektrum vonA selbst aussieht. Das Spektrum symmetrischer Matrizen z.B. liegt inR, die mögliche Ausdehnung ist damit schon mal „um eine Dimension“ verkleiner. Für positiv definite Matrizen steht sogar nur noch eine Hälfte der reellen Achse zu Verfügung, weshalb mit ihnen besonders leicht umzugehen ist.
Betrachten wir z.B. das Richardson–Verfahren B =τ I, xn+1 = (I− 1
τA)xn+ 1
τb, (272)
so sieht man, dass das Spektrum von Gτ =I− 1τA gerade durch spec(Gτ) =
{ 1−λ
τ :λ ∈spec(A) }
gegeben ist. Die optimale Wahl erhält man, indem manτ so wählt, dass die Funktion gτ(λ) = 1− λ
τ
in der Mitte des Intervalles [λmin, λmin] ihre Nullstelle hat, wobeiλmin und λmax den betragskleinsten bzw. -größten Eigenwert von A bezeichnen. Das ist offenbar für
τopt = λmax+λmin 2 der Fall, und dort wird
ρ(G) = λmax−λmin
λmax+λmin = K −1 K + 1 mit der spektralen Konditionszahl
K = λmax λmin
Wie die Zeichnung 97 zeigt, nimmt die mit ρ(K)verbundenen Konvergenzrate mit zuneh-mender Kondition ab.
−5 0 5 10 15 20 25 30
−0.2 0 0.2 0.4 0.6 0.8 1 1.2
Kondition der Matrix
ρ(G)
Abbildung 97: Konvergenzrate und Kondition
Wie wir unten sehen werden, ist dies Verhalten typisch für alle Iterationsverfahren.
Zunächst wollen wir aber einige Konvergenzresultate für positiv definite Systemmatrix A festhalten. Die Beweise aller Resultate lassen sich eleganter führen, wenn man es gewohnt ist, mit durch SPD-Matrizen definierten inneren Produkten zu rechnen196.
Wir schreiben die Beweise hier so auf, wie sie mit der Linearen Algebra der ersten zwei Semester der TUHH verstehbar sind.
Relativ häufig anwendbar ist der folgende Satz.
196Vgl. dazu etwas [DW] S. 156-162
Satz 5.16 (Konvergenz bei hinreichend großer Splitting-Matrix)
Sind sowohl die SystemmatrixA ∈R(n,n) als auch die SplittingmatrixB ∈R(n,n) SPD und ist die Splitting-Matrix in folgendem Sinne wenigstens „halb so groß“ wie A:
vTAv <2vTBv für alle v ∈Rn (273) so ist ρ(I−B−1A)<1.
Beweis: Wir wollen das Spektrum von G = I −B−1A untersuchen. G ist keine symme-trische Matrix, was die Untersuchungen schwieriger zu machen droht. G ist aber wegen der SPD-Eigenschaft von B eine sogenannte „symmetrisierbare Matrix“ (vgl. hierzu gege-benenfalls [HY]).
Durch Ähnlichkeitstransformation wird ausG eine symmetrische Matrix G−→G˜ :=B1/2GB−1/2 =I −B−1/2AB−1/2, die dasselbe Spektrum hat wie G.
Sei nun (λ, x) ein Eigenpaar vonG˜ mit xTx= 1. Dann ist
λ=λxTx=xTGx˜ =xTx−xTB−1/2AB−1/2x=: 1−α. (274) Setzen wirv :=B−1/2x, so ergibt sich wegen der positiven Definitheit vonAundB erstens
α=xTB−1/2AB−1/2x >0, und zweitens wegen (273) die obere Abschätzung
α=xTB−1/2AB−1/2x <2xTB−1/2BB−1/2x= 2xTx= 2.
Mit (274) schließen wir aus α∈(0,2), dass
λ∈(−1,1).
2 Die erste Anwendung formulieren wir als
Aufgabe 5.17 Zeigen Sie:
Sei A ∈ R(n,n) SPD und sei λmax der größte Eigenwert von A. Dann konvergiert das Richardsonverfahren (272), wenn
τ > λmax/2 ist.
Auch die Konvergenz des Jacobi-Verfahrens bei starker Diagonaldominanz ist im SPD-Fall schnell mit Satz 5.16 zeigbar:
Satz 5.18 (Jacobi konvergiert für stark-diagonaldominante SPD-Matrix)
Sei A∈R(n,n) SPD und stark diagonaldominant197. Dann konvergiert die Jacobi-Iteration zur Lösung von Ax=b.
Beweis: Es ist B = D mit der Diagonalmatrix D von A, die selbst auch positiv definit ist. Dann schätzt man ab
xTAx=xTD1/2D−1/2AD−1/2D1/2x≤ρ(D−1/2AD−1/2)·xTD1/2D1/2x=ρ(D−1A)xTDx.
Da
ρ(D−1A)≤ ∥D−1A∥∞<2, haben wir
xTAx <2xTBx
gezeigt. 2
Satz 5.19 (GS konvergiert für alle SPD-Systeme)
Das Gauss-Seidel-Verfahren konvergiert für jedes Gleichungssystem mit SPD Systemmatrix A.
Beweis: Die Iterationsmatrix zum Gauss-Seidelverfahren war G = (D−L)−1LT bei der Zerlegung A=D−L−LT. Sei (x, λ) ein Eigenpaar198 von G. Dann ist
[Gx=λx]⇐⇒[LTx=λ(D−L)x].
Durch Multiplikation mit x¯T und Auflösen nach λ erhält man199 λ= x¯TLTx
¯
xTDx−x¯TLx
Seix¯TLTx=a+ibmita, b∈Rundi2 =−1. Dann istx¯TLx=a−ib. Weiter seid= ¯xTDx.
Da D SPD ist, ist d >0. mit den eingeführten Größen ist λ = a+ib
d−a+ib und daher
|λ|2 = a2+b2 (d−a)2 +b2
Der Beweis ist erbracht, wenn wir (d−a)2 > a2 zeigen können.
Aus der positiven Definitheit vonA selbst erhalten wir
0<x¯TAx= ¯xTDx−x¯TLTx−x¯TLx=d−2a.
Wenn a negativ ist, so ist sicher(d−a)> a > 0. Ista andererseits positiv, so folgt aus d >2a
dass
(d−a)>2a−a =a
ist, und der Beweis ist beendet. 2
Mit dem letzten Satz kommt man schnell zu einem konvergenten Iterationsverfahren für jedes lineare Gleichungssystem
Ax=b
198möglicherweise komplex
199Dass der Nenner im nächsten Ausdruck immer von Null verschieden ist, ergibt sich gleich aus der weiteren Rechnung.
dessen Systemmatrix regulär ist. Man setzt nämlich einfach x=ATy
und gelangt zu dem Gleichungssystem
AATy=b (275)
mit SPD-Systemmatrix A := AAT, für welches das Gauss-Seidel-Verfahren - wie gerade bewiesen - konvergiert.
Man kann nun das Gauss-Seidel mit der MatrixAund einer Folge vony-Werteny0, y1, y2, . . . ausführen und am Ende aus dem letztenykeine Approximationxk :=ATykgewinnen. Man kann die ganze Iteration aber auch direkt in denx-Werten schreiben. wie dies in der näch-sten Aufgabe geschieht:
Aufgabe 5.20
Zeigen Sie: Für i ∈ {1, . . . , n} sei ai := ATei das Transponierte der i-ten Zeile von A.
Dann lässt sich ein Schritt des Gauss-Seidel-Verfahrens für (275) in der folgenden Form schreiben, wenn man während der Iteration bei jedem Update einer Komponente von y gleich wieder alles in die zugehörigen x-Vektoren umrechnet.
For i=1:n
x=x−(aTi ai)−1(
aTi x−bi)
·ai; end
Anmerkung: Man beachte, dass im i-ten Unterschritt jeweils der ganze x-Vektor ange-passt wird.
Wir wenden uns jetzt noch einigen Aussagen über die SOR-Verfahren zu und haben für diese deren Iterationsmatrix
GSOR =I −ω(D−ωL)−1A zu untersuchen.
Allgemeiner betrachten wir200 gleich die allgemeineren Iterationsmatrizen G(ω, τ) :=I−τ(D−ωL)−1A.
Je nach Parameterwahl sind in dieser Klasse enthalten G(0,1) Jacobi-Iteration
G(0, τ) τ-extrapolierte Jacobi-Iteration G(1,1) Gauss-Seidel-Iteration
G(1, τ) τ-extrapoliertes Gauss-Seidel G(ω, ω) ω-Sor
Wir merken an, dass sich die G(ω, τ)-Methode wegen G(ω, τ) = τ
ωG(ω, ω) + (1− τ ω)I
allgemein als gedämpfte Overrelaxation interpretieren läßt201.
Wir interessieren uns hier vorwiegend für die SOR-Iteration und ihre gedämpfte Variante-Fall202.
Für den SOR-Fall beweisen wir zunächst den folgenden einfachen Satz, der die bei Defini-tion des SOR-Verfahrens getroffene Einschränkung ω∈(0,2)begründet.
Satz 5.21 (Sinnvoller Overrelaxationsparameter-Bereich) Für beliebige Matrizen A mit regulärer Diagonale gilt
ρ(G(ω, ω)) ≥ |ω−1| ∀ω,
so dass als Relaxationsparameter ω nur Werte aus (0,2) in Frage kommen.
Beweis: Es ist det(I −ωL) = 1, so dass für das charakteristische Polynom φ(λ) von G(ω, ω)gilt
φ(λ) : =det(λI−G(ω, ω))
=det((I−ωL)(λI−G(ω, ω)))
=det((I−ωL) [λI−(I−ωL)−1{(1−ω)I+ωR}])
=det((I−ωL)λ−(1−ω)I−ωR). Die Matrix
V(λ, ω) := (I−ωL)λ−(1−ω)I−ωR
wird für λ= 0 zur unteren Dreiecksmatrix, so dass sich einerseits ergibt
|φ(0)|=|ω−1|n. (276)
Auf der anderen Seite ist φ(0) aber auch das Produkt der Eigenwerte von G(ω, ω) φ(0) =
∏n i=1
λi(G(ω, ω)). (277)
Aus (276) und (277) folgt nun sofort die Behauptung.
2 Für die in diesem Abschnitt im Vordergrund stehenden SPD-Matrizen erhält man im ganzen Parameterbereich (0,2) Konvergenz. Etwas allgemeiner haben wir den folgenden Satz.
Satz 5.22 (Konvergenzsatz für SOR im SPD-Fall) Sei A∈Cn×n hermitesche und positiv definit203. Dann ist
ρ(G(ω, τ)) <1, ∀0< τ ≤ω < 2.
201Gedämpfte Relaxationsverfahren sind später bei sogenannten Mehrgitter-Verfahren von Interesse
202Für die anderen Fälle haben wir ja teilweise schon gesonderte Konvergenzresultate geliefert.
Beweis:
Sicherheitshalber wiederholen wir, was hermitesch und positiv definit für komplexe Matri-zen bedeuten: A=AH := ¯AT, xHAx >0,∀x̸= 0.
Um Übersicht über die Teilschlüsse zu behalten, teilen wir den Beweis in drei Unterab-schnitte auf:
Abschnitt A. Die Splitting-Matrix B =B(ω, τ)von G(ω, τ) ist B = 1
τD−ω τE.
Dafür rechnet man aus
B+BH −A = 1
τD−ω τE+ 1
τD−ω
τF −(D−E−F)
= 2−ω τ D+
(ω τ −1
) A Da A und D positiv definit sind und
2−ω τ >0,
(ω τ −1
)≥0,
istB +BH −A positiv definit.
Abschnitt B. Wir zeigen hierRe(spec(A−1(2B−A)))>0.
Sei dazu xein Eigenvektor vonA−1(2B−A)mit Eigenwertλ. Dann hat man (2B−A)x= λAx, woraus
xH(2B−A)x=λxHAx (278)
und durch Übergang zum konjugiert Komplexen
-xH(2BH −A)x= ¯λxHAx (279)
folgen. Durch Addition von (278) und (279) ergibt sich xH(B+BH −A)
| {z }
pos.def.
x=Re(λ)xH (A)
|{z}
pos.def.
x,
was Re(λ)>0 bedeutet.
Abschnitt C.Zu Q:=A−1(2B−A) = 2A−1B−I existiert (Q+I)−1 = (2A−1B)−1, und es ist
(Q−I)(Q+I)−1 =I−B−1A=G(ω, τ).
Ist µein Eigenwert von G(ω, τ)und xzugehöriger Eigenvektor, so ist (Q−I) (Q+I)−1x
| {z }
y̸=0
=G(ω, τ)x=µx,
so dass
(Q−I)y=µ(Q+I)y bzw.
(1−µ)Qy = (1 +µ)y ist. Wegen y̸= 0 ist sicher µ̸= 1, und wir haben daher weiter
Also ist
λ = 1 +µ 1−µ
Eigenwert von Q=A−1(2B−A) und hat damit nach B. positiven Realteil.
Umgekehrt ist nun µ= λλ+1−1 und
|µ|2 = |λ|2+ 1−2Re(λ)
|λ|2+ 1 + 2Re(λ).
und daher |µ|<1. 2
5.1.3 Konsistente Ordnung
Das grundsätzliche Wissen, dass das SOR-Verfahren konvergiert, ist natürlich schon sehr schön. Allein hätte man wohl gern näher gewusst, wie man denn den Overrelaxationspa-rameter bestmöglich wählen sollte.
Eine geschlossene Theorie existiert hier nur für eine recht enge Klasse von Matrizen, die sogenannten „konsistent geordneten“ Matrizen. Dies sind mit der Standard-Zerlegung A=D(I−E−F)diejenigen Matrizen A, bei denen die Eigenwerte der Matrix
J(α) := αE+α−1F nicht von α̸= 0 abhängen.
Solche Matrizen gibt es, wie die folgende Aufgabe zeigt.
Aufgabe 5.23
Alle Block-Tridiagonal-Matrizen, deren Blöcke wiederum Diagonalmatrizen sind, sind kon-sistent geordnet.
Tatsächlich ist diese Klasse relativ klein. In den Anwendungen gehören zu ihnen u.a. die Finite-Differenzen-Diskretisierungen einfacher elliptischer Differentialgleichungen, wenn man die Variablen auch noch in bestimmter Weise sortiert. Färbt man das Gitter der Abbildung 6 z.B. wie ein Schachbrett ein und wählt man analog der Reihenfolge der Standardnumerie-rung zeilenweise von unten nach oben und in den Zeilen von links nach rechts die Variablen zunächst auf den schwarzen Feldern und dann auf den weißen Feldern, so gelangt man zu einer 2×2-Blockmatrix mit Diagonalmatrizen als Diagonalblöcke.
Für die Klasse der konsistent geordneten Matrizen gibt es zwei sehr hübsche Resultate, die wir hier ohne Beweis angeben (Beweise hierzu findet man etwa bei [STBU] ).
Satz 5.24 (Vergleichssatz für Jacobi- und Gauss-Seidel-Verfahren) Wenn A konsistent geordnet ist, so gilt
ρ(GGS) = ρ(GJ)2.
Für konsistent geordnete Matrizen konvergiert das Gauss-Seidel-Verfahren also asympto-tisch doppelt so schnell wie das Jacobi-Verfahren, wenn letzteres denn konvergiert.