Grundlagen der Numerischen MathematikIWenn das Integral R ju 0j2dx in J mit R ju 0jdx ersetzt wird,...

Grundlagen der Numerischen Mathematik

a.o.Univ.Prof. Mag.Dr. Stephen Keelinghttp://math.uni-graz.at/keeling/

Literatur:Numerical Analysis von R.L. Burden und J.D. Faires

Unterlagen:http://math.uni-graz.at/keeling/teaching.html

1

Inhaltsverzeichnis IEinführung

Motivierendes Beispiel: BildverarbeitungApproximation von Ableitungen und IntegralenMittelwertsatzOptimalitätsbedingung fürs BildverarbeitungsbeispielRandwertproblem fürs BildverarbeitungsbeispielFehleranalyseSatz von TaylorO-NotationZwischenwertsatzIterative Bestimmung einer NullstelleKonvergenzgeschwindigkeitStabilitätDarstellung von Zahlen im ComputerAufrundung und AbrundungSignifikante ZiffernFehler der FließkommadarstellungAuslöschungHorner Algorithmus

Lineare Gleichungssysteme: Direkte MethodenGaußsche EliminationRückwärts SubstitutionDiskretisierung einer PDGAutomatisierung der Gaußschen EliminationPseudo-Kode der Gaußschen Eliminationflops der Gaußschen EliminationRekursive Lösungen: EvolutionsgleichungLU ZerlegungRekursive Lösungen mit LUPivot StrategienPseudo-Kode für Gaußsche Elimination mit Pivot-SucheLU Zerlegung mit PermutationSkalierte Pivotsuche2

Inhaltsverzeichnis IIPseudo-Kode mit Skalierter PivotsucheTotale PivotsucheInverse und DeterminanteBesondere Eigenschaften einer MatrixGerschgorin SatzPseudo-Kode für den Cholesky AlgorithmusBandmatrizenSparse Speicherformat

Lineare Gleichungssysteme: Iterative MethodenVektornormenÄquivalenz von NormenQualitative Eigenschaften der NormenMatrixnormenCharakterisierung bekannter MatrixnormenSpektralradiusKonvergente MatrizenZerlegungen für IterationenJacobi MethodePseudo-Kode der Jacobi MethodeGauß-Seidel MethodeSymmetrische Gauß-Seidel MethodeApproximierte InversenSOR MethodeKonvergenz von Iterativen VerfahrenKonjugierte GradientenPräkonditionierungFehler Abschätzungen

Eigenwerte und EigenvektorenVektoriterationInverse VektoriterationBerechnung aller EigenwerteDie Householder MethodeHessenberg Matrizen3

Inhaltsverzeichnis IIIPseudo-Kode zur Transformation auf Hessenberg oder Tridiagonale FormQR AlgorithmusGivens TransformationenKonvergenz des QR-AlgorithmusPseudo-Kode für den QR-Algorithmus

AusgleichsproblemeLineare RegressionPolynome RegressionSingulärwert ZerlegungEigenschaften der SWZPseudoinverseBerechnung der SWZBerechnung der SWZVereinfachung auf Bidiagonale FormPseudo-Kode zur Transformation auf Bidiagonale FormDiagonalisierende IterationFolge der Bidiagonalen MatrizenPseudo-Kode für den QR Algorithmus zur Singulärwert-Zerlegung

InterpolationLagrange Interpolierende PolynomeGenauigkeit globaler InterpolationIterierte InterpolationNevilles AlgorithmusDividierte DifferenzenHermite InterpolationStückweise Polynome InterpolationStückweise Kubisch Hermite InterpolationSplinesDie Kanonischen SplinesGlattheitsbedingungenInterpolation mit Kubischen SplinesTensor Produkte für 2D Interpolation

Nicht Lineare Gleichungen4

Inhaltsverzeichnis IVPseudo-Kode für das BisektionsverfahrenFixpunktiterationPseudo-Kode für eine FixpunktiterationVergleich zweichen Bisektions- und FixpunktiterationNewton VerfahrenPseudo-Kode für das Newton VerfahrenPseudo-Kode für Sekant VerfahrenVergleich zweichen Newton und BisektionsiterationAsymptotische KonvergenzrateSysteme von Nicht Linearen GleichungenAbstiegsverfahrenImplizite FixpunktiterationNewton Verfahren für SystemeRichtungsfelder der LösungsverfahrenEntrauschen Ergebnisse

Numerisches Differenzieren und IntegrierenApproximation der AbleitungenRichardson ExtrapolationNumerische IntegrationTrapez-RegelSimpson-RegelGenauigkeit einer geschlossenen Newton-Cotes FormelGenauigkeit einer offenen Newton-Cotes FormelMittelpunkt-RegelZusammengesetzte Mittelpunkt-RegelZusammengesetzte Simpsons-RegelnZusammengesetzte Trapez-RegelnRomberg IntegrationGauß Quadratur

Gewöhnliche DifferentialgleichungenDas Eulersche VerfahrenWärmegleichungKonservative Verfahren5

Inhaltsverzeichnis VWellengleichungCrank-Nicholson VerfahrenRunge-Kutta MethodenRunge-Kutta-Fehlberg MethodenKonsequenz, Stabilität, KonvergenzAbsolute StabilitätSteife SystemeRandwertproblemeDiskretisierung eines Sturm-Liouville RandwertproblemsRayleigh-Ritz Methode

6

Was ist Numerische Mathematik?

I Experiment zu Hause:I Im Taschenrechner: 2 eingeben, dann 10×

√· und 10×x2.

I Subtrahiere 2. Ergebnis 6= 0. Warum?I Zahlen werden nicht genau gespeichert!

I Wie soll man Ax = b, f (x) = 0, u′′(x) = f (x) lösen?Antwort: Näherungsweise.

I Ein gutes numerisches Verfahren soll nicht erlauben, dassFehler sich schlecht aufbauen!

7

Motivierendes Thema: Signal- und Bildverarbeitung

I Dieses Beispiel wird mehrmals in verschiedenenKontexten analysiert.

I v(x) ist ein gegebenes verrauschtes Signal.I u?(x) ist das gewünschte Signal.I u(x) ist eine Abschätzung von u?(x).

8


I Wie bestimmt man u ≈ u??I Ansatz: Ω = (0,1), u? ≈ u = argminwJ(w),

J(u) =∫

Ω(u − v)2dx + µ

∫Ω|u′|2dx

I Qualitativ:I µ→∞⇒ u = Mittelwert von v (übergeglättet)I µ→ 0⇒ u = v (untergeglättet)

Das beste µ? ist irgendwo in der Mitte.I Optimalitätsbedingung zur Bestimmung von u?

Diskreter Ansatz...

9

Motivierendes Thema: Signal- und BildverarbeitungI Die Daten werden üblicherweise sowieso diskret gegeben:

vi = v(xi) auf einem Gitter

xi = ih, i = 0, . . . ,N, h = 1/N

x0

0

xN

1

xi

h

I Fragestellung: Es gelten

hN∑

i=0

[u(xi)− v(xi)]2h→0−→

∫Ω|u(x)− v(x)|2dx

µhN∑

i=1

[u(xi)− u(xi−1)

h

]2h→0−→ µ

∫Ω|u′(x)|2dx

unter welchen Bedingungen?10

Motivierendes Thema: Signal- und BildverarbeitungSatz (Mittelwertsatz): Wenn f ∈ C[a,b] differenzierbar in (a,b)ist, dann ∃c ∈ (a,b) sodass f ′(c) = [f (b)− f (a)]/[b − a].I Also ∃ξi ∈ (xi−1, xi) sodass

u(xi)− u(xi−1)h

=u(xi)− u(xi−1)

xi − xi−1= u′(ξi)

und

µhN∑

i=1

[u(xi)− u(xi−1)

h

]2= µ

N∑i=1

[u′(ξi)]2h ≈ µ∫

Ω|u′(x)|2dx

Satz (Mittelwertsatz für Integrale): Wenn f ∈ C[a,b] undg(x) ≥ 0, ∀x ∈ [a,b], dann ∃c ∈ (a,b) sodass∫ b

af (x)g(x)dx = f (c)

∫ ba

g(x)dx

I Themen: Approximation von Ableitungen und Integralen.11

Motivierendes Thema: Signal- und BildverarbeitungI Für Daten {vi} und Approximationen {ui ≈ u(xi)} seien

u = 〈u0, . . . ,uN〉T und v = 〈v0, . . . , vN〉T. Wir minimieren:

Jh(u) = hN∑

i=0

(ui − vi)2 + µhN∑

i=1

[ui − ui−1

h

]2I Optimalität: ∇Jh(u) = ∇uJh(u) = 0.∂Jh∂uk︸︷︷︸

0


I oder ∇Jh(u) = 2h(u − v) + 2(µ/h)Du wobei

D =

1 −1 0−1 2 −1

. . .. . .

. . .

−1 2 −10 0 −1 1

I Optimalität:

0 =1

2h∇Jh(u) = u − v +

µ

h2Du

oder das lineare Gleichungssystem:[ µh2

D + I]

u = v

13


I Optimalitätssystem komponentenweise:

−µui+1 − 2ui + ui−1h2

+ ui = vi , i = 1, . . . ,N − 1

−µu1 − u0h2

+ u0 = v0, i = 0

−µ −uN + uN−1h2

+ uN = vN , i = N

wobei in der zweiten und dritten Zeilen,“u0 = u−1” beziehungsweise “uN = uN+1”

(virtuelle Randbedingungen u′ = 0) zu sehen sind.I Kontinuum-Ansatz für Optimalität: Richtungableitung sind

Null für alle Störungen w :

0 =δJδu

(u; w) = lim�→0

dd�

J(u + �w)

14


I Gilt wenn u das Randwertproblem erfüllt:{−µu′′ + u = v , Ω

u′ = 0, ∂Ω

I Vergleich mit dem diskreten Ergebnis zeigt:

u′′(xi) ≈ui+1 − 2ui + ui−1

h2

Approximation der zweiten Ableitung.I Hausaufgabe: Schreibe eine Approximation Jh(u) für

J(u) =∫

Ω(u − v)2dx + µ

∫Ω

√|u′|2 + ε2dx

und leite die Optimalitätsbedingung ∇Jh(u) = 0 her.

15

Zusammenfassung unserer Fragestellungen

I Wie leitet man Formeln zur Approximation eine Integralsher?

I Formeln für Ableitungen?I Wann wissen wir, dass diese Formeln konvergieren wenn

h→ 0? Und zwar schnell?I Wie löst man die entstehenden Gleichungssysteme?I Wenn das Integral

∫Ω |u

′|2dx in J mit∫

Ω |u′|dx ersetzt

wird,I ist das Entrauschen-Ergebnis viel besser, aberI das Gleichungssystem ist nicht linear.I Wie löst man das nicht lineare Problem?

Dies werden wir später mit dem Newtonschen Verfahrenlösen.

16

FehleranalyseI Fehler heissen Abbruchsfehler, wenn sie durch eine

Diskretisierung und trotz exakter Arithmetik entstehen, z.B.der Rest in

u′′(xi)−1h2

[u(xi+1)− 2u(xi) + u(xi−1)]

wird nicht berechnet sondern abgebrochen.I Fehler heissen Rundungsfehler, wenn Zahlen einer

Rechnung nicht exakt gespeichert werden können.I Diese können üblicherweise so grafisch dargestellt

werden:

17

FehleranalyseI Wie können wir die Abbruchsfehler abschätzen?

Satz (Taylor): Sei u ∈ Ck [a,b] und u(k+1) existiert in [a,b]. Seix0 ∈ [a,b]. ∀x ∈ [a,b], ∃ξ zwischen x und x0 sodass

u(x) =k∑

m=0

u(m)(x0)m!

(x − x0)m +u(k+1)(ξ)(k + 1)!

(x − x0)k+1

I Anwendung für die zweite Ableitung: xi+1 − xi = h

u(xi+1) = u(xi) + u(1)(xi)h + u(2)(xi)h2

2+ u(3)(xi)

h3

6+ u(4)(ξi+1)

h4

24

u(xi−1) = u(xi)− u(1)(xi)h + u(2)(xi)h2

2− u(3)(xi)

h3

6+ u(4)(ξi−1)

h4

24

Summe:

u(xi+1)− 2u(xi) + u(xi−1)h2

= u′′(xi) + F (h)

18

FehleranalyseI Der Abbruchsfehler ist:

F (h) = [u(4)(ξi+1) + u(4)(ξi−1)]h2

24

I Man sagt, F (h) = O(h2).Def: F (x) = O(G(x)) für x → x0 (x0 ∈ R ∪ {±∞}) wennlimx→x0 |F (x)/G(x)|

KonvergenzBemerkung: Wir haben gezeigt, der Abbruchsfehler dieserApproximation der zweiten Ableitung konvergiert zu Null mitVerfeinerung. Wir haben aber noch nicht gezeigt,[ µ

h2D + I

]u = v . . .

{−µu′′ + u = v , Ω

u′ = 0, ∂Ω

dass die Lösung des linearen Gleichungssystems zur exaktenLösung des Randwertproblems mit Verfeinerung konvergiert,

maxi |ui − u(xi)| → 0, h→ 0aber der Taylorsatz gehört zum Werkzeug.

Bemerkung: Die Optimalitätsbedingung ∇Jh(u) = 0 für

J(u) =∫

Ω(u − v)2dx +

∫Ω

√|u′|2 + ε2dx

ist nicht linear. Wir brauchen ein Verfahren, um die Nullstelle zubestimmen...20

Nullstellen

I Einfaches Beispiel: f : R→ R, finde die Nullstelle x0. Wirhaben durch Versuch und Irrtum x1 und x2 gefunden,sodass f (x1)f (x2) < 0. Gibt es eine Nullstelle zwischen x1und x2?

Satz (Zwischenwert): Wenn f ∈ C[a,b] und k zwischen f (a)und f (b) liegt, dann existiert c ∈ (a,b) sodass f (c) = k gilt.I Gegenbeispiele:

21

Nullstellen

I Bekannte Verfahren zur Bestimmung einer Nullstelle:Bisektion

f (a)f (b) < 0 . . . (a + b)/2 = c →{

a, f (c)f (b) < 0b, f (a)f (c) < 0

und Newton

xk+1 = xk − f (xk )/f ′(xk ), k = 1,2, . . .

I Unsere Fragen:I Ist ein Verfahren konvergent?I Was soll das Abbruchskriterium sein?I Was ist die Geschwindigkeit der Konvergenz?

22

Konvergenz und Stabilität

Def: Die Folge {xk} konvergiert zu x∗ mitKonvergenzgeschwindigkeit O(yk ) wenn xk = x∗ +O(yk ).

I Beispiel: xk = (k + 1)/k2 = 0 +O(1/k):limk→∞ |xk/(1/k)| = limk→∞(k + 1)/k = 1.

Ähnlich yk = (k + 3)/k3 = 0 +O(1/k2).

I Falls (Rundungs)fehler in einer Iteration auftauchen,können sie sich schlecht aufbauen?

I D.h. ist das Verfahren stabil?I Wir stellen die Fehler als Anfangsstörung E0 dar...

23

Konvergenz und StabilitätDef: Sei Ek der Fehler eines Verfahrens nach k Operationen,wenn der Anfangsfehler E0 ist. Wenn gilt

|Ek | ≤ ck |E0| 0 < c 6= c(k)sagt man, das Fehlerwachstum ist linear. Wenn nur gezeigtwerden kann,

|Ek | ≤ cak |E0| 0 < c 6= c(k),a > 1sagt man, das Fehlerwachstum ist exponentiell. DaFehlerwachstum in der Regel unvermeidlich ist, sagen wir, einVerfahren mit höchstens linearem Fehlerwachstum ist stabil.I Beispiel: Wir berechnen die Folge 3−k .

I Methode 1: p0 = 1, pk = 13 pk−1. Probiere sie! Sie ist stabil!I Methode 2: p0 = 1, p1 = 13 , pk =

103 pk−1 − pk−2, z.B.

p2 = 10313 − 1 =

132 (mit exakter Arithmetik). Probiere sie!

Sie ist nicht stabil!I Die Fehler in diesem Beispiel sind Rundungsfehler. Wie

entstehen sie genau?

24

Darstellung von Zahlen im Computer

I Beispiel: IEEE Fließkommazahl einfacher Genauigkeit(32 bits, 4 bytes).

e0· · · →e7 m1 · · · · · · m23

8bits︷︸︸︷ 8bits︷︸︸︷ 8bits︷︸︸︷ 8bits=1byte︷︸︸︷σ︸︷︷︸

e︸︷︷︸

m

σ →Vorzeichen, e→Exponent, m→Mantissa

e = 〈e0, . . . ,e7〉, m = 〈m1, . . . ,m23〉, ek ,mk ∈ {0,1}

V = (−1)σ, E = −127 +7∑

k=0

ek2k , M = 1 +23∑

k=1

mk2−k

Fließkommazahl: x = V · 2E ·M

25


I Beispiel: x = 1

x = (1)10 = (1.0)2 × 20 (M ∈ [1,2))M = (1)2 ⇒ mk = 0, k = 1, . . . ,230 = E = −(127)10 + e⇒ e = (127)10

Wie bringt man (127)10 auf die binäre Darstellung?

Algorithmus für die dezimalen Ziffer von 127:12710 = 12 +

710 ,

1210 = 1 +

210 ,

110 = 0 +

110

Algorithmus für die binären Ziffer von 127:127

2 = 63 +12 ,

632 = 31 +

12 ,

312 = 15 +

12 ,

152 = 7 +

12 ,

72 = 3 +

12 ,

32 = 1 +

12 ,

12 = 0 +

12

Also e = (127)10 = (1111111)2

x (IEEE): 00111111 10000000 00000000 0000000026


I Beispiel: x = 0.375

Algorithmus für die dezimalen Ziffer von 0.375:0.375× 10 = 3 + 0.75,0.750× 10 = 7 + 0.50, 0.50× 10 = 5 + 0Algorithmus für die binären Ziffer von 0.375:0.375× 2 = 0 + 0.75,0.750× 2 = 1 + 0.50, 0.50× 2 = 1 + 0Alsox = (0.375)10 = (0.011)2 = (1.1)2 × 2−2 (M ∈ [1,2))M = (1.1)2 ⇒ m1 = 1,mk = 0, k = 2, . . . ,23e − (127)10 = E = −(2)10 ⇒ e = (125)10 = (01111101)2

x (IEEE): 00111110 11000000 00000000 00000000

27


I Beispiel: x = 12.375.

x = (12)10 + (0.375)10x = (1100)2 + (0.011)2 = (1100.011)2

x = (1.100011)2 × 23, (M ∈ [1,2))M = (1.100011)2 ⇒ m1 = 1,m5 = 1,m6 = 1, sonst mk = 0e − (127)10 = E = (3)10 ⇒ e = (130)10 = (10000010)2

x (IEEE): 01000001 01000110 00000000 00000000

28

Darstellung von Zahlen im ComputerI Beispiel: x = −0.1

Algorithmus für die binären Ziffer von 0.1:0.1× 2 = 0 + 0.2, 0.2× 2 = 0 + 0.4,0.4× 2 = 0 + 0.8, 0.8× 2 = 1 + 0.6,0.6× 2 = 1 + 0.2→ zurück zum Anfang

Also x = −(0.1)10 = −(0.00011)2 = −(1.10011)2 × 2−4

Probe: x = −∑∞

k=1 2−4k − 12

∑∞k=1 2

−4k = − 110

M = (1.10011)2 ⇒ m1 = 1, {m4k = m4k+1 = 1}5k=1,sonst mk = 0?

e − (127)10 = E = −(4)10 ⇒ e = (123)10 = (1111011)2

Mit Aufrundung (m23 = 1)x (IEEE): 10111101 11001100 11001100 11001101Mit Abrundung (m23 = 0)x (IEEE): 10111101 11001100 11001100 11001100

29

Darstellung von Zahlen im ComputerI Fortsetzung des Beispiels: x = −0.1

Mit Aufrundung (m23 = 1)x (IEEE): 10111101 11001100 11001100 11001101

→ (0.100000001490116119384765625)10Mit Abrundung (m23 = 0)x (IEEE): 10111101 11001100 11001100 11001100

→ (0.0999999940395355224609375)10Mit einfacher Genauigkeit gibt es keine speicherbarenZahlen zwischen diesen!

I IEEE Fließkommazahl doppelter Genauigkeit (64 bits, 8bytes)

x = V · 2E ·M, E = −1023 +10∑

k=0

ek2k

V = (−1)σ, M = 1 +52∑

k=1

mk2−k

30


Fehler in der Darstellung?

Def: Wenn p̃ eine Approximation zu p ist, ist der absoluteFehler |p − p̃|. Der relative Fehler ist |p − p̃|/|p| wenn p 6= 0.

Def: p̃ approximiert p zu t signifikanten Ziffern wenn t ∈ N diegrößte Zahl ist, bei der der relative Fehler kleiner als 5× 10−tist.

I Beispiel: x = −0.1, x̃ = −0.10000000149 · · · (Aufrundung)

5× 10−9 < |x − x̃ ||x |

=1.49 · · · × 10−9

10−1< 5× 10−8 ⇒ t = 8

31

Darstellung von Zahlen im ComputerI Relativer Fehler der Fließkommazahl-Darstellung fl(x), mit

Abrundung und p-Ziffer Genauigkeit M =∑p−1

k=0 mk2−k ,

|x − fl(x)||x |

=|V · 2E (

∑∞k=0)mk2

−k − V · 2E (∑p−1

k=0)mk2−k |

|V · 2E (∑∞

k=0)mk2−k |

=

∣∣∣∣∣∣∞∑

k=p

mk2−k

∣∣∣∣∣∣/ ∣∣∣∣∣

∞∑k=0

mk2−k∣∣∣∣∣≥1

≤∞∑

k=p

mk2−k

≤∞∑

k=p

2−k = 2−p∞∑

k=0

2−k =2−p

1− 2−1= 2−(p−1)(< tol)

I IEEE einfache Genauigkeit, p = 24,2−p+1 ≈ 1.19× 10−7, t = 7

I IEEE doppelte Genauigkeit, p = 53,2−p+1 ≈ 2.22× 10−16, t = 16

32


I Verlust von signifikanten Ziffern: Auslöschung.Bildlich gesehen:

x1 = (−1)σ2EM1x2 = (−1)σ2EM2

gleich anders

x2 − x1 E #######000 · · · 000

(−1)σ2Ẽ (M2 −M1) Ẽ #######000 · · · 000

Nur die Ziffer ####### überleben Subtraktion. In exakterArithmetik wäre der letzte Bereich 000 · · · 000 mit richtigenZiffern besetzt.

33


I Beispiel: Berechne die Nullstellen vonp(x) = x2 + 62.10x + 1

mit 4-Dezimal-zifferiger Genauigkeit.

−62.10 + [(62.10)2 − 4.000 · 1.000 · 1.000]12

2.000

→ −62.10 + [3856− 4.000]12

2.000→ −0.04000

2.000= −0.02000

6= x1 = −0.01610723

I Besserer Weg: A− B = (A2 − B2)/(A + B)⇒

−b +√

b2 − 4ac2a

=−b2 + (b2 − 4ac)2a[b +

√b2 − 4ac]

→ −0.01610

I Auch anfällig für Auslöschung: Skalarprodukt

34


I Auch anfällig für Auslöschung: PolynomauswertungI Beispiel: An der Stelle x = 4.71 das Polynom

p(x) = x3 − 6x2 + 3x − 0.149mit 3-Dezimal-zifferiger Genauigkeit auswerten.exakt:104.487111− 133.1046 + 14.13− 0.149 = −14.6364893-zifferig:105− 133 + 14.1− 0.149→ −14.0, 4% relativer Fehler

I Horner Algorithmus:p(x) = x [x2 + 6x + 3]− 0.149 = x [x(x − 6) + 3]− 0.149

p ← x − 6, p ← x · p, p ← p + 3, p ← x · p, p ← p − 0.149

Ergebnis: -14.5, 0.25% relativer Fehler

35

Lineare Gleichungssysteme: Direkte MethodenGaußsche Elimination zur Lösung von Ax = b. Beispiel:I Das System:

E1E2E3

3x1 + 6x2 9x3 = −122x1 + 5x2 − 8x3 = −11

x1 − 4x2 − 7x3 = −10

I E2 ← E2 − 23E1, E3 ← E3 −13E1,

3x1 + 6x2 + 9x3 = −12x2 − 14x3 = −3

− 6x2 − 10x3 = −6

I E3 ← E3 + 61E2,

3x1 + 6x2 + 9x3 = −12x2 − 14x3 = −3− 94x3 = −24

36

Gaußsche Elimination und Rückwärts SubstitutionI Rückwärts Substitution:

3x1 + 6x2 + 9x3 = −12 ⇒ x1 = 13(−12− 6x2 − 9x3) = −27847

x2 − 14x3 = −3 ⇒ x2 = −3 + 14x3 = 2747 ↑− 94x3 = −24 ⇒ x3 = 1247 ↑

I Zur Lösung von Ax = b kann man die Schritte auf dieerweiterte Matrix [A|b] durchführen:

[A|b] =

1 1 0 | 42 1 −1 | 13 −1 −1 | −13

∈ R3×4 (Rn×(n+1))I Wenn es mehrere rechte Seiten gibt,

Ax i = bi , i = 1,2, . . . ,m oderAX = B, X = [x1, . . . ,xm], B = [b1, . . . ,bm]

kann man die Schritte gleichzeitig auf die erweiterte Matrix[A|B] durchführen. Aber nicht wenn bi von x i−1 abhängt!

37

Diskretisierung einer PDGI Die Diskretisierung [ µh2 D + I]u = v des Randwertproblems,{

−µu′′ + u = v , Ωu′ = 0, ∂Ω

wird sehr groß mit Verfeinerung, aber sie ist trotzdemimmer tridiagonal; deswegen ist es ziemlich billig,Gaußsche Elimination für dieses System durchzuführen.

I Bessere Motivation: Ω = (0,1)2, ∆u = uxx + uyy{−µ∆u + u = v , Ω

∂u/∂n = 0, ∂Ω

1 · · · N+1N+2 · · · 2N+2· · ·· · ·· · ·

(N+1)2

uxx + uyy ≈ui+1,j − 2uij + ui−1,j

h2+

ui,j+1 − 2uij + ui,j−1h2

38

Diskretisierung einer PDG

I Gleichungssystem:

. . .

. . .. . .

. . .. . . 0

− µh2

− µh2

4µh2

+ 1 − µh2

− µh2

0. . .

. . .. . .

. . .. . .

u0,0

...uN,0

...u0,N

...uN,N

=

v0,0

...vN,0

...v0,N

...vN,N

� -O(N) � -O(N)I Obwohl die Matrix dünn besetzt ist (sparse, d.h. ganz

wenig nicht triviale Elemente), führt Gaußsche Eliminationzur Auffüllung (fill-in) der Zwischendiagonalen.

39

Automatisierung der Gaußschen Elimination

I Um die Lösungsschritte der Gaußschen Elimination zuautomatisieren, müssen wir diese symbolisch schreiben:

A ← [A|b]

2. Zeile:{

A21 = A21/A11 (Multiplikator: m21)A2j = A2j − A21 · A1j , j = 2, . . . ,n(+1)

3. Zeile:{


...

i . Zeile:{

Ai1 = Ai1/A11 (Multiplikator: mi1)Aij = Aij − Ai1 · A1j , j = 2, . . . ,n(+1)

i = 2, . . . ,n

I Ergebnis:

A11m21

mn1

...

� Pivot

� neues b

� neue Zielspalte

-

sonst Null40

Automatisierung der Gaußschen Elimination

I Mit der nächsten Zielspalte:

3. Zeile:{


...

i . Zeile:{

Ai2 = Ai2/A22 (Multiplikator: mi2)Aij = Aij − Ai2 · A2j , j = 3, . . . ,n(+1)

i = 3, . . . ,n

I Ergebnis:

A11A22m21

mn1

...

m32

mn2

...

...

� Pivot� neues b

� neue Zielspalte

-

sonst Null

I Wird mit Zielspalten fortgesetzt, bis man ein Dreieck für dieMultiplikatoren hat.

41

Automatisierung der Gaußschen EliminationI Pseudo-Kode zur Gaußschen Elimination:A

flops der Gaußschen EliminationI Wie viele Operationen kostet Gaußsche Elimination?

I Divisionenk = 1, ...,n− 1

i = k + 1, ...,n → (n − k)

} n−1∑k=1

(n − k)

I Multiplikationen

k = 1, ...,n− 1i = k + 1, ...,n → ×(n − k)

j = k + 1, ...,n(+1) → (n + 1− k)

n−1∑k=1

(n−k)(n+1−k)

I Subtraktionen

...ebenso :n−1∑k=1

(n − k)(n + 1− k)

I Nützliche Summen: (vgl.∫ m

0 kp−1dk = mp/p)

m∑k=1

1 = m,m∑

k=1

k =m(m + 1)

2,

m∑k=1

k2 =m(m + 1)(2m + 1)

643

flops der Gaußschen EliminationI Kosten für Gaußsche Elimination:

I Divisionen

n−1∑k=1

(n−k) = nn−1∑k=1

1−n−1∑k=1

k = n(n−1)− (n − 1)n2

=n(n − 1)

2

I Multiplikationen

n−1∑k=1

(n − k)(n + 1− k)︸︷︷︸n(n+1)−(2n+1)k+k2

= (n2 + n)n−1∑k=1

1− (2n + 1)n−1∑k=1

k +n−1∑k=1

k2 =

(n2 + n)(n − 1)− (2n + 1) (n − 1)n2

+(n − 1)n[2(n − 1) + 1]

6=

n(n2 − 1)3

I Subtraktionen: auch n(n2 − 1)/3I Gesamte Operationen:

n(n − 1)/2 + 2n(n2 − 1)/3 = n(n − 1)(4n + 7)/6 = O(n3)44

flops der Gaußschen Elimination

I Kosten für Rückwärts Substitutionen:I Divisionen: nI Subtraktionen: n − 1I Additionen:

i = n− 1, ...,1j = i + 1, ...,n → (n − i)

} n−1∑i=1

(n − i) = n(n − 1)2

I Multiplikationen: auch n(n − 1)/2I Gesamte Operationen:

2n − 1 + 2n(n − 1)/2 = n2 + n − 1 = O(n2)I Botschaft: Rückwärts Substitution ist viel billiger!I Hausaufgabe: Einen Pseudo-Kode für Vorwärts

Substitution schreiben und die gesamten Operationenausrechnen.

45

Rekursive Lösungen: EvolutionsgleichungI Wenn man sehr viele Systeme lösen muss,

Ax i = bi , i = 1,2, . . . ,mist es vorteilhaft wenn man die folgenden Schrittedurchführen kann:I Gaußsche Elimination nur einmal mit Kosten O(n3)

machen,I das Ergebnis irgendwie speichern undI dann m-Mal Rückwärts und Vorwärts Substitutionen mit

Kosten O(mn2) machen.Sonst sind die Kosten O(mn3)!

I Beispiel: Evolutionsgleichung,{u′ = Bu

u(0) = u0

zeitlich

−→diskretisieren

u(tk+1)− u(tk )tk+1 − tk

= Bu(tk+1)

oder mit τ = tk+1 − tk ,

Auk+1 = [I − τB]uk+1 = uk , k = 0,1, . . .46

LU Zerlegung

I Wie kann man die O(n3) Arbeit für Gaußsche Eliminationgünstig speichern?

I Definiere A(1) = A. Ergebnis nach der ersten Zielspalte istA(2) = M(1)A(1) wobei:

A(2) =0

0

...M(1) =

0· · ·01−m21

−mn1

... I(n−1)

. . .

. . .

Probe: Sei A =

z1...zn

, M(1)A = z1z2 −m21z1

...

47

LU Zerlegung

I Ergebnis nach der zweiten Zielspalte ist A(3) = M(2)A(2)

wobei:

A(3) =0

0

...0

0

...M(2) =

0· · ·00· · ·0

1

10

0

...

−m32

−mn2

... I(n−2). . .

. . .

I Am Ende von Gaußscher Elimination:A(n) = M(n−1)A(n−1) = M(n−1) · · ·M(1)A

wobei A(n) eine obere Dreiecksmatrix ist.I Definiere U = A(n). Bemerke:

A = M(1)−1 · · ·M(n−1)−1UI Was sind die Inversen von M(i) explizit?

48

LU ZerlegungI Bemerke:

1 0 · · · 0. . .

. . .

I(n−1)

m21

mn1

...∗

0· · ·01-m21

-mn1

... I(n−1)

. . .

. . .

= M(1)−1M(1) = I(n)

I Es gelten

M(i)−1 =

001. . .0 1

mi+1,i

mn,i

0 ... I(n−i)

. . .

. . .

i = 1, . . . ,n − 1

I Bemerke:1 0 · · · 0. . .

. . .

I(n−1)

m21

mn1

...∗

0· · ·011 0 · · · 0

m32

mn2

... I(n−2)

. . .

. . .

=

0· · ·011 0 · · · 0m21

mn1

m32

mn2

... ... I(n−2)

. . .

. . .49

LU ZerlegungI Es gilt:

0· · ·011

1

01mn,n−1

0 · · · 0m21

mn1

m32

mn2 · · ·

. . .. . .

. . ....... ...

M(1)−1 · · ·M(n−1)−1 = = L

I Wir haben gezeigt:

Satz: Wenn alle Pivot-Elemente a(k)kk nicht Null sind und alleMultiplikatoren mjk = a

(k)jk /a

(k)kk wohl definiert sind, kann die

Matrix mit A = LU faktorisiert werden, wobei

U = {uij}, uij =

{a(n)ij , 1 ≤ i ≤ n, i ≤ j ≤ n

0, 2 ≤ i ≤ n, 1 ≤ j < i

L = {`ij}, `ij =

mij , 2 ≤ i ≤ n, 1 ≤ j < i

1, 1 ≤ i ≤ n, j = i0, 1 ≤ i ≤ n, i < j ≤ n

50

Rekursive Lösungen mit LU

I Mit dieser Zerlegung kann man mehrere sequentielleSysteme günstig lösen: Ax i = bi ⇔ Ly i = bi ,Ux i = y i .

I Beispiel: Evolutionsgleichung,

{u′ = Bu

u(0) = u0

zeitlich

−→diskretisieren

u(tk+1)− u(tk )tk+1 − tk

= Bu(tk+1)

mit τ = tk+1 − tk ,

Auk+1 = [I − τB]uk+1 = uk , k = 0,1, . . .

I Durch die LU Zerlegung A = LU, O(n3),

Vorwärts Substitution: Luk+12 = uk , O(n2)

Rückwärts Substitution: Uuk+1 = uk+12 , O(n2)

k = 0,1, . . .

51

Pivot StrategienPivot Strategien:I Beispiel:[

0.003000 59.145.291 −6.130

] [x1x2

]=

[59.1746.78

]I Exakte Lösung:

x∗1 = 10.00, x∗2 = 1.000

I Mit 4-zifferiger Arithmetik,

m21 =5.291

0.003000= 1763.6̄6→ 1764

I Erster Schritt von Gaußscher Elimination,[0.003000 59.14 59.17

(1764) −104300 −104400

]I In exakter Arithmetik wäre es gewesen,[

0.003000 59.14 59.17(1763.6̄6) −104309.376̄6 −104309.376̄6

]52

Pivot Strategien

I Nun Rückwärts Substitution mit 4-zifferiger Arithmetik

x2 =−104400−104300

→ 1.001 ≈ x∗2

x1 =59.17− (59.14)(1.001)

0.003000→ −10.00 6= 10.00 = x∗1

I Hausaufgabe: Zeige, wenn die Zeilen getauscht werden,[5.291 −6.130 46.78

0.003000 59.14 59.17

]bekommt man

x2 → 1.000, x1 → 10.00

53

Pivot Strategien

Tauschmit I(·)

k →p →

I Pseudo-Kode zur Gaußschen Elimination mit Pivot-Suche:A

LU Zerlegung mit PermutationI Wie kann man mit Pivot-Suche die Matrix schön zerlegen?

Def: Eine Permutationsmatrix P hat genau einen nicht trivialenEintrag (mit dem Wert 1) in jeder Spalte und in jeder Zeile, undes gilt PPT = I.I Beispiel:

P =[

0 11 0

]A =

[0 13 2

]PA =

[3 20 1

]Bemerkung: {I(i) : i = 1, . . . ,n} ist eine gegebenePermutation von den Zahlen {1, . . . ,n}. Die Permutationsmatrix

Pi,j = δI(i),j erfüllt Ã = PA wobei Ãij = AI(i),j .

I Wegen dem obigen Satz haben wir Ã = LU von demAlgorithmus mit Pivot-Suche, und PA = LU wobei

Pi,j = δI(i),j Lij =

A(n)I(i),j , i > j1, i = j0, i < j

Uij =

{A(n)I(i),j , i ≤ j0, i > j

55

Skalierte Pivotsuche

Weitere Pivotsuche-StrategienI Das 2× 2 Beispiel hätten wir so umschreiben können:[

30.00 5914005.291 −6.130

] [x1x2

]=

[59170046.78

](← ×104)

I Erster Schritt von Gaußscher Elimination ohne Tauschen,

m21 =5.29130.00

→ 0.1764[

30.00 591400 591700(0.1764) −104300 −104400

]führt zur gleichen Lösung wie vorher ohne Pivotsuche:

x2 → 1.001, x1 → −10.00 6= 10.00 = x∗1

I Neue Strategie: Man sollte die Zeilen skalieren.

56

Skalierte PivotsucheI Pseudo-Kode zur G.E. mit skalierter Pivotsuche:

A

Totale PivotsucheBemerkung: Kosten für Tauschen und Skalieren sind O(n2)

Weitere Pivotsuche-Strategien: Totale PivotsucheI Virtueller Zeilentausch mit einem Zeiger I(·).I Virtueller Spaltentausch mit einem Zeiger J(·).I Bestimmung des Pivot-Elements:

(p,q)=argmax (k≤i,j≤n) |A I(i),J(j)|I Zusätzliche Arbeit ist O(n3).I Geeignet wenn das Problem sehr verschiedene Skalen

hat, d.h. wenn die Matrix A sehr steif ist.I Verwendung:

Zeilentausch: Ax = b → P1Ax = LUx = P1bSpaltentausch: Ax = b → AP2P2x = LUP2x = b

Beides: (P1AP2)P2x = LUP2x = P1b

Ly = P1b, Uz = y , x = PT2 z58

Inverse und DeterminanteBemerkung: Falls die Inverse notwendig ist, kann man siegünstig durch die Lösung des folgenden Systems berechnen:

AX = I ⇒ X = A−1

Die direkte Berechnung ist aber viel teurer.

Bemerkung: Die Determinante,

δk = det

a11 · · · a1k... ...ak1 · · · akk

kann man mit Gaußscher Elimination so berechnen,

δk = a(1)11 · a

(2)22 · · · a

(k)11

falls keine Pivotsuche notwendig ist. Die direkte Berechnung istaber viel teurer.59

Diagonal Dominante MatrizenWann wird keine Pivotsuche notwendig?

Def: Eine Matrix {aij} = A ∈ Rn×n ist streng diagonal dominantwenn

|aii | >n∑

i 6=j=1|aij |, i = 1, . . . ,n

Bei Gleichheit ist sie schwach diagonal dominant.

I Beispiel: Wir haben die folgende Diskretisierung,{−µu′′ + u = v , Ω

u′ = 0, ∂Ω−→ Au =

[ µh2

D + I]

u = v

wobei

A =

µh2 + 1 −

µh2

− µh22µh2 + 1 −

µh2

. . .. . .

. . .

− µh22µh2 + 1 −

µh2

− µh2µh2 + 1

60

SPD MatrizenI A ist streng diagonal dominant:

(x0, xN) i = 0,N :µh2 + 1 > | −

µh2 |

√

(xi) 0 < i < N :2µh2 + 1 > | −

µh2 |+ | −

µh2 |

√

Satz: Für eine streng diagonal dominante Matrix, kannGaußsche Elimination ohne Pivotsuche stabil durchgeführtwerden.

Def: Eine symmetrische Matrix A ∈ Rn×n ist positiv definit(SPD) wenn xTAx > 0, ∀x ∈ Rn, x 6= 0.

Bemerkung: Wenn A ∈ Rn×n symmetrisch ist, sind dieEigenwerte reel und es gilt

min{λ ∈ σ(A)} ≤ xTAx

xTx≤ max{λ ∈ σ(A)}

Wenn min{λ ∈ σ(A)} > 0 gezeigt werden kann, ist A SPD.61

Gerschgorin SatzI Dies kann man für das obige Beispiel mit dem folgenden

Satz zeigen.

Satz (Gerschgorin): Für {aij} = A ∈ Rn×n gilt

σ(A) ⊂n⋃

i=1

Ri wobei Ri =

z ∈ C : |z − aii | ≤n∑

i 6=j=1|aij |

Weiters liegen genau k Eigenwerte in ∪kj=1Rij wenn dieseMenge einen leeren Schnitt mit den anderen Scheiben hat.I Beispiel: Wo liegen die Eigenwerte für das obige A?

R0 = B(µh2 + 1,

µh2 )

Ri = B(2µh2 + 1,

2µh2 )

RN = B(µh2 + 1,

µh2 )

C

&%'$

1

µ

h2+ 1

2µh2

+ 1

4µh2

+ 1

@

�

�

@62

Notwendigkeit einer PivotsucheI Das obige A ist symmetrisch und es gilt

1 ≤ min{λ ∈ σ(A)} ≤ xTAx

xTxalso ist A SPD.

Satz: Für eine SPD Matrix kann Gaußsche Elimination ohnePivotsuche stabil durchgeführt werden.

Botschaft: Für viele Matrizen in Anwendungen ist einePivotsuche nicht notwendig, aber es ist sicher nicht immer so.

Folgesatz: Eine Matrix A ∈ Rn×n ist SPD genau dann wenn∃L ∈ Rn×n eine untere Dreiecksmatrix sodass A = LLT.I Mit der Gleichung für LLT = A = {aij}, L = {`ij},

aij =min(i,j)∑

k=1

`ik`jk

a11 = `211a21 = `21`11a12 = `11`21

. . .

landet man auf dem folgenden Algorithmus:63

Pseudo-Kode für den Cholesky AlgorithmusL 11 =

√A 11

for i=2,...,nL i1 = A i1 / L 11

endfor j=2,...,n-1

s = 0for k=1,...,j-1

s = s + L jk * L jkendL jj =

√A jj− s

for i=j+1,...,ns = 0for k=1,...,j-1

s = s + L ik * L jkendL ij = [A ij - s]/L jj

endends = 0for k=1,...,n-1

s = s + L nk * L nkendL nn =

√A nn− s

I Hausaufgabe: Mit demobigen Gleichungssystemleite diese Methode her.

I Bemerke: Der CholeskyAlgorithmus ist nur für SPDMatrizen.

I Hausaufgabe: Schätze dieflops des Kodes ab.

64

BandmatrizenI Die Bandbreite einer solchen Matrix,

A =

a1,1 · · · a1,q+1 0 · · · 0...

. . .. . .

. . .. . .

...

ap+1,1. . .

. . .. . .

. . . 0

0. . .

. . .. . .

. . . an−q,n...

. . .. . .

. . .. . .

...0 · · · 0 an,n−p · · · an,n

∈ Rn×n

is p + q + 1, die kleinste Anzahl der benachbartenDiagonalen, zu denen die nicht trivialen Elementeeingeschränkt sind.

I Bemerke:

p ≤ n − 1, q ≤ n − 1, Bandbreite ≤ 2n − 1.

Wenn p + q + 1

Bandmatrizen

I Um Speicherplatz zu reduzieren, wird A als Liste derDiagonalen gespeichert:

B =

0 · · · 0 a1,1 · · · · · · a1,q+1... . .

.. .. ...

......

...

0 . .. ...

......

... an−q,n

ap+1,1...

......

... . ..

0...

......

... . ..

. .. ...

an,n−p · · · · · · an,n 0 · · · 0

∈ Rn×(p+q+1)

I Bemerke: Die nicht trivialen Elemente von A = {ai,j}lassen sich bezüglich der Elemente von B = {bα,β} sodarstellen:

ai,j = bi,j−i+p+1

66

Sparse SpeicherformatI Hausaufgabe: Schreibe einen Pseudo-Kode zur

Implementierung des Gauß Algorithmus für A bezüblich Bmit möglichst wenig arithmetischen Operationen.

I Es gibt auch sparse format:

A = sparse(i,j,w,m,n) ⇒ A ∈ Rm×n, Ai(k),j(k) = w(k)

I Beispiel:

A =

µh2 + 1 −

µh2

− µh22µh2 + 1 −

µh2

. . .. . .

. . .

− µh22µh2 + 1 −

µh2

− µh2µh2 + 1

A = spdiags([ µh2 + 1; (

2µh2 + 1)ones(N − 1,1);

µh2 + 1],0,N + 1,N + 1)

− spdiags( µh2 ones(N + 1,1),+1,N + 1,N + 1)− spdiags( µh2 ones(N + 1,1),−1,N + 1,N + 1)

67

Lineare Gleichungssysteme: Iterative MethodenIterative Verfahren zur Lösung von Ax = b.I Um Ax = b zu lösen, werden wir solche Iterationen

formulieren,

xk+1 = T xk + c, k = 0,1, . . .

I Um Konvergenz zu untersuchen, müssen wir eine gewisseInfrastruktur aufbauen.

Def: Eine Vektornorm ist eine Funktion ‖ · ‖ : Rn → R mit denEigenschaften:

1. ‖x‖ ≥ 0, ∀x ∈ Rn (∑n

i=1 xi keine Norm)2. ‖x‖ = 0⇔ x = 0, (|x1| keine Norm für n > 1)3. ‖αx‖ = |α|‖x‖, ∀x ∈ Rn, ∀α ∈ R (

∑ni=1 x

2i keine Norm)

4. ‖x + y‖ ≤ ‖x‖+ ‖y‖, ∀x ,y ∈ Rn

(f (√∑n

i=1 x2i ) keine Norm wenn f nicht konvex)

68

VektornormenSatz: Die folgenden sind Normen:

‖x‖`p =

[n∑

i=1

|xi |p] 1

p

,1 ≤ p

Äquivalenz von Normen

Satz: In Rn sind alle Normen äquivalent, d.h. für 2 beliebigeNormen ‖x‖A und ‖x‖B auf Rn ∃c1, c2 3

c1‖x‖A ≤ ‖x‖B ≤ c2‖x‖A ∀x ∈ Rn

Insbesondere:

Satz: ∀x ∈ Rn,

‖x‖∞ ≤ ‖x‖2 ≤√

n‖x‖∞

Beweis: Für x ∈ Rn sei |xk | = ‖x‖∞. Dann gilt

‖x‖2∞ = |xk |2 ≤n∑

i=1

|xi |2 = ‖x‖22 ≤n∑

i=1

|xk |2 = n|xk |2 = n‖x‖2∞

Hausaufgabe: Zeige Äquivalenz der Normen ‖ · ‖1 und ‖ · ‖2.70

Qualitative Eigenschaften der NormenI Soll man einen Fehler ‖x − x∗‖ mit der `1-, `2- oder`∞-Norm messen? Mit einer anderen Norm?

I Antwort: Es hängt vom Kontext ab!I Beispiel: Die gemessenen pH-Werte einer chemischen

Lösung sind f = 〈a,b, . . . ,b〉 ∈ Rn+1.

Dieses Beispiel übertreibt die Situation, in der ein einziger Wert einstatistischer Ausreißer ist, während die anderen Werte näher beieinander liegen.

I In einem Blick ist b die beste Abschätzung vom pH.I Hausaufgabe: Für e = 〈1, . . . ,1〉 zeige

a + nb1 + n

= arg minc∈R‖f−ce‖22 während b = arg minc∈R ‖f−ce‖1

I Also läßt die ‖ · ‖2-Norm sich von Ausreißern beeinflussen.I Botschaft: Die ‖ · ‖1-Norm ist statistisch robuster.

71

MatrixnormenAuch müssen wir den Abstand zwischen Matrizen messen.

Def: Eine Matrixnorm ist eine Funktion ‖ · ‖ : Rn×n → R mit denEigenschaften:

1. ‖A‖ ≥ 0, ∀A ∈ Rn×n

2. ‖A‖ = 0⇔ A = 03. ‖αA‖ = |α|‖A‖, ∀A ∈ Rn×n, ∀α ∈ R4. ‖A + B‖ ≤ ‖A‖+ ‖B‖, ∀A,B ∈ Rn×n

5. ‖AB‖ ≤ ‖A‖‖B‖, ∀A,B ∈ Rn×n

Def: Für eine Vektornorm ‖ · ‖V ist die entsprechendeinduzierte (oder natürliche) Matrixnorm ‖ · ‖M definiert durch:

‖A‖M = max‖x‖V=1

‖Ax‖V

Hausaufgabe: Zeige dass diese Funktion ‖ · ‖M eineMatrixnorm ist.

72

Charakterisierung bekannter MatrixnormenSatz: Die Matrixnormen ‖ · ‖1, ‖ · ‖2 und ‖ · ‖∞ sind für eineMatrix {Aij} = A ∈ Rn×n explizit so gegeben:

‖A‖1 = max‖x‖1=1

‖Ax‖1 = max1≤j≤n

n∑i=1

|aij |

‖A‖2 = max‖x‖2=1

‖Ax‖2 =[max{λ ∈ σ(ATA)}

] 12

‖A‖∞ = max‖x‖∞=1

‖Ax‖∞ = max1≤i≤n

n∑j=1

|aij |

Def: Eine Matrixnorm ‖ · ‖M ist mit einer Vektornorm ‖ · ‖Vkompatibel (oder verträglich) wenn gilt

‖Ax‖V ≤ ‖A‖M‖x‖V, ∀x ∈ Rn, ∀A ∈ Rn×n

Hausaufgabe: Zeige für eine Matrix {Aij} = A ∈ Rn×n, dieFrobenius Matrixnorm,

‖A‖F =

n∑i,j=1

|Aij |2 12

ist mit der Vektornorm ‖x‖2 kompatibel.73

SpektralradiusDef: Sei Rn mit einer gewissen Vektornorm ‖ · ‖V versehen. Fürxk ,x ∈ Rn gilt limk→∞ xk = x genau dann wenn ∀� > 0, ∃K , 3∀k ≥ K , ‖x − xk‖V < �.

Satz: Für xk = 〈(x1)k , . . . , (xn)k 〉,x = 〈x1, . . . , xn〉 ∈ Rn giltlimk→∞ xk = x genau dann wenn limk→∞(xi)k = xi .

Beweis: Hausaufgabe.

Def: Die Eigenwerte einer Matrix A ∈ Rn×n sind die Nullstellendes charakteristischen Polynoms

p(λ) = det(A− λI), σ(A) = {λ : p(λ) = 0}

und die entsprechenden Eigenvektoren sind die nicht trivialenLösungen (x 6= 0) zur Gleichung Ax = λx .

Def: Der Spektralradius einer Matrix A ∈ Rn×n ist

ρ(A) = max{|λ| : λ ∈ σ(A)}74

SpektralradiusI ρ(A) ist der größte Vergrößerungsfaktor.

Beispiel: A = SΛS−1, S = [x1,x2],

Λ =

[1 00 2

],

Λê1 = ê1, Ax1 = x1Λê2 = 2ê2, Ax2 = 2x2 ←

Satz: Sei ‖ · ‖? eine beliebige Matrixnorm auf Rn×n, die von derVektornorm ‖ · ‖? auf Rn induziert ist. Dann gilt

ρ(A) ≤ ‖A‖?, ∀A ∈ Rn×n

Beweis: Wähle λ, x so aus, dass Ax = λx , |λ| = ρ(A) und‖x‖? = 1 gelten. Dann gilt

ρ(A) = |λ| = |λ|‖x‖? = ‖λx‖? = ‖Ax‖? ≤ ‖A‖?‖x‖? = ‖A‖?.

I Wie kann man eine Lösungsiteration

xk+1 = T xk + c

von Ax = b herleiten? Zerlegung: A = M + N, wobei75

Konvergente MatrizenSysteme mit M leicht zu lösen sind:

Mxk+1 = b − Nxk

so c = M−1b und T = −M−1N. Die exakte Lösung x? = A−1berfüllt

Mx? = b − Nx?,

Subtraktion,(xk+1 − x?) = T (xk − x?)

Wenn gilt ‖T‖ < 1 für eine induzierte (oder verträgliche)Matrixnorm, dann gilt

‖xk+1−x?‖ ≤ ‖T‖‖xk−x?‖, ‖xk+1−x?‖ ≤ ‖T‖k‖x0−x?‖k→∞−→ 0

Def: T ∈ Rn×n ist konvergent wenn gilt

limk→∞

(T k )ij = 0, ∀1 ≤ i , j ≤ n76

Konvergente MatrizenI Beispiele:

A =[

2 00 1

], Ak =

[2k 00 1

]−→×

[0 00 0

]B =

[ 12 00 13

], Bk =

[ 12k 00 13k

]−→

[0 00 0

]C =

[ 12 014

12

], Ck =

[ 12k 0k

2k+112k

]−→

[0 00 0

]D =

[−32 −

165

65

52

], Dk =

[1 23 4

]−1 [( 910)

k 00 ( 110)

k

] [1 23 4

]−→

[0 00 0

]I ‖D‖1 = 5.7, ‖D‖2 = 4.4922, ‖D‖∞ = 4.7 und ρ(D) = 0.9.I Gibt es eine Matrixnorm ‖ · ‖M wobei ‖D‖M < 1?

Satz: Sei A ∈ Rn×n gegeben. ∀� > 0 gibt es eine induzierteMatrixnorm ‖ · ‖M� wobei gilt ρ(A) + � > ‖A‖M� .I Bemerkung: Diese Matrixnorm ‖ · ‖M� hängt von A ab.77

Zerlegungen für IterationenSatz: Die folgenden sind äquivalent:

1. T ∈ Rn×n ist konvergent2. limk→∞ ‖T k‖ = 0 für eine induzierte Matrixnorm3. limk→∞ ‖T k‖ = 0 für alle induzierten Matrixnormen4. ρ(T ) < 15. limk→∞ T kx = 0, ∀x ∈ Rn

I In der Zerlegung A = M + N,

Mxk+1 = b − Nxk

welches M ist geeignet? Besonders wenn A diagonaldominant ist, ist die folgende Zerlegung natürlich:

A = D + L + U, M = D,N = L + U

wobei D eine diagonale Matrix, L ist eine streng untereDreiecksmatrix und U ist eine streng obere Dreiecksmatrix.

78

Jacobi Methode

I Die Jacobi Methode: TJ = −D−1(L + U)

xk+1 = D−1[b − (L + U)xk ]= D−1[b − (L + U + D)= Axk + Dxk ]= xk + D−1[b − Axk ]= xk + D−1rk

wobei r = b − Ax das Residuum für gegebenes x ist.I Die Jacobi Methode, komponentenweise,

∑ji

aijxkj =n∑

j=1

aijxj = bi , 1 ≤ i ≤ n

I Im folgenden Pseudo-Kode wird xi sequentiell in einerSchleife i = 1, . . . ,n berechnet:

79

Pseudo-Kode der Jacobi Methoded i = A ii, i=1,...,nx i = b i, i=1,...,nfor it=1,...,itmax

for i=1,...,nr i = b ifor j=1,...,n

r i = r i - A ij * x jenddx i = r i / d i

endfor i=1,...,n

x i = x i + dx iendif (||dx|| < tol * ||x||)

breakend

end

% Abbruchsparameter: itmax� O(n)% tol = relative Fehlertoleranz

80

Jacobi MethodeI Wenn ein neuer Wert x i+dx i im Pseudo-Kode verfügbar

ist, warum nicht in x i sofort speichern? Somit wären dieZeilen im Pseudo-Kode

dx i = r i / d ix i = x i + dx i

neben einander innerhalb einer einzigen i-Schleife.I Da im Pseudo-Kode die Schleife i=1,...,n so läuft,

würde die neue Iteration so aussehen:∑ji

aijx(k)j =

n∑j=1

aijxj = bi , 1 ≤ i ≤ n

oder

x (k+1)i =1aii

bi −∑ji

aijx(k)j

I In Matrixnotation,

81

Gauß-Seidel MethodeI In Matrixnotation,

xk+1 = D−1[b − Lxk+1 − Uxk ]

oder die Gauß-Seidel Methode,

(D + L)xk+1 = b − Uxk

I Die Zerlegung A = M + N ist M = D + L, N = U undTGS = −M−1N = −(D + L)−1U.

I Umgeschrieben,

xk+1 = (D + L)−1︸︷︷︸Vorwärts Substitution

[b − Uxk ]

= (D + L)−1[b − (D + L + U)︸︷︷︸=A

xk + (D + L)xk ]

= xk + (D + L)−1[b − Axk ] = xk + (D + L)−1rk

82

Symmetrische Gauß-Seidel Methode

I Es gibt eine gewisse Asymmetrie bei der Gauß-SeidelMethode.

I Wir hätten die Schleife ebensogut i=n,...,1 laufenlassen können.

I Mit der Schleife in dieser entgegengesetzten Richtungbekommt man xk+1 = (D + U)−1[b − Lxk ].

I Mit beiden Schleifen bekommt man die SymmetrischeGauß-Seidel Methode,{

xk+12 = (D + L)−1[b − Uxk ]

xk+1 = (D + U)−1[b − Lxk+12 ]

I Hausaufgabe: Finde die Iterationsmatrix TSGS für dieSymmetrische Gauß-Seidel Methode.

83

Approximierte Inversen

Bemerkung: Die Zerlegung A = M + N kann soumgeschrieben werden:

xk+1 = −M−1Nxk + M−1b= [−M−1(M + N) + I]xk + M−1b= [I −M−1A]xk + M−1b

und so wird M−1 eine Approximierte Inverse benannt, da‖I −M−1A‖ klein sein sollte.

Hausaufgabe: Oben sind genannt worden: MJ = D undMGS = (D + L).I Leite die Approximierte Inverse MSGS für die Symmetrische

Gauß-Seidel Methode.I Zeige, wenn A symmetrisch ist, ist MSGS symmetrisch,

während MGS im allgemeinen nicht symmetrisch ist.

84

SOR MethodeI Gauß-Seidel noch einmal komponentenweise,

x (k+1)i =1aii

bi −∑j

aijx(k)j

ω+x (k)i (1−ω)Man könnte das alte und das neue so gewichten.

I In Matrix-Notation,

xk+1 = D−1[b − Lxk+1 − Uxk ]ω + xk (1− ω)Dxk+1 = [b − Lxk+1 − Uxk ]ω + Dxk (1− ω)

(D + ωL)xk+1 = [b − Uxk ]ω + Dxk (1− ω)xk+1 = (D + ωL)−1[(1− ω)D − ωU]xk + ω(D + ωL)−1b

I Diese Methode heisst Successive Over-Relaxation (SOR)I 0 < ω < 1⇒ under relaxed, d.h. sie konvergiert wegen

Dämpfung, wenn sie bei ω = 1 möglicherweise nichtkonvergiert.

I ω > 1⇒ over relaxed, d.h. sie konvergiert schneller wegenExtrapolation, wenn sie bei ω = 1 schon zuverlässigkonvergiert.

85

Konvergenz von Iterativen VerfahrenI Die Symmetrische Successive Over-Relaxation (SSOR)

Methode:{xk+

12 = (D + ωL)−1[(1− ω)D − ωU]xk + ω(D + ωL)−1b

xk+1 = (D + ωU)−1[(1− ω)D − ωL]xk+12 + ω(D + ωU)−1b

Konvergenz

Satz: Wenn für das System Ax? = b gilt x? = T x? + c, dannerfüllt das iterative Verfahren xk+1 = T xk + c die Konvergenzlimk→∞ xk = x?, genau dann wenn ρ(T ) < 1.

Bemerkung: Wenn für eine induzierte Norm ‖T‖ < 1 gezeigtwerden kann, folgt ρ(T ) < 1.

Satz: Wenn A streng diagonal dominant ist, giltlimk→∞ xk = x?, ∀x0 ∈ Rn, für die Jacobi und Gauß-SeidelMethoden.86

Konvergenz von Iterativen VerfahrenI Wie wählt man ω für SOR aus? Sodass ρ(TSOR) < 1,

wobei TSOR = (D + ωL)−1[(1− ω)D − ωU].Satz: Wenn für {aij} = A ∈ Rn gilt aii 6= 0, 1 ≤ i ≤ n, folgtρ(TSOR) ≥ |ω − 1|.

Deswegen gilt ρ(TSOR) < 1 nur für ω ∈ (0,2) gilt. Umgekehrt:

Satz (Ostrowski-Reich): Wenn A SPD ist und 0 < ω < 2 gilt,folgt für SOR die Konvergenz limk→∞ xk = x?, ∀x0 ∈ Rn.

I Beispiel: Für die Diskretisierung Au = [ µh2 D + I]u = v desRandwertproblems,{

−µu′′ + u = v , x ∈ (0,1)u′ = 0, x ∈ {0,1}

haben wir schon gezeigt, A ist streng diagonal dominantsowohl SPD. Deswegen konvergieren Jacobi undGauß-Seidel sowohl SOR, ∀ω ∈ (0,2).

87

Konjugierte GradientenDie Methode der Konjugierten Gradienten

Satz: Sei f ∈ C2(Ω). Wenn x? ∈ Ω erfüllt ∇f (x?) = 0 und∇2f (x?) SPD ist, ist x? ein lokales Minimum für f .

Explizite Bespiele:I Für f (x) = bTx gelten ∇f (x) = b und ∇2f (x) = 0.

∂

∂xk

n∑i=1

bixi =n∑

i=1

bi∂xi∂xk

=n∑

i=1

biδik = bk

I Für f (x) = xTMx , M ∈ Rn×n, gelten ∇f (x) = (M + MT)xund ∇2f (x) = M + MT.

∂

∂xk

n∑i=1

n∑j=1

xiMijxj =n∑

i=1

n∑j=1

(∂xi∂xk

Mijxj + xiMij∂xj∂xk

)=

n∑i=1

n∑j=1

(δikMijxj + xiMijδjk

)=

n∑j=1

Mkjxj +n∑

i=1

Mikxi

88

Konjugierte GradientenSatz: Sei A ∈ Rn×n SPD. Dann löst x? das System Ax = bgenau dann wenn x? die Funktion g(x) = 12x

TAx − bTxminimiert.

Beweis: (⇒) Aus Ax? = b folgt ∇g(x?) = Ax? − b = 0. Weil ASPD ist, hat g genau einen kritischen Punkt. Da ∇2g(x) = ASPD ∀x ∈ Rn ist, ist g überall konvex, und x? ist ein lokalesMinimum. Aus der konvexen Analysis ist x? ein globalesMinimum. (⇐) Mit einem Minimum g(x?) von g ∈ C∞(Rn) giltnotwendigerweise ∇g(x?) = 0. Da A SPD ist, ∃!x? 3∇g(x?) = Ax? − b = 0.

Vorausgesetzt ist A SPD, und Ax = b soll gelöst werden.

Entwicklung der Methode der Konjugierten Gradienten:

I Sei x0 eine Approximation der Lösung x?.I Sei v1 eine gegebene Suchrichtung zur Minimierung von g.I Wo ist das Minimum von g entlang v1 von x0 weg?89

Konjugierte Gradienten

I Antwort durch die Ketten-Regel,

ddt

g(x0+tv1) = ∇g(x0+tv1)·v1 = [A(x0+tv1)−b]·v1 = 0

und das Minimum ist in x1 = x0 + t1v1 wobei

t1 =(−Ax0 + b) · v1

v1 · Av1=

r0 · v1

v1 · Av1

I Eigenschaft des Gradienten:∇g(x1) ⊥ {x : g(x) = c}, c = g(x1)

I Orthogonalität −r1 · v1 = ∇g(x1) · v1= [A(x0 + t1v1)− b] · v1 = 0

bedeutet v1 ‖ {x : g(x) = g(x1)}.I Wie soll man die Suchrichtung v2 auswählen?I Antwort: Minimiere g über span{x1 + r1,x1 + v1}.

90

Konjugierte Gradienten

I Die Funktion (t , σ) 7→ g(x1 + tr1 + σv1) wird minimiert:

ddσ

g(x1 + tr1 + σv1) = ∇g(x1 + tr1 + σv1) · v1 =σ=ts· · · =

[A(x1 + tr1 + tsv1)− b] · v1 = (Ax1 − b) · v1︸︷︷︸−r1·v1=0

+t [A (r1 + sv1)︸︷︷︸=v2

] · v1 = 0

I Nimm v2 = r1 + s1v1 wobei

s1 = − r1 · Av1

v1 · Av1

I Bemerke: v2 · Av1 = (r1 + s1v1) · Av1 = 0,d.h. v1 und v2 sind A-orthogonal oder A-konjugiert.

Iddt

g(x1 + tr1 + ts1v1) =ddt

g(x1 + tv2) = 0 führt zu t = t2.

I Die obigen Formeln sind allgemein für jede Iteration...91

Konjugierte GradientenZusammenfassung der Formeln:

x0 sei gegeben Für k = 1,2, . . .v1 = r0 = b − Ax0 tk = rk−1 · vk/vk · Avk

xk = xk−1 + tkvk

rk · vk = 0 rk = b − Axkvk+1 · Avk = 0 sk = −rk · Avk/vk · Avk

vk+1 = rk + skvkI Vereinfachungen:

rk · vk = 0⇒ tk = rk−1 · [rk−1 + sk−1vk−1]

vk · Avk=

rk−1 · rk−1

vk · AvkAlso

rk = b − Axk = b − A(xk−1 + tkvk ) = rk−1 − tkAvk⇒ rk · rk−1 = rk−1 · rk−1 − tk (vk − sk−1vk−1) · Avk = 0⇒ rk · rk = rk · rk−1 − tk rk · Avk = −tk rk · Avk

undsk = − r

k · Avk

vk · Avk=

rk · rk/tk

rk−1 · rk−1/tk=

rk · rk

rk−1 · rk−192

Konjugierte GradientenDie Methode der Konjugierten Gradienten:

x0sei gegeben Für k = 1,2, . . .v1 = r0 = b − Ax0 tk = rk−1 · rk−1/vk · Avk

xk = xk−1 + tkvk

rk = rk−1 − tkAvksk = rk · rk/rk−1 · rk−1

vk+1 = rk + skvk

I Hausaufgabe: Schreibe einen effizienten Pseude-Kode zurImplementierung dieser Methode.

Vorausgesetzt dass A SPD ist:

Satz: Es gelten rk · r l = 0, 0 ≤ l < k und vk · Av l = 0, 1 ≤ l < k .

Satz: Mit exakter Arithmetik gilt Axk = b für ein k ≤ n.93

PräkonditionierungI Angenommen gilt C2 ≈ A aber Systeme Cy = c mit SPD

C sind leichter zu lösen als Systeme Ax = b.I Das Problem Ax = b kann so umgeschrieben werden:

(C−1AC−1)(Cx) = C−1b

I Die Methode der Konjugierten Gradienten führt schnellerzur Konvergenz, wenn sie auf das hergeleitete Systemangewendet wird, falls die SPD Matrix (C−1AC−1) besserkonditioniert ist als die Matrix A.

Def: Sei ‖ · ‖? eine induzierte Matrixnorm. Die Konditionszahleiner Matrix A bezüglich dieser Norm ist κ?(A) = ‖A‖?‖A−1‖?.

Bemerkung: 1 = ‖I‖? = ‖AA−1‖? ≤ ‖A‖?‖A−1‖? = κ?(A).Beispiel:

A =[

1 1− �1− � 1

]�∈(0,1)

σ(A) = {�,2− �}κ2(A) = (2− �)/�

94

Fehler Abschätzungen

Satz: Eine gegebene Vektornorm und die entsprechendeinduzierte Matrixnorm sei mit ‖ · ‖ bezeichnet, und κ ist dieentsprechende Konditionszahl. Dann für A ∈ Rn×n undb,x? = A−1b, x̃ ,∈ Rn gilt:

‖x? − x̃‖‖x?‖

≤ κ(A)‖b − Ax̃‖‖b‖

Beweis: Aus Ax? = b folgt

‖b‖ = ‖Ax?‖ ≤ ‖A‖‖x?‖ ⇒ 1‖x?‖

≤ ‖A‖‖b‖

Mit e = x? − x̃ und r = b − Ax̃ = Ax? − Ax̃ = Ae gilt

‖e‖ = ‖A−1Ae‖ ≤ ‖A−1‖‖Ae‖ = ‖A−1‖‖r‖

und daher‖e‖‖x?‖

≤ ‖A−1‖‖A‖ ‖r‖‖b‖

95

Fehler AbschätzungenI Ist (‖b−Axk‖ < τ · ‖b‖) ein gutes Abbruchskriterium einer

Iteration xk k→∞−→ x??I Wegen der obigen Abschätzung bleibt es offen, dass‖x? − xk‖/‖x?‖ sehr groß sein kann. Pessimistisch?

I Beispiel:

A =[

1 1− �1− � 1

]�∈(0,1)

σ(A) = {�,2− �}κ2(A) = (2− �)/�

x? =[

cos(θ)sin(θ)

]θ=π4

b = Ax? =[

cos(θ) + (1− �) sin(θ)sin(θ) + (1− �) cos(θ)

]‖x?‖2 = 1, ‖b‖22 = 1+(1−�)2+4 cos(θ) sin(θ)︸︷︷︸

2 sin(2θ)=2

(1−�) = (2−�)2

x̃ = x? +zδ, z =

1√2

[1−1

], Az = �z, ‖z‖2 = 1

‖x? − x̃‖2 =‖z‖2δ

=1δ, ‖b − Ax̃‖2 =

‖Az‖2δ

=�

δ96

Fehler Abschätzungen

I Zusammenfassung: mit δ =√�,

‖x? − x̃‖2‖x?‖2︸︷︷︸

1/δδ→0−→∞

= κ2(A)︸︷︷︸(2−δ2)/δ2δ→0−→∞

‖b − Ax̃‖2‖b‖2︸︷︷︸

δ/(2−δ2)δ→0−→0

d.h. es kann doch sein, dass das Residuum sehr kleinwird, während der Fehler groß bleibt.

I Abbruchskriterium? Zu empfehlen ist:

(‖xk+1 − xk‖ < τ · ‖xk+1‖)

wobei die Toleranz τ rein von der Genauigkeit derFließkommazahl-Darstellung gewählt werden kann.

97

Eigenwerte und EigenvektorenZuerst eine gewisse Infrastruktur:

Satz: Hat A ∈ Rn×n die verschiedenen Eigenwerte {λi}ki=1,(k ≤ n), mit den entsprechenden Eigenvektoren {x i}ki=1, dannsind diese Eigenvektoren linear unabhängig. Wenn k = n gilt,dann bilden die Eigenvektoren eine Basis für Rn.

Def: P ∈ Rn×n ist orthogonal wenn P−1 = PT gilt.Def: A,B ∈ Rn×n sind ähnlich wenn ∃S ∈ Rn×n 3 A = S−1BS .Satz: A = S−1BS ⇒ σ(A) = σ(B) undAx = λx ⇒ BSx = λSx .Satz: Ist A ∈ Rn×n symmetrisch und D = diag({λi}ni=1),{λi}ni=1 = σ(A), dann ∃P ∈ R

n×n orthogonal wobei D = PTAP.

Satz: Ist A ∈ Rn×n symmetrisch, gilt σ(A) ⊂ R.Satz: Ist A ∈ Rn×n symmetrisch, dann bilden die Eigenvektoreneine orthonormale Basis für Rn.

Satz: A ist SPD⇔ A ist symmetrisch und σ(A) ⊂ R+.98

VektoriterationZur Bestimmung eines dominanten EigenwertsI Der Plan: Für A ∈ Rn×n und einen Startvektor x0, soll Akx0

einigermaßen zu einem Eigenvektor konvergieren, derzum dominanten Eigenwert gehört.

I Angenommen hat A ∈ Rn×n die Eigenwerte {λi}ni=1 mitden entsprechenden Eigenvektoren {v i}ni=1, die eine Basisfür Rn bilden.

I Zusätzlich gelten

|λ1| > |λ2| ≥ · · · ≥ |λn| ≥ 0.

d.h. λ1 ist dominant.I Ein beliebiger Startvektor für die Vektoriteration kann so

dargestellt werden:x0 =

n∑i=1

αiv i

I Angenommen gilt α1 6= 0 und x0 ist normalisiert sodass‖x0‖∞ = 1 gilt.

99

VektoriterationI Die Vektoriteration: Für k = 1,2, . . .

yk = Axk−1

xk = yk/ykpkµk = ykpk−1

wobei{|ykpk | = ‖y

k‖∞d.h. ‖xk‖∞ = 1

I Behauptungen:

µkk→∞−→ λ1 und skxk

k→∞−→ v1

‖v1‖∞wobei s2k = 1

I Zuerst eine Darstellung von xk ,

xk =yk

ykpk=

Axk−1

ykpk=

Ayk−1

ykpk yk−1pk−1

= · · · = Ak−1y1

ykpk · · · y1p1

=Akx0∏ki=1 y

ipi

I Dann eine Darstellung für µk ,

µk =ykpk−1

xk−1pk−1(= 1)=

(Axk−1)pk−1xk−1pk−1

=(Akx0)pk−1/

∏k−1i=1 y

ipi

(Ak−1x0)pk−1/∏k−1

i=1 yipi

=(Akx0)pk−1

(Ak−1x0)pk−1100

VektoriterationI Weil λ1 dominant ist, gilt

µk =

λk1

[α1v1 +

n∑i=2

(λiλ1

)k↓0

v i]

pk−1

λk−11

[α1v1 +

n∑i=2

(λiλ1

)k−1↓0

v i]

pk−1

k→∞−→ λ1

I Wegen der obigen Darstellung von xk gilt

yk = ykpk xk = ykpk

Akx0∏ki=1 y

ipi

=λk1∏k−1

i=1 yipi︸︷︷︸

Q=1/(α1νk )

[α1v1 +

n∑i=2

αi

(λiλ1

)k↓0

v i]

Mit |v1p | = ‖v1‖∞ gilt

|ykpk ||ykp |

=‖yk‖∞|ykp |

=‖νkyk‖∞|νkykp |

k→∞−→ ‖v1‖∞|v1p |

= 1

101

Inverse VektoriterationI Also konvergieren die Vektoren xk wie behauptet:

ykpkykp

∣∣∣∣∣→±1

xk =yk

ykp=

Axk−1

(Axk−1)p=

Akx0/∏k−1

i=1 yipi

(Akx0)p/∏k−1

i=1 yipi

=Akx0

(Akx0)pk

=

λk1

[α1v1 +

n∑i=2

αi

(λiλ1

)k↓0

v i]

λk1

[α1v1 +

n∑i=2

αi

(λiλ1

)k↓0

v i]

p

k→∞−→ v1

v1p= sign(v1p )

v1

‖v1‖∞

Inverse VektoriterationI Um λk zu berechnen, wird die Vektoriteration auf

(A− qI)−1 angewendet, wobei |λk − q| < |λi − q|,∀i 6= k .I Hausaufgabe: Es gilt µk

k→∞−→ (λk − q)−1 für die Iteration,(A− qI)yk = xk−1

xk = yk/ykpkµk = ykpk−1

wobei{|ykpk | = ‖y

k‖∞d.h. ‖xk‖∞ = 1

102

Berechnung aller EigenwerteI Die Vektoriteration und die Inverse Vektoriteration sind nur

zur Berechnung bestimmter Eigenwerte geeignet.I Nun suchen wir Verfahren, die alle Eigenwerte und ihre

entsprechenden Eigenvektoren liefern.I Nach dem Abel-Ruffini Satz gibt es für n ≥ 5 keine

allgemeine endliche Formel mit Arithmetik und Wurzeln zurBestimmung der Nullstellen des charakteristischenPolynoms.

I Daher müssen Eigenwerte iterativ berechnet werden.I Der Plan vom QR-Algorithmus: A(1) = A (SPD)

Für k = 1,2 . . . A(k) = Q(k)R(k), A(k+1) = R(k)Q(k) wobei

Q(k) orthogonal ist und R(k) eine obere Dreiecksmatrix ist.I Aber diese Iteration ist viel billiger, wenn zuerst so

transformiert wird: Ã = PAPT, wobei P orthogonal und ÃSPD tridiagonal sind.

103

Die Householder MethodeI Zur Bestimmung der Transformation Ã = PAPT, A SPD,

wobei P orthogonal und Ã SPD tridiagonal sind:I Mit P = P(n−2) · · ·P(1), A(1) = A, soll gelten

A(2) = P(1)A(1)P(1) =

� � 0 · · · 0� + − − +0 | |... | |0 + − − +

I und im nächsten Schritt

A(3) = P(2)A(2)P(2) =

� � 0 · · · 0� � � 0 · · ·0 � + − +... 0 | |

0... + − +

I usw bis Ã = A(n−1).

104

Householder TransformationenDef: Sei ŵ ∈ Rn mit ŵ · ŵ = 1. Dann heisst P = I − 2ŵŵT

eine Householder Transformation.

Satz: Eine Householder Transformation ist symmetrisch undorthogonal.

PTP = P2 = (I−2ŵŵT)(I−2ŵŵT) = I−4ŵŵT+4ŵŵTŵŵT = I

Satz: Seien 0 6= u ∈ Rn und θ = 12‖u‖22. Dann ist

P = I − uTu/θ eine Householder Transformation.

ŵ = u/‖u‖2 ⇒ P = I − 2ŵTŵ

Satz: Seien x ∈ Rn und σ = ±‖x‖2. Angenommen x 6= −σê(1).Seien u = x + σê(1) und θ = 12‖u‖

22. Dann ist P = I − uTu/θ

eine Householder Transformation und es gilt Px = −σê(1).

Beweis: x 6= −σê(1) ⇒ u = x + σê(1) 6= 0. Nach dem obigenSatz ist P eine Householder Transformation. Dann gilt105

Householder Transformationenθ = 12‖u‖

22 =

12(x + σê

(1))T(x + σê(1))

= 12(xTx + 2σxTê(1) + σ2) (σ2 = ‖x‖22)

= 12(σ2 + 2σxTê(1) + σ2) = σ2 + σx1

Daher gilt(I − uuT/θ)x = x − uuTx/θ

= x − (x + σê(1))(x + σê(1)

)Txσ2 + σx1

= x − (x + σê(1))(xTx + σx1)σ2 + σx1

= −σê(1)

Bemerkung: Um Auslöschung zu vermeiden, soll σ mit demgleichen Vorzeichen wie x1 ausgewählt werden. P wird auchnicht verändert, wenn x zur Sicherheit skaliert wird.

Algorithmus: Hausaufgabe: θ = 12‖u‖22

v = x/‖x‖∞ θ = ρu1ρ = s(v1)‖v‖2 σ = ρ‖x‖∞u = v + ρê(1) P = I − uTu/θ

s(t) ={

1, t ≥ 0−1, sonst

Px = −σê(1)106

Transformation auf Tridiagonale FormI Zur Bestimmung von P(1) sei

a(1) =

a(1)11a(1)21...

a(1)n1

=

a(1)11

ã(1)

P̃ = I − ũTũ/θ, ũ ∈ Rn−1

wobei P̃ã(1) = −σê(1) ∈ Rn−1

I Dann

P = I − uTu/θ, u =[

0ũ

]∈ Rn

erfüllt

Pa(1) =

1 0 · · · 00 + − +... | P̃ |0 + − +

a(1)11

ã(1)

=

a(1)11−σ0...

107

Transformation auf Tridiagonale FormI Zur Bestimmung von P(2) sei

a =

a(2)12a(2)22...

a(2)n2

=

a(2)12a(2)22

ã(2)

P̃ = I − ũTũ/θ, ũ ∈ Rn−2

wobei P̃ã(2) = −σê(1) ∈ Rn−2

I Dann

P = I − uTu/θ, u =

00ũ

∈ Rnerfüllt

Pa(2) =

1 0 0 · · · 00 1 0 · · · 00 0 + − +...

... | P̃ |0 0 + − +

a(2)12a(2)22

ã(2)

=

a(2)12a(2)22−σ0...

108

Hessenberg MatrizenI Dann sehen die Matrizen P(i) · · ·P(1)A so aus:

P(1)A =

� + − − +� | |0 | |... | |0 + − − +

P(2)P(1)A =

� � · · · · · · �� + − +0 � | |... 0 | |

0... + − +

I Und P(n−2) · · ·P(1)A ist eine obere Hessenberg Matrix:

P(n−2) · · ·P(1)A =

� · · · · · · · · · �

�. . .

...

0. . .

. . ....

.... . .

. . .. . . �

0 · · · 0 � �

109

Transformation auf Tridiagonale FormI Wenn A(1) = A SPD ist, gilt A(2) = P(1)A(1)P(1) =

1 0 · · · 00 + − +... | P̃ |0 + − +

a(1)11 ã(1)T

+ − +ã(1) | Ã |

+ − +

1 0 · · · 00 + − +... | P̃ |0 + − +

=

a(1)11 −σ 0 · · ·−σ + − +0 | P̃ÃP̃ |... + − +

I P(2) wird von A(2) = P(1)A(1)P(1) wie vorher bestimmt, und

ähnlich ist A(3) = P(2)A(2)P(2) tridiagonal.I Dann mit

P = P(n−2) · · ·P(1)

ist PAPT tridiagonal.110

Pseudo-Kode zur Transformation auf Hessenbergoder Tridiagonale Form

H = A, Q = I, ASym = (A = A’)for k=1,...,n-2

h i = H k+i,k, i=1,...,n-k(theta,sigma,u) = Householder(h) % (I-uu’/theta)h=-sigma e1for l=k,...,n % H

QR AlgorithmusI Für die Iteration, A(k) = Q(k)R(k), A(k+1) = R(k)Q(k), gilt

A(k+1) = R(k)Q(k) = Q(k)TA(k)Q(k) = · · · =Q(k)T · · ·Q(1)T︸︷︷︸

QT

A Q(1) · · ·Q(k)︸︷︷︸Q

also sind die Matrizen {A(k)} alle ähnlich zu A.I Unter gewissen Bedingungen gilt

A(k) k−→∞−→ Λ = diag(σ(A)).I Der erste Bedarf ist, ein Algorithmus zur QR-Zerlegung.I Householder Transformationen funktionieren: Für Q(1)

a(1) =

a(1)11a(1)21...

a(1)n1

Q(1) = I − uTu/θ, u ∈ Rn

wobei Q(1)a(1) = −σê(1) ∈ Rn

112

QR AlgorithmusI Für Q(2)

a(2) =

a(2)12a(2)22...

a(2)n2

=

a(2)12

ã(2)

Q̃(2) = I − ũTũ/θ, ũ ∈ Rn−1

wobei Q̃(2)ã(2) = −σê(1) ∈ Rn−1

Q(2) = I − uTu/θ, u =[

0ũ

]I Es gelten

Q(1)A=

−σ1 + − − +

0 | |... | |0 + − − +

Q(2)Q(1)A=−σ1 � + − +

0 −σ2 | |... 0 | |

0... + − +

I Nach n − 1 solchen Transformationen ergibt sich eine

obere Dreiecksmatrix,

QTA = Q(n−1) · · ·Q(1)A = R113

Givens TransformationenI Aber wenn A schon tridiagonal ist, kann Q billiger

konstruiert werden:

A =

a1 b1 0 · · · 0

b1 a2 b2. . .

...

0. . .

. . .. . . 0

.... . . bn−2 an−1 bn−1

0 · · · 0 bn−1 an

Def: Eine Givens Transformation G = Gij(θ) hat die Form

GGT = I, G =

1. . . 1qii qij

1. . . 1qji qjj

1. . . 1

qii = cos(θ) = qjjqij = sin(θ) = −qjifür irgendwelche

i , j und θ

114

Givens TransformationenI Für x1 = a1, y1 = b1, soll θ so ausgewählt werden, dass

cos(θ) sin(θ)− sin(θ) cos(θ)

11. . . 1

x1 y1b1 a2 b2

. . .. . .

. . .

bn−1 an

=

g1 d1 e10 x2 y2

b2 a3 b3. . .

. . .. . .

bn−1 an

I Mit s1 = sin(θ1), c1 = cos(θ1),

−s1x1 + c1b1 = 0s21 + c

21 = 1

} s1 = b1/√b21 + x21c1 = x1/

√b21 + x

21

I Hausaufgabe: Für g1 > 0 ist Q eindeutig bestimmt.115

Die QR-ZerlegungI Initialisierung: R = A, Q = I.I Für k = 1, . . . ,n − 1,

Q = Gk ,k+1(θk )Q, R = Gk ,k+1(θk )R

wobei sin(θk ) = bk/√

b2k + x2k , cos(θk ) = xk/

√b2k + x

2k .

I Es gilt1. . . 1

cos(θk ) sin(θk )− sin(θk ) cos(θk )

1. . . 1

. . .. . .

. . .

gk−1 dk−1 ek−10 xk yk 0

bk ak+1 bk+1. . .

. . .. . .

. . .. . .

. . .

gk dk ek0 xk+1 yk+1 0

bk+1 ak+2 bk+2. . .

. . .. . .

I Am Ende setze Q = QT und es gilt A = QR. Eindeutig?

116

QR AlgorithmusSatz: Sei A ∈ Rn×n mit Rang(A) = n. Dann gibt es eineZerlegung A = QR, wobei Q orthogonal ist und R eine obereDreiecksmatrix ist. Weiters ist die Zerlegung eindeutig, wenndie diagonalen Einträge von R positiv sind.

Bemerkung: Für die Konvergenz des QR-Algorithmus zurBestimmung der Eigenräume von A ist diese eindeutigeZerlegung nicht notwendig.

I Nun haben wir eine QR-Zerlegung A(k) = Q(k)R(k).I Zu zeigen ist, dass

A(k+1) = R(k)Q(k) = Q(k)TA(k)Q(k)

wieder tridiagonal ist.I Zuerst wird gezeigt, Q(k) ist obere Hessenberg.

117

Tridiagonale Form Bleibt ErhaltenI Bemerke:

1c2 s2−s2 c2

1. . . 1

c1 s1−s1 c1

11. . . 1

=

c1 s1� � ��

11. . . 1

I Induktiv:

1. . . 1

ck sk−sk ck

1. . . 1

� �

�. . .

. . ....

. . .. . .

. . ....

. . .. . . �

� · · · · · · � �1

1. . . 1

=

� �

�. . .

. . ....

. . .. . .

. . ....

. . .. . . �

.... . . � �

� · · · · · · · · · � �1. . . 1

I Also ist Gk ,k+1(θk ) · · ·G1,2(θ1) untere Hessenberg.118

Tridiagonale Form Bleibt ErhaltenI Daher haben wir folgende Struktur für A(k+1) = R(k)Q(k) =

� � �. . .

. . .. . .

. . .. . . �. . . �

�

� � · · · · · · �

�. . .

. . ....

. . .. . .

. . ....

. . .. . . �

� �

und A(k+1) ist

eine obereHessenberg

Matrix. =

� � · · · · · · �

�. . .

. . ....

. . .. . .

. . ....

. . .. . . �

� �

I Aber A(k+1) = Q(k)TA(k)Q(k) ist symmetrisch und

deswegen tridiagonal!119

Konvergenz des QR-AlgorithmusSatz: Sei A ∈ Rn×n symmetrisch mit Eigenwerten {λi}ni=1, dieerfüllen |λ1| > |λ2| > · · · > |λn| > 0. Seien {v i}ni=1 dieentsprechenden Eigenvektoren, die eine orthonormale Basisbilden. Für die Matrix Q = [v1, . . . ,vn] soll es eineLU-Zerlegung QT = LU geben, ohne dass QT mit einerPermutationsmatrix multipliziert werden muss. Dann mitA(1) = A und einer (beliebigen) QR-Zerlegung Ak = QkRk ,konvergiert die Folge A(k+1) = R(k)Q(k) zur diagonalen MatrixΛ = diag({λi}ni=1).

Bemerkung: Mit A(k) = {a(k)ij } gilt

|a(k)i+1,i | = O

(∣∣∣∣λi+1λi∣∣∣∣k)

oder |a(k)i+1,i | = O

(∣∣∣∣λi+1 − qλi − q∣∣∣∣k)

wenn man eine Verschiebung macht,

A(k) − qI = Q(k)R(k), A(k+1) = R(k)Q(k) + qI120

Pseudo-Kode für den QR-Algorithmusa i = H ii, i=1,...,n, b i = H i,i+1 = H i+1,i, i=1,...,n-1for k=1,...,kmax

x 1 = a 1, y 1 = b 1for j=1,...,n-1 % compute R, Q (storing c’s & s’s)

g j =√x jˆ2 + b jˆ2

c j = x j / g j, s j = b j / g jd j = c j * y j + s j * a j+1x j+1 = c j * a j+1 - s j * y jif j 6= (n-1)

e j = s j * b j+1, y j+1 = c j * b j+1end

endg n = x na 1 = s 1 * d 1 + c 1 * g 1, b 1 = s 1 * g 2for j=2,...,n-1 % compute R * Q, overwriting H

a j = s j * d j + c j * c j-1 * g jb j = s j * g j+1

enda n = c n-1 * g nif (||b|| < tol * ||a||) then break

end% a i, i=1,...,n ≈ Eigenwerte von H121

AusgleichsproblemeBeispiel (Lineare Regression): Gegeben sind {(xi , yi)}m−1i=0 .Was ist die beste Gerade durch diese Daten?

I Es soll minimiert werden:

E(a,b) =m−1∑i=0

[yi − (a + bxi)]2 V = V (x) = [e,x ]

= ‖y − Vc‖22 c = 〈a,b〉T (ei = 1)

= cTV TVc − 2cTV Ty + yTy

I E(c?) = min⇒ 0 = ∇cE(c?) = 2V TVc? − 2V TyI Hausaufgabe: {xi}m−1i=0 verschieden⇒ V

TV ist SPD.I Wenn c? die normalen Gleichungen löst,

V TVc = V Ty

folgt E(c?) = min.122

Lineare und Polynome RegressionI Hausaufgabe: Zeige, für lineare Regression ist die Lösung

explizit so gegeben:

a? = ȳ−b?x̄, b? = x̄y − x̄ ȳx̄2 − x̄2

wobei z.B. x̄y =1m

m−1∑i=0

xiyi

I Polynome Regression: (V=Vandermonde Matrix)

E(c) =m−1∑i=0

[yi − P(xi ; c)]2, P(x ; c) =n−1∑k=0

ckxk

= ‖y − Vc‖22 V = [e,x , . . . ,xn−1]

= cTV TVc − 2cTV Ty + yTy c = 〈c0, . . . , cn−1〉T

I Lösung gegeben durch die normalen Gleichungen,

V TVc = V Ty

aber V TV kann für hohes n schlecht konditioniert sein.123

Singulärwert ZerlegungI Ein alternativer Ansatz ist durch die Singulärwert

Zerlegung gegeben.

Satz (SWZ): Für A ∈ Rm×n ∃ orthogonale Matrizen U ∈ Rm×m,V ∈ Rn×n sodass UTAV = Σ = diag{σ1, . . . , σp} ∈ Rm×n wobeip = min{m,n} und σ1 ≥ σ2 ≥ · · · ≥ σp ≥ 0.

I Also UTAV = Σ⇒ A = UUTAVV T = UΣV T.

Def: A = UΣV T ist die Singulärwert Zerlegung (SWZ) von A,und {σi}pi=1 sind die Singulärwerte von A.

Satz: Sei A = UΣV T die SWZ für A ∈ Rm×n, wobei m > n undσ1 ≥ · · · ≥ σr > σr+1 = · · · = σn = 0 gelten. Dann wird‖Ax − b‖2 durch x? = V Σ†UTb ∈ Rn minimiert, wobeiΣ† = diag{σ−11 , . . . , σ

−1r ,0, . . . ,0} ∈ Rn×m. Wenn x̃ ∈ Rn auch

minimierend ist und x̃ 6= x? gilt, dann gilt ‖x?‖2 < ‖x̃‖2.124

Eigenschaften der SWZ

Satz (Eigenschaften der SWZ): Sei A = UΣV T die SWZ fürA ∈ Rm×n mit p = min{m,n}. Es gelten:I U = {u1, . . . ,um}, V = {v1, . . . ,vn} ⇒

Av i = σiui , ATui = σiv i , i = 1, . . . ,p.I σ1 ≥ · · · ≥ σr > σr+1 = · · · = σn = 0⇒N (A) = span{v i}ni=r+1, R(A) = span{u

i}ri=1.I ‖A‖2 = σ1.I ‖A‖2F = σ21 + · · ·+ σ2p.

I κ?(A) = max‖x‖?=1

‖Ax‖?/ min‖x‖?=1

‖Ax‖?m≥n⇒ κ2(A) = σ1/σn.

I ATA = V ΣTΣV T und AAT = UΣΣTUT.I {σ2i }

pi=1 sind die eindeutigen Nullstellen von det(λI − A

TA)für p = n oder die eindeutigen Nullstellen vondet(λI − AAT) für p = m.

125

PseudoinverseDef: Sei A = UΣV T die SWZ für A ∈ Rm×n. Dann ist istA† = V Σ†UT die Pseudo-Inverse von A.

Satz: Für A ∈ Rm×n, N (A) = 0⇒ A† = (ATA)−1AT.

Beweis: Es gilt xTATAx = ‖Ax‖22 ≥ 0 und zusätzlich folgt‖Ax‖22 = 0⇔ x = 0 sowohl m ≥ n aus N (A) = 0. Daher istATA SPD. Aus ATA = (UΣV T)T(UΣV T) = V ΣTUTUΣV T =V ΣTΣV T folgt, dass die positiven Eigenwerte durch 0 < ΣTΣ =diag{σ21, . . . , σ2n} ∈ Rn×n gegeben sind. Es folgt (ATA)−1AT =V [ΣTΣ]−1V TV ΣTUT = V [ΣTΣ]−1ΣTUT = V Σ†UT = A†.

Satz: Für A ∈ Rm×n und orthogonale Matrizen P ∈ Rm×m undQ ∈ Rn×n haben A und PAQ die gleichen Singulärwerte.

Beweis: Für m > n sind die Eigenwerte {σ2i }ni=1 von A

TA =V ΣTΣV T die Eigenwerte von (PAQ)T(PAQ) = (QTV )ΣTΣ(V TQ).Für n > m sind die Eigenwerte {σ2i }

mi=1 von AA

T = UΣΣTUT dieEigenwerte von (PAQ)T(PAQ) = (QTU)ΣΣT(UTQ).126

Konditionierung und die SWZI Für A ∈ Rm×n, m > n, mit SWZ, A = UTΣV ,

Σ =

(Σ̂0

)∈ Rm×n, Σ̂ = diag{σi}ni=1, σi > 0

sei die Projektion P = ΣΣ† von Rm in Rn definiert.I Für ein gegebenes b sei die Zerlegung definiert:

Ub = 〈c,d〉T, c ∈ Rr , d ∈ Rm−rPUb = 〈c,0〉T, (I − P)Ub = 〈0,d〉T

I Das Ausgleichsproblem kann so umformuliert werden:‖Ax − b‖22 = ‖UTΣVx − b‖22 = ‖ΣVx − Ub‖22= ‖P(ΣVx − Ub)‖22 + ‖(I − P)(ΣVx − Ub)‖22

+2〈P(· · · ), (I − P)(· · · )〉=0 = ‖Σ̂Vx − c‖22 + ‖d‖22I Da κ2(Σ̂) = κ2(Σ̂V ) gilt, erfüllen die Lösungen der

Systeme Σ̂Vx = c und Σ̂V x̃ = c̃ die Fehlerabschätzung‖x − x̃‖2‖x‖2

≤ κ2(Σ̂)‖c − c̃‖2‖c‖2

I Vergleiche mit κ2(ATA) = κ2(Σ̂)2 für die normalenGleichungen.

127

Berechnung der SWZBerechnung der SWZ (zur Lösung eines Ausgleichsproblems)I Da die Singulärwerte durch die Eigenwerte von ATA

(m > n) gegeben sind, könnten wir den QR-Algorithmusauf ATA anwenden, um ΣTΣ zu bestimmen.

I Der folgende Algorithmus mit bidiagonalen Matrizen isteffizienter.

I Der Plan:I Vereinfachung PAQ = 〈BT,0〉T, A ∈ Rm×n, P ∈ Rm×m und

Q ∈ Rn×n orthogonal und B ∈ Rn×n obere bidiagonal.I Mit B(1) = B ist B(l)TB(l) tridiagonal. (warum?)I Setze B(l)TB(l) = Q(l)R(l). (QR-Zerlegung)I Dann ist B(l)Q(l) untere bidiagonal. (warum?)I Berechne P(l) sodass B(l+1) = P(l)B(l)Q(l)

obere bidiagonal ist. (wie?)I Zeige, es gilt B(l+1)TB(l+1) = R(l)Q(l)!I Im Endeffekt wird ein Schritt eines QR-Algorithmus effizient

durchgeführt, und es folgt |B(l)| → diag{σi}ni=1, l →∞.

128

Vereinfachung auf Bidiagonale FormI Für die Vereinfachung PAQ = (B,0)T setze

A = A(1) = (d (1)|C(1)), d (1) ∈ Rm×1, C(1) ∈ Rm×(n−1)I Für k = 2, . . . ,n − 1,

A(k) = P(k−1)A(k−1)Q(k−1) =

(B(k) 0

a(k)0

0 d (k) C(k)

)wobei

k > 2

B(k) ∈ R(k−1)×(k−1) a(k) ∈ R1

C(k) ∈ R(m−k+1)×(n−k) d (k) ∈ R(m−k+1)×1

I Im nächsten Schritt,

P(k) =(

I(k−1) 00 P̃(k)

)Q(k) =

(I(k) 00 Q̃(k)

)wobei durch Householder Transformationen

P̃(k)d (k) = −σê(1) ∈ Rm−k+1 c(k) = 1.Zeile von P̃(k)C(k)

Q̃(k)Tc(k)T = −τ ê(1) ∈ Rn−k a(k+1) = −τ129

Pseudo-Kode zur Transformation auf BidiagonaleForm

d i = A i,1, C i,j = A i,j+1, i=1,...,m, j=1,...,n-1P = I (m×m), Q = I (n×n)for k=1,...,n

% d is (m-k+1)×1(u,sigma,theta) = Householder(d) % (I-uu’/theta)d=-sigma e1B k,k = -sigma % B growsif (k>1) then B k-1,k = afor l=1,...,n-k % C̃

Pseudo-Kode zur Transformation auf BidiagonaleForm

(u,tau,theta) = Householder(c) % (I-uu’/theta)c=-tau e1a = -taufor l=2,...,m-k+1 % C̃

Diagonalisierende Iteration

I Bei k = n − 1,A(n) = P(n−1)A(n−1)

und Q(n−1) ist nicht notwendig, weil C(n) nur eine einzigeSpalte hat, und c(n−1) ∈ R1 ist schon a(n).

I Bei k = n,A(n+1) = P(n)A(n)

und Q(n) ist nicht notwendig, weil C(n) nicht existiert. Dad (n) ∈ R(m−n+1)×1, ist P(n) aber notwendig.

I Dann ist B = B(n+1) bidiagonal.I Für die Diagonalisierung von B setze B(1) = B.I Sei die obere Hessenberg Matrix

Q(l) = G(1,2)(θ1) · · ·G(n−1,n)(θn−1)

durch Givens Transformationen definiert, damit

132

Diagonalisierende Iteration

B(l)Q(l) =

� �

. . .. . .

. . . ��

Q(l) =� ↙0

�. . .

. . .

. . .. . . ↙0� �

eine untere bidiagonale Matrix ist.

I Dann durch Multiplikation links mit B(l)T

R(l)T = B(l)T(B(l)Q(l)) =

�

�. . .

. . .. . .

� �

�

�. . .

. . .. . .

� �

ist R(l) eine obere tridiagonale Matrix.

133

Folge der Bidiagonalen MatrizenI Daher ist

B(l)TB(l) = Q(l)R(l)

eine QR Zerlegung von der tridiagonale Matrix B(l)TB(l),� �

�. . .

. . .

. . .. . . �

� �

= B(l)TB(l) = Q(l)R(l) =

� � · · · · · · �

�. . .

. . ....

. . .. . .

. . ....

. . .. . . �

� �

� � �. . .

. . .. . .

. . .. . . �. . . �

�

134

Folge der Bidiagonalen MatrizenI Sei die untere Hessenberg Matrix

P(l) = G(n−1,n)(θn−1) · · ·G(1,2)(θ1)

durch Givens Transformationen definiert, damit

B(l+1) =P(l)(B(l)Q(l))=P(l)

�

�. . .

. . .. . .

� �

=� �

0↗. . .

. . .

. . .. . . �0↗ �

eine obere bidiagonale Matrix ist.

I Schliesslich,

B(l+1)TB(l+1) = (P(l)B(l)Q(l))T (P(l)B(l)Q(l))= (Q(l)TB(l)TP(l)T) (P(l)B(l)Q(l))= Q(l)T(B(l)TB(l))Q(l) = Q(l)T(Q(l)R(l))Q(l)

= R(l)Q(l)

und ein Schritt des QR Verfahrens für BTB ist vollständig.135

Pseudo-Kode für den QR Algorithmus zurSingulärwert-Zerlegung

Q = I, P = I, Bk = Bfor k=1,...,kmax

a i = Bk i,i, i=1,...,nb i = Bk i,i+1, i=1,...,n-1x 1 = a 1for j=1...,(n-1) % compute Qk by storing c’s & s’s


c j = x j / g js j = -b j / g jd j = -s j * a j+1x j+1 = c j * a j+1

endg n = x n % Bk*Qk has (lower) diags {d,g}

136


Qk = Ifor j=1,...,(n-1)

q1 i = Qk i,j * c j - Qk i,j+1 * s j, i=1,...,nq2 i = Qk i,j+1 * c j + Qk i,j * s j, i=1,...,nQk i,j = q1 i, Qk i,j+1 = q2 i, i=1,...,n

end % (Bk*Qk)’*Bk = RkQ = Q*Qka i = (Bk*Qk) i,i = g i, i=1,...,nb i = (Bk*Qk) i+1,i = d i, i=1,...,n-1x 1 = a 1for j=1...,(n-1)


c j = x j / g js j = b j / g jd j = s j * a j+1x j+1 = c j * a j+1

endg n = x n % Tk has (upper) diags {g,d}

137


Pk = Ifor i=1,...,(n-1)

p1 j = Pk i,j * c i + Pk i+1,j * s i, j=1,...,np2 j = Pk i+1,j * c i - Pk i,j * s i, j=1,...,nPk i,j = p1 j, Pk i+1,j = p2 j, j=1,...,n

end % Pk*(Bk*Qk) = TkP = Pk*PBk i,i = g i, i=1,...,nBk i,i+1 = d i, i=1,...,(n-1) % Bk’*Bk = Rk*Qkif (||d|| < tol * ||g||) then break

end% |g i|, i=1,...,n ≈ Singulärwerte von BSn i,i = sign(g i), i=1,...,nS i,i = |g i|, i=1,...,nP=Sn*P% S = P*B*Q

138

InterpolationI Für Polynome Regression haben wir das

Ausgleichsproblem gehabt:

E(c) =m−1∑i=0

[yi − P(xi ; c)]2, P(x ; c) =n−1∑k=0

ckxk

mit m > n. Nun n = m.I Die explizite Lösung ist gegeben durch die Lagrange

interpolierenden Polynome.I Beispiel: Gegeben sind (x0, y0) und (x1, y1). Die Lagrange

interpolierenden Polynome sind:

L0(x) =x − x1x0 − x1

L1(x) =x − x0x1 − x0

Bemerke:L0(x0) = 1 L1(x0) = 0L0(x1) = 0 L1(x1) = 1

}oder Li(xj) = δij

und

P(x) =1∑

i=0

yiLi(x) erfüllt P(xj) =1∑

i=0

yiLi(xj) =1∑

i=0

yiδij = yj139

Lagrange Interpolierende PolynomeI Beispiel: Gegeben sind (x0, y0), (x1, y1) und (x2, y2). Die

Lagrange interpolierenden Polynome sind:

L0(x) =(x−x1)(x−x2)

(x0−x1)(x0−x2)L1(x) =

(x−x0)(x−x2)(x1−x0)(x1−x2)

L2(x) =(x−x0)(x−x1)

(x2−x0)(x2−x1)und

P(x) =2∑

i=0

yiLi(x) erfüllt P(xj) =2∑

i=0

yiLi(xj) =2∑

i=0

yiδij = yj

I Im allgemeinen werden diese Polynome so definiert:

LN,i(x) =N∏

i 6=j=0

x − xjxi − xj

PN(x) =N∑

i=0

yiLN,i(x)

PN(xj) =N∑

i=0

yiLN,i(xj) =N∑

i=0

yiδij = yj

I Was ist die Qualität dieser Interpolation?140

Genauigkeit globaler InterpolationI Beispiel: Die Funktion f (x) = x2/(1 + 20x2) wird mit

n = 10 oder m = 100 Punkten folgendermasseninterpoliert:I P(xi ) = f (xi ), xi = −1 + 2i/n, i = 0, . . . ,nI Q(tj ) = f (tj ), tj = cos(π(j + 1/2)/(N + 1)), j = 0, . . . ,nI∑m

k=0 |R(yk )− f (yk )|2 = min, yk = −1 + 2k/m, k = 0, . . . ,m

(Gleichmäßig)

(Tchebychev)

(Ausgleich)

I Die Ergebnisse sehen so aus:

wobei die theoretischen Fehler rechts gezeigt werden.141

Genauigkeit globaler InterpolationSatz: Wenn {xi}Ni=0 ⊂ [a,b] und f ∈ C

N+1([a,b]) dann∀x ∈ [a,b], ∃ξ(x) ∈ [a,b] 3

f (x) = PN(x) +f (N+1)(ξ(x))

(n + 1)!ψN+1(x ; x0, . . . , xN)

wobei

PN(x) =N∑

i=0

f (xi)LN,i(x) ψN+1(x ; x0, . . . , xN) =N∏

i=0

(x − xi)

Insbesondere gilt

|f (x)− PN(x)| ≤ M|ψN+1(x ; x0, . . . , xN)|wenn |f (N+1)(ξ)| ≤ M(N + 1)!, ∀ξ ∈ [a,b].

I Wie kann man |ψN+1(x ; x0, . . . , xN)| minimieren?I Diese Funktion ist oben grafisch dargestellt für die Fälle:

I Gleichmäßig verteilte Stützstellen undI Tchebychev Stützstellen.

I Das zweite Ergebnis ist klar besser, aber man kann nichtimmer die Stützstellen auswählen.

142

Iterierte InterpolationDef: Gegeben sind die Daten {(xi , fi)}Ni=0. Mit verschiedenen{mj}kj=1, 0 ≤ mj ≤ N, wird das Polynom hier mit Pm1,...,mkbezeichnet, das die Daten {(xmj , fmj )}kj=1 interpoliert.

I Beispiel: Für xi = i , i = 0,1,2,

P0,2(x) =(x − x2)(x0 − x2)

f0 +(x − x0)(x2 − x0)

f2

Satz: Gegeben sind die Daten {(xi , fi)}Ni=0, P0,...,j−1,j+1,...,N undP0,...,i−1,i+1,...,N . Dann gilt

P0,...,N(x) =(x − xj)P0,...,j−1,j+1,...,N(x)− (x − xi)P0,...,i−1,i+1,...,N(x)

xi − xj

I Dieser Satz führt zu einem Algorithmus zur Berechnungdes Polynoms P0,...,N .

143

Iterierte Interpolation

I Beispiel: Gegeben seien {(xi , fi)}2i=0 und eineAuswertungsstelle x für P012. Man bekommt P012(x) ausder Tabelle:

x0 P0(x) = f0x1 P1(x) = f1 P01(x) =

(x−x1)P0−(x−x0)P1x0−x1

x2 P2(x) = f2 P12(x) =(x−x2)P1−(x−x1)P2

x1−x2 P012(x) =(x−x2)P01−(x−x0)P12

x0−x2

I Wenn ein neuer Datenpunkt verfügbar wird, dann kann dieTabelle einfach ergänzt werden:

x0 P0(x) = f0

x1 P1(x) = f1 P01(x) =(x−x1)P0−(x−x0)P1

x0−x1

x2 P2(x) = f2 P12(x) =(x−x2)P1−(x−x1)P2

x1−x2P012(x) =

(x−x2)P01−(x−x0)P12x0−x2

x3 P3(x) = f3 P23(x) =(x−x3)P2−(x−x2)P3

x2−x3P123(x) =

(x−x3)P12−(x−x1)P23x1−x3

P0123(x) =(x−x3)P012−(x−x0)P123

x0−x3

144

Nevilles AlgorithmusAlgorithmus (Neville): Eingaben sind Daten {(xi , fi)}Ni=0 unddie Auswertungsstelle x , Ausgabe ist der Wert P0,...,N(x).

Qi,0 = fi , i = 0, . . . ,Nfor i = 1, . . . ,N

for j = 1, . . . , i

Qij =(x − xi−j)Qi,j−1 − (x − xi)Qi−1,j−1

xi − xi−jend

endreturn QNN

Bildlich gesehen:

x0 Q00 Qij = Pi−j,...,i...

... Q11 QNN = PN−N,...,N...

.... . .

xN QN0 QN1 · · · QNN145

Dividierte DifferenzenDef: PN(x) bezeichnet hier das Polynom, das die Daten{(xi , fi)}Ni=0 interpoliert.I Nun wollen wir nicht nur den Wert PN(x) an einer

Auswertungsstelle x berechnen, sondern auchKoeffizienten eines Polynoms bestimmen.

I PN(x) kann so dargestellt werden:PN(x) = a0 + a1(x−x0) + · · ·+ aN(x−x0)×· · ·× (x−xN−1)

I Die Koeffizienten {ai}Ni=0 finden wir so:f0 = PN(x0) = a0 =: f [x0]f1 = PN(x1) = a0 + a1(x1 − x0) = f0 + a1(x1 − x0)

⇒ a1 =f1 − f0x1 − x0

=: f [x0, x1]

I Die restlichen Koeffizienten werden induktiv durch ähnlicheDifferenzen gegeben:

ak = f [x0, . . .

Grundlagen der Numerischen MathematikIWenn das Integral R ju 0j2dx in J mit R ju 0jdx ersetzt wird,...

Documents

Transcript of Grundlagen der Numerischen MathematikIWenn das Integral R ju 0j2dx in J mit R ju 0jdx ersetzt wird,...