7.1 Motivation
Lineare Gleichungssysteme sehen auf den ersten Blick trivial einfach aus, ihre L¨osung ist aber Kernst¨uck sehr vieler Numerikan-wendungen, weil sowohl Diskretisierungen (Gitter, finite Elemente, usw.) als auch Basisfunktionsdarstellungen vieler Differenti-algleichungen letztlich lineare Gleichungssysteme liefern. Das Paradebeispiel ist numerische Str¨omungsmechanik (computational fluid dynamics, CFD) mit unz¨ahligen Anwendungsgebieten in den Natur- und Ingenieurwissenschaften (Umstr¨omung von Ro-torbl¨attern, Tragfl¨ugeln, Schiffsr¨umpfen, usw.; Wetter-, Klima- und Ozeansimulationen; Verbrennungsvorg¨ange in Automotoren, Flugzeug- und Raketentriebwerken, usw.). Gelegentlich wird das auch in der Chemie bzw. Molek¨ulphysik wichtig:
Chemische Reaktionen A + B→C + D k¨onnen als “Streuprozeß” beschrieben werden: In der Reaktanden-Asymptote (Entfernung A-B → ∞ und Zeit t → −∞) sind die Zust¨ande des Systems
”A und B“ durch einen Hamiltonoperator ˆH0 beschreibbar. Bei vorw¨artslaufender Zeit ttreffen Teilchen A und B zur Zeit t = 0 aufeinander und reagieren mit einer gewissen Wahrscheinlichkeit zu C und D. Dabei ist der reagierende Komplex durch einen Hamiltonoperator ˆH gegeben. Bei t → ∞ sind die Produkte C und D wieder
”unendlich“ weit voneinander entfernt. Im Experiment kann man typischerweise in der Asymptote t → −∞ die Reaktanden mit einer gewissen Verteilung ¨uber wenige Quantenzust¨ande kontrolliert pr¨aparieren und dann in der Asymptote t→ ∞ die Verteilung ¨uber die Produkt-Quantenzust¨ande detektieren, hat aber keinen direkten Zugang zum Reaktionsgeschehen um t= 0.
In der Streutheorie betrachten wir daher die zu l¨osende Schr¨odingergleichung
HˆΨ = EΨ (137)
f¨ur die Gesamtsituation, sowie eine bereits gel¨oste Schr¨odingergleichung
Hˆ0Ψ0 = EΨ0 (138)
f¨ur die Eduktasymptote t → −∞. Dabei sei
Hˆ = ˆH0 +V , Ψ = Ψ0 +χ (139)
Einsetzen von Gl. 139 in Gl. 137 liefert
( ˆH0 +V)(Ψ0 +χ) = E(Ψ0 +χ) (140)
Unter Verwendung von Gl. 138 wird daraus:
( ˆH −Eˆ1)χ = −VΨ0 (141)
F¨ur eine gegebene Streu-Gesamtenergie E ist der Operator auf der linken Seite bekannt; aufgrund der L¨osung von Gl. 138 ist die rechte Seite auch bekannt. Gesucht ist der Wellenfunktionsanteil χ, der den eigentlichen Streuvorgang im Wechselwirkungsgebiet wiedergibt und diesen Bereich mit den Edukt- und Produktasymptoten verbindet. Aus dessen Projektion auf die Produktzust¨ande kann man die Reaktionswahrscheinlichkeiten einzelner Eduktzust¨ande in einzelne Produktzust¨ande berechnen, und daraus wieder durch geeignete Summationen/Integrationen die Reaktionsgeschwindigkeitskonstante k(T).
Gl. 141 kann abgek¨urzt geschrieben werden als
Aχˆ = −VΨ0 (142)
Verwenden wir einen geeigneten Satz von Basisfunktionen {φi}, k¨onnen wir sowohl Ψ0 als auch χ in diese Basis entwickeln χ = X
i
xiφi , Ψ0 = X
i
aiφi (143)
und diese Entwicklungen in Gl. 142 einsetzen:
X
i
xiAφˆ i = −X
i
aiV φi (144)
Multiplikation dieser Gleichung von links mit hφj| liefert X
i
xihφj|A|φˆ ii = −X
i
aihφj|V|φii (145)
Alle Integrale k¨onnen berechnet werden und liefern jeweils eine Zahl f¨ur jede Indexkombination i und j. Auf der rechten Seite die Entwicklungskoeffizienten {ai} bekannt, also kann dort die Summe ¨uber i ausgef¨uhrt werden und wir erhalten:
X
i
Ajixi = bj (146)
Dies ist die j-te Zeile der folgenden Matrix-Vektor-Gleichung
A~x =~b (147)
Dabei ist die Matrix A auf der linken Seite und der Vektor~b auf der rechten Seite bekannt; gesucht wird der Vektor ~x. Es handelt sich um ein lineares Gleichungssystem.
7.2 Allgemeines
Ein lineares Gleichungssystem besteht aus M Gleichungen mit linearen Termen aus N Unbekannten und N ×M Koeffizienten.
Einzelne Terme ohne Unbekannte kann man auf der rechten Seite zusammenfassen:
a11x1 +a12x2 +· · ·+a1NxN = b1 (148)
• mehrere rechte Seiten~bi bei gleicher Matrix A:
A~x1 =~b1 , A~x2 =~b2 , . . . AX = B (154)
Z.B. ein Streuproblem mit mehreren Anfangsbedingungen.
• Wenn diese~bi die Spalten der Einheitsmatrix1 sind, dann sind wegen AA−1 = 1 die zu jedem igeh¨orenden L¨osungsvektoren
~
xi die Spalten von A−1 ⇒ Matrixinversion
Formal w¨urde Multiplikation von Gl. 153 von links mit A−1 direkt die L¨osung liefern:
A−1A~x = A−1~b ⇒ ~x = A−1~b (155)
Dies wird in der numerischen Praxis nie so gemacht, weil die Matrixinversion nach obigem Schema so aufwendig ist wie eine N-malige L¨osung des Gleichungssystems Gl. 153 und nach anderen Schemata die numerisch ung¨unstige Berechnung von Deter-minanten ben¨otigt.
Tats¨achlich kann man in der Praxis eine scheinbar n¨otige Matrixinversion ~x = A−1~b h¨aufig durch eine L¨osung des entsprechenden linearen Gleichungssystems A~x =~b ersetzen. Alternativ kann man die obige Idee verwenden (Gl. 154 mit B = 1).
7.3 Kleines Matrix-Nomenklatur-Lexikon: AMoN quadratisch: M = N
rechteckig: M 6= N
diagonal, tridiagonal, band-diagonal, block-diagonal
schwach besetzt, d¨unn besetzt (sparse):
Dreiecksmatrix: Hier abgek¨urzt: Obere/untere 4-Matrix
symmetrisch (hermitesch/selbstadjungiert): A= AT (A = A† = (AT)∗)
• Alle Eigenwerte reell
• Eigenvektoren bilden vollst¨andiges Orthonormalsystem positiv definit:
• Symmetrische (hermitesche) Matrizen mit positiven Eigenwerten
• Haben immer ein Inverses positiv semi-definit:
• Neben positiven Eigenwerten auch der Eigenwert 0 erlaubt
orthogonal (unit¨ar): A† = A−1
• Die Spaltenvektoren bilden ein Orthonormalsystem
• Bei A~v bleibt die Norm von ~v erhalten: kA~vk = k~vk normal: AA† = A†A
• Die Eigenvektoren bilden ein vollst¨andiges Orthonormalsystem
• unit¨ar diagonalisierbar: A = UDU†, mit D = diag{λi}
• Bei reellen Matrizen sind alle symmetrischen und alle orthogonalen Matrizen normal singul¨ar:
• Linear abh¨angige Zeilen (oder Spalten), det(A) = 0 regul¨ar: det(A) 6= 0, invertierbar.
7.4 L¨osbarkeit
AMoN~xN =~bM. (156)
Inhomogene Systeme:~b 6=~0
• M = N
eindeutige L¨osung, wenn A nicht singul¨ar: det(A) 6= 0.
Wenn det(A) = 0, dann – Keine L¨osung, oder
– Ein- oder mehrfach unendliche L¨osungsmannigfaltigkeit, d.h. eine oder mehrere Unbekannte xi bleiben unbestimmt.
• M < N: Unterbestimmt
(effektiv auch der Fall, wenn det(A) = 0)
Singularit¨aten von A analysierbar und L¨osungsmannigfaltigkeit bestimmbar mit Singul¨arwertzerlegung (singular value de-composition, SVD, s.u.)
• M > N: ¨Uberbestimmt
Und letztlich widerspr¨uchlich! (wenn nicht, dann r¨uckf¨uhrbar auf M ≤ N)
M¨ogliche “Kompromiß-L¨osung” durch lineare Regression finden (Minimierung der Fehler mit least squares); ggf. auch wieder mit SVD machbar.
7.5 Konditionierung
(Herleitungen z.T. in Henrik Larssons Numerikskript WS15/16)
• Der relative Fehler in ~x ist bei der L¨osung von A~x =~b proportional zur Kondition der Matrix A.
cond(A) ≡ kA−1kkAk (157)
Vorsicht: Durch die n¨otige Matrixinversion ist die Konditionsberechnung aufwendiger als die L¨osung des linearen Gleichungs-systems mit nachfolgender Berechnung des sog. Residuums ~r =~b−A~x
• Mit der euklidischen (2er-)Norm gilt:
cond(A)2 = kA−1k2kAk2 = s
λmax λmin
. (158)
Je gr¨oßer also das Verh¨altnis zwischen gr¨oßtem und kleinstem Eigenwert λ, desto gr¨oßer auch die Kondition.
• F¨ur eine gute Kondition auf dem Computer muss mit der Maschinengenauigkeit mach gelten:
cond(A)mach 1 (159)
Mit t-stelligen real-Zahlen und einer Kondition von cond(A) ≈ 10α kann ein Gleichungssystem daher nur mit einer Genau-igkeit von t−α −1 Dezimalstellen gel¨ost werden (bezogen auf die betragsgr¨oßte Komponente).
7.6 Gaußsches Eliminationsverfahren
• Wichtig: Der L¨osungsvektor ~x bleibt unver¨andert, wenn beliebige Zeilen des Gleichungssystems vertauscht werden, oder wenn von einer Zeile das Vielfache einer anderen Zeile subtrahiert wird.
• Idee: Bringe das System auf 4-Form
A11x1 +A12x2 +· · ·+A1NxN = b1 (160) 0 + ˜A22x2 +· · ·+ ˜A2NxN = ˜b2 (161)
...
0 + ˜AM NxN = ˜bM (162)
• Beispiel: Um Ai1x1 in Zeile 2 bis M zu eliminieren, muss von allen Eintr¨agen ab Spalte 2 AAi1
11A1j subtrahiert werden. A11 nennt man dann Pivotelement und AAi1
11 Multiplikator.
• Wichtig f¨ur numerische Stabilit¨at: Zeilen werden so vertauscht, dass immer der Eintrag mit dem gr¨oßten absoluten Wert Pivotelement wird. Dies nennt man teilweise Pivotierung (Spaltenpivotierung).
• Algorithmus f¨ur ANoN~xN =~bN:
p ← 1, . . . , N B Speicherung der Permutationen
f¨ur Spalten k = 1, . . . , N −1 :
imax ←max(|Aik|, i= k, . . . , N) wenn |Aimaxk| < dann
Matrix ist Singul¨ar; Abbruch.
wenn ende
Vertausche kte Zeile mit imaxter Zeile B Spaltenpivotierung
Vertausche pk und pimax f¨ur Zeilen i = k + 1, . . . , N :
d ← AAik
kk
f¨ur Spalten j = k+ 1, . . . , N : Aij ← Aij −Akj ×d f¨ur ende
Aik ←0 B Unteres 4 wird Null
bi ←bi −bk ×d f¨ur ende
f¨ur ende
• Vollst¨andige Pivotierung (Zeilen- und Spaltentausch) auch m¨oglich, aber aufw¨andiger und in der Praxis oft nicht n¨otig.
• Beachte: A und~b werden ¨uberschrieben.
• Anschließend R¨uckw¨artssubstitution zur Berechnung von ~x:
f¨ur i = N, . . . ,1 : s ← bi
f¨ur j = i+ 1, . . . , N : s ←s−Aijxj f¨ur ende
xi ← As
ii
f¨ur ende
• Zuletzt Beachtung der Permutation: ~x ←~x(~p).
• Algorithmus skaliert mit O(N3), R¨uckw¨artssubstitution nur mit O(N2).
(Deswegen ist die Erweiterung auf den sog. Gauß-Jordan-Algorithmus, bei dem man A mit ganz analogen Operationen nicht auf Dreiecks-, sondern auf Diagonalform bringt und sich dadurch die R¨uckw¨artssubstitution spart, zwar didaktisch interessant, aber f¨ur die Praxis irrelevant.)
7.6.1 Beispiel: Notwendigkeit von Pivotierung
• Gegeben sei
10−4x1 +x2 = 1, (163)
x1 +x2 = 2, (164)
⇒ x1 = 104
9999, x2 = 9998
9999, cond(A) = 2.6 (165)
• Mit einer Mantissenl¨ange von 3 Dezimalstellen ergibt sich:
x1 ≈ 0.100×101, (166)
x2 ≈ 0.100×101. (167)
• Ohne Pivotierung berechnet sich das Gleichungssystem nach
0.100×10−3x1 + 0.100×101x2 = 0.100×101 1.Zeile
0.100×101x1 + 0.100×101x2 = 0.200×101 2.Zeile (168) MultiplikatorA21
A11
= 0.100×101
0.100×10−3 = 0.100×105 (169)
Anderung der zweiten Zeile:¨
y(0.100×101 −0.100×105 ×0.100×10−3)x1 + (0.100×101 −0.100×105 ×0.100×101)x2
= 0.200×101 −0.100×105 ×0.100×101 (170)
Ver¨andertes Gleichungssystem (mit Rundungsfehlern in A22 und b2):
⇒ 0.100×10−3x1 + 0.100×101x2 = 0.100×101 1.Zeile
−0.100×105x2 = −0.100×105 2.Zeile (171)
L¨osung nach R¨uckw¨artssubstitution:
x2 = 0.100×101 (172)
x1 = 0.000. (173)
• x1 ist offensichtlich falsch! F¨ur eine Mantissenl¨ange von drei ist das L¨osen schlecht konditioniert. (Pivotierung ist aber fast immer n¨otig, auch bei sehr viel gr¨oßeren Mantissenl¨angen, da A durchaus auch Nullen enthalten kann, ohne dass die Matrix singul¨ar ist.)
• Fehlerursache: In Gl. 169 wird aus einem kleinen Pivotelement ein großer Multiplikator, der bei der Berechnung der modi-fizierten 2.Zeile (von Gl. 170 zu Gl. 171) zu Rundungsfehlern/Weghebeph¨anomenen f¨uhrt (bei der Differenzbildung in A22 und b2 sind 105 und 101 wegen der nur dreistelligen Mantisse nicht in derselben Zehnerpotenz darstellbar).
• Mit Pivotierung (Vertauschung der Reihenfolge der Gleichungen):
0.100×101x1 + 0.100×101x2 = 0.200×101 1.Zeile
0.100×10−3x1 + 0.100×101x2 = 0.100×101 2.Zeile (174) Multiplikator 0.100×10−3
0.100×101 = 0.100×10−3 (175)
⇒ 0.100×101x1 + 0.100×101x2 = 0.200×101 1.Zeile
0.100×101x2 = 0.100×101 2.Zeile (176)
⇒ x2 = 0.100×101 (177)
x1 = 0.100×101 (178)
• Die Pivotierung ist also essentiell. Damit ergeben sich kleinere Multiplikatoren und damit im weiteren Verlauf Zahlen mit einer ¨ahnlichen Gr¨oßenordnung ⇒ kein/vermindertes Weghebeph¨anomen.
7.7 LU-Zerlegung
• Problem beim Gauß-Verfahren: Es muß bei gleicher MatrixA aber neuen rechten Seiten~bi erneut angewandt werden. Oftmals l¨ost man aber mehrere Gleichungssysteme mit derselben Koeffizientenmatrix A, f¨ur zahlreiche rechte Seiten~bi.
• Daher besser: Finde ein Verfahren, welches A so variiert, dass Gleichungssysteme einfacher gel¨ost werden k¨onnen.
• Ein Gauß-Eliminationsschritt s mit Multiplikatoren lik = AAik
kk ergibt sich letztlich durch Multiplikation mit der sogenannten Frobenius-Matrix
Das Inverse von Ls ergibt sich durch Vorzeichenwechsel der lik.
• Wiederholte Anwendung f¨uhrt zur oberen 4-Matrix, was dem Endprodukt der Gauß-Elimination entspricht:
LN−1. . .L1A =
• A muß entsprechend erhalten werden, wenn man in umgekehrter Reihenfolge L−1s auf Uanwendet. Das Produkt der Inversen von Ls f¨uhrt zu einer unteren Dreiecksmatrix mit Einsen auf der Diagonalen:
L−11 . . .L−1N−1 =
• Dies ist die LU-Zerlegung (lower/upper triangular matrix; im Deutschen: LR-Zerlegung).
• L¨osung von Gleichungssystemen mittels Vorw¨arts- und R¨uckw¨artssubstitution:
A~x = b, (183)
⇔L U~x
|{z}
~z
=~b, (184)
L~z =~b, L¨osung per Vorw¨artssubstitution, (185) U~x = ~z, L¨osung per R¨uckw¨artssubstitution. (186)
• Wenn die Diagonale von L nicht abgespeichert wird, kann sowohl L als auch U in A gespeichert werden (in place).
• Mit teilweiser Pivotierung l¨asst sich schreiben:
PA = LU, (187)
A = P−1LU = ˜LU. (188)
P ist die Permutationsmatrix (als Vektor speicherbar!). ˜L ist keine 4-Matrix.
• Crout-Algorithmus mit teilweiser Pivotierung und in-place-Speicherung:
p ← 1, . . . , N B Speicherung der Permutationen
f¨ur Spalten k = 1, . . . , N −1 :
imax ←max(|Aik|, i= k, . . . , N) wenn |Aimaxk| < dann
Matrix ist Singul¨ar; Abbruch.
wenn ende
Vertausche kte Spalte mit imaxter Spalte Vertausche pk und pimax
f¨ur Zeilen i = k + 1, . . . , N : Aik ← AAik
kk B Uberschreibe¨ Aik mit lik
f¨ur Spalten j = k+ 1, . . . , N : Aij ← Aij −Akj ×Aik f¨ur ende
f¨ur ende f¨ur ende
• Wie Gauß-Verfahren, nur daß die Multiplikatoren im unteren 4 gespeichert werden. (Und wie beim Gauß-Verfahren ist dies nur der (hier: LU-)Zerlegungsteil; dazu kommen noch die Substitutionsschritte Gl. 185 und 186.)
7.8 Cholesky-Zerlegung
• F¨ur symmetrische, positiv definite Matrizen A gilt 1. ∀~x ∈ RN : ~xTA~x > 0,
2. Alle Eigenwerte sind positiv,
3. Alle Diagonalelemente sind positiv,
4. maxij |Aij| = maxi|Aii| (nur Diagonalpivotierung n¨otig!),
5. Bei jeder Gauß-Elimination ohne Pivotsuche ist die Restmatrix wiederum positiv definit.
• Letzte Eigenschaft l¨asst sich wie folgt zeigen:
A =
Multiplikation von A(1) von rechts mit LT1 eliminiert ~zT: L1ALT1 =
A11 ~0T
~0 B(1)
. (191)
Eine ¨Ahnlichkeitstransformation L1ALT1 ¨andert nicht die Positiv-Definitheit. Also muss B(1) auch positiv definit sein.
• Fortsetzung der Multiplikation mit Ls wie in Gl. 191 f¨uhrt zur LDL-Zerlegung (die auch f¨ur nicht positiv-definite, aber symmetrische Matrizen funktioniert):
• Wurzelziehung der Diagonalmatrix f¨uhrt zur Cholesky-Zerlegung
A = LD12D12LT = CCT. (194)
(Dies ist wegen der Wurzelziehung nur f¨ur positiv definite Matrizen m¨oglich (D ist positiv definit)).
• Der Algorithmus ergibt sich ¨uber:
• Also lautet der Algorithmus:
f¨ur Spalten k= 1, . . . , N :
wenn Akk < dann B Beachte: Kein Betrag
Matrix ist nicht positiv-definit; Abbruch.
wenn ende Ckk =√
Akk
f¨ur Zeilen i=k+ 1, . . . , N : Cik = CAik
f¨ur Spaltenkk j =k+ 1, . . . , i : B Beachte: Schleife geht nur bisi
Aij ←Aij −CikCjk f¨ur ende
f¨ur ende f¨ur ende
• Da die letzte Schleife nur bis i geht und keine Pivotierung ben¨otigt wird, ist der Algorithmus mindestens doppelt so schnell wie LU-Zerlegung.
• Da A symmetrisch und C eine 4-Matrix ist, gibt es kompakte Speicherformen, die die Matrix jeweils als Vektor der L¨ange N(N + 1)/2 speichert.
7.8.1 Anwendungsbeispiele
Die Cholesky-Zerlegung ist von sehr großer Bedeutung und wird f¨ur viele Probleme angewandt. Zwei sehr h¨aufige werden kurz vorgestellt.
Verallgemeinerte Eigenwertprobleme
• Neben dem typischen Eigenwertproblem (L¨osung per Diagonalisierung von H)
HX = diag(~e)X, H† = H (198)
kommt sehr h¨aufig das verallgemeinerte Eigenwertproblem
HX = diag(~e)SX, (199)
mit einer symmetrischen, positiv definiten Matrix S vor.
• folgende L¨osungsidee nicht sinnvoll, da S−1H nicht mehr hermitesch ist:
[S−1H]X = diag(~e)X (200)
• Besser: per Cholesky-Zerlegung von S in ein normales Eigenwertproblem ¨uberf¨uhren:
S = CC†,(Cholesky) (201)
y HX = diag(~e)CC†X, (202)
⇒ C−1HX = diag(~e)C†X, (203)
⇔C−1H([C†]−1
| {z }
H˜
C†)X = diag(~e)C†X, (204)
y H˜X˜ = diag(~e) ˜X, (205)
H˜ = C−1H[C−1]†, (206)
X˜ = C†X. (207)
• H˜ ist eine ¨Ahnlichkeitstransformation von H und damit erhalten sich sowohl die Hermitizit¨at als auch die Eigenwerte.
• Am Ende: R¨ucktransformation per X = [C†]−1X.˜
• Cholesky-Zerlegung und Invertierung von C problemlos m¨oglich, sofern S numerisch positiv definit ist (keine zu große Kondition). (Die Invertierung kann aber auch durch L¨osung von CA = B umgangen werden.)
Orthogonalisierung
• Finde mit positiv definiter Matrix S eine Matrix Y, sodaß
X†X = S, (208)
Xorth = XY, (209)
⇒X†orthXorth = 1, (210)
• L¨osung wieder per Cholesky-Zerlegung:
S = C†C,(Cholesky) (211)
y X†X = C†C, (212)
⇒[C†]−1X†X = C, (213)
⇒ [C†]−1X†XC−1 = 1, (214)
⇔ [XC−1]†XC−1 = 1, (215)
⇒Xorth = XC−1. (216)
• Diese (
”kanonische“) Orthogonalisierung ist identisch zur Gram-Schmidt-Orthogonalisierung.
7.9 Nachiteration
• Das Erreichen von hoher Genauigkeit ist f¨ur lineare Gleichungssysteme mit schlecht konditionierten Matrizen schwierig. Eine Nachbesserung der erhaltenen L¨osung ist aber durch Nachiteration (iterative refinement) sehr einfach m¨oglich, wennA oder
~b nicht schon mit Fehlern behaftet sind und die Kondition von A nicht zu groß ist.
• Gegeben sei A~x =~b mit der mit Rundungsfehlern behafteten L¨osung ~x(0).
• Der Residuenvektor/das Residuum ist definiert durch
~r(0) =~b−A~x(0) = A(~x−~x(0)) =A∆~x(0). (217)
• ∆~x(0) ist die Korrektur zu ~x(0):
~x(1) = ~x(0) + ∆~x(0), (218)
mit ~x(1) als verbesserte L¨osung.
• ∆~x(0) kann durch L¨osen des Gleichungssystems in Gl. 217 berechnet werden:
A∆~x(0) = ~r(0). (219)
• Dieses Verfahren kann iterativ bis zur gew¨unschten Genauigkeitsverbesserung durchgef¨uhrt werden.
• Da die Faktorisierung von A schon gegeben ist, hat dieses Verfahren nur einen Aufwand von O(N2) und ist damit im Verh¨altnis zum Gesamtaufwand der Faktorisierung deutlich geringer.
7.10 Iterative Verfahren
• In vielen Anwendungen (z.B. L¨osung von partiellen Differentialgleichungen) sind die erhaltenen Matrizen sehr groß (N ∼ 106– 109), aber schwach besetzt/sparse.
• Die gesamte Matrix kann dann nicht abgespeichert werden, sehr wohl aber alle Elemente, die nichtverschwindend sind.
• Klassische Algorithmen sind f¨ur solche Matrizen nicht zu gebrauchen, da sich dort das Muster der Besetztheit ¨andert und die Matrix andauernd ver¨andert wird.
• Besser: Iterative Verfahren, welche nur Matrix-Vektor-Multiplikationen verwenden. Um ein Matrix-Vektor-Produkt A~b =~c effizient durchf¨uhren zu k¨onnen, reicht es zur Berechnung jedes einzelnen Elements ci, nur die eine, jeweils zugeh¨orige Zeile von A im Speicher zu halten.
• Es gibt hierf¨ur eine F¨ulle von Algorithmen (Gauß-Seidel, SOR, etc.). Hier nur Besprechung von conjugate gradient.
7.10.1 Conjugate Gradient
• Siehe Abschnitt 6.3.3 zur Minimierung.
Dort wurde f(~x) =~cT~x+ 12~xTH~x minimiert, was effektiv die L¨osung des Gleichungssystems H~x = −~c ist.
• vgl. Gl. 129 mit Gl. 217: statt den Gradienten ∇f(~x) zu minimieren, bis er Null wird, minimieren wir hier das Residuum
~
r =~b −A~x, mit einem weitestgehend analogen Algorithmus, der wieder iterativ ist.
• Pro Iteration braucht es nur eine Matrix-Vektor-Multiplikation.
• Theoretisch braucht der Algorithmus N Schritte f¨ur eine N oN-Matrix.
• Allerdings ist die numerische Konvergenz schon nach wenigen ∼ 10−100 Schritten erreicht. Dies h¨angt aber sehr von der Kondition von A ab!
• Erweiterung auf nicht-symmetrisch-positiv-definite Matrizen ergibt das BiCG-Verfahren (biconjugate gradient). Weitere Variationen sind vorhanden.
7.10.2 Pr¨akonditionierung
• Iterative Verfahren konvergieren nur schnell, wenn die Kondition der Matrix klein ist.
• Verringerung der Konvergenz durch einen Vorkonditionierer M−1, der n¨aherungsweise dem Inversen von A entspricht (oder eine entsprechende Funktion, welche Gleichungssysteme f¨ur M l¨ost).
• Anstelle von A~x =~b wird also gel¨ost:
M−1A~x = M−1~b, (220)
⇔M−1(A~x−~b) = 0 (221)
• Im obigen Algorithmus wird M−1 auf die Residuen ~r angewandt.
• Da Rundungsfehler auftreten k¨onnen, ist die Polak-Ribi`ere-Korrektur wieder zu empfehlen (Gl. 132).
• Die Wahl geeigneter Pr¨akonditionierer ist schwierig: Das Gleichungssystem M~y =~b soll schnell gel¨ost werden und M nicht viel Arbeitsspeicher belegen, aber dennoch die Kondition sehr stark verringern.
• M¨ogliche Wahlen
– Jacobi-Konditionierung: M ist eine Diagonalmatrix mit den Diagonaleintr¨agen von A (sinnvoll f¨ur diagonal dominante Matrizen).
– Sog. unvollst¨andige Cholesky-Zerlegung f¨ur schwachbesetzte Matrizen, welche dasselbe Besetzungsmuster von A f¨ur die Dreiecksmatrix annimmt.
– L¨osung einer einfacheren Differentialgleichung, etc.
• Pr¨akonditionierung und die Suche nach einem guten Konditionierer lohnt sich h¨aufig, ist aber keine exakte Wissenschaft. . .