n 1) n 1) Dahmen-Reusken Kapitel 13 1 - igpm.rwth-aachen.de · Fill-in bei einem direkten L¨oser...

KAPITEL 13. Große dunnbesetzte LGS, iterative Loser

Erstes Beispielproblem: die diskretisierte Poisson-Gleichung

A1x = b, A1 ∈ Rm×m, m := (n− 1)2 , nh = 1,

wobei

A1 =1

h2

T −I

−I T −I ∅. . . . . . . . .

∅ −I T −I

−I T

, T =

4 −1−1 4 −1 ∅

. . . . . . . . .

∅ −1 4 −1−1 4

.

T ∈ R(n−1)×(n−1) und I ist die (n− 1) × (n− 1)-Identitatsmatrix.

Dahmen-Reusken Kapitel 13 1

Zweites Beispielproblem:

die diskretisierte Konvektions-Diffusionsgleichung

A2x = b, A2 ∈ Rm×m, m := (n− 1)2 ,

wobei

A2 =1

h2

T −nI−sI T −nI ∅

. . . . . . . . .

∅ −sI T −nI−sI T

, T =

z −o−w z −o ∅

. . . . . . . . .

∅ −w z −o−w z

,

mit

n = ε− hs−, s = ε+ hs+ ,

o = ε− hc−, w = ε+ hc+ ,

z = n+ s+ o+ w = 4ε+ h(|c| + |s|) ,und

c := cos(β), c+ := max(c,0), c− := min(c,0),

s := sin(β), s+ := max(s,0), s− := min(s,0).


13.2 Eigenschaften von Steifigkeitsmatrizen

Hochdimensionalitat:

Anzahl der Unbekannten ist bei d Ortsvariablen proportional zu h−d.

Dunnbesetztheit:

Steifigkeitsmatrizen sind immer dunnbesetzt, d. h. nur eine gleichmaßig

beschrankte Anzahl von Eintragen pro Zeile ist von Null verschieden.

0 50 100 150 200

0

50

100

150

200

nz = 1065


Blockstruktur:

Bei Diskretisierungen auf regelmaßigen Rechteckgittern haben die sich

ergebenden Steifigkeitsmatrizen oft eine regelmaßige Blockstruktur.

Schlechte Kondition:

Die Konditionszahl von Steifigkeitsmatrizen nimmt oft rasch zu, wenn

die Schrittweite des Gitters kleiner wird.

Symmetrie, Positiv-Definitheit:

Fur ein symmetrisches elliptisches Randwertproblem ist die entpre-

chende Steifigkeitsmatrix oft symmetrisch positiv definit.

Diagonaldominanz:

Die Matrizen A1 und A2 sind beide irreduzibel diagonaldominant.

Die durch Finite-Elemente-Methoden und Finite-Volumen-Verfahren

erzeugten Steifigkeitsmatrizen sind im allgemeinen nicht irreduzibel

diagonaldominant.


Fill-in bei einem direkten Loser

Bemerkung 13.2. Als direkter Loser fur ein Gleichungssystem mit der

Matrix A1 wurde sich die Cholesky-Zerlegung A1 = LDLT anbieten.

Nichtnulleintrage der Matrizen A1 und L:

h 1/16 1/32 1/64 1/128

nz(A1) 1065 4681 19593 80137 ≈ 5h−2

nz(L) 3389 29821 250109 2048509 ≈ h−3

Das Verhalten nz(L) ≈ h−3 gegenuber nz(A1) ≈ 5h−2 zeigt den ungunsti-

gen Effekt des”fill-in“.


Nichtnullelemente von L bei der Cholesky-Zerlegung der Matrix A1:

0 50 100 150 200

0

50

100

150

200

nz = 3389


13.3 Lineare Iterationsverfahren

Aufgabe:

Fur A ∈ Rn×n (nichtsingular) und b ∈ Rn ist das Gleichungssy-

stem

Ax = b

zu losen. Wir nehmen an, daß n groß (n > 10000) und A dunn-

besetzt ist.

Typische Beispiele sind die Matrizen A1 und A2.


Einfacher Ansatz: das Gleichungssystem als Fixpunktgleichung

x = x + C(b − Ax) =: Φ(x),

wobei C ∈ Rn×n eine geeignet zu wahlende nichtsingulare Matrix ist.

Dies fuhrt auf die Fixpunktiteration

xk+1 = Φ(xk) = xk + C(b − Axk)

= (I − CA)xk + Cb, k = 0,1,2, . . . .

Fur den Fehler ek := xk − x∗ ergibt sich

ek+1 = xk+1 − x∗ = Φ(xk) − Φ(x∗) = (I − CA)ek,

also

ek = (I − CA)ke0, k = 0,1,2, . . . .


Satz 13.3. Dieses Verfahren konvergiert fur jeden Startwert x0

gegen die Losung x∗ des Gleichungssystems genau dann, wenn

ρ(I − CA) < 1

gilt, wobei ρ(I − CA) der Spektralradius, also der betragsmaßig

großte Eigenwert, von I − CA ist.

Folgerung 13.5. Fur eine beliebige Vektornorm ‖·‖ mit zugehori-

ger Operatornorm gilt:

‖xk − x∗‖ ≤ ‖I − CA‖k ‖x0 − x∗‖, k = 0,1,2, . . .

Die Konvergenz ist folglich fur jeden Startwert x0 gesichert, falls

fur irgendeine Norm die Bedingung

‖I − CA‖ < 1

erfullt ist.

Die Große ‖I−CA‖ heißt die Kontraktionszahl der Fixpunktiteration.


Konvergenzrate

ρ(I − CA) ist ein sinnvolles Maß fur die Konvergenzgeschwindigkeit.

Je kleiner ρ(I − CA) ist, desto schneller wird der Fehler reduziert. Es

gilt:

σk =

(‖ek‖‖e0‖

)1k → |λ1| = ρ(I − CA) fur k → ∞.

Fur eine hinreichend große Anzahl k von ausgefuhrten Iterationsschrit-

ten ist ρ(I − CA) etwa die mittlere Fehlerreduktion σk.

Um den Anfangsfehler ‖e0‖ um den Faktor 1e zu reduzieren, sind

k = (− lnσk)−1 Iterationsschritte notig.

Die sogenannte asymptotische Konvergenzrate

− ln(ρ(I − CA))

Ist ein sinnvolles Maß fur die asymptotische Konvergenzgeschwindig-

keit.


Wie sollte man C wahlen?

Es geht um folgenden Balanceakt:

Zu A ∈ Rn×n wahle C ∈ Rn×n so, daß

• A−1 durch C genugend gut in dem Sinne approximiert wird, daß

ρ(I − CA) moglichst klein ist,

• die Operation

y 7→ Cy

mit moglichst geringem Aufwand durchfuhrbar ist.

Beachte:

Bei den meisten iterativen Verfahren wird die Operation y 7→ Cy durch-

gefuhrt, ohne daß die Matrix C explizit berechnet wird.


Komplexitat eines iterativen Verfahrens

Um unterschiedliche Verfahren miteinander vergleichen zu konnen, ge-

hen wir davon aus, daß

• eine Klasse von Gleichungssystemen Ax = b vorliegt

(z.B. A1x = b wie oben mit n = 1h ∈ N beliebig),

• vorgegeben ist, mit welchem Faktor R ein (beliebiger) Startfehler

reduziert werden soll (‖ek‖ ≤ ‖e0‖R ).

Komplexitat eines iterativen Verfahrens: die Großenordnung der An-

zahl arithmetischer Operationen, die benotigt werden, um fur das vor-

liegende Problem (aus der Problemklasse) eine Fehlerreduktion um

den Faktor R zu bewirken.

Liegt dem Gleichungssystem eine Diskretisierung mit Gitterweite h und

einem Diskretisierungsfehler von der Ordnung O(hℓ) zugrunde, ware

R−1 ≈ hℓ sinnvoll.


13.3.2 Das Jacobi-Verfahren

Beim Jacobi- oder Gesamtschrittverfahren:

C = (diag(A))−1.

Die Iterationsvorschrift lautet

xk+1 = (I − D−1A)xk + D−1b.

In Komponentenschreibweise ergibt sich:

Algorithmus 13.7 (Jacobi-Verfahren).

Gegeben: Startvektor x0 ∈ Rn. Fur k = 0,1, . . . berechne:

ai,ixk+1i = −

n∑

j=1j 6=i

ai,jxkj + bi, i = 1,2, . . . , n.

Beim Jacobi-Verfahren wird die i-te Gleichung nach der i-ten Unbe-

kannten xi aufgelost, wobei fur die ubrigen Unbekannten (xj, j 6= i)

Werte aus dem vorherigen Iterationsschritt verwendet werden.


Rechenaufwand. Der Rechenaufwand pro Iterationsschritt ist beim

Jacobi-Verfahren angewandt auf eine dunnbesetzte Matrix A ∈ Rn×n

vergleichbar mit einer Matrix-Vektor-Multiplikation Ax, beansprucht

also O(n) Operationen.

Konvergenz

Satz 13.9. Fur das Jacobi-Verfahren gelten folgende Konvergenz-

kriterien:

• Falls sowohl A als auch 2D−A symmetrisch positiv definit sind,

folgt

ρ(I − D−1A) < 1.

• Falls A irreduzibel diagonaldominant ist, gilt

ρ(I − D−1A) ≤ ‖I − D−1A‖∞ < 1.


Beispiel 13.10.

Das Diffusionsproblem mit der Matrix A1. Es gilt

ρ(I − CA)= ρ(I − D−1A1)

= max{ |1 − 14h

2λ| | λ Eigenwert von A1}= 1 − 2 sin2(1

2πh) = cos(πh) ≈ 1 − 12π

2h2.

Fur die asymptotische Konvergenzrate gilt somit

− ln(ρ(I − D−1A1)) ≈ − ln(1 − 12π

2h2) ≈ 12π

2h2.

Um einen Startfehler um einen Faktor R zu reduzieren, sind etwa

K =− lnR

ln ρ(I − D−1A1)≈ 2

π2h2lnR

Iterationsschritte erforderlich.

Beachte: Die geschatzte Anzahl der notwendigen Iterationsschritte

wachst in Abhangigkeit der Schrittweite h wie die Konditionszahlen

der betreffenden Matrizen.


Sei # die Anzahl von Iterationsschritten, die zur Reduktion des Start-

fehlers um einen Faktor R = 103 benotigt werden, und K die theore-

tische Schatzung von # aus:

K =− ln 103

ln(cosπh)≈ 2

π2h2ln 103.

Ergebnisse:

h 1/40 1/80 1/160 1/320# 2092 8345 33332 133227

K 2237 8956 35833 143338

Komplexitat. Dam ∼ h−2, folgt, daß K ∼ m Iterationsschritte benotigt

werden. Da jeder Schritt einen zu m proportionalen Aufwand erfordert,

ergibt sich, daß die Komplexitat des Jacobi-Verfahrens (fur diese Pro-

blemstellung) etwa cm2 ist.


Beispiel 13.11.

Das Konvektions-Diffusionsproblem mit der Matrix A2.

Wir fuhren fur β = π6 ein Experiment wie in Beispiel 13.10 durch.

Ergebnisse (#):

h 1/40 1/80 1/160 1/320ε = 1 2155 8611 34415 137597

ε = 10−2 194 587 1967 7099

ε = 10−4 75 148 293 597

Fur den Fall ε = 1 (dominante Diffusion) gilt κ2(A2) ∼ 1h2 und fur

ε = 10−4 (dominante Konvektion) gilt κ2(A2) ∼ 1h.

Diese Beziehungen entsprechen den Wachstumsfaktoren 4 und 2 in

den Zeilen fur ε = 1 bzw fur ε = 10−4.


13.3.3 Das Gauß-Seidel-Verfahren

Zerlegung:

A = D − L − U.

Mit C := (D − L)−1 ergibt sich:

xk+1 = (I − (D − L)−1A)xk + (D − L)−1b.

Aquivalent:

(D − L)xk+1 = Uxk + b.

In Komponentenschreibweise heißt dies

Algorithmus 13.12 (Gauß-Seidel-Verfahren).

Gegeben: Startvektor x0 ∈ Rn. Fur k = 0,1, . . . berechne:

xk+1i = a−1

i,i

(

bi −i−1∑

j=1

ai,jxk+1j −

n∑

j=i+1

ai,jxkj

)

, i = 1,2, . . . , n.

Bei der Berechnung der i-ten Komponente von xk+1 verwendet man

die bereits vorher berechneten Komponenten der neuen Annaherung.


Rechenaufwand. Fur eine dunnbesetzte Matrix A ∈ Rn×n ist

der Rechenaufwand pro Iterationsschritt bei der Gauß-Seidel-

Methode vergleichbar mit dem Aufwand beim Jacobi-Verfahren,

betragt also O(n) Operationen.

Konvergenz

Satz 13.13. Fur das Gauß-Seidel-Verfahren gelten folgende

Konvergenzkriterien:

• Falls A symmetrisch positiv definit ist, folgt

ρ(I − (D − L)−1A) < 1.

• Falls A irreduzibel diagonaldominant ist, gilt

ρ(I − (D − L)−1A) ≤ ‖I − (D − L)−1A‖∞ < 1.


Beispiel 13.15.

Das Diffusionsproblem mit der Matrix A1. Dafur gilt:

ρ(I − CA1) = ρ(I − (D − L)−1A1) = (ρ(I − D−1A1))2

Hieraus ergibt sich:

ρ(I − CA1) = cos2(πh) ≈ 1 − π2h2 .

− ln(ρ(I − (D − L)−1A1)) ≈ − ln(1 − π2h2) ≈ π2h2.

Um einen Startfehler um einen Faktor R zu reduzieren sind (asympto-

tisch) etwa K Iterationsschritte erforderlich, wobei

K =− lnR

ln(ρ(I − (D − L)−1A1))≈ 1

π2h2lnR.

Ergebnisse:

h 1/40 1/80 1/160 1/320# 1056 4193 16706 66694

K 1119 4478 17916 71669


Komplexitat

Fur das Gauß-Seidel-Verfahren angewandt auf die diskrete Poisson-

Gleichung

A1x = b, A1 ∈ Rm×m, m = (n− 1)2, n =

1

h

ist die Komplexitat cm2 Operationen.


Beispiel 13.16.

Die Resultate des Gauß-Seidel-Verfahrens hangen von der Anordnung

der Unbekannten ab.

Beispiel: Die diskrete Konvektions-Diffusionsgleichung mit β = π6.

h = 1160, R = 103 und x0, b wie in Beispiel 13.11.

Gauß-Seidel-Verfahren mit Standardanordnung:

ε 1 10−2 10−4

# 17197 856 14

Gitterpunkte in umgekehrter Reihenfolge:

ε 1 10−2 10−4

# 17220 1115 285

Fur ein Problem, bei dem die Konvektion dominant ist, ist es fur

das Gauß-Seidel-Verfahren vorteilhaft, bei der Anordnung der Unbe-

kannten (Gitterpunkte) die Stromungsrichtung zu berucksichtigen.


13.3.4. SOR-Verfahren

Der Iterationsschritt beim Gauß-Seidel Verfahren laßt sich auch als

xk+1i = xki − a−1

i,i

( i−1∑

j=1

ai,jxk+1j +

n∑

j=i

ai,jxkj − bi

)

, i = 1,2, . . . , n,

darstellen.

Algorithmus 13.17 (SOR-Verfahren).

Gegeben: Startvektor x0 ∈ Rn, Parameter ω ∈ (0,2).

Fur k = 0,1, . . . berechne:

xk+1i = xki − ωa−1

i,i

( i−1∑

j=1

ai,jxk+1j +

n∑

j=i

ai,jxkj − bi

)

, i = 1,2, . . . , n.

Das SOR-Verfahren hat folgende Matrix-Darstellung:

(D − ωL)xk+1 = Dxk − ω(D − U)xk + ωb


Rechenaufwand. Beim SOR-Verfahren ist der Rechenaufwand

pro Iterationsschritt vergleichbar mit dem Aufwand beim Gauß-

Seidel-Verfahren, also O(n) Operationen fur ein dunnbesetztes

A ∈ Rn×n.

Konvergenz

Satz 13.18. Sei Mω := I− (1ωD−L)−1A die Iterationsmatrix des

SOR-Verfahrens. Es gilt:

• ρ(Mω) ≥ |ω − 1| .• Fur A s.p.d. gilt:

ρ(Mω) < 1 fur alle ω ∈ (0,2).

• Sei A irreduzibel diagonaldominant mit ai,j ≤ 0 fur alle i 6= j

und ai,i > 0 fur alle i. Dann gilt:

ρ(Mω) < 1 fur alle ω ∈ (0,1].


Beispiel 13.19.

Die Diffusionsgleichung wie in den Beispielen 13.10 und 13.15 fur den

Fall h = 1160.

Anzahl der Iterationen in Abhangigkeit von ω:

1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 210

2

103

104

105


Die Konvektions-Diffusionsgleichung wie in den Beispielen 13.11 und

13.16 fur den Fall ε = 10−2, h = 1160.

Anzahl der Iterationen in Abhangigkeit von ω:

0.9 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.910

1

102

103

104


Optimaler Wert von ω

Der optimale Wert von ω ist stark problemabhangig.

Fur die meisten Probleme ist dieser optimale Wert nicht bekannt.

Eine Ausnahme bildet die diskretisierte Poisson-Gleichung (Matrix A1):

Satz 13.20.

Wir betrachten die diskretisierte Poisson-Gleichung A1x = b.

Sei µ := ρ(I−D−1A1) < 1 der Spektralradius der Iterationsmatrix

des Jacobi-Verfahrens.

Sei Mω die Iterationsmatrix des SOR-Verfahrens.

Dann ist ρ(Mω) fur den Relaxationsparameter

ωopt :=2

1 +√

1 − µ2= 1 +

µ

1 +√

1 − µ2

2

optimal und

ρ(Mωopt) = ωopt − 1.


Komplexitat.

Fur die diskretisierte Poisson-Gleichung A1x = b ergibt sich beim

SOR-Verfahren mit dem optimalen ω-Wert eine sehr große Komple-

xitatsverbesserung im Vergleich zum Jacobi- und zum Gauß-Seidel-

Verfahren.

ρ(Mωopt) = ωopt − 1 =

(

cos(πh)

1 + sin(πh)

)2

=1 − sin(πh)

1 + sin(πh)≈ 1 − 2πh.

K = − lnR

ln ρ(Mωopt)≈ 1

2πhlnR ≈ lnR

2π

√m.

Da der Aufwand pro Iteration proportional zu m ist, hat dieses Ver-

fahren eine Komplexitat (fur die vorliegende Problemstellung) von

cm1.5 Operationen.


Beispiel 13.21.

Poisson-Gleichung A1x = b mit x0, b und R wie in Beispiel 13.10.

Wir verwenden das SOR-Verfahren mit ω = ωopt.

Ergebnisse:

h 1/40 1/80 1/160 1/320# 73 146 292 585

K 44 88 176 352


SSOR-Verfahren

Eine SSOR-Iteration sieht also wie folgt aus:

xk+1

2i = xki − ω a−1

i,i

( i−1∑

j=1

ai,jxk+1

2j +

n∑

j=i

ai,jxkj − bi

)

, i = 1,2, . . . , n

xk+1i = x

k+12

i − ω a−1i,i

( n∑

j=i+1

ai,jxk+1j +

i∑

j=1

ai,jxk+1

2j − bi

)

, n ≥ i ≥ 1

Die Iterationsmatrix dieser Methode ist

Mω = I − CωA, Cω := ω(2 − ω)(D − ωU)−1D(D − ωL)−1 .

Wenn die Matrix A symmetrisch positiv definit ist, ist fur ω ∈ (0,2)

die Matrix Cω in der Iterationsmatrix der SSOR-Methode auch sym-

metrisch positiv definit.


13.4 Die Methode der konjugierten Gradienten

Wenn A ∈ Rn×n s.p.d. ist, gelten die folgenden zwei grundlegenden

Eigenschaften:

• Fur x,y ∈ Rn definiert 〈x,y〉A := xTAy ein Skalarprodukt.

• Sei

f(x) := 12xTAx − bTx, b,x ∈ R

n.

Dann gilt, daß f ein eindeutiges Minimum hat und

Ax∗ = b ⇔ f(x∗) = minx∈Rn

f(x).

Lemma 13.22. Sei f wie oben. Die Richtung des steilsten Abstiegs

von f an der Stelle x, d. h. s ∈ Rn so, daß die Richtungsableitung

d

dtf

(

x + ts

‖s‖2

)

∣

∣

∣ t=0= (∇f(x))T

(

s

‖s‖2

)

minimal ist, wird durch s = −∇f(x) = b − Ax gegeben.


Sei U ein Unterraum von Rn. Es gilt:

f(x) = minx∈U

f(x) ⇐⇒ ‖x − x∗‖A = minx∈U

‖x − x∗‖A

Das folgende Lemma zeigt, wie man eine solche Aufgabe lost.

Lemma 13.24. Sei Uk ein k-dimensionaler Teilraum von Rn,

und p0,p1, . . . ,pk−1 eine A-orthogonale Basis dieses Teilraumes:⟨

pi,pj⟩

A= 0 fur i 6= j. Sei v ∈ Rn, dann gilt fur uk ∈ Uk :

‖uk − v‖A = minu∈Uk

‖u − v‖A

genau dann, wenn uk die A-orthogonale Projektion von v auf Uk ist.

Außerdem hat uk die Darstellung

uk =k−1∑

j=0

⟨

v,pj⟩

A⟨

pj,pj⟩

A

pj.


Herleitung der CG-Methode

Der Einfachheit halber: x0 = 0.

Die folgenden Teilschritte definieren die beim CG-Verfahren erzeugten

Naherungen x1,x2, . . . der Losung x∗:

U1 := span{r0},Fur k = 1,2,3, . . . , falls rk−1 = b − Axk−1 6= 0 :

CGa : Bestimme A-orthogonale Basis p0, . . . ,pk−1 von Uk.

CGb : Bestimme xk ∈ Uk, so daß

‖xk − x∗‖A = minx∈Uk

‖x − x∗‖A.

CGc : Erweiterung des Teilraumes:

Uk+1 := span{p0, . . . ,pk−1, rk}, wobei rk := b − Axk.


Erlauterung dieser Teilschritte:

CGa: Die A-orthogonale Basis von Uk wird benotigt, um die Annahe-

rung xk in Teilschritt CGb effizient und stabil zu berechnen.

CGb: xk ist die (bezuglich ‖ · ‖A) optimale Approximation von x∗ in Uk.

CGc: Uk+1 = Uk ⊕ span{rk}:Die Erweiterung des Raumes Uk ist optimal in dem Sinne, daß rk

an der Stelle xk die Richtung des steilsten Abstiegs von f ist


Verfahren der konjugierten Gradienten (CG)

Es ergibt sich der beruhmte

Algorithmus 13.27 (Verfahren der konjugierten Gradienten).

Gegeben: A ∈ Rn×n s.p.d., b ∈ Rn, Startvektor x0 ∈ Rn,

β−1 := 0. Berechne r0 = b − Ax0. Fur k = 1,2, . . . , falls rk−1 6= 0 :

pk−1 = rk−1 + βk−2pk−2, βk−2 =

⟨

rk−1, rk−1⟩

⟨

rk−2, rk−2⟩ (k ≥ 2),

xk = xk−1 + αk−1pk−1, αk−1 =

⟨

rk−1, rk−1⟩

⟨

pk−1,Apk−1⟩,

rk = rk−1 − αk−1Apk−1.

Die CG-Methode ist ein nichtlineares Verfahren:

ek+1 = ψ(ek),

wobei ψ eine nichtlineare Funktion ist.


Rechenaufwand. Im Algorithmus 13.27 werden pro Iterations-

schritt nur eine Matrix-Vektor-Multiplikation, zwei Skalarpro-

dukte und drei Vektor-Operationen der Form x + αy benotigt.

Der Aufwand pro Iterationsschritt ist vergleichbar mit dem des

Jacobi- oder Gauß-Seidel-Verfahrens, also O(n) Operationen fur

ein dunnbesetztes A ∈ Rn×n.

Konvergenz

Satz 13.29. Gegeben sei Ax = b mit A ∈ Rn×n s.p.d. Sei x∗ die

Losung von Ax = b. Dann gilt fur die durch Algorithmus 13.27

gelieferten Naherungen xk

‖xk − x∗‖A ≤ 2

√

κ2(A) − 1√

κ2(A) + 1

k

‖x0 − x∗‖A, k = 0,1,2, . . . .


Beispiel 13.30.

Das Diffusionsproblem A1x = b.

Fur die Konditionszahl der Matrix A1 gilt

κ2(A1) =cos2(1

2πh)

sin2(12πh)

≈(

2

πh

)2,

also

ρ :=

√

κ2(A) − 1√

κ2(A) + 1≈

2πh − 12πh + 1

=1 − 1

2πh

1 + 12πh

≈ 1 − πh.

K := − ln(2R)

ln(1 − πh)≈ 1

πhln(2R)

Ergebnisse:

h 1/40 1/80 1/160 1/320# 65 130 262 525


Superlineare Konvergenz

Fur den Fall h = 180 wird die Fehlerreduktion in der Energie-Norm

ρk :=‖xk − x∗‖A

‖xk−1 − x∗‖A

,

in den ersten 250 Iterationsschritten dargestellt.

0 50 100 150 200 2500.65

0.7

0.75

0.8

0.85

0.9

0.95

1


Komplexitat. Das CG-Verfahren fur die diskrete Poissongleichung:

A1x = b, A1 ∈ Rm×m, m = (n− 1)2, n =

1

h.

Als Fehlerreduktionsfaktor wird R = 103 gewahlt.

Die Komplexitat des CG-Verfahrens (fur diese Problemstellung) ist

cm1.5 Operationen,

also von der gleichen Großenordnung wie beim SOR-Verfahren mit

optimalem ω-Wert.

Ein großer Vorteil der CG-Methode im Vergleich zum SOR-Verfahren

ist, daß sie keinen vom Benutzer zu wahlenden Parameter enthalt.


13.5 Vorkonditionierung

Idee: Die Matrix A wird mit einer geeigneten zu wahlenden symme-

trisch positiv definiten Matrix W vorkonditioniert.

Da die Matrix W symmetrisch positiv definit ist, existiert die Cholesky-

Zerlegung

W = LDLT =: L1LT1 ,

Wir betrachten das transformierte System

Ax = b mit A := L−11 AL−T

1 , x := LT1x, b = L−11 b .

Die Matrix A ist symmetrisch positiv definit.

Der CG-Algorithmus ist also auf dieses Problem anwendbar.


Wenn man den CG-Algorithmus fur das transformierte System in die

ursprunglichen Variablen x = L−T1 x umschreibt, ergibt sich folgender

Algorithmus 13.32 (PCG-Verfahren).

Gegeben: A,W ∈ Rn×n s.p.d., b ∈ Rn, Startvektor x0 ∈ Rn,

β−1 := 0. Berechne r0 = b − Ax0, z0 = W−1r0 (lose Wz0 = r0).

Fur k = 1,2,3, . . . , falls rk−1 6= 0 :

pk−1 = zk−1 + βk−2pk−2, βk−2 =

⟨

zk−1, rk−1⟩

⟨

zk−2, rk−2⟩ (k ≥ 2),

xk = xk−1 + αk−1pk−1, αk−1 =

⟨

zk−1, rk−1⟩

⟨

pk−1,Apk−1⟩,

rk = rk−1 − αk−1Apk−1,

zk = W−1rk (lose Wzk = rk).

Beachte: nur die Matrizen A und W treten auf und nicht A, L1.


Die Konvergenz des PCG-Verfahrens wird durch die Konditionszahl

κ2(A) = κ2(L−11 AL−T

1 ) = κ2(W−1A)

bedingt.

Um eine gute Effizienz des PCG-Verfahrens zu gewahrleisten, ist W

so zu wahlen, daß

• κ2(W−1A) moglichst klein ist,

• die Losung von

Wzk = rk

mit moglichst geringem Aufwand (O(n) Operationen) berechnet

werden kann.

Manchmal ist es schon hilfreich, W = diag(A) zu nehmen.


Unvollstandige Cholesky-Zerlegung als Vorkonditionierung

A laßt sich in der Form

A = L1LT1

mit einer unteren Dreiecksmatrix L1 zerlegen (Cholesky).

Um Dunnbesetztheit zu erhalten greift man auf eine unvollstandige

Cholesky-Zerlegung zuruck. Sei

E := {(i, j) | ai,j 6= 0}Es wird eine naherungsweise Faktorisierung konstruiert:

A = L1LT1 + R,

mit:

• L1 eine untere Dreiecksmatrix,

• (L1)i,j = 0 fur (i, j) /∈ E.

Die Konstruktion einer naherungsweisen Faktorisierung A ≈ L1LT1 be-

ruht auf der unvollstandigen Durchfuhrung einer Methode zur Berech-

nung der vollstandigen Cholesky-Zerlegung.


Algorithmus 13.34 (Unvollstandige Cholesky-Methode).

Gegeben: A ∈ Rn×n s.p.d., Muster E.

Berechne fur k = 1,2, . . . , n :

lk,k =

(

ak,k −∑

j

′l2k,j

)12

;

fur i = k+ 1, . . . , n : falls (i, k) ∈ E,

li,k =

(

ai,k −∑

j

′′li,j lk,j

)/

lk,k ;

Hierbei sind die Summen nur uber Indizes aus dem Muster zu fuhren:

∑

j

′=

k−1∑

j=1,(k,j)∈E,

∑

j

′′=

k−1∑

j=1,(i,j)∈E,(k,j)∈E.

Im PCG-Algorithmus: W = L1LT1 .

Fur viele Probleme: κ2(W−1A) = κ(L−1

1 AL−T1 ) ≪ κ(A).

Die Systeme L1LT1zk = rk sind einfach losbar.


Algorithmus 13.36 (Modifiziertes unvollstandiges Cholesky)

Gegeben: A ∈ Rn×n s.p.d., Muster E. Berechne fur k = 1,2, . . . , n :

lk,k =

(

ak,k −∑

j

′lk,j

)12

;

fur i = k+ 1, . . . n :

falls (i, k) ∈ E :

li,k =

(

ai,k −∑

j

′′li,j lk,j

)/

lk,k ;

sonst:

ai,i = ai,i −∑

j

′′li,j lk,j ;

Man kann zeigen, daß fur L1 := (li,j)1≤i,j≤n:

(L1LT1)i,j = ai,j fur (i, j) ∈ E, i 6= j,

n∑

j=1

(L1LT1)i,j =

n∑

j=1

ai,j fur alle i.


Beispiel 13.37.

Die diskretisierte Poisson-Gleichung.

PCG-Verfahren:

ICCG: CG + unvollstandige Cholesky Vorkonditionierung.

MICCG: CG + mod. unvollstandige Cholesky Vorkonditionierung.

#: Anzahl der Iterationsschritte, die zur Reduktion des Startfehlers

‖x0 − x∗‖2 um einen Faktor R = 103 benotigt werden (x0 := 0).

Ergebnisse:

h 1/40 1/80 1/160 1/320

W = L1LT1 , # 20 40 79 157

W = L1LT1 , # 8 11 14 20


Komplexitat. PCG-Verfahren fur die diskrete Poissongleichung.

Als Fehlerreduktionsfaktor wird R = 103 gewahlt.

Die Konvergenz der ICCG-Methode ist schneller als die der Methode

ohne Vorkonditionierung. Die Komplexitat ist aber fur beide Methoden

cm1.5 Operationen,

wobei die Konstante c beim PCG-Verfahren kleiner ist.

Man kann zeigen, daß die MICCG-Methode die Komplexitat

cm1.25 Operationen

hat.


SSOR-Verfahren als Vorkonditionierung

Sei

xk+1 = xk + C(b − Axk) , A symmetrisch positiv definit.

Man kann die Matrix C als eine Naherung fur A−1 interpretieren.

Deshalb liegt der Ansatz

W = C−1

als Vorkonditionierung fur das CG-Verfahren nahe.

Allerdings muß die Matrix W symmetrisch positiv definit sein.

Beispiele:

- die Jacobi-Methode, mit C−1 = diag(A).

- die SSOR-Methode.

In PCG muß zk = W−1rk also zk = Crk berechnet werden.

zk ist das Resultat einer Iteration des linearen iterativen Verfahrens

angewandt auf Az = rk mit Startvektor 0:

b = rk, x0 = 0 ⇒ x1 = Crk.


Beispiel 13.39.

Die diskretisierte Poisson-Gleichung.

Das PCG-Verfahren mit der SSOR-Methode als Vorkonditionierung.

Als Fehlerreduktionsfaktor wird R = 103 gewahlt; x0 = 0.

Ergebnisse:

h 1/40 1/80 1/160 1/320ω = 2/(1 + sin(πh)) 1.854 1.924 1.961 1.981

# 11 16 22 32


n 1) n 1) Dahmen-Reusken Kapitel 13 1 - igpm.rwth-aachen.de · Fill-in bei einem direkten L¨oser...

Documents

Transcript of n 1) n 1) Dahmen-Reusken Kapitel 13 1 - igpm.rwth-aachen.de · Fill-in bei einem direkten L¨oser...