7. Iterative L¨osung linearer...

7. Iterative Losung

linearer Gleichungssysteme

1

Grundlagen und Wiederholung (1)

Die Grundlagen decken sich mit dem Stoff, der einen Teil des Kapitels

2 - Numerik ausmacht und bereits in Mathematik behandelt wurde.

Eine zentrale Rolle bei numerischen Berechnungen spielen lineare Glei-

chungssysteme

• Es sind die am haufigsten auftretenden numerischen Probleme

• Anwendungsgebiete sind z.B.

* fast alle naturwissenschaftlich-technischen Problemstellungen

vom Wetterbericht bis zur Warmeentwicklung auf einer Koch-

platte oder der Planung der Leiterbahnen auf Mikrochips,

* Bildverarbeitung oder z.B Beleuchtungsprobleme in der Com-

putergrafik,

* wirtschaftlichen Fragestellungen wie Versicherungskosten oder

Borsenkursvorhersage.

2


Losungsverfahren

• Die direkten Verfahren liefern eine mit Rundungsfehlern behaftete

Losung nach endlich vielen Schritten.

• Die iterativen Verfahren beginnen mit einer Anfangsnaherung und

produzieren eine verbesserte Naherungslosung nach endlich vielen

Schritten.

• Falls moglich wird das Problem mit einem direkten Verfahren be-

rechnet und anschließend werden die Rundungsfehler mit einem

iterativen Verfahren verringert.

3


Problemstellung: Berechne den Vektor x = (x1, x2, . . . xn) aus

a1,1x1 + a1,2x2 + · · · a1,nxn = b1

a2,1x1 + a2,2x2 + · · · a2,nxn = b2

·

·

an,1x1 + an,2x2 + · · · an,nxn = bn

oder in Matrix-Schreibweise

Ax = b

Bemerkung: Vektoren werden hier ohne Vektorpfeil geschrieben

4

Wiederholung Gauß-Verfahren (1)

Schulbeispiel:

E1 : x1 + x2 + 3x4 = 4E2 : − x2 − x3 − 5x4 = −7E3 : − 4x2 − x3 − 7x4 = −15E4 : 3x2 + 3x3 + 2x4 = 8

Mit Matrix:

A =

1 1 0 3

0 −1 −1 −5

0 −4 −1 −7

0 3 3 2

, x =

x1

x2

x3

x4

, b =

4

−7

−15

8

5


Erlaubte Transformationen zur Losung des Gleichungssystems

• Multiplizieren einer Zeile (Gleichung) mit einer Zahl verschieden

von Null

• Addieren eines Vielfachen einer Zeile zu einer anderen Zeile

• Vertauschen von Zeilen (Gleichungen) bzw. Spalten (Unbekann-

ten, entspricht Umnummerierung)

Mit Hilfe dieser Transformationen reduziere Gleichungssystem auf ein

Dreieckssystem.

6


Zweite Spalte des Schulbeispiels:

1 1 0 3

0 −1 −1 −5

0 −4 −1 −7

0 3 3 2

·

x1

x2

x3

x4

=

4

−7

−15

8

·

a3,i = a3,i − a3,2/a2,2 · a2,i b3 = b3 − a3,2/a2,2 · b2

a4,i = a4,i − a4,2/a2,2 · a2,i b4 = b4 − a4,2/a2,2 · b2

7


Schulbeispiel wird zu

E1 : x1 + x2 + 3x4 = 4E2 : − x2 − x3 − 5x4 = −7E3 : 3x3 + 13x4 = 13E4 : − 13x4 = −13

oder allgemein

a1,1 a1,2 · · · a1,n

a2,2 · · · a2,n

. . . ...

an,n

·

x1

x2

...

xn

=

b1

b2

...

bn

·

8


Ein Dreieckssystem ist leicht zu losen. Aus E4:

x4 = 1

x4 einsetzten in E3:

3x3 +13 = 13 → x3 = 0

x3, x4 einsetzten in E2:

−x2 − 5 = −7 → x2 = 2

x2, x3, x4 einsetzten in E1:

x1 +2+ 3 = 4 → x1 = −1

Der Gauß-Algorithmus ist von der Ordnung O(n3).

9

Norm von Vektoren und Matrizen

Unterschiedliche Definitionen von Normen sind moglich, hier nur die

sogenannte 2-Norm:

• 2-Norm oder euklidische Norm von Vektoren (Lange eines Vek-

tors):

||x||2 =

√

√

√

√

n∑

i=1

x2i

• 2-Norm oder Spektralnorm von Matrizen (Wurzel aus der Summe

der Quadrate der Diagonalelemente bei einer Diagonalmatrix):

||A||2 = maxx6=0||Ax||

||x||=

√

ρ(ATA)

10

Iterative Verfahren (1)

Da das Gauß-Verfahren

• O(n3) ist und damit oft zu lange dauert,

• viele Matrizen nur dunn besetzt sind, das Gauß-Verfahren aber die

Matrix fullt

• und das Gauß-Verfahren zu viele Rundungsfehler hat,

werden meist iterative Gleichungsloser verwendet. Beispiele:

• alle Programme fur Computergraphiken, wie z.B. PovRay,

• alle Programme fur Computersimulationen wie im IMH.

11


Idee: Fixpunktgleichung

Eine Gleichung der Form F(x) = x heißt Fixpunktgleichung. Ihre

Losungen, also der Werte, der die Gleichung F(x) = x erfullen, heißen

Fixpunkte

Definition:

Gegeben sei F : [a, b] → R,

x0 ∈ [a, b]. Die rekursive Folge

xn+1 := F(xn), n = 0,1, . . .

heißt Fixpunktiteration von F

zum Startwert x0.

X0

X1

X2

X3

X1

X2

X3

12


Anwendungsbeispiel

Bestimme die Nullstelle von p(x) = x3−x+0.3 oder lose die Gleichung

x = x3 +0.3.

Methode: Fuhre die Fixpunktiteration xn+1 = F(xn) = x3n+0.3, durch

Ergebnis:

• Startwert x0 = 0 : konvergiert gegen x = 0.3389 . . .

• Startwert x0 = 1: divergiert

Z.B. kann die Gleichung x = 2sin(x) + 1 nur iterativ gelost werden!

13


Iterative Losungsverfahren konstruieren ganz allgemein eine Losung

aus einer Startlosung:

x(0) → x(1) → x(2) → . . .

Angewendet auf die Berechnung der Losung eines linearen Gleichungs-

systems ist die Konvergenz abhangig von den Eigenschaften der Matrix

und den Details des iterativen Losungsverfahrens.

Es gibt zwei Gruppen von Verfahren

• die “klassischen” Verfahren, das Jacobi- und das Gauß-Seidel-

Iterationsverfahren und darauf aufbauend Relaxationsverfah-

ren, entwickelt im spaten 18. Jahrhundert, werden heute noch

angewandt.

• Krylov Unterraum-Methoden, die allgemeiner und oft schneller

sind, jedoch ist in vielen Fallen die Konvergenz unklar.

14

Richardson Iterationsverfahren (1)

Richardson Iterationsverfahren:

Idee: Formuliere ein passendes Fixpunktproblem

b = Ax = (A− I)x+ x ⇔ x = b+ (I− A)x := F(x)

Iterationsvorschrift:

x(t+1) = b+ (I− A)x(t) = x(t) + (b−Ax(t)) = x(t) + r(t)

Fur x(t+1) = x(t) gilt r(t) = 0 oder b− Ax(t) = 0

Der gesuchte Vektor x ist der Fixpunkt der Gleichung und wir oft auch

als x bezeichnet.

Der Vektor r(t) heißt Rest- oder Residuenvektor (manchmal Residu-

umsvektor oder kurz Residuum) und ist ein Maß fur den Fehler.

15


Abbruchbedingung der Iterationsverfahren: Der Betrag des Re-

siduenvektors ||r(t)|| = ||b − Ax(t)|| = ||A(x − x(t))|| oder besser der

relative Betrag ||r(t)||/||x(t)|| muss klein sein.

||x|| steht fur die 2-Norm oder der Lange des Vektors x

||x|| =

√

√

√

√

n∑

i=1

x2i

Fehlerbetrachtung im t-ten Schritt

x− x(t) = x− x(t−1) − (b− Ax(t−1))

= x− x(t−1) − (Ax− Ax(t−1))

= (I− A)(x− x(t−1))

16


Mit der 2-Norm

||x− x(t)||2 = ||(I− A)(x− x(t−1))||2

≤ ||(I− A)||2 · ||(x− x(t−1))||2

≤ ||(I− A)||22 · ||(x− x(t−2))||2

· · ·

≤ ||(I− A)||t2 · ||(x− x(0))||2

Konvergenz liegt sicher vor fur

||(I− A)|| < 1

oder A ist fast eine Einheitsmatrix.

(Die 2-Norm der Matrix ist die “Spektralnorm” und wird hier nicht

naher betrachtet, siehe Kapitel 2 zu Matrixnorm)

17

Jacobi Iterationsverfahren (1)

Verbesserung:

In vielen praktischen Problemen sind die Matrixeintrage auf der Dia-

gonalen groß. Betrachte dann das modifizierte Problem

D−1Ax = D−1b mit D = diag(A)

Diese Anderung andert die Losung nicht!

D−1A ≈ I bzw. ||(I−D−1A)|| < 1 ist aber eher erfullt als ||(I−A)|| < 1

Beispiel:

A =

5 2 12 7 4−1 2 8

, D−1 =

1/5 0 00 1/7 00 0 1/8

, I−D−1A =

0 2/7 1/82/5 0 4/8−1/5 2/7 0

,

Matrizen, deren großte Eintrage auf der Diagonalen sind, kommen z.B.

bei Beleuchtungsproblemen (Kapitel 2) oder bei sogenannten partiel-

len Differentialgleichungen vor (Kapitel 10).

18


Aus Ax = b wurde D−1Ax = D−1b, also, ersetze im Richardson-

Verfahren A → D−1A und b → D−1b:

x(t+1) = x(t) + (D−1b−D−1Ax(t)) oder

Dx(t+1) = Dx(t) + (b−Ax(t)) = b+ (D − A)x(t)

Das ist das Jacobiverfahren

Historisch wurde das Verfahren anders hergeleitet, so wie es auch meist

in der Literatur zu finden ist.

Ersetze Ax durch (D + A−D)x und forme den Ausdruck in eine Fix-

punktgleichung um

Ax = Dx+ (A−D)x = b

Das fuhrt zur selben Iterationsvorschrift:

Dx(t+1) = b+ (D −A)x(t)

19


Elementweise geschrieben:

ajjx(t+1)j = bj + ajjx

(t)j −

∑

k

ajkx(t)k

x(t+1)j =

1

ajj

bj −∑

k 6=j

ajkx(t)k .

• Multipliziere die Matrix ohne Diagonalelemente mit x(t)

• Ziehe das Ergebnis vom Vektor b ab

• Teile jedes Element durch das entsprechende Diagonalelement

• Fuhre das solange durch, bis der Residuenvektor klein ist

20

Gauß-Seidel Verfahren (1)

Weitere Verbesserung: Zerlege die Matrix in 3 Teile:

• D: Diagonalteil von A

• −L: Unterdiagonalteil von A

• −U : Oberdiagonalteil von A

Lose

b = Ax = (D − L− U)x = (D − L)x− Ux

uber die Iterationsvorschrift

(D − L)x(t+1) = b+ Ux(t) = b− (A− (D − L))x(t)

= (D − L)x(t) + (b−Ax(t)) = (D − L)x(t) + r(t)

21


Das Gauß-Seidel Verfahren lautet

x(t+1) = x(t) + (D − L)−1(b− Ax(t))

und ist gleich einer Richardson Iteration angewandt auf

(D − L)−1Ax = (D − L)−1b

Das Verfahren ist konvergent, falls ||I− (D − L)−1A||2 < 1

In vielen Anwendungen, basierend auf der diskretisierten Form soge-

nannter partieller Differentialgleichungen, kommt eine Variante dieser

Methode zum Einsatz (Kapitel 10).

22


Elementweise geschrieben:

(D − L)x(t+1) = b+ Ux(t)

ajjx(t+1)j +

∑

k<j

ajkx(t+1)k = bj −

∑

k>j

ajkx(t)k .

Berechnet x(t+1)1 , und damit x

(t+1)2 usw.

a11x(t+1)1 = b1 −

∑

k>1

a1kx(t)k

a22x(t+1)2 + a21x

(t+1)1 = b2 −

∑

k>2

a2kx(t)k

a33x(t+1)3 + a31x

(t+1)1 + a32x

(t+1)2 = b3 −

∑

k>3

a3kx(t)k

... = ...

23

Relaxationsverfahren

Die Wahl der Matrix, mit der die ursprunglich Fixpunktgleichung mo-

difiziert wird, ist im Prinzip beliebig (siehe Prakonditionierung)

Das SOR (successive over-relaxation) Verfahren modifiziert das Gauß-

Seidel Verfahren durch

(D − L) → (D

ω− L)

Eingesetzt und multipliziert mit ω folgt

(D − ωL) x(t+1) = [(1− ω)D + ωU ]x(t) + ωb

Die richtige Wahl vom ω kann den Algorithmus stark beschleunigen.

24

Prakonditionierung

Verallgemeinerung der obigen Methoden: Anstatt das Originalproblem

zu losen, betrachte

M−1Ax = M−1b

M sollte leicht zu invertieren sein und in der Nahe von A liegen, damit

M−1A leicht zu berechnen ist und in der Nahe der Einheitsmatrix liegt.

Die Richardson Iterationsvorschrift fur das prakonditionierte System

lautet

x(t+1) = (1−M−1A)x(t)+M−1b = x(t)+M−1(b−Ax(t)) = x(t)+M−1r(t)

mit hoffentlich kleiner Iterationsmatrix C = 1 − M−1A. Jacobi bzw.

Gauß-Seidel Verfahren entsprechen den Prakonditionierungsmatrizen

M = D bzw. M = D − L.

25

Krylov Unterraum-Methoden (1)

Problem:

Die bis jetzt vorgestellten Methoden konvergieren nur, wenn die Itera-

tionsmatrix ||C|| = ||1−M−1A|| klein ist.

Fur das Jacobi-, Gauß-Seidel- und SOR-Verfahren bedeutet das, dass

die Nebendiagonalelemente von A im Vergleich zu den Diagonalele-

menten klein sein mussen.

Krylov Unterraum-Methoden beruhen auf einer Verallgemeinerung

der bisher vorgestellten Fixpunktgleichungen. Aus dem prakonditio-

niertem Richardson Iterationsverfahren wird

x(t+1) = x(t) +M−1r(t) ⇒ x(t+1) = x(t) + α(t)p(t).

Je nach Wahl von α(t) und p(t) ergeben sich unterschiedliche Metho-

den, wobei nur sicher gestellt werden muss, dass x(t) jede mogliche

Richtung annehmen kann, so dass die Losung gefunden werden kann,

und dass der Betrag des Residuums ||b− Axk|| immer kleiner wird.

26


Was ist ein Krylov Unterraum?

Beispiel: Fur die einfachste Iterationsvorschrift: Die Richardson Itera-

tion

x(t+1) = x(t) + r(t).

gilt

r(t) = b−Ax(t)

= b−A(x(t−1) + r(t−1))

= b−A(x(t−1) + b−Ax(t−1))

= (I−A)(b− Ax(t−1))

= (I−A)r(t−1)

= (I−A)2r(t−2)

= . . . = (I− A)tr(0)

27


Mit dem Startvektor x(0) = 0 und damit r(0) = b−Ax0 = b folgt

x(t+1) = x(t) + r(t)

= x(t−1) +t∑

j=t−1

r(j)

= x(0) +t∑

j=0

r(j)

=t∑

j=0

(I−A)j b

28


Die Iterationslosung x(t) =∑t−1

j=1(I− A)j b ist somit Element des Un-

terraums

x(t) ∈ {b, Ab, . . . , At−1b} = Kt(A,b)

Dieser Raum wird als Krylov Unterraum der Dimension t, Kt(A,b)

bezeichnet.

Die Iterationslosung xt liegt bei einem Krylov-Unterraumverfahren in

Kt(A,b).

Je nachdem, wie die Iterationslosung innerhalb des Unterraums

bestimmt wird, gibt es verschiedene Iterations-Strategien.

29


Die bekanntesten sind

• Ritz-Galerkin Ansatz: r(t) orthogonal zu Kt(A,b) (z.B. CG).

Das bedeutet genauer:

xt ∈ Kt und rt = (b−Axt) ∈ K⊥t

• Petrov-Galerkin Ansatz: r(t) orthogonal zu einem allgemeinen Un-

terraum Lt (z.B. BI-CG oder QMR))

• Minimierung von ||r(t)||2 (z.B. MINRES oder GMRES)

• Hybride Verfahren (z.B. Bi-CGSTAB)

• Minimierung des Fehlers ||x− x(t)||2 (z.B. GMERR)

Es werden laufend neue Verfahren, angepasst an spezielle Probleme

entwickelt.

30

Einschub: Schreibweisen fur Vektoren und Matrizen

• Vektoren werden mit einem Vektorpfeil nur geschrieben, wenn es

sich um Vektoren im Raum handelt, ansonsten meist durch “fetten

Schriftsatz” gekennzeichnet.

• Ein Skalarprodukt zweier Vektoren wird geschrieben als

* a◦b oder a·b, haufig bei Ingenieuren und in der Schule verwendet

* 〈a,b〉 oder (a,b), Bra-Ket Schreibweise, beliebt bei z.B. Physi-

kern und angewandten Mathematikern

* aTb, Matrix-Schreibweise, beliebt bei Mathematikern

• Dementsprechend gilt z.B.

a ◦ (Mb) ≡ aTMb ≡ (a,Mb)

31

CG (1)

Das konjugierte Gradientenverfahren (CG) ist Grundlage vieler moder-

ner Iterationsverfahren (entwickelt 1952 von Stiefel und Hestens).

• Es konvergiert (falls keine Rundungsfehler vorliegen) fur positiv

definite und symmetrische Matrizen

a(x) = xTAx = (x, Ax) > 0; Ai,j = Aj,i

der Große n × n in n Schritten und ist somit eine schnelle und

sichere Methode zur Losung des Gleichungssystems.

• Die Rundungsfehler fuhren bei großen Systemen jedoch dazu, dass

man in vielen Fallen nicht n sondern ca. 3n Schritte anwendet.

• Es ist das einzige iterative Verfahren, fur das die Konvergenz im

Allgemeinen bewiesen wurde.

• Da jeder Schritt einer Matrixmultiplikation entspricht, lohnt sich

das Verfahren besonders fur dunn besetzte Matrizen.

32

CG (2)

Grundlage ist ein Optimierungsproblem. Die Losung von

Ax = b

ist ein Minimum der Funktion

f(x) =1

2(x, Ax)− (b,x)

und umgekehrt. Begrundung: Sie x der Losungsvektor von Ax = b, so

gilt fur einen beliebigen Vektor x+ p

f(x+ p) =1

2(x+ p, A(x+ p))− (b,x+ p)

= f(x) + (p, (Ax− b)) +1

2(p, Ap) = f(x) +

1

2(p, Ap) > f(x)

Das Verfahren setzt sich aus 2 Teilen zusammen: die Methode des

“steilsten Abstiegs” und die Methode der “konjugierten Richtung”.

33

CG (3)

Die Methode des “steepest descent” der Gradientenverfahren:

• Richtung: Beginnt man mit einem Vektor x+ p und mochte ent-lang des steilsten Abstiegs das Minimum erreichen, so muss manentlang des Negativen der Ableitung der Funktion f(x) nach x

gehen, also entlang

−f′(x) = −(Ax− b) = r.

Das Residuum gibt die Richtung des steilsten Abstiegs an.

• Relaxationsfaktor: Betrachte die Funktion f(x+αp) und bestimmtfur einen Vektor p das Minimum im Parameter α:

df(x+ αp)

dα=

d

dα(f(x) + (αp, (Ax− b)) +

1

2α2(p, Ap))

= (p, (Ax− b)) + α(p, Ap) = 0

oder

αopt =(p, r)

(p, Ap)

34

CG (4)

Algorithmus: Ausgangspunkt ist eine verallgemeinerte Richardson Ite-

rationsvorschrift x(i+1) = x(i) + αip(i).

• Verwende die Richtung des steilsten Abstiegs p = r und αopt.

• Setze x(0) und berechne r(0) = b−Ax(0)

• Fuhre folgende Schritte durch:

αi =(r(i), r(i))

(r(i), Ar(i))

x(i+1) = x(i) + αir(i)

r(i+1) = b−Ax(i+1) = r(i) − αiAr(i)

Das Verfahren konvergiert immer fur positiv definite und symmetri-

sche Matrizen, aber meist sehr langsam, da die Richtung des steilsten

Abstiegs zu einer oszillierenden Bewegung der Losung fuhren kann.

35

CG (5)

Beispiel aus G.Opfer:

Bestimme die Losung des Gleichungssystem(

1 00 10

)

·

(

x1x2

)

=

(

00

)

mit der Methode des steilsten Abstiegs. Das aquivalente System ist:

Bestimme das Minimum der Funktion

f(x) =1

2(x, Ax)− (b,x)

=1

2

(

x1 x2)

·

(

1 00 10

)

·

(

x1x2

)

−(

0 0)

·

(

x1x2

)

=1

2(x21 +10x22)

36

CG (6)

Wahle als Startvektor

x(0)T= (1.0,1.0)

Der Losungsvektor

oszilliert zur Losung

-0.1

-0.08

-0.06

-0.04

-0.02

0

0.02

0.04

0.06

0.08

0.1

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

min(x_1^2 + 10 x_2^2)

37

CG (7)

Die Methode der “konjugierten Richtung”:

Ausgangspunkt fur Methode des “steepest descent” war die Funktion

f . Als Richtung wurde der Gradient festgelegt und die Lange des

Gradientenvektors durch das Minimum der Funktion f

f(x+ αp) = f(x) + α(p, (Ax− b)) +1

2α2 (p, Ap)

= f(x) + α(p, Ax) +1

2α2 (p, Ap)− α(p,b)

als Funktion von α bestimmt. Verfahren zur Vermeidung der Oszilla-

tionen (nur zum Verstandnis, ohne Beweise):

• Die Richtungen, in der sich x(i) entwickelt, seien in einem Unter-

raum P (i) =< p(0), p(1), . . . p(i−1) >

38

CG (8)

• Der Vektor x(i) lasst sich dann schreiben als

x(i) = x(i−1) + α(i−1)p(i−1) = x(0) +i−1∑

j=0

α(j)p(j)

• Wahle als nachste Richtung zur Berechnung von x(i+1) nicht den

Gradienten, sondern eine Richtung p(i), so dass der 2. Summand

in f(x+ αp) fur alle Betrage von x in Richtung x(i) wegfallt

(p(i), Ax(i)) = 0

• Da x(i) eine Linearkombination der x(j), j = 1, . . . , i−1 ist, bedeutet

das

(p(i), Ap(j)) = 0 fur j ≤ i

39

CG (9)

• Die Vektoren p(0),p(1), . . . ,p(i) heißen paarweise zu A-konjugierte

Richtungen, da gilt

(p(k), Ap(j)) = 0 fur k 6= j

• Bestimme nun, wie gehabt, das optimale α.

αopt =(p(i), r(i))

(p(i), Ap(i))

• und berechne damit

x(i+1) = x(i) + αp(i)

40

CG (10)

Algorithmus: (bis jetzt)

Verwende einen Satz von konjugierten Suchrichtungen pi.

Setze x(0) und r(0) = b− Ax(0)

Dann fuhre folgende Schritte durch:

αi =(r(i),p(i))

(p(i), Ap(i))

x(i+1) = x(i) + αip(i)

r(i+1) = b−Ax(i+1) = r(i) − αiAp(i)

Die verbleibende Aufgabe ist es, die Suchrichtungen p(i) gunstig zu

wahlen, dass sich schnell eine gute Naherung ergibt.

41

CG (11)

Die “konjugierte Gradienten-Methode” von Hestens und Stiefel

Hier werden die konjugierten Richtungen aus den Restvektoren kon-

struiert werden.

p(0) = r(0)

p(i) = r(i) + β(i)p(i−1).

Ist β(i) = 0, ergibt sich die Methode des Gradientenverfahren. Die

Koeffizienten β(i) werden so gewahlt, dass die Richtungen p(i) wie

gefordert zueinander konjugiert sind.

0 = (p(m), Ap(i)) = (r(m), Ap(i))+β(m)(p(m−1), Ap(i)) fur i = 0, . . . ,m−1

42

CG (12)

Da

(p(m−1), Ap(i)) = 0 fur i < m− 1

folgt

β(m−1) = −(r(m), Ap(m−1))

(p(m−1), Ap(m−1))

Ohne Beweis: Jetzt gilt fur den neuen Residuenvektor

(r(i+1),p(j)) = 0 fur j = 0, . . . , i,

d.h. der Vektor steht senkrecht auf dem Raum P (i). Da sich der Vek-

torraum bei jedem Iterationsschritt um eine Dimension vergroßert, ist

sicher gestellt, dass das Verfahren nach n Schritten eine Losung findet.

Nach einigen Umrechnungen ergibt sich CG-Algorithmus.

43

CG (13)

Endgultiger Algorithmus (ohne Beweis!):

Starte mit p(0) = r(0) = b−Ax(0). Dann fuhre folgende Schritte durch:

α(i−1) =|r(i−1)|2

(p(i−1), Ap(i−1))

x(i) = x(i−1) + αi−1p(i−1)

r(i) = r(i−1) − αi−1Ap(i−1)

β(i) =|r(i)|2

|r(i−1)|2

p(i) = r(i) + β(i)pi−1

44

CG (14)

Nochmal das Beispiel

aus G.Opfer(

1 00 10

)

·

(

x1x2

)

=

(

00

)

Wahle als Startvektor

wieder x(0)T= (1.0,1.0)

Der Losungsvektor

oszilliert nicht mehr und

die Losung wird in 2

Schritten erreicht.

-0.1

-0.08

-0.06

-0.04

-0.02

0

0.02

0.04

0.06

0.08

0.1

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

min(x_1^2 + 10 x_2^2)

45

CG (15)

Rechenaufwand pro Iteration

1 Matrixmultiplikation A,

2 Skalarprodukte,

3 AXPYs

Es werden (theoretisch) n Iterationen benotigt.

Einfache Erweiterung fur nicht-symmetrische Matrizen

Ist die Matrix nicht symmetrisch und positiv definite, betrachte

ATAx = AT b

und konstruiere die Losung in Ki(ATA,AT b) unter Verdoppelung der

Anzahl der Matrix-Vektor Multiplikationen (geht besser durch Erwei-

terung des Krylov-Unterraums).

46

7. Iterative L¨osung linearer...

Documents

Transcript of 7. Iterative L¨osung linearer...