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).
Dahmen-Reusken Kapitel 13 2
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
Dahmen-Reusken Kapitel 13 3
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.
Dahmen-Reusken Kapitel 13 4
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“.
Dahmen-Reusken Kapitel 13 5
Nichtnullelemente von L bei der Cholesky-Zerlegung der Matrix A1:
0 50 100 150 200
0
50
100
150
200
nz = 3389
Dahmen-Reusken Kapitel 13 6
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.
Dahmen-Reusken Kapitel 13 7
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, . . . .
Dahmen-Reusken Kapitel 13 8
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.
Dahmen-Reusken Kapitel 13 9
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.
Dahmen-Reusken Kapitel 13 10
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.
Dahmen-Reusken Kapitel 13 11
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.
Dahmen-Reusken Kapitel 13 12
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.
Dahmen-Reusken Kapitel 13 13
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.
Dahmen-Reusken Kapitel 13 14
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.
Dahmen-Reusken Kapitel 13 15
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.
Dahmen-Reusken Kapitel 13 16
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.
Dahmen-Reusken Kapitel 13 17
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.
Dahmen-Reusken Kapitel 13 18
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.
Dahmen-Reusken Kapitel 13 19
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
Dahmen-Reusken Kapitel 13 20
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.
Dahmen-Reusken Kapitel 13 21
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.
Dahmen-Reusken Kapitel 13 22
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
Dahmen-Reusken Kapitel 13 23
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].
Dahmen-Reusken Kapitel 13 24
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
Dahmen-Reusken Kapitel 13 25
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
Dahmen-Reusken Kapitel 13 26
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.
Dahmen-Reusken Kapitel 13 27
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.
Dahmen-Reusken Kapitel 13 28
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
Dahmen-Reusken Kapitel 13 29
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.
Dahmen-Reusken Kapitel 13 30
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.
Dahmen-Reusken Kapitel 13 31
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.
Dahmen-Reusken Kapitel 13 32
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.
Dahmen-Reusken Kapitel 13 33
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
Dahmen-Reusken Kapitel 13 34
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.
Dahmen-Reusken Kapitel 13 35
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, . . . .
Dahmen-Reusken Kapitel 13 36
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
Dahmen-Reusken Kapitel 13 37
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
Dahmen-Reusken Kapitel 13 38
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.
Dahmen-Reusken Kapitel 13 39
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.
Dahmen-Reusken Kapitel 13 40
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.
Dahmen-Reusken Kapitel 13 41
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.
Dahmen-Reusken Kapitel 13 42
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.
Dahmen-Reusken Kapitel 13 43
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.
Dahmen-Reusken Kapitel 13 44
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.
Dahmen-Reusken Kapitel 13 45
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
Dahmen-Reusken Kapitel 13 46
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.
Dahmen-Reusken Kapitel 13 47
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.
Dahmen-Reusken Kapitel 13 48
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
Dahmen-Reusken Kapitel 13 49
Top Related