Lineare Iterationsverfahren für M-Matrizen

Lineare Iterationsverfahren für M-Matrizen

Simon Renner

22.06.2010

INHALTSVERZEICHNIS 2

Inhaltsverzeichnis

1 Einleitung 3

2 M-Matrizen 3

2.1 Vorbemerkungen der linearen Algebra . . . 3 2.2 Kriterien für M-Matrizen . . . 5

3 Reguläre Aufspaltungen 8

4 Anwendungen 10

4.1 Blockweise Iterationsverfahren . . . 10 4.2 Konvergenz für M-Matrizen . . . 11

1 EINLEITUNG 3

1 Einleitung

Zur Lösung der Gleichung A· x = b können entweder direkte oder iterative Verfahren verwendet werden. In direkten Verfahren wird die Gleichung in endlich vielen Schritten bis auf Rundungsfehler exakt gelöst. Ein iteratives Verfahren generiert in jedem Schritt ein neuesx, das mit jedem Schritt näher an die exakte Lösung heranrückt. Wie viele Schritte für eine bestimmte Genauigkeit benötigt werden ist u.a ein Maßstab für die Qualität des Verfahrens. Ist die MatrixAsehr groß und dünn besetzt, so sind iterative Verfahren aufgrund des geringeren Rechen- und Speicheraufwandes die bessere Wahl.

Wird die Matrix A in D+ L+ R mit D = diag(a_αβ),R = (a_αβ) α < β,L = (a_αβ) α > β aufgespalten, so entsteht im FalleDregulär folgendes Iterationsverfahren:

x^t =−D⁻¹(L+R)x^t−1+D⁻¹b

Dies ist das sog. Jacobi Verfahren. Ein weiteres wichtiges Verfahren ist das Gauß-Seidel Verfahren:

x^t =−(D+L)⁻¹Rx^t−1+(D+L)⁻¹b

Es stellt sich nun die Frage, wann ein solches Verfahren konvergiert. Eine Antwort hierauf liefert der folgende Satz.

Satz 1. Die durch

x^t =Bx^t−1+b

erzeugten Iterierten x^t ∈ R,t = 1,2... konvergieren genau dann für jeden Startwert x⁰ gegen die Lösung x∈Rder Fixpunktgleichung, wenn

ρ(B)<1 ρ(B)bezeichnet den Spektralradius der Matrix.

Die Eigenschaften der IterationsmatrixBhängen offensichtlich von der AusgangsmatrixA ab. Wir wollen nun die Eigenschaften solcher Verfahren bei der Anwendung auf M-Matrizen untersuchen.

2 M-Matrizen

2.1 Vorbemerkungen der linearen Algebra

Zunächst einige Hilfssätze und Definitionen, die in den nächsten Kapiteln benötigt werden:

Hilfssatz 1. Für eine Matrix A∈K^I×Igilt, falls A≥0:

1. ρ(A)∈σ(A)

2. λ=ρ(A)hat einen positiven Eigenvektor x 3. ρ(B)≥ρ(A)für alle B≥A

Man beachte, dass die Ordnungsrelation≥für jedes Element definiert ist:A≥B⇔ a_αβ ≥b_αβ für alleα, β∈I.

2 M-MATRIZEN 4 Hilfssatz 2. Für eine Matrix A∈K^I×I,A≥0gilt:

ρ(A)<1⇔

∞

X

ν=0

A^ν =(1−A)⁻¹

Definition 1. Sei A∈K^I×Ieine Matrix. Als Graph G(A)der Matrix wird folgende Menge bezeichnet:

G(A)={(α, β)∈I×I :a_αβ ,0}

Definition 2. Eine Matrix A ∈ K^I×I heißt irreduzibel, wenn in G(A)jedesα ∈ I mit jedemβ ∈ I verbunden ist. Andernfalls heißt sie reduzibel.

Diese Formulierung ist sehr nützlich, da sie am Rechner leicht überprüft werden kann.

Eine andere Formulierung der Irreduzibilität, die wir später noch brauchen werden, liefert die folgende Bemerkung:

Bemerkung 1. Eine Matrix A∈K^I×I ist genau dann reduzibel, wenn man die Indizes so anordnen kann, dass A die Blockgestalt

A= A11 A12

0 A22

!

wobei A11, A22nichtleere quadratische Blöcke sind.

Definition 3. Eine Matrix A∈K^I×I heißt stark diagonaldominant, wenn:

|a_αα|>X

β∈I β,α

a_αβ

∀α∈I

Definition 4. Eine Matrix A∈K^I×I heißt schwach diagonaldominant, wenn:

|a_αα| ≥X

β∈I β,α

a_αβ

∀α∈I

Definition 5. Eine Matrix A∈K^I×Iheißt irreduzibel diagonaldominant, falls A irreduzibel, schwach diagonaldominant ist und außerdem gilt:

|a_αα|>X

β∈I β,α

a_αβ

für mindestens einα∈I

Definition 6. Sei A∈ K^I×Ieine Matrix und fürγ∈I sei

G_γ B {β∈ I :γmitβim Matrixgraphen G(A)verbunden}. A heißt im wesentlichen diagonaldomi- nant, wenn A schwach diagonaldominant ist und für alleγ ∈I gilt:

|a_αα|>X

β∈I β,α

a_αβ

für mindestens einα∈G_γ

Hilfssatz 3. Für eine Matrix B∈K^I×Igilt, falls|B| ≤A∈K^I×I: ρ(B)≤ρ(A)

Hilfssatz 4. Für eine Matrix A∈K^I×Igelten die Implikationen:

stark diagonaldominant⇒irreduzibel diagonaldominant⇒im wesentlichen diagonaldominant

2 M-MATRIZEN 5

2.2 Kriterien für M-Matrizen

Definition 7. Eine Matrix A∈R^I×I heißt M-Matrix, falls 1. a_αα >0für alleα∈I

2. a_αβ ≤0für alleα,β 3. A regulär und A⁻¹≥0

Die ersten zwei Kriterien dieser Definition sind leicht nach zu prüfen. Da Matrizeninversion jedoch ein numerisch recht aufwändiger Prozess ist, werden wir nach anderen Kriterien suchen, um die Regularität vonAund Positivität vonA⁻¹fest zu stellen. Einen ersten Schritt in diese Richtung tut der folgende Satz.

Satz 2. A∈R^I×Ierfülle: a_αβ ≤0für alleα,β. D= diag(a_αα : α∈I)bezeichne die Diagonale von A.

Dann sind die folgenden Aussagen äquivalent:

1. A regulär und A⁻¹≥0

2. a_αα >0, M:=1−D⁻¹A≥0, ρ(M)<1

Im Fall 1 oder 2 ist A eine M-Matrix. Umgekehrt gilt für jede M-Matrix Aussage 2.

Beweis. ⇒:

a_αα >0: Seis^γdieγ ∈Ientsprechende Spalte vonA. DaA⁻¹A=1

⇒A⁻¹s^γ =e^γ(Einheitsvektor) Daa_αβ ≤0 ist, gilt, fallsa_αα ≤0 :a^γ ≤0

⇒A⁻¹a^γ ≤0⇒e^γ ≤0

⇒a_αα >0 M=1−D⁻¹A≥0: Daa_αα >0⇒D≥0 und regulär.

⇒ A⁰ :=D⁻¹A⇒A⁰⁻¹ =A⁻¹D≥0

M=1−A⁰hat die DiagonaleinträgeM_αα =1−1=0. Die restlichen Einträge sind M_αβ =0−a⁻_αα¹a_αβ ≥0

⇒M≥0

ρ(M)<1: Nach Hilfssatz 1 existiert zuλ:=ρ(M) ein positiver Eigenvektorx≥0:

Mx=λx⇔(1−A⁻¹D)x=A⁻¹Dλx

⇔(A⁻¹D−A⁻¹Dλ)x=x

⇔A⁰⁻¹(1−λ)x=x Da sowohlA⁰⁻¹ ≥0 als auchx≥0gilt:

1−λ≥0⇒0≤ρ(M)=λ <1

⇐:

Aus Hilfssatz 2 folgt

(1−M)⁻¹≥0

⇒ 0≤(1−M)⁻¹D⁻¹ =(D⁻¹A)⁻¹D⁻¹ =A⁻¹DD⁻¹ =A⁻¹

2 M-MATRIZEN 6 Man kann den Satz sogar noch ein wenig verschärfen:

Satz 3. Für eine Matrix A∈R^I×Igelte: a_αβ ≤0für alleα,β. D bezeichne die Diagonale der Matrix, MB1−D⁻¹A. Dann sind die folgenden Aussagen äquivalent:

1. A regulär und A⁻¹>0

2. a_αα >0, M≥0, ρ(M)<1, M irreduzibel Beweis. ⇒:

WäreAreduzibel, so gäbe es eine Blockstruktur wie in Bemerkung 1:

A= A11 A12

0 A22

!

Wobei die inverse MatrixCBA⁻¹ folgende Struktur hätte:

C= A⁻₁₁¹ −A⁻₁₁¹A12A⁻₂₂¹ 0 A⁻₂₂¹

!

wobeiC21 =0, was im Widerspruch steht zuA⁻¹ >0. Also istAirreduzibel. Da sichMundA nur auf ihrer Diagonalen unterscheiden gilt:G(A)=G(M), somit ist auchMirreduzibel.

⇐:

Aus dem Beweis von Satz 2 entnehmen wir:

A⁻¹ =(

∞

X

ν=0

M^ν)D⁻¹ DaP^∞

ν=0M^νundD⁻¹ >0⇒A⁻¹>0

Satz 4. Die Matrix A∈K^I×I sei stark diagonaldominant, im wesentlichen diagonaldominant oder irreduzibel diagonaldominant. Dann gilt:

ρ(M)<1

für die Jacobi Iterationsmatrix M B 1−D⁻¹A, wobei D die Diagonale von A ist. Gilt außerdem a_αα >0und a_αβ ≤0∀α,β, so ist A eine M-Matrix.

Beweis. Für stark diagonaldominante Matrizen istD offensichtlich regulär. Dies gilt auch für irreduzibel diagonaldominante Matrizen, da in jeder Spalte ein Eintrag außerhalb der Diagonale ungleich 0 ist. Bei im wesentlichen diagonaldominanten Matrizen lässt sich für jedesa_αα einβ∈G_γfinden, sodass auch hierDregulär ist.

Wir definieren nun eine neue MatrixM⁰ B|M|mit den ElementenM⁰_αα =0 für alleα=βund M⁰_αβ = |a_αβ/aαα|fürα,β. Nach Hilfssatz 3 giltρ(M)≤ρ(M⁰) weshalb es reichtρ(M⁰)<1 zu zeigen.

Da M⁰ ≥ 0 ⇒ λ B ρ(M⁰) ∈ σ(M⁰) und es gibt einen zugehörigen positiven normierten Eigenvektorx: kxk_∞=1.

Seiαein Index mitx_α =1

⇒λ=λxα =(M⁰x)_α

2 M-MATRIZEN 7

=









 X

β,α

|a_αβ|x_β











/|a_αα| ≤









 X

β,α

|a_αβ|











/|a_αα| ≤1

Hier wurde benutzt, dassM⁰_αα =0 um die Summe zu reduzieren,kxk_∞ =1 schafft die erste Umgleichung, die schwache Diagonaldominanz die zweite.

Fallsλ < 1⇒ρ(M⁰)<1.

Fallsλ =1⇒x_γ=1 für alleγ∈G_α

Da stark diagonaldominant⇒irreduzibel diagonaldominant⇒im wesentlichen diagonaldo- minant (Hilfssatz 4) gibt es einen Indexγ ∈G_αmit|a_γγ>P

β,γ|a_βγ|, so dass wegenx_γ =x_β =1 fürγ, β∈G_αgilt:

λ=λx_γ =









 X

β,γ

|a_γβ|x_β











/|a_γγ|=









 X

β,γ











/|a_γγ|<1

Wenden wir Satz 2 auf diesen Satz an, so erhalten wir:

Satz 5. Eine Matrix A ∈R^I×I erfülle a_αα >0und a_αβ ≤0∀α,βund sei irreduzibel diagonaldomi- nant. Dann ist A eine M-Matrix mit A⁻¹ >0.

Wir haben nun also mehrere Kriterien, die am Rechner leicht zu überprüfen sind, gewonnen um fest zu stellen, ob eine gegebene MatrixAeine M-Matrix ist. Für einen späteren Beweis benötigen wir aber noch folgenden Hilfssatz.

Hilfssatz 5. Sei A∈R^I×Ieine M-Matrix und für B∈R^I×Igelte: B≥A sowie b_αβ ≤0für alleα,β.

Dann ist auch B eine M-Matrix und es gilt

0≤B⁻¹ ≤A⁻¹

Beweis. Da B ≥ A gilt b_αα ≥ a_αα für alle α ∈ I. Somit sind die Elemente der Inversen Diagonalmatrix vonB1/bαα ≤1/aαα. Es gilt also

0≤D⁻_B¹ ≤D⁻¹ sowie, daM_αα =0

0≤MB ≤M Daρ(A)<1 gilt nach Hilfssatz 2

∞

X

ν=0

M^ν =(1−M)⁻¹ DaMB durchMmajorisiert wird, gilt:

∞

X

ν=0

M^ν_B konvergiert und erneute Anwendung von Hilfssatz 2 ergibt

∞

X

ν=0

M^ν_B =(1−MB)⁻¹ ⇔ρ(MB)<1

3 REGULÄRE AUFSPALTUNGEN 8

3 Reguläre Aufspaltungen

Die AufspaltungA=W−Reiner Matrix induziert das IterationsverfahrenWx^t+1 =Rx^t+b fallsWregulär ist.

Definition 8. Die Matrix W∈R^I×Ibeschreibt eine reguläre Aufspaltung von A ∈R^I×I, falls W regulär, W⁻¹ ≥0, W ≥A

Die Iterationsmatrix des induzierten Iterationsverfahrens ist M = W⁻¹R, zudem wird durch die Definition impliziert:M≥0 für reguläre Aufspaltungen. Mit Hilfe dieser Bedingung lässt sich die Definition etwas abschwächen:

Definition 9. Die Matrix W ∈R^I×Ibeschreibt eine schwach reguläre Aufspaltung von A∈R^I×I, falls W regulär, W⁻¹≥0, M=W⁻¹R≥0

Die Konvergenzeigenschaften eines solchen induzierten Iterationsverfahrens für M- Matrizen werden im folgenden Satz deutlich.

Satz 6. Die Matrix A∈R^I×I sei eine M-Matrix. W beschreibe eine schwach reguläre Aufspaltung von A. Dann konvergiert das induzierte Iterationsverfahren und es gilt:

ρ(M)=ρ(W⁻¹R)= ρ(A⁻¹R) 1+ρ(A⁻¹R) <1

Beweis. MitCBA⁻¹Rreicht es, die Gleichheit fürρ(W⁻¹R)= ρ(C)/(1+ρ(C)) zu zeigen. Da M≥0 gilt:

0≤M=W⁻¹R

⇔[A⁻¹W]⁻¹A⁻¹WM=[A⁻¹W]⁻¹A⁻¹R

=[A⁻¹(A+R)]⁻¹A⁻¹R=[1+C]⁻¹C

Wir verwenden nun wieder Hilfssatz 1: Zuλ=ρ(M)∈σ(M) gehört wegenM≥0 ein positiver Eigenvektorx.

⇒λx=Mx=(1+C)⁻¹Cx

⇔ λx+λCx=Cx Angenommenλ=1⇒x=0.

⇒Cx= λ 1−λx FallsC≥0⇒ ^λ

1−λ ≥0 also 0 ≤λ=ρ(M)<1. Daher zeigen wir nun, dassC≥0:

DaM≥0 undW⁻¹ ≥0 gilt 1. 0≤Pm−1

ν=0 M^ν W⁻¹ 2. W⁻¹=(1−M)A⁻¹ ≤A 3. Pm−1

ν=0 M^ν(1−M)=1−M^m

3 REGULÄRE AUFSPALTUNGEN 9 Einsetzen von zwei in eins liefert:

0≤









m−1

X

ν=0

M^ν









(1−M)A⁻¹ Kombination mit der dritten Gleichung ergibt:

0≤(1−M^m)A⁻¹ ≤A⁻¹ ⇒0≤M^mA⁻¹ ≤A⁻¹

M ist also beschränkt ⇒ λ = ρ(M) ≤ 1. Da λ = 1 jedoch zum Widerspruch führt, muss λ=ρ(M)<1 gelten. Dies impliziert unter Verwendung von Hilfssatz 2:

C=A⁻¹R=[W(1−M)]⁻¹R

=(1−M)⁻¹W⁻¹R=









∞

X

ν=0

M^ν







M≥0 Nun zuρ(M)=ρ(C)/(1−ρ(C)):

Die Gleichung

Cx= λ 1−λx bedeutet, dass

λist Eigenwert vonM⇔µ= λ

1−λ ist Eigenwert vonC

Daµ= ₁^λ−λ monoton mitλwächst, ist|µ|=µmaximal fürλ=ρ(M)∈σ(M). DaC≥0 liefert Hilfssatz 1, dass der größte Eigenwert gleich dem Spektralradius ist.

⇒ρ(C)= ρ(M) 1−ρ(m)

⇔ρ(M)= ρ(C) 1+ρ(C)

Der folgende Satz erlaubt einen Vergleich der Konvergenzgeschwindigkeiten zweier regulärer Aufspaltungen.

Satz 7. Für die Matrix A∈R^I×I gelte A≥0. Durch W1und W2 seien zwei reguläre Aufspaltungen gegeben. Wenn W1und W2in der Form

A≤W1 ≤W2

vergleichbar sind, lassen sich auch die zugehörigen Konvergenzraten vergleichen:

0≤ρ(M1)≤ρ(M2)<1 wobei Mi BW⁻_i¹Ri, Ri BWi−A

Beweis. Die MatrizenBBA⁻¹R1undCBA⁻¹R2erfüllen 0≤B≤Cund damit nach Hilfssatz 3: 0≤ρ(B)≤ρ(C). Anwendung von Satz 6 liefert:

0≤ρ(M1)= ρ(B)

1+ρ(B) ≤ ρ(C)

1+ρ(C) =ρ(M2)<1

4 ANWENDUNGEN 10

4 Anwendungen

4.1 Blockweise Iterationsverfahren

Sowohl das Jacobi- als auch das Gauß-Seidel-Verfahren können in einer Blockversion durch- geführt werden:

Definition 10. Für eine Matrix A∈R^I×Isei eine Zerlegung der Indexmenge I in disjunkte, nichtleere Teilmengen gegeben: I = S

κ∈BI_κ, wobei B die Indexmenge der Blöcke ist. Mit D wird nun die Blockdiagonale der Matrix bezeichnet:

DBblockdiag(A)=blockdiag(A^κκ:κ ∈B)

Dabei sind A^κκ die Diagonalblöcke der Matrix. Das blockweise Jacobi-Verfahren wird durch die Iteration mit

WBD,RBD−A beschrieben.

Bemerkung 2. Offenbar ist dieses Verfahren genau dann wohldefiniert, wenn D regulär ist, in diesem Fall also, wenn

A^κκregulär für alleκ ∈B

Definition 11. Für eine Matrix A∈R^I×Isei eine Zerlegung der Indexmenge I in disjunkte, nichtleere Teilmengen gegeben: I = S

κ∈BI_κ, wobei B die Indexmenge der Blöcke ist. Mit D wird nun die Blockdiagonale der Matrix bezeichnet:

DBblockdiag(A)=blockdiag(A^κκ:κ ∈B)

Dabei sind A^κκ die Diagonalblöcke der Matrix. Mit L B A^κφ für κ > φ wird eine strikte untere Blockdreiecksmatrix, mit R⁰ B A^κφfürκ < φ eine strikte obere Blockdreiecksmatrix definiert. Das blockweise Jacobi-Verfahren wird durch die Iteration mit

WBD−L,RBR⁰ beschrieben.

Bemerkung 3. Auch dieses Verfahren ist genau dann wohldefiniert, wenn D regulär ist, denn dann ist auch

D−L=A^κφregulär für alleκ ≥φ

Da bei den blockweisen Verfahren in jedem Schritt Gleichungssysteme zu lösen sind, haben sie einen höheren Rechenaufwand (Faktor 1,4). Wie wir im nächsten Kapitel sehen werden konvergieren sie jedoch schneller als die punktweisen Verfahren, da sich die Ausgangsmatrix Aund die AufspaltungsmatrixWweniger unterscheiden.

4 ANWENDUNGEN 11

4.2 Konvergenz für M-Matrizen

Wir wollen nun die Konvergenzeigenschaften des Jacobi- und des Gauß-Seidel-Verfahrens für M-Matrizen untersuchen.

Satz 8. Sei A ∈R^I×Ieine M-Matrix. Dann konvergiert sowohl das punktweise als auch das blockweise Jacobi-Verfahren, wobei letzteres schneller konvergiert:

ρ(M^BlockJac)≤ρ(M^Jac)<1 Ist D die Diagonale D^pkt bzw. D^block von A, so gilt:

D beschreibt eine reguläre Aufspaltung Beweis. DaAeine M-Matrix ist, gilt

D≥0

Nach Hilfssatz 5 ist D eine M-Matrix und somit D⁻¹ ≥ 0.Daraus folgt: D beschreibt eine reguläre Aufspaltung. DaD^pkt ≥D^blockfolgt mit Satz 7

0≤ρ(M^BlockJac)≤ρ(M^Jac)<1

Satz 9. Sei A ∈R^I×Ieine M-Matrix. Dann konvergiert sowohl das punktweise als auch das blockweise Gauß-Seidel-Verfahren, wobei letzteres schneller konvergiert:

ρ(M^BlockGS)≤ρ(M^GS)<1

Ist D die Diagonale D^pkt bzw. D^block von A und L die untere Dreiecks- bzw. Blockdreiecksmatrix, die Durch A =D+L+R entsteht, so gilt:

D−L beschreibt eine reguläre Aufspaltung Beweis. DaAeine M-Matrix ist, gilt

D≥0

undL≤ 0. Nach Hilfssatz 5 istD−Leine M-Matrix und somit (D−L)⁻¹ ≥ 0.Daraus folgt:

D−Lbeschreibt eine reguläre Aufspaltung. Somit istD−Leine reguläre Aufspaltung. Da (D−L)^pkt ≥(D−L)^block folgt mit Satz 7

0≤ρ(M^BlockGS)≤ρ(M^GS)<1

Nachdem nun die Konvergenz für block-und punktweise Jacobi- und Gauß-Seidel- Verfahren gezeigt ist, liefert der folgende Satz noch einen Vergleich der beiden Verfahren.

Satz 10. Für eine M-Matrix A∈R^I×Igilt:

1. ρ(M^GS)≤ρ(M^Jac)<1

2. ρ(M^BlockGS)≤ρ(M^BlockJac)<1

Beweis. Für die beiden Aufspaltungen gilt im punktweisen wie im blockweisen Fall:

W^GS =D−L≤D=W^Jac

Satz 8,9 (Konvergenz der Verfahren) und 7 (Ungleichung) liefern dann die Aussage des

Satzes.

LITERATUR 12

Literatur

[1] Wolfgang Hackbusch,Iterative Lösung großer schwachbesetzter Gleichungssysteme, Teubner- Verlag, Stuttgart, 2. überarb. und erw. Aufl., 1994

[2] Rolf Rannacher,Vorlesungsskriptum Einführung in die Numerische Mathematik, 2006

[3] Richard S. Varga,Matrix iterative Analysis, Prentice-Hall-Verlag, Englewood Cliffs NJ, 1962