Direkte Loser und ihre Nachteile III.1 III.2 Klassische ... · L¨osung des linearen...
Transcript of Direkte Loser und ihre Nachteile III.1 III.2 Klassische ... · L¨osung des linearen...
Prof. Dr. Barbara Wohlmuth
Lehrstuhl fur Numerische Mathematik
III. Iterative Loser
III.1 Direkte Loser und ihre Nachteile
III.2 Klassische lineare Iterationsverfahren
– Typeset by FoilTEX – 1
Prof. Dr. Barbara Wohlmuth
Lehrstuhl fur Numerische Mathematik
Erinnerung: Lineares Gleichungssystem bei FDM
Diskretisierung einer linearen PDGL 2. Ordnung mit Finiten Differenzen/FEMfuhrt zu einem linearen Gleichungssystem Ahuh = h2fh, welches
• groß ist (typischerweise: Dimension = O(h−2d)),
• dunn besetzt ist (O(h−d) nicht-Null Eintrage)
• schlecht konditioniert ist (κ(Ah) = O(h−2)),
• haufig symmetrisch, positiv definit ist, und
• bei geeigneter Nummerierung eine Bandstruktur besitzt.
– Typeset by FoilTEX – 2
Prof. Dr. Barbara Wohlmuth
Lehrstuhl fur Numerische Mathematik
Verschiedene Loser: MATLAB
Aus der Matlab Dokumentation: Direkte Loser
– Typeset by FoilTEX – 3
Prof. Dr. Barbara Wohlmuth
Lehrstuhl fur Numerische Mathematik
Verschiedene Loser: MATLAB
Aus der Matlab Dokumentation: Iterative Loser
– Typeset by FoilTEX – 4
Prof. Dr. Barbara Wohlmuth
Lehrstuhl fur Numerische Mathematik
Gauß-Elimination
Zerlege die invertierbare Matrix A ∈ Rn×n in das Produkt zweier Dreiecksmatrizen
L (untere Dreiecksmatrix mit 1 als Diagonaleintrage) und U (obereDreiecksmatrix), so dass
PA = LU,
mit einer Permutationsmatrix P . Damit lasst sich das Gleichungssystem durchSubstitution losen:
Lz = P T b, Ux = z
Aufwand bei vollbesetzten Matrizen:
• Berechnung der LU–Zerlegung: O(n3)
• Ruckwarts/Vorwarts–Substitution: O(n2)
– Typeset by FoilTEX – 5
Prof. Dr. Barbara Wohlmuth
Lehrstuhl fur Numerische Mathematik
Einfluss der Bandstruktur auf Fill-In
• Permutationsmatrizen P erlauben die Vertauschung der Zeilen (beiMultiplikation von links) und der Spalten (bei Multiplikation von rechts) einerMatrix A. Fur diese gelten P−1 = P⊤.
Ax = b ⇔ (PA)x = Pb,
AP⊤Px = b ⇔ (AP⊤)x = b mit x = Px.
• Die Bandstruktur der Matrix A bleibt bei der LU-Zerlegung erhalten, d.h.außerhalb der Bandstruktur treten keine weiteren Eintrage auf, jedoch konnendie Nulleintrage innerhalb der Bandstruktur verschwinden (Fill-In).
• Ziel: Finden von geeigneten Permutationen der Matrix A, so dass moglichstgeringe Bandbreite entsteht. Hierzu gibt es Minimierungsalgorithmen z.B. vonCuthill-McKee.
– Typeset by FoilTEX – 6
Prof. Dr. Barbara Wohlmuth
Lehrstuhl fur Numerische Mathematik
LU–Zerlegung: Matrix-Struktur
Beispiel: 5-Punkte-Stern aus
Finite-Differenzen-Diskretisierung (lexikographische
Nummerierung)
nz = 288
n = 64
nz = 519
nz = 1216
n = 256
nz = 4111
nz = 4992
n = 1024
nz = 32799
A
L
– Typeset by FoilTEX – 7
Prof. Dr. Barbara Wohlmuth
Lehrstuhl fur Numerische Mathematik
LU–Zerlegung: Matrix-Struktur
Beispiel: 5-Punkte-Stern aus
Finite-Differenzen-Diskretisierung (Zufalls-Nummerierung)
nz = 288
n = 64
nz = 622
nz = 1216
n = 256
nz = 5900
nz = 4992
n = 1024
nz = 79321
A
L
– Typeset by FoilTEX – 8
Prof. Dr. Barbara Wohlmuth
Lehrstuhl fur Numerische Mathematik
LU–Zerlegung: Rechenzeit
n0 1000 2000 3000 4000
0
10
20
30
40
ZufallLexikographisch
Rechenzeit stark (asymptotisch) von Nummerierung abhangig!
– Typeset by FoilTEX – 9
Prof. Dr. Barbara Wohlmuth
Lehrstuhl fur Numerische Mathematik
III. Iterative Loser
III.1 Direkte Loser und ihre Nachteile
III.2 Klassische lineare Iterationsverfahren
– Typeset by FoilTEX – 10
Prof. Dr. Barbara Wohlmuth
Lehrstuhl fur Numerische Mathematik
Lineare Iterationsverfahren
Voruberlegung: Eigentlich ist man gar nicht an einer exakten Losung desGleichungssystems interessiert, sondern an einer guten Approximation der PDGL.Sei u∗
h eine Naherung von uh mit Ahuh = fh, so gilt fur den Fehler
‖u− u∗
h‖h≥≤‖u− uh‖h
−+‖uh − u∗
h‖h
d. h. es genugt, wenn ‖uh − u∗
h‖h ∼ ‖u− uh‖h.
Idee: Berechne uh nicht direkt, sondern wiederum uber ein Naherungsverfahren.
– Typeset by FoilTEX – 11
Prof. Dr. Barbara Wohlmuth
Lehrstuhl fur Numerische Mathematik
Lineare Iterationsverfahren: Definitionen
1. Ein Losungsverfahren zur Berechnung von Ax = b heißt iterativ, falls ausgehendvon einem Startwert x0 eine Folge xk von Iterierten bestimmt wird.
2. Ein Iterationsverfahren heißt konvergent, falls unabhangig vom Startwert gilt
limk → ∞
xk = x,
wobei x die exakte Losung ist.
3. Es heißt konsistent, falls aus xk = x folgt, dass xk+1 = x.
4. Es heißt linear, falls xk+1 linear von xk und b abhangt, d.h. es gibt zweiMatrizen M,N ∈ R
n×n mit xk+1 = Mxk +Nb
– Typeset by FoilTEX – 12
Prof. Dr. Barbara Wohlmuth
Lehrstuhl fur Numerische Mathematik
Lineare IterationsverfahrenEin allgemeines konsistentes lineares Iterationsverfahren zur Losung von Ax = b
hat die Form
xk+1 = Mxk +Nb mit M = (Id−NA).
Mit der Zerlegung A = L+D + U erhalten wir folgende Verfahren:
• Jacobi-Verfahren: N = D−1
xk+1 = −D−1(L+ U)xk +D−1b
• vorwartiger Gauß–Seidel: N = (D + L)−1
xk+1 = −(D + L)−1Uxk + (D + L)−1b
• ruckwartiger Gauß–Seidel: N = (D + U)−1
xk+1 = −(D + U)−1Lxk + (D + U)−1b
• SOR-Verfahren fur ω ∈ (0, 2): N = ω(D + ωL)−1
(D + ωL)xk+1 =(
(1− ω)D − ωU)
xk + ωb
– Typeset by FoilTEX – 13
Prof. Dr. Barbara Wohlmuth
Lehrstuhl fur Numerische Mathematik
Definition: Konvergenzrate und Fehlerreduktion
Sei ρ(M) := max{|λj| : λj EW von M} der Spektralradius der Matrix M .
1. ρ(M) heißt Konvergenzrate
2. It(M) := − 1ln(ρ(M)) definiert ein Maß fur die Fehlerreduktion
Bemerkung:
1. Je kleiner ρ(M), desto weniger Iterationen werden benotigt, um den Fehler umeinen vorgegebenen Faktor zu reduzieren.
2. It(M) gibt an, wieviele Iterationen notwendig sind, um den Fehler um denFaktor 1
ezu reduzieren.
Folgerung: ρ(M1) = ρ(M2)2 ⇒ It(M1) =
12It(M2)
– Typeset by FoilTEX – 14
Prof. Dr. Barbara Wohlmuth
Lehrstuhl fur Numerische Mathematik
Historische Bemerkungen
C.F. Gauß in einem Brief vom 26.12.1823 an Gerling:Ich empfehle Ihnen diesen Modus zur Nachahmung. Schwerlich werden Sie jewieder direct eliminieren, wenigstens nicht, wenn Sie mehr als 2 Unbekanntehaben. Das indirecte Verfahren lasst sich halb im Schlafe ausfuhren, oder mankann wahrend desselben an andere Dinge denken.
[C.F. Gauß: Werke Bd. 9, S. 280f, Gottingen 1903]
Block Gauß–Seidel Verfahren: (1819-1822)Supplementum theoriae combinationis observationum erroribus minime obnoxiae
C.G. Jacobi: 1845Uber eine neue Auflosungsart der bei der Methode der kleinsten Quadratevorkommenden linearen Gleichungen
– Typeset by FoilTEX – 15
Prof. Dr. Barbara Wohlmuth
Lehrstuhl fur Numerische Mathematik
Vergleich iterativer Losungsverfahren
Poisson-Matrix
Ah ∈ RN2
h×N2h, Ah = AT
h ,
Ah :=
B −I
−I . . . . . .. . . . . . −I
−I B
mit I,B ∈ RNh×Nh,
B :=
4 −1−1 . . . . . .
. . . . . . −1−1 4
.
– Typeset by FoilTEX – 16
Prof. Dr. Barbara Wohlmuth
Lehrstuhl fur Numerische Mathematik
Iterative Loser: Einfuhrung
Aufwandsabschatzung
• Bandbreite: ω = O(Nh)
• Die Anzahl der Nicht-Nulleintrage wachst mit O(Nh)
• Anwendung von z.B. LU-Zerlegung erfordert Aufwand von O(N2h)
• Matrix-Vektor-Multiplikation hat Aufwand von O(Nh)
– Typeset by FoilTEX – 17
Prof. Dr. Barbara Wohlmuth
Lehrstuhl fur Numerische Mathematik
Iterative Loser: Einfuhrung
Poisson-Matrix
0 200 400 600 800 10000
0.2
0.4
0.6Zeit, bis Fehler < 10
−8
n2
Ze
it
LU (voll)GSCG
Gut geeignet um große, gut konditionierte, dunnbesetzte Matrizen zu losen
Das Verfahren der konjugierten Gradienten (CG) ist ein nichtlinearesKrylovraum-Verfahren, welches fur spd Matrizen besonders schnell konvergiert.
– Typeset by FoilTEX – 18
Prof. Dr. Barbara Wohlmuth
Lehrstuhl fur Numerische Mathematik
SOR fur verschiedene Dampfungsparameter
Poisson-Matrix
0 50 100 150 200
10−8
10−4
100
Anzahl der Iterationen
Fe
hle
rn
orm
w = 0.2
w = 0.6
w = 1
w = 1.4
w = 1.8
Figure 1: Dampfungsparameter imBereich [0.2, 1.8]
0 10 20 30
10−8
10−4
100
Anzahl der Iterationen
Feh
lern
orm
w = 0.8
w = 0.9
w = 1
w = 1.1
w = 1.2
Figure 2: Dampfungsparameter imBereich [0.8, 1.2]
– Typeset by FoilTEX – 19
Prof. Dr. Barbara Wohlmuth
Lehrstuhl fur Numerische Mathematik
Konvergenz fur SOR-Verfahren
Fur die Matrix Mω des SOR-Verfahrens gilt
Mω = (D + ωL)−1(
(1− ω)D − ωU)
.
Konvergenz liegt vor, falls fur den Spektralradius ρ(Mω) < 1 ⇔ 0 < ω < 2gilt. Der Spektralradius ρ(Mω) nimmt sein Minimum fur den optimalenDampfungsparameter
ωopt =2
1 +√
1− ρ2J
an, wobei ρJ den Spektralradius der Iterationsmatrix MJ = −D−1(L+ U) desJacobi-Verfahrens bezeichnet. Dann gilt fur die Konvergenzrate
ρ(Mωopt) =1−
√
1− ρ2J
1 +√
1− ρ2J.
Allgemein gilt
ρ(Mω) =
{
ω − 1 fur ωopt ≤ ω ≤ 2,
1− ω + 12ω
2ρ2J + ωρJ
√
1− ω + 14ω
2ρ2J fur 0 ≤ ω ≤ ωopt.
– Typeset by FoilTEX – 20
Prof. Dr. Barbara Wohlmuth
Lehrstuhl fur Numerische Mathematik
Asymptotische Konvergenzrate fur SOR-Verfahren
0 0.5 1 1.5 20
0.2
0.4
0.6
0.8
1
Dämpfungsparameter ω
Konverg
enzra
te ρ
ρJ = 0.3
ρJ = 0.5
ρJ = 0.7
ρJ = 0.9
– Typeset by FoilTEX – 21
Prof. Dr. Barbara Wohlmuth
Lehrstuhl fur Numerische Mathematik
Jacobi, Gauß-Seidel und optimales SOR im Vergleich
Poisson-Matrix
0 20 40 600
1000
2000
3000
4000
5000
6000
7000
Number of unknowns n2
Num
ber
of itera
tions
opt. SOR Gauss−SeidelJacobi
0 10 20 30 40 50 600.4
0.5
0.6
0.7
0.8
0.9
1
Anzahl d. Iterationen η
Faustregel:
ηJ ≈ 2ηGS
Konvergenzraten ρ
Faustregel:
ρGS ≈ ρ2J
ρJ ≈ 1− cJh2
ρSOR ≈ 1− cSORh
– Typeset by FoilTEX – 22
Prof. Dr. Barbara Wohlmuth
Lehrstuhl fur Numerische Mathematik
Vergleich iterativer Losungsverfahren
Poisson-Matrix, Fehler gegen Anzahl der Iterationen,
n = p = 100
Losung des linearen Gleichungsystems Ahuh = fh mit Jacobi, Gauß-Seidel,symmetrischer Gauß-Seidel Methode und mit dem Verfahren der konjugierten
Gradienten (CG)
0 20 40 60 8010
−8
10−4
100
JacobiGSSGSCG
– Typeset by FoilTEX – 23