• Keine Ergebnisse gefunden

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

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

xij|A|φˆ ii = −X

i

aij|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 = AA

• 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

xiAs

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 : AikAAik

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:

ur Spalten k= 1, . . . , N :

wenn Akk < dann B Beachte: Kein Betrag

Matrix ist nicht positiv-definit; Abbruch.

wenn ende Ckk =

Akk

ur Zeilen i=k+ 1, . . . , N : Cik = CAik

ur Spaltenkk j =k+ 1, . . . , i : B Beachte: Schleife geht nur bisi

Aij Aij CikCjk ur ende

ur ende 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)CCX, (202)

⇒ C−1HX = diag(~e)CX, (203)

⇔C−1H([C]−1

| {z }

H˜

C)X = diag(~e)CX, (204)

y H˜X˜ = diag(~e) ˜X, (205)

H˜ = C−1H[C−1], (206)

X˜ = CX. (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ß

XX = S, (208)

Xorth = XY, (209)

⇒XorthXorth = 1, (210)

• L¨osung wieder per Cholesky-Zerlegung:

S = CC,(Cholesky) (211)

y XX = CC, (212)

⇒[C]−1XX = C, (213)

⇒ [C]−1XXC−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. . .