Numerische Mathematik 1 - uni-due.de · Kapitel 1 Einleitung Die Aufgabe der Numerik ist die...

Prof. Dr. Irwin Yousept Universität Duisburg-Essen

Numerische Mathematik 1Typeset und Layout: Roman HändlerFassung vom 4. März 2017

Inhaltsverzeichnis

1 Einleitung 1

2 Grundlagen 32.1 Zahlendarstellung und Rundungsfehler . . . . . . . . . . . . . . . . . . . 32.2 Vektor- und Matrixnormen . . . . . . . . . . . . . . . . . . . . . . . . . . 62.3 Kondition und Fehlerabschätzung . . . . . . . . . . . . . . . . . . . . . . 14

3 Direkte Verfahren 253.1 Lineare Gleichungssysteme einfacher Strukturen . . . . . . . . . . . . . . 253.2 LR-Zerlegung ohne Pivotisierung . . . . . . . . . . . . . . . . . . . . . . . 283.3 LR-Zerlegung mit Pivotisierung . . . . . . . . . . . . . . . . . . . . . . . 393.4 Cholesky-Zerlegung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

4 Lineares Ausgleichsproblem 534.1 QR-Zerlegung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 584.2 Gram-Schmidt-Orthogonalisierung . . . . . . . . . . . . . . . . . . . . . . 624.3 Householder-Spiegelungen . . . . . . . . . . . . . . . . . . . . . . . . . . 66

5 CG-Verfahren 755.1 Herleitung des CG-Verfahrens . . . . . . . . . . . . . . . . . . . . . . . . 755.2 Charakterisierung des CG-Verfahrens . . . . . . . . . . . . . . . . . . . . 825.3 Konvergenzanalyse des CG-Verfahrens . . . . . . . . . . . . . . . . . . . 845.4 Das präkonditionierte CG-Verfahren . . . . . . . . . . . . . . . . . . . . . 92

6 Nichtlineares Gleichungssystem 956.1 Differentialrechnung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 956.2 Newton-Verfahren . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1026.3 Lokale Konvergenz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1046.4 Das vereinfachte Newton-Verfahren . . . . . . . . . . . . . . . . . . . . . 109

7 Interpolation 1137.1 Polynominterpolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1147.2 Hermite-Interpolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1217.3 Splines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122

iii

8 Numerische Integration 1298.1 Newton-Cotes-Formel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1298.2 Fehler von allgemeinen Quadraturformeln . . . . . . . . . . . . . . . . . 133

Kapitel 1Einleitung

Die Aufgabe der Numerik ist die Analyse und Konstruktion von Lösungsverfahrenmathematischer Aufgaben. Das folgende Diagramm gibt einen kurzen Einblick indie Arbeitsweise der Numerik und zeigt den Bezug zu anderen wissenschaftlichenBereichen auf.

Problemstellungenaus Natur- oderWirtschaftswissenschaften

Praktische Informatik(Computer)

Analysis und lineare AlgebraFunktionalanalysis

Numerik: Konstruktion vonAlgorithmen, Fehler- undKonvergenzanalyse

Daten

LösungenProblem

Theo

rie

Theorie

Problem

Problem

Alg

orit

hmen

1

Kapitel 2Grundlagen

2.1 Zahlendarstellung und Rundungsfehler

Von einem Rechner bzw. Taschenrechner können nur endlich viele verschiedene Zahlendargestellt und gespeichert werden, diese heißen Maschinenzahlen. Die Menge allerMaschinenzahlen wird mit M bezeichnet, also

M = Menge aller Maschinenzahlen.

Zur Speicherung von Maschinenzahlen verwenden die Rechner die Gleitpunktdarstel-lung:

± a0.a1a2 . . . at−1 · be,

wobei

t ∈N : Mantissenlänge,ai ∈N∪ {0} : 0 ≤ ai ≤ b− 1, a0 6= 0b ∈N : die Basis des Zahlensystems,e ∈ Z : der Exponent mit e ∈ [−n, n] mit n ∈N die Exponentenlänge.

Eine Ausnahme ist die Zahl Null, die als Maschinenzahl angesehen wird und für diealle ai = 0 sind.

Bemerkung 2.1. Im Falle b = 2 spricht man vom Binärsystem. Im Falle b = 10 sprichtman vom Dezimalsystem.

Beispiel 2.2. Betrachte b = 2, t = 3, e ∈ [−4, 4]. Es gilt

Zahl Gleitpunktdarstellung auf dem Rechner4 +1.00 · 22

−0.0625 −1.00 · 2−4

3.5 +1.11 · 21

24 +1.10 · 24

3


Beispiel 2.3. Betrachte b = 10, t = 6, e ∈ [−9, 9]. Es gilt

Zahl Gleitpunktdarstellung auf dem Rechner3 +3.00000 · 100

−26.4 −2.64000 · 101

0.005 +5.00000 · 10−3

1234567 +1.23457 · 106 (Nicht exakt, wird also gerundet)13 = 0.33 . . . +3.33333 · 10−1(gerundet)

Ist x ∈ R eine beliebige reelle Zahl, so wird x auf dem Rechner durch die nächstgele-gene Maschinenzahl ersetzt (gerundet, bei Doppeldeutigkeit wird die betragsmäßiggrößere Zahl gewählt), die wir mit

rd : R→ M, R 3 x 7→ rd(x) ∈ M

bezeichnen. Die Rundung rd : R→ M erfüllt

1. x ∈ M ⇒ rd(x) ∈ M und

2. |x− rd(x)| = miny∈M|x− y| ist die Rundung zur nächstgelegenen Maschinenzahl.

Es gilt also|x− rd(x)| ≤ |x− y| für alle y ∈ M.

Bemerkung 2.4. In der Gleitpunktarithmetik gelten die üblichen Rechenregeln, wiedas Assoziativ- oder auch das Distributivgesetz, nicht mehr.

Beispiel 2.5. Bei Maschinenzahlen mit b = 10, t = 3, e ∈ [−1, 1] ergeben sich

(1.17 · 101 ⊕ 1.84 · 100)⊕ 2.43 · 100 = 1.35 · 101 ⊕ 2.43 · 100 = 1.59 · 101,

1.17 · 101 ⊕ (1.84 · 100 ⊕ 2.43 · 100) = 1.17 · 101 ⊕ 4.27 · 100 = 1.60 · 101,

wobei wir mit dem Symbol ⊕ die Addition auf dem Rechner kennzeichnen. Der exakteWert ist 15.97. Die zweite Rechnung liefert eine bessere Approximation.

Definition 2.6 (Exakte und relative Fehler). Sei x ∈ R eine reelle Zahl und rd(x) ∈ Mdie zugehörige Maschinenzahl. Dann definieren wir

eabs(x) := x− rd(x), erel(x) :=x− rd(x)

x.

Hier heißt eabs(x) exakter Fehler und erel(x) relativer Fehler von x.

Definition 2.7 (Maschinengenauigkeit). Mit

eps := min {y ∈ M | 1⊕ y > 1 und y > 0}

bezeichnen wir die Maschinengenauigkeit.

4


Bemerkung 2.8. Für rd(x) = x(1 + ε) gilt nach Definition |ε| ≤ eps und somit

|erel(x)| =∣∣∣∣x− rd(x)

x

∣∣∣∣ = |ε| ≤ eps.

Im Folgenden untersuchen wir, welchen Einfluss der relative Fehler bei der Additionund der Multiplikation hat.Lemma 2.9. Seien x, y ∈ R \ {0} zwei reelle Zahlen, rd(x) und rd(y) die zugehörigenMaschinenzahlen sowie erel(x) und erel(y) die relativen Fehler. Dann gilt:

(i) x+y−(rd(x)+rd(y))x+y = x

x+y · erel(x) + yx+y · erel(y).

(ii) xy−rd(x)rd(y)xy = erel(x) + erel(y)− erel(x)erel(y).

Beweis.

Zu (i):x

x + y· erel(x) +

yx + y

· erel(y)Def.=

xx + y

x− rd(x)x

+y

x + yy− rd(y)

y

=x− rd(x) + y− rd(y)

x + y.

Zu (ii): erel(x) + erel(y)− erel(x)erel(y) =x− rd(x)

x+

y− rd(y)y

− (x− rd(x))(y− rd(y))xy

=xy− rd(x)rd(y)

xy.

Bemerkung 2.10. Da erel(x) und erel(y) sehr klein sind, ist das Produkt erel(x)erel(y)im Vergleich zu erel(x) + erel(y) vernachlässigbar klein, so dass

xy− rd(x)rd(y)xy

= erel(x) + erel(y)− erel(x)erel(y) ≈ erel(x) + erel(y)

gilt. Der relative Fehler bei der Multiplikation ist also klein!Hingegen kann es bei der Addition zu einem großen relativen Fehler kommen, fallsx + y ≈ 0. In diesem Fall spricht man von Auslöschung.

Beispiel 2.11. Seien x = 9970 ≈ 1.41428571 und y =

√2 ≈ 1.41421356. Bei Ma-

schinenzahlen mit b = 10, t = 6, e ∈ [−1, 1] gilt dann rd(x) = 1.41429 · 100 undrd(y) = 1.41421 · 100. Dann ist

eabs(x) = x− rd(x) ≈ −4.3 · 10−6,

eabs(y) = y− rd(y) ≈ 3.6 · 10−6

5

2.2 Vektor- und Matrixnormen

sowie

erel(x) =x− rd(x)

x≈ −3.0 · 10−6,

erel(y) =y− rd(y)

y≈ 2.5 · 10−6.

Also folgt für die Addition

x + y− (rd(x) + rd(y))x + y

Lemma=

xx + y

· erel(x) +y

x + y· erel(y)

≈ −1.5 · 10−6 + 1.25 · 10−6.

Für die Multiplikation gilt

xy− rd(x)rd(y)xy

Lemma= erel(x) + erel(y)− erel(x)erel(y) ≈ −3.0 · 10−6 + 2.5 · 10−6.

Bei der Subtraktion jedoch folgt

x− y− (rd(x)− rd(y))x− y

≈ 0.11,

wir können also den Effekt der Auslöschung feststellen.

Fazit 2.12. Vermeide die Subtraktion annähernd gleichgroßer Zahlen in der Gleitpunk-tarithmetik.


Definition 2.13 (Norm). Sei X ein linearer Raum (Vektorraum) über K (= R oder= C). Eine Abbildung ‖·‖ : X → R heißt Norm, falls sie den folgenden Bedingungengenügt:

(N1) ‖x‖ ≥ 0 ∀x ∈ X und ‖x‖ = 0 ⇔ x = 0,

(N2) ‖λx‖ = |λ| ‖x‖ ∀λ ∈ K ∀x ∈ X (positive Homogenität),

(N3) ‖x + y‖ ≤ ‖x‖+ ‖y‖ ∀x, y ∈ X (Dreiecksungleichung).

Bemerkung 2.14. Da x = x− 0 gilt, kann man ‖x‖ als Abstand von x zum Nullpunktin X interpretieren.(N1) besagt, dass alle Abstände positiv sind, sofern es sich nicht um den Nullpunkthandelt.(N2) besagt, dass die Länge eines Vielfachen gleich das Vielfache der Länge ist.(N3) besagt, dass die Strecke die kürzeste Verbindung zweier Punkte ist.

6


In dieser Vorlesung betrachten wir hauptsächlich endlich-dimensionale Vektorräumeüber K = R.Bemerkung 2.15. Eine Norm auf Rn heißt auch Vektornorm. Im Falle X = Rm×n sprichtman auch von einer Matrixnorm.

Definition 2.16 (lp-Norm). Die lp-Norm auf Rn ist definiert durch

‖x‖p := (n

∑i=1|xi|p)

1p für 1 ≤ p < ∞

‖x‖∞ := maxi∈{1,...,n}

|xi| für p = ∞.

Drei wichtige Spezialfälle sind:

p = 1 : ‖x‖1 =n

∑i=1|xi| (Betragssummennorm)

p = 2 : ‖x‖2 = (n

∑i=1|xi|2)

12 (euklidische Norm)

p = ∞ : ‖x‖∞ = maxi∈{1,...,n}

|xi| (Maximumnorm).

Definition 2.17 (Äquivalenz von Normen). Sei X ein reeller Vektorraum und seien‖·‖α : X → R sowie ‖·‖β : X → R zwei Normen auf X. Dann heißen ‖·‖α und ‖·‖β

äquivalent, falls es positive reelle Konstanten 0 < c1 ≤ c2 existieren, so dass

c1 ‖x‖α ≤ ‖x‖β ≤ c2 ‖x‖α für alle x ∈ X

gilt.

Lemma 2.18. Alle Normen auf Rn sind äquivalent.

Beweis. Sei ‖·‖ : Rn → R eine beliebig aber fest gewählte Norm auf Rn. Wir zeigen,dass ‖·‖ und ‖·‖2 äquivalent sind.Die endlich-dimensionale Einheitssphäre

K = {x ∈ Rn | ‖x‖2 = 1}ist kompakt (Heine-Borel). Da jede Norm stetig ist, nimmt ‖·‖ : Rn → R ihr Minimumund Maximum auf dem Kompaktum K an (Weierstraß). Somit existieren die Konstanten

c1 := minx∈K‖x‖ > 0 und c2 := max

x∈K‖x‖ ≥ c1.

Für x ∈ Rn \ {0} ist y := x‖x‖2∈ K. Dann gilt

c1 = minx∈K‖x‖

y∈K≤ ‖y‖

y∈K≤ max

x∈K‖x‖ = c2 ⇔ c1 ≤

∥∥∥∥x‖x‖2

∥∥∥∥ ≤ c2

⇔ c1 ≤1‖x‖2

‖x‖ ≤ c2

⇔ c1 ‖x‖2 ≤ ‖x‖ ≤ c2 ‖x‖2 .

7


Somit haben wir gezeigt, dass ‖·‖ und ‖·‖2 äquivalent sind. Folglich sind alle Normenauf Rn äquivalent.

Bemerkung 2.19.

(i) Alle Normen auf einem endlich-dimensionalen Vektorraum sind äquivalent.

(ii) Auf einem unendlich-dimensionalen Vektorraum sind zwei beliebige Normen imAllgemeinen nicht mehr äquivalent.

Im folgenden Satz zeigen wir, wie man eine Matrixnorm aus Vektornormen konstruie-ren kann.Satz 2.20. Seien ‖·‖X : Rn → R und ‖·‖Y : Rm → R Normen auf Rn und Rm. Dann wirddurch

‖A‖YX : = supx 6=0

‖Ax‖Y‖x‖X

= sup‖x‖X=1

‖Ax‖Y, A ∈ Rm×n

eine Matrixnorm auf Rm×n definiert.

Beweis. Da die Einheitssphäre K = {x ∈ Rn | ‖x‖X = 1} kompakt ist und die Abbil-dung f : Rn → R, x 7→ ‖Ax‖Y stetig ist, existiert

‖A‖YX = max‖x‖X=1

‖Ax‖Y = maxx 6=0

‖Ax‖Y‖x‖X

∈ R

für alle A ∈ Rm×n. Nun zeigen wir, dass ‖·‖YX : Rm×n → R eine Norm ist:

Zu (N1): Auf Grund der Konstruktion gilt ‖A‖YX ≥ 0 für alle A ∈ Rm×n und

0 = ‖A‖YXDef.= sup

x 6=0

‖Ax‖Y‖x‖X

⇔ ‖Ax‖Y = 0 ∀x ∈ Rn

⇔ Ax = 0 ∀x ∈ Rn

⇔ A = 0.

Zu (N2): Seien λ ∈ R und A ∈ Rm×n. Dann gilt laut Definition

‖λA‖YX = sup‖x‖X=1

‖λAx‖Y = |λ| sup‖x‖X=1

‖Ax‖Y = |λ| ‖A‖YX .

Zu (N3): Seien A, B ∈ Rm×n. Dann gilt

‖A + B‖YX = sup‖x‖X=1

‖(A + B)x‖Y = sup‖x‖X=1

‖Ax + Bx‖Y

≤ sup‖x‖X=1

(‖Ax‖Y + ‖Bx‖Y)

≤ sup‖x‖X=1

‖Ax‖Y + sup‖x‖X=1

‖Bx‖Y

= ‖A‖YX + ‖B‖YX .

8


Somit ist ‖·‖YX : Rm×n → R eine Matrixnorm.

Lemma 2.21. Seien ‖·‖X, ‖·‖Y und ‖·‖Z Vektornormen auf Rn, Rm und Rp. Dann gilt:

(i) ‖Ax‖Y ≤ ‖A‖YX ‖x‖X für alle A ∈ Rm×n und alle x ∈ Rn,

(ii) ‖AB‖YZ ≤ ‖A‖YX ‖B‖XZ für alle A ∈ Rm×n, B ∈ Rn×p.

Beweis.

Zu (i): Nach Definition gilt für A ∈ Rm×n

‖A‖YXDef.= sup

x 6=0

‖Ax‖Y‖x‖X

≥ ‖Ax‖Y‖x‖X

∀x ∈ Rn \ {0}.

Also ist‖Ax‖Y ≤ ‖A‖YX ‖x‖X ∀x ∈ Rn \ {0}

und daraus folgt die Aussage (i).

Zu (ii): Seien A ∈ Rm×n und B ∈ Rn×p gegeben. Dann gilt AB ∈ Rm×p und

‖AB‖YZDef.= sup

x 6=0x∈Rp

‖A(Bx)‖Y‖x‖Z

≤ supx 6=0

x∈Rp

‖A‖YX ‖Bx‖X‖x‖Z

= ‖A‖YX supx 6=0

x∈Rp

‖Bx‖X‖x‖Z

= ‖A‖YX ‖B‖XZ ,

wodurch die Aussage (ii) gezeigt ist.

Bemerkung 2.22. Die Eigenschaft (i) bezeichnet man als Verträglichkeit von Vektor-und induzierter Matrixnorm, die Eigenschaft (ii) als Submultiplikativität von induziertenMatrixnormen.

Einige wichtige Matrixnormen auf Rm×n sind gegeben durch:

‖A‖1 := maxj=1,...,n

m

∑i=1|aij| (Spaltensummennorm)

‖A‖2 :=√

λmax(AT A) (Spektralnorm)

‖A‖∞ := maxi=1,...,m

n

∑j=1|aij| (Zeilensummennorm).

Hierbei bezeichnet λmax(AT A) den größten Eigenwert von der Matrix AT A.

9


Notation 2.23. Für eine Matrix A führen wir folgende Bezeichnung ein:

Rm×n 3 A = (aij) =

a11 · · · a1n... . . . ...

am1 · · · amn

.

Beispiel 2.24. Betrachte

A =

−1 20 21 2

∈ R3×2.

Dann ist

‖A‖1 = max{| − 1|+ 0 + 1, 2 + 2 + 2} = 6,‖A‖∞ = max{| − 1|+ 2, 0 + 2, 1 + 2} = 3.

Für die Spektralnorm berechnen wir zuerst AT A:

AT A =

(−1 0 12 2 2

)−1 20 21 2

=

(2 00 12

).

Auf Grund der Diagonalgestalt der Matrix können wir die Eigenwerte der Matrix AT Adirekt ablesen: λ1 = 2, λ2 = 12. Somit ist

‖A‖2 =√

λmax(AT A) =√

12.

Definition 2.25. Sei A = (aij) ∈ Rn×n eine quadratische Matrix. Dann heißt A

• symmetrisch, falls gilt: AT = A, d.h. aij = aji für alle i, j = 1, . . . , n

• positiv semidefinit, falls gilt: xT Ax ≥ 0 für alle x ∈ Rn

• positiv definit, falls gilt: xT Ax > 0 für alle x ∈ Rn \ {0}

• orthogonal, falls gilt: AAT = I, wobei I ∈ Rn×n die Einheitsmatrix bezeichnet.

Aus der linearen Algebra verwenden wir folgenden Hilfssatz ohne Beweis.Hilfssatz 2.26. Sei A ∈ Rn×n symmetrisch. Dann existiert eine orthogonale Matrix Q ∈Rn×n, so dass

QT AQ = D =

λ1 0. . .

0 λn

= diag(λ1, . . . , λn) ∈ Rn×n

gilt. Die Diagonalelemente λ1, . . . , λn von D sind gerade die Eigenwerte von A.

10


Korollar 2.27. Jede symmetrische Matrix ist genau dann positiv definit (positiv semidefinit),wenn alle Eigenwerte positiv (nichtnegativ) sind.

Korollar 2.28. Sei B ∈ Rn×n eine symmetrische Matrix. Dann gilt

λminxTx ≤ xTBx ≤ λmaxxTx ∀x ∈ Rn,

wobei λmin bzw. λmax den kleinsten bzw. den größten Eigenwert von B bezeichnen.

Beweis. Da B ∈ Rn×n symmetrisch ist, existiert eine orthogonale Matrix Q ∈ Rn×n mit

QTBQ = D = diag(λ1, . . . , λn).

Sei x ∈ Rn. Dann gilt mit y := QTx:

xTBxQQT=I= xTQ QTBQ︸︷︷︸

D

QTx = xTQDQTx = yTDy =n

∑i=1

λiy2i .

Insgesamt folgt:

xTBx =n

∑i=1

λiy2i ≥ λminyTy

y=QT x= λminxTQTQx = λminxTx.

Analog gilt:

xTBx =n

∑i=1

λiy2i ≤ λmaxyTy

y=QT x= λmaxxTQTQx = λmaxxTx.

Somit gilt die Behauptung.

Satz 2.29. Wählen wir sowohl im Raum Rn als auch im Raum Rm die Betragssummennorm‖·‖1, die euklidische Norm ‖·‖2 und die Maximumnorm ‖·‖∞, so sind die jeweils induziertenMatrixnormen gegeben durch:

(i) ‖A‖1 := max1≤j≤n

m

∑i=1|aij| = sup

x 6=0

‖Ax‖1‖x‖1

∀A ∈ Rm×n (Spaltensummennorm)

(ii) ‖A‖2 :=√

λmax(AT A) = supx 6=0

‖Ax‖2‖x‖2

∀A ∈ Rm×n (Spektralnorm)

(iii) ‖A‖∞ := max1≤i≤m

n

∑j=1|aij| = sup

x 6=0

‖Ax‖∞‖x‖∞

∀A ∈ Rm×n (Zeilensummennorm)

11


Beweis.

Zu (i): Für alle x ∈ Rn gilt

‖Ax‖1 =m

∑i=1

∣∣∣∣∣n

∑j=1

aijxj

∣∣∣∣∣ ≤m

∑i=1

n

∑j=1|aij||xj| =

n

∑j=1

m

∑i=1|aij||xj| ≤

n

∑j=1

max1≤j≤n

m

∑i=1|aij||xj|

= max1≤j≤n

m

∑i=1|aij|

n

∑j=1|xj|

= ‖A‖1 ‖x‖1 .

Daraus folgt

supx 6=0

‖Ax‖1‖x‖1

≤ ‖A‖1 .

Umgekehrt wähle ein j mit ‖A‖1 = ∑mi=1 |aij|, x = ej. Dann ist

‖Ax‖1 =∥∥Aej

∥∥1 =

m

∑i=1|aij| = ‖A‖1 ‖x‖1 .

Also folgt für x = ej‖Ax‖1‖x‖1

= ‖A‖1 .

Somit ist

supx 6=0

‖Ax‖1‖x‖1

≥ ‖A‖1

und insgesamt ergibt sich

‖A‖1 ≤ supx 6=0

‖Ax‖1‖x‖1

≤ ‖A‖1

und daraus ergibt sich die Behauptung.

Zu (ii): Für alle x ∈ Rn gilt

‖Ax‖22

Def.= (Ax)T Ax = xT AT Ax

Kor.≤ λmax(AT A)xTx = λmax(AT A) ‖x‖2

2 .

Daraus folgt

‖Ax‖2 ≤√

λmax(AT A) ‖x‖2

und somit

supx 6=0

‖Ax‖2‖x‖2

≤√

λmax(AT A) = ‖A‖2 .

12


Umgekehrt sei v ∈ Rn ein Eigenvektor von AT A zum Eigenwert λmax(AT A).Dann gilt

‖Av‖22 = vT AT Av = λmax(AT A)vTv = λmax(AT A) ‖v‖2

2 .

Folglich ist

‖Av‖2 =√

λmax(AT A) ‖v‖2 = ‖A‖2 ‖v‖2

und somit

‖A‖2 =‖Av‖2‖v‖2

≤ supx 6=0

‖Ax‖2‖x‖2

.

Insgesamt gilt also

‖A‖2 ≤ supx 6=0

‖Ax‖2‖x‖2

≤ ‖A‖2

und daraus folgt die Behauptung.

Zu (iii): Es gilt für alle x ∈ Rn

‖Ax‖∞ = max1≤i≤m

|n

∑j=1

aijxj| ≤ max1≤i≤m

n

∑j=1|aij| max

1≤j≤n|xj| = ‖A‖∞ ‖x‖∞ .

Daraus folgt

supx 6=0

‖Ax‖∞‖x‖∞

≤ ‖A‖∞ .

Sei nun k ∈ {1, . . . , m} beliebig aber fest. Wir definieren den Vektor v ∈ Rn mit

vj :=

|akj|akj

falls akj 6= 0

1 falls akj = 0.

Offenbar ist ‖v‖∞ = 1. Somit gilt

supx 6=0

‖Ax‖∞‖x‖∞

≥ ‖Av‖∞‖v‖∞

= ‖Av‖∞ = max1≤i≤m

|(Av)i| ≥ |(Av)k|

= |n

∑j=1

akjvj|

Def.= |

n

∑j=1|akj||.

Daher ist

supx 6=0

‖Ax‖∞‖x‖∞

≥n

∑j=1|akj|.

13

2.3 Kondition und Fehlerabschätzung

Da aber der Zeilenindex k frei gewählt wurde, so folgt

supx 6=0

‖Ax‖∞‖x‖∞

≥ max1≤i≤m

n

∑j=1|aij| = ‖A‖∞ .

Insgesamt haben wir also gezeigt, dass

‖A‖∞ ≤ supx 6=0

‖Ax‖∞‖x‖∞

≤ ‖A‖∞

gilt und daraus folgt die Behauptung.

Folgerung 2.30. Es gilt für µ ∈ {1, 2, ∞}:(i): ‖Ax‖µ ≤ ‖A‖µ ‖x‖µ (Verträglichkeit),

(ii): ‖AB‖µ ≤ ‖A‖µ ‖B‖µ (Submultiplikativität).


Definition 2.31. Sei A ∈ Rm×n eine Matrix. Dann definieren wir:

Kern(A) := {x ∈ Rn | Ax = 0}Bild(A) := {y ∈ Rm | ∃x ∈ Rn : y = Ax}Rang(A) := dim(Bild(A)).

Aus der linearen Algebra verwenden wir ohne Beweis den folgenden Satz.Satz 2.32 (Dimensionssatz). Sei A ∈ Rm×n eine Matrix. Dann gilt

n = dim(Kern(A)) + dim(Bild(A)).

Lemma 2.33 (Lösbarkeitskriterium für Ax = b). Sei A ∈ Rn×n und b ∈ Rn. Dann sinddie folgenden Aussagen äquivalent:

(i) Das lineare Gleichungssystem Ax = b hat mindestens eine Lösung.

(ii) Es gilt b ∈ Bild(A).

(iii) Rang(A) = Rang(A, b), wobei (A, b) ∈ Rn×(n+1) diejenige Matrix ist, die aus Adurch Hinzunahme von b als weitere Spalte entsteht.

Beweis. Die Äquivalenz zwischen (i) und (ii) folgt direkt aus der Definition des Bildes.Aus (ii) folgt wieder mit der Definition des Bildes Bild(A) = Bild(A, b) und darauswieder (iii). Umgekehrt folgt aus (iii), dass b als Linearkombination von A darstellbarist. Dies liefert b ∈ Bild(A), woraus sich (ii) ergibt. Damit ist insgesamt die Äquivalenzaller Aussagen gezeigt.

14


Lemma 2.34 (Lösungsmenge). Seien A ∈ Rn×n und b ∈ Rn. Es existiere x∗ ∈ Rn, so dassAx∗ = b ist. Dann gilt für die Lösungsmenge:

L := {x ∈ Rn | Ax = b} = x∗ + Kern(A).

Beweis.

⊆: Sei x ∈ L. Dann gilt:

A(x− x∗) = Ax− Ax∗ = b− b = 0 Def.⇒ x− x∗ ∈ Kern(A).

Damit ergibt sich:x = x∗ + x− x∗︸︷︷︸

∈Kern(A)

∈ x∗ + Kern(A).

⊇: Sei x = x∗ + z mit z ∈ Kern(A). Dann gilt:

Ax = A(x∗ + z) = Ax∗ + Az = b + 0.

Damit ergibt sich x ∈ L.

Definition 2.35. Sei A ∈ Rn×n eine quadratische Matrix. Dann heißt A regulär, falls

Kern(A) = {0}.

Folgerung 2.36. Sei A ∈ Rn×n eine quadratische Matrix. Dann sind die folgendenAussagen äquivalent:

(i) A ist regulär (Kern(A) = {0}).

(ii) Bild(A) = Rn.

(iii) Das lineare Gleichungssystem Ax = b hat für jedes b ∈ Rn genau eine Lösungx ∈ Rn.

Bemerkung 2.37. Aus der linearen Algebra ist auch bekannt:

A ist regulär ⇔ det(A) 6= 0.

Jede reguläre Matrix A ∈ Rn×n besitzt eine inverse Matrix A−1 ∈ Rn×n, die durch

A−1A = I ⇔ AA−1 = I

charakterisiert ist.Das Produkt zweier regulärer Matrizen A, B ∈ Rn×n ist regulär mit

(AB)−1 = B−1A−1,

15


denn es ist(B−1A−1)AB = B−1(A−1A)B = B−1 IB = I.

Ist A ∈ Rn×n regulär, so ist die transponierte Matrix AT auch regulär mit

(AT)−1 = (A−1)T,

denn es gilt(A−1)T AT = (A−1A)T = IT = I.

Definition 2.38 (Orthogonalraum). Sei U ⊆ Rn nichtleer. Mit

U⊥ := {y ∈ Rn | 〈y, u〉 = 0 ∀u ∈ U}

bezeichnen wir den Orthogonalraum von U. Hierbei ist

〈·, ·〉 : Rn ×Rn → R

das euklidische Skalarprodukt 〈x, y〉 = xTy.

Bemerkung 2.39. Ist U ⊆ Rn nichtleer, so gilt:

(i) U⊥ ⊆ Rn ist ein Untervektorraum.

(ii) U ⊆ (U⊥)⊥.

Satz 2.40. Sei U ⊆ Rn ein Untervektorraum (Unterraum) und U⊥ ⊆ Rn der Orthogonal-raum von U. Dann gilt

Rn = U ⊕U⊥.

Mit anderen Worten gibt es zu jedem x ∈ Rn genau ein u ∈ U und genau ein u⊥ ∈ U⊥, sodass

x = u + u⊥.

Beweis. Sei {u1, . . . , um} eine Orthonormalbasis von U. Definiere

P : Rn → U, P(x) =m

∑i=1〈x, ui〉ui.

Sei x ∈ Rn. Wir zeigen:

x− P(x) ∈ U⊥ Def.⇔ 〈x− P(x), u〉 = 0 ∀u ∈ U.

Sei also u ∈ U beliebig aber fest. Dann hat u die Darstellung

u =m

∑i=1

λiui, λi ∈ R,

16


da {u1, . . . , um} eine Basis ist. Folglich gilt

〈x− P(x), u〉 = 〈x, u〉 − 〈P(x), u〉

=m

∑i=1

(λi〈x, ui〉 − λi〈P(x), ui〉)

=m

∑i=1

(λi〈x, ui〉 − λi

m

∑j=1〈x, uj〉〈uj, ui〉

)

〈uj,ui〉=δij=

m

∑i=1

(λi〈x, ui〉 − λi〈x, ui〉)

= 0.

Daher ist 〈x− P(x), u〉 = 0 für alle u ∈ U. Mit anderen Worten: x− P(x) ∈ U⊥.Setzen wir

u := P(x) ∈ U und u⊥ := x− P(x) ∈ U⊥,

so hat x die Darstellung

x = P(x) + x− P(x) = u + u⊥.

Zur Eindeutigkeit: Gibt es ein u ∈ U und ein u⊥ ∈ U⊥ mit

x = u + u⊥,

so müssen wir nun zeigen, dass u = u und u⊥ = u⊥ ist. Es gilt

0 = u− u + u⊥ − u⊥,

woraus0 = 〈u− u, u− u〉+ 〈u⊥ − u⊥, u− u〉

folgt. Jetzt ist aber per Definition 〈u⊥ − u⊥, u− u〉 = 0, denn es ist u⊥ − u⊥ ∈ U⊥ undu− u ∈ U. Also muss auch

〈u− u, u− u〉 = 0

gelten. Aus der Definition der Norm folgt aus 0 = ‖u− u‖22, dass u = u gilt und somit

auch u⊥ = u⊥.

Folgerung 2.41. Ist U ⊆ Rn ein Untervektorraum, dann gilt

U = (U⊥)⊥.

Satz 2.42. Sei A ∈ Rm×n eine Matrix. Dann gilt:

(i) Kern(A) = Bild(AT)⊥ und Kern(A)⊥ = Bild(AT).

(ii) Kern(A) = Kern(AT A).

17


(iii) Bild(AT) = Bild(AT A).

Beweis.

Zu (i):

⊆: Ist v ∈ Kern(A), so ist per Definition Av = 0. Daraus folgt 〈v, ATx〉 =〈Av, x〉 = 0 für alle x ∈ Rm. Mit anderen Worten ist also v ∈ Bild(AT)⊥.

⊇: Ist v ∈ Bild(AT)⊥, so ist per Definition 〈v, ATx〉 = 0 für alle x ∈ Rm. Diesist gleichbedeutend mit 〈Av, x〉 = 0 für alle x ∈ Rm. Mit der Substitutionx = Av folgt 〈Av, Av〉 = 0. Daher ist Av = 0 bzw. v ∈ Kern(A).

Zu (ii):

⊆: Ist v ∈ Kern(A), so ist Av = 0. Daraus folgt AT Av = 0. Das bedeutetv ∈ Kern(AT A).

⊇: Ist v ∈ Kern(AT A), so ist AT Av = 0. Daraus folgt vT AT Av = 0 bzw.(Av)T Av = 0, woraus sich Av = 0 ergibt. Dies bedeutet v ∈ Kern(A).

Zu (iii): Es gilt:

Bild(AT) = (Bild(AT)⊥)⊥(i)= Kern(A)⊥ = Kern(AT A)⊥

(i)= (Bild(AT A)⊥)⊥

= Bild(AT A).

Wir wollen nun die Stabilität der Lösung des linearen Gleichungssystems

Ax = b

mit einer regulären Matrix A ∈ Rn×n und einem vorgegebenen Vektor b ∈ Rn untersu-chen. Diese Aufgabe hat genau eine Lösung x ∈ Rn, da A regulär ist. Betrachten wirden Vektor b mit einer „kleinen“ Störung, also b ∈ Rn, dann hat

Ax = b

genau eine Lösung x ∈ Rn. Wir stellen also in diesem Zusammenhang die Frage, wiesich

‖x− x‖verhält.

Satz 2.43. Es sei A ∈ Rn×n eine reguläre Matrix sowie b, b ∈ Rn Vektoren mit b 6= 0. Weiterseien x und x Lösungen der linearen Gleichungssysteme

Ax = b und Ax = b

18


(b ist die gestörte rechte Seite). Dann gilt

‖x− x‖‖x‖ ≤ ‖A‖ ‖A−1‖

∥∥b− b∥∥

‖b‖ .

Hierbei bezeichnet ‖·‖ eine (beliebige) Vektornorm bzw. die hierdurch induzierte Matrixnorm.

Beweis. Sei ‖·‖ : Rn → R eine Vektornorm und die induzierte Matrixnorm ist gegebendurch

‖B‖ = supx 6=0

‖Bx‖‖x‖ ∀B ∈ Rn×n.

Da A regulär ist, gilt

‖x− x‖ = ‖A−1(b− b)‖ ≤ ‖A−1‖∥∥b− b

∥∥ .

Außerdem gilt

‖b‖ = ‖Ax‖ ≤ ‖A‖ ‖x‖ b 6=0, x 6=0⇒ 1‖x‖ ≤

‖A‖‖b‖ .

Insgesamt ergibt sich also

‖x− x‖‖x‖ ≤ ‖A‖ ‖A−1‖

∥∥b− b∥∥

‖b‖ .

Bemerkung 2.44. Der Satz besagt

‖x− x‖‖x‖︸︷︷︸

Relativer Fehler in x

≤ ‖A‖ ‖A−1‖︸︷︷︸Verstärkungsfaktor

‖b− b‖‖b‖︸︷︷︸

Relativer Fehler in b

.

Ist der Verstärkungsfaktor klein, so ist der relative Fehler in x auch klein. Problematischist, wenn ‖A‖ ‖A−1‖ groß ist. In diesem Fall führt eine kleine Störung der rechtenSeite b auf eine große Änderung der Lösung.

Definition 2.45 (Kondition). Sei A ∈ Rn×n regulär und ‖·‖ : Rn×n → R eine Ma-trixnorm. Dann heißt

cond(A) := ‖A‖ ‖A−1‖Kondition von A. Im Falle der Spektralnorm

‖A‖2 =√

λmax(AT A)

wird die Konditioncond2(A) := ‖A‖2 ‖A−1‖2

auch als Spektralkondition bezeichnet.

19


Bemerkung 2.46. Ist ‖·‖ : Rn×n → R eine induzierte Matrixnorm, so lässt sich derSatz wie folgt formulieren:

‖x− x‖‖x‖ ≤ cond(A)

‖b− b‖‖b‖ .

Es gilt für cond(A)

1 = ‖I‖ = ‖AA−1‖Submultiplikativität

≤ ‖A‖ ‖A−1‖ = cond(A).

Definition 2.47. Eine reguläre Matrix A ∈ Rn×n heißt schlecht konditioniert, falls dieKondition cond(A) groß ist (cond(A) ≫ 1).

Lemma 2.48. Ist Q ∈ Rn×n orthogonal (d.h. QTQ = QQT = I), so gilt für die Spektralkon-dition von Q:

cond2(Q) = 1.

Beweis. Sei Q ∈ Rn×n orthogonal. Dann gilt:

‖Qx‖22 = (Qx)TQx = xTQTQx = xTx = ‖x‖2

2 .

Daraus folgt:

‖Q‖2 = supx 6=0

‖Qx‖2‖x‖2

= 1.

Analog gilt:

‖Q−1x‖22 = (Q−1x)TQ−1x = xT(Q−1)TQ−1x = xT(QT)−1Q−1x = xT(QQT)−1x

= ‖x‖22 .

Daraus folgt

‖Q−1‖2 = supx 6=0

∥∥Q−1x∥∥

2‖x‖2

= 1

und insgesamt giltcond2(Q) = ‖Q‖2 ‖Q−1‖2 = 1.

Lemma 2.49. Sei A ∈ Rn×n symmetrisch und positiv semidefinit. Dann gilt

λmax(AA) = λmax(A)2.

Beweis. Da A ∈ Rn×n symmetrisch ist, existiert eine orthogonale Matrix Q ∈ Rn×n, sodass

QT AQ = D = diag(λ1, . . . , λn) ∈ Rn×n,

20


wobei λ1, . . . , λn die Eigenwerte von A sind. Weiter sind

λ1, . . . , λn ≥ 0,

weil A positiv semidefinit ist. Wir zeigen, dass

λ21, . . . , λ2

n

die Eigenwerte von AA sind.Es gilt

AA = QQT AQQT AQQT = QD2QT

mitD2 = diag(λ2

1, . . . , λ2n).

Folglich gilt für jedes i = 1, . . . , n

AA Qei︸︷︷︸=w

= QD2QTQei = QD2ei = Qλ2i ei = λ2

i Qei = λ2i w,

was bedeutet, dass λ2i Eigenwert von AA mit Eigenvektor w = Qei 6= 0 ist. Somit sind

λ21, . . . , λ2

n Eigenwerte von AA. Folglich ist die Aussage richtig.

Beispiel 2.50. Die Aussage gilt nicht für allgemeine Matrizen.

A =

(1 00 −10

), AA =

(1 00 100

).

Hier ist λmax(AA) = 100 6= 1 = λmax(A)2.

Lemma 2.51. Es gilt‖BTB‖2 = ‖B‖2

2 ∀B ∈ Rn×n.

Beweis. Es ist

‖BTB‖2 =√

λmax((BTB)T(BTB)) =√

λmax((BTB)(BTB)).

Die Matrix BTB ∈ Rn×n ist symmetrisch und positiv semidefinit. Somit liefert dasLemma 2.49

λmax((BTB)T(BTB)) = λmax(BTB)2.

Daher gilt

‖BTB‖2 =√

λmax((BTB)T(BTB)) = λmax(BTB) = ‖B‖22 .

Satz 2.52. Sei A ∈ Rn×n regulär. Dann gilt für die Spektralkondition:

cond2(AT A) = cond2(A)2.

21


Beweis. Zunächst zeigen wir, dass die Eigenwerte von AT A und AAT gleich sind:

λ ist ein Eigenwert von AT A mit Eigenvektor v ∈ Rn \ {0}⇔ AT Av = λv

A regulär⇔ AAT Av︸︷︷︸w

= λAv

⇔ AATw = λw

⇔ λ ist ein Eigenwert von AAT mit Eigenvektor w ∈ Rn \ {0}.

Daraus folgtλmax(AT A) = λmax(AAT).

Daher ist‖A‖2 =

√λmax(AT A) =

√λmax(AAT) = ‖AT‖2.

Die gleiche Aussage gilt auch für A−1

‖A−1‖2 = ‖(A−1)T‖2.

Das obige Lemma liefert‖AT A‖2 = ‖A‖2

2

und

‖A−1(A−1)T‖2B=(A−1)T

= ‖(A−1)T‖22

s.o.= ‖A−1‖2

2.

Daraus folgt

cond2(AT A) = ‖AT A‖2‖(AT A)−1‖2 = ‖A‖22‖(AT A)−1‖2 = ‖A‖2

2 ‖A−1(A−1)T‖2

= ‖A‖22 ‖A−1‖2

2

= cond2(A)2.

Bemerkung 2.53. Ist A schlecht konditioniert, dann ist AT A noch schlechter konditio-niert.

Satz 2.54. Ist A ∈ Rn×n symmetrisch und positiv definit, so gilt

cond2(A) =λmax(A)

λmin(A).

Beweis. Es gilt

‖A‖22 = λmax(AT A)

A symm.= λmax(AA)

Lemma + A s.p.d.= λmax(A)2.

22


Ist λ Eigenwert von A, so ist λ−1 Eigenwert von A−1. Folglich ist

‖A−1‖22 = λmax((A−1)T A−1) = λmax(A−1A−1) = λmax(A−1)2 =

1λmin(A)2 .

Daraus folgt

cond2(A) = ‖A‖2 ‖A−1‖2 =λmax(A)

λmin(A).

23

Kapitel 3Direkte Verfahren

Dieses Kapitel befasst sich mit direkten Verfahren zur Lösung von linearen Gleichungs-systemen der Gestalt Ax = b.

3.1 Lineare Gleichungssysteme einfacher Strukturen

Wir betrachten zunächst ein lineares Gleichungssystem der Form

Dx = b

mit einer Diagonalmatrix

D = diag(d1, . . . , dn) =

d1. . .

dn

∈ Rn×n.

Die Lösung x ist offenbar gegeben durch

dixi = bi ∀i = 1, . . . , n.

Hieraus ergibt sich:

Algorithmus 3.1 Lösung von Dx = b mit einer Diagonalmatrix D ∈ Rn×n

1: for i = 1 : n do2: xi := bi/di3: end for

Der Algorithmus 3.1 ist genau dann durchführbar, wenn di 6= 0 für alle i = 1, . . . , n,also wenn die Matrix D regulär ist.Sei nun L ∈ Rn×n eine untere Dreiecksmatrix

L =

l11 0... . . .

ln1 · · · lnn

.

25


Wir betrachten das Gleichungssystem Lx = b. Das ist äquivalent zu:

l11x1 = b1

l21x1 + l22x2 = b2

......

ln1x1 + ln2x2 + · · ·+ lnnxn = bn.

Aus der ersten Zeile lässt sich x1 berechnen, danach ergibt sich x2 aus der zweitenZeile, usw. Diese Vorgehensweise heißt Vorwärtseinsetzen oder Vorwärtssubstitution.

Algorithmus 3.2 Lösung von Lx = b mit einer unteren Dreiecksmatrix L ∈ Rn×n

1: x1 := b1/l112: for i = 2 : n do3: xi := (bi −∑i−1

j=1 lijxj)/lii4: end for

Der Algorithmus 3.2 ist genau dann durchführbar, wenn lii 6= 0 für alle i = 1, . . . , n,also wenn L regulär ist. Zur Berechnung von xi benötigen wir:

• Für i = 1: Eine Division.

• Für i > 1: Eine Division, i− 1 Multiplikationen und i− 1 Subtraktionen.

Zur Durchführung des Algorithmus 3.2 benötigen wir also

1 +n

∑i=2

(2(i− 1) + 1) = n2

Rechenoperationen.Sei nun R ∈ Rn×n eine obere Dreiecksmatrix:

R =

r11 · · · r1n. . . ...

0 rnn

.

Das lineare Gleichungssystem lautet

Rx = b,

was äquivalent ist zu

r11x1 + r12x2 + · · ·+r1nxn = b1

r22x2 + · · ·+r2nxn = b2

. . . ...rnnxn = bn.

26


Dieses System lässt sich explizit lösen, indem man zunächst xn aus der n-ten Gleichungbestimmt, danach ergibt sich xn−1 aus der vorletzten Gleichung.

Algorithmus 3.3 Rückwärtssubstitution für Rx = b mit einer oberen DreiecksmatrixR ∈ Rn×n

1: xn := bn/rnn2: for i = (n− 1) : −1 : 1 do3: xi := (bi −∑n

j=i+1 rijxj)/rii4: end for

Zur Durchführung des Algorithmus 3.3 benötigen wir ebenso n2 Rechenoperationen.Dieser Algorithmus ist durchführbar, wenn rii 6= 0 ∀i = 1, . . . , n.Zum Abschluss betrachten wir ein lineares Gleichungssystem mit einer Permutations-matrix.

Definition 3.1 (Permutationsmatrix). Mit ei ∈ Rn bezeichnen wir den i-ten Einheits-vektor. Eine Matrix P ∈ Rn×n heißt Permutationsmatrix, falls P die folgende Gestalthat:

P =

(eπ(1))T

(eπ(2))T

...(eπ(n))

T

mit einer bijektiven Abbildung (Permutation)

π : {1, . . . , n} → {1, . . . , n}.

Notation: P = P(π).

Beispiel 3.2. Die folgenden zwei Matrizen sind Beispiele für Permutationsmatrizenim R3×3:

P =

1 0 00 1 00 0 1

oder P =

0 1 01 0 00 0 1

.

Bemerkung 3.3. Jede Permutationsmatrix ist orthogonal, also PTP = PPT = I.

Betrachte das GleichungssystemPx = b

mit einer Permutationsmatrix Rn×n 3 P = P(π) und

π : {1, . . . , n} → {1, . . . , n}

27

3.2 LR-Zerlegung ohne Pivotisierung

ist eine bijektive Abbildung. Dieses System ist äquivalent zu

(eπ(1))T

(eπ(2))T

...(eπ(n))

T

x = b.

Äquivalent lässt sich das schreiben als:

eTπ(i)x = bi

xπ(i) = bi ∀i = 1, . . . , n.

Hieraus erhalten wir:

Algorithmus 3.4 Lösung von Px = b mit einer Permutationsmatrix P = P(π)

1: for i = 1 : n do2: xπ(i) := bi3: end for

Analog lässt sich die Lösung vonPTx = b

wie folgt berechnen:

PTx = b ⇒ PPTx = Pb ⇒ x = Pb.

Algorithmus 3.5 Lösung von PTx = b mit einer Permutationsmatrix P = P(π)

1: for i = 1 : n do2: xi := bπ(i)3: end for


Im Folgenden sei A ∈ Rn×n regulär und b ∈ Rn. Gesucht sei die eindeutige Lösungvon

Ax = b,

mit Hilfe des Eliminationsverfahrens von Gauß. Diese Methode führt auf eine Zerle-gung von A der Form

A = LR

28


mit einer normierten unteren Dreiecksmatrix

L =

1 0. . .

∗ 1

∈ Rn×n

und einer oberen Dreiecksmatrix

R =

∗ · · · ∗

. . . ...0 ∗

∈ Rn×n.

Definition 3.4 (LR-Zerlegung). Unter einer LR-Zerlegung einer Matrix versteht maneine Zerlegung der Form

A = LR

mit einer normierten unteren Dreiecksmatrix L ∈ Rn×n und einer oberen Dreiecksma-trix R ∈ Rn×n.

Mit Hilfe der LR-Zerlegung lässt sich Ax = b sehr einfach lösen. Denn:

Ax = b ⇔ L(Rx) = b.

Dazu löst manLy = b

mittels Vorwärtseinsetzen und dann

Rx = y

durch Rückwärtssubstitution.

Algorithmus 3.6 Lösung von Ax = b mittels LR-Zerlegung

(S1) Bestimme LR-Zerlegung(S2) Löse Ly = b durch Vorwärtseinsetzen (Algorithmus 2.2)(S3) Löse Rx = y durch Rückwärtssubstitution (Algorithmus 2.3)

Definition 3.5 (Frobenius-Matrix). Eine Matrix Li ∈ Rn×n heißt Frobenius-Matrix, wennLi die folgende Form hat:

Li =

1 0. . .

1

−li+1,i. . .

...−ln,i 1

.

29


Beispiel 3.6. Betrachte das lineare Gleichungssystem Ax = b mit

2 1 1 04 3 3 18 7 9 56 7 9 8

und b ∈ Rn.

Das System Ax = b ist äquivalent zu

2x1 + 1x2 + 1x3 + 0x4 = b1

4x1 + 3x2 + 3x3 + 1x4 = b2

8x1 + 7x2 + 9x3 + 5x4 = b3

6x1 + 7x2 + 9x3 + 8x4 = b4.

Definiere

L1 :=

1−2 1−4 1−3 1

.

Dann ist

L1A =

1−2 1−4 1−3 1

2 1 1 04 3 3 18 7 9 56 7 9 8

=

2 1 1 01 1 13 5 54 6 8

=: A(2)

L2A(2) =

11−3 1−4 1

︸︷︷︸=:L2

2 1 1 01 1 13 5 54 6 8

=

2 1 1 01 1 1

2 22 4

=: A(3)

L3A(3) =

11

1−1 1

︸︷︷︸=:L3

2 1 1 01 1 1

2 22 4

=

2 1 1 01 1 1

2 22

= R.

Daraus folgt A = L−11 L−1

2 L−13 R mit Frobeniusmatrizen L1, L2 und L3. Die Inversen von

Li lauten:

L−11 =

12 14 13 1

, L−1

2 =

113 14 1

, L−1

3 =

11

11 1

.

30


Damit ergibt sich:

L = L−11 L−1

2 L−13 =

12 14 3 13 4 1 1

undA = LR.

Lemma 3.7. Es sei li := (0, . . . , 0, li+1,i, . . . , ln,i) ∈ Rn. Definiere die Frobenius-Matrix wiefolgt:

Li := I − lieTi =

1. . .

1

−

0...0

li+1,i...

ln,i

(0, . . . , 1, . . . , 0

)

=

1. . .

1

−li+1,i. . .

...−ln,i 1

.

Dann gilt:

(i) L−1i = I + lieT

i =

1. . .

1

li+1,i. . .

...ln,i 1

,

(ii) L−11 · L−1

2 · . . . · L−1k = I + ∑k

i=1 lieTi =

1

li,1. . .

... 1

lk+1,1 ··· lk+1,k. . .

...... . . .

ln,1 ··· ln,k 1

für k = 1, . . . , n.

Beweis.

Zu (i): (I − lieTi )︸︷︷︸

=Li

(I + lieTi ) = I − lieT

i + lieTi − li eT

i li︸︷︷︸=0

eTi = I.

31


Zu (ii): Beweis durch vollständige Induktion.

Induktionsanfang k=1: Wurde bereits in (i) gezeigt.

Induktionsannahme: Die Aussage gelte für ein festes aber beliebiges k mit 1 ≤k ≤ n− 2. Wir zeigen, dass die Aussage für k + 1 gilt.

L−11 · L−1

2 · . . . · L−1k+1 = (I +

k

∑i=1

lieTi )L−1

k+1

(i)= (I +

k

∑i=1

lieTi )(I + lk+1eT

k+1)

= I + lk+1eTk+1 +

k

∑i=1

lieTi +

k

∑i=1

li eTi lk+1︸︷︷︸=0

eTk+1

= I +k+1

∑i=1

lieTi .

Algorithmus 3.7 Grundversion einer LR-Zerlegung

1: Setze A(1) := A2: for i = 1 : (n− 1) do

3: Definiere eine Frobeniusmatrix Li =

1. . .

1

−li+1,i. . .

...−ln,i 1

mit lj,i := a(i)j,i /a(i)i,i für j = i + 1, . . . , n.

Setze A(i+1) := Li A(i)

4: end for5: Setze R := A(n), L := L−1

1 · . . . · L−1n−1

Bemerkung 3.8. Der Algorithmus 3.7 liefert eine LR-Zerlegung, wenn er durchführbarist, denn

Ln−1 · . . . · L1A = R ⇔ A = L−11 ·

. . . · L−1n−1R = LR.

Dabei ist L eine untere normierte Dreiecksmatrix und R eine obere Dreiecksmatrix.Der Algorithmus ist durchführbar, wenn

a(i)ii 6= 0 ∀i = 1, . . . , n− 1.

32


Definition 3.9. Das i-te Diagonalelement a(i)ii der Matrix

A(i) = Li−1 · . . . · L1A

heißt Pivot-Element.

Definition 3.10. Für eine quadratische Matrix A ∈ Rn×n heißt

A[k] :=

a11 · · · a1k... . . . ...

ak1 · · · akk

∈ Rk×k, k ∈ {1, . . . , n}

die führende k− k−Hauptachsenabschnittsmatrix (Hauptabschnittsmatrix) von A. DieDeterminante det(A[k]) heißt die führende k−te Hauptachsenabschnittsdeterminantevon A.

Beispiel 3.11. Es sei

A =

1 2 34 5 67 8 9

.

Dann ist

A[1] = (1), A[2] =(

1 24 5

), A[3] = A.

Lemma 3.12. Es sei A ∈ Rn×n und L ∈ Rn×n eine untere (nicht notwendigerweise normierte)Dreiecksmatrix. Dann gilt

(LA)[k] = L[k]A[k] ∀k = 1, . . . , n.

Beweis. Sei k ∈ {1, . . . , n} beliebig aber fest gewählt. Für i, j ∈ {1, . . . , k} gilt

((LA)[k])i,j =n

∑m=1

li,mam,j =k

∑m=1

li,mam,j +n

∑m=k+1

li,m︸︷︷︸=0

am,j,

denn L ∈ Rn×n ist eine untere Dreiecksmatrix, also li,j = 0 ∀j > i. Also folgt

((LA)[k])i,j =k

∑m=1

li,mam,j = (L[k]A[k])i,j.

Satz 3.13 (Existenz einer LR-Zerlegung). Es sei A ∈ Rn×n eine reguläre Matrix. Dannbesitzt A genau dann eine LR-Zerlegung, wenn

det(A[k]) 6= 0 ∀k = 1, . . . , n.

33


Beweis.

„⇒“: A habe eine LR-Zerlegung, das heißt

A = LR

mit einer normierten unteren Dreiecksmatrix L ∈ Rn×n und einer oberen Drei-ecksmatrix R ∈ Rn×n. Dann gilt

det(A) = det(L)det(R) (Determinantenmultiplikationssatz).

Da A regulär ist, gilt

0 6= det(A) ⇒ det(L) 6= 0 und det(R) 6= 0.

Hieraus folgt

det(L[k]) 6= 0 und det(R[k]) 6= 0 ∀k = 1, . . . , n,

denn

det

λ1 ∗. . .

λn

= λ1 · . . . · λn und det

λ1. . .

∗ λn

= λ1 · . . . · λn.

Das Lemma 3.12 liefert

det(A[k]) Lemma 3.12= det(L[k]R[k]) = det(L[k])det(R[k]) 6= 0 ∀k = 1, . . . , n.

„⇐“: Es geltedet(A[k]) 6= 0 ∀k = 1, . . . , n.

Wir zeigen nun, dass A eine LR-Zerlegung besitzt. Dazu zeigen wir, dass derAlgorithmus 3.7 durchführbar ist, also dass

a(i)ii 6= 0 ∀i = 1, . . . , n

gilt. Beweis durch vollständige Induktion:

Induktionsanfang i=1: a(1)11De f= det(A[1]) 6= 0 (Voraussetzung)

Induktionsannahme: Es gelte a(1)11 · . . . · a(i)ii 6= 0 für ein festes aber beliebiges imit 1 ≤ i ≤ n− 2. Nun zeigen wir:

a(i+1)i+1,i+1 6= 0.

Aufgrund der Induktionsannahme existiert

A(i+1) = Li · . . . · L1A (siehe Algorithmus 3.7).

34


Hieraus folgt

A(i+1)[i + 1] Lemma 3.12= Li[i + 1] · . . . · L1[i + 1]A[i + 1].

Daher ist

det(A(i+1)[i + 1]) = det(Li[i + 1])︸︷︷︸=1

· . . . · det(L1[i + 1])︸︷︷︸=1

det(A[i + 1])︸︷︷︸6=0

6= 0.

Laut Konstruktion (Algorithmus 3.7) ist

Ai+1[i + 1] =

a(i+1)11 ∗

. . .

a(i+1)i+1,i+1

.

Damit folgt

0 6= det(A(i+1)[i + 1]) = a(i+1)11 · . . . · a(i+1)

i+1,i+1.

Insgesamt ergibt sicha(i+1)

i+1,i+1 6= 0

und daraus folgt die Behauptung.

Wir zeigen nun, dass es genau eine LR-Zerlegung gibt, falls diese existiert. Dazuverwenden wir folgende Lemmata.Lemma 3.14. Es gelten die folgenden Aussagen:

(i) Sind L(1), L(2) ∈ Rn×n zwei untere (normierte) Dreiecksmatrizen, so ist das ProduktL := L(1)L(2) eine untere (normierte) Dreiecksmatrix.

(ii) Sind R(1), R(2) ∈ Rn×n zwei obere Dreiecksmatrizen, so ist das Produkt R := R(1)R(2)

eine obere Dreiecksmatrix.

Beweis. Seien L(1), L(2) ∈ Rn×n untere Dreiecksmatrizen, setze L := L(1)L(2). Sei i ∈{1, . . . , n− 1} beliebig aber fest. Dann gilt für j > i:

l(1)im = 0 ∀m = j, . . . , n

l(2)mj = 0 ∀m = 1, . . . , j− 1,

da L(1) und L(2) untere Dreiecksmatrizen sind. Somit gilt für j > i:

lij =n

∑m=1

l(1)im l(2)mj =j−1

∑m=1

l(1)im l(2)mj︸︷︷︸=0

+n

∑m=j

l(1)im︸︷︷︸=0

l(2)mj .

35


Daraus folgtlij = 0 ∀j > i.

Also ist L eine untere Dreiecksmatrix.Sind L(1) und L(2) zusätzlich normiert, d.h. l(1)ii = l(2)ii = 1 ∀i = 1, . . . , n, so folgt

lii =n

∑m=1

l(1)im l(2)mj =i−1

∑m=1

l(1)im l(2)mj︸︷︷︸=0

+ l(1)ii l(2)ii +n

∑m=i+1

l(1)im︸︷︷︸=0

l(2)mj

= l(1)ii l(2)ii = 1 ∀i = 1, . . . , n.

Der Beweis der Aussage (ii) erfolgt analog.

Lemma 3.15. Es gelten die folgenden Aussagen:

(i) Ist L ∈ Rn×n eine reguläre untere Dreiecksmatrix, so ist L−1 auch eine untere Dreiecks-matrix. Ist L zusätzlich normiert, so ist auch L−1 normiert.

(ii) Ist R ∈ Rn×n eine reguläre obere Dreiecksmatrix, so ist R−1 auch eine obere Dreiecksma-trix. Ist R zusätzlich normiert, so ist auch R−1 normiert.

Beweis. Da L regulär ist, existiert die Inverse

L−1 := X =

| | |

x1 x2 . . . xn| | |

.

Dabei ist xi ∈ Rn die i-te Spalte von X. Laut Definition gilt

LL−1 = LX = I ⇔ Lxi = ei ∀i = 1, . . . , n. (3.1)

Wir zerlegen die Matrix L und die Vektoren xi und ei wie folgt:

L =

(L(i)

11 0L(i)

21 L(i)22

), xi =

(x(i)1

x(i)2

), ei =

(0

e(i)1

),

mitL(i)

11 ∈ R(i−1)×(i−1), L(i)21 ∈ R(n−i+1)×(i−1), L(i)

22 ∈ R(n−i+1)×(n−i+1)

sowiex(i)1 ∈ Ri−1, x(i)2 ∈ Rn−i+1 und ei ∈ Rn−i+1.

Dabei ist also ei der 1-te Einheitsvektor des Rn−i+1. Auf diese Weise lässt sich (3.1) wiefolgt darstellen: (

L(i)11 0

L(i)21 L(i)

22

)(x(i)1

x(i)2

)=

(0

e(i)1

).

36


Hieraus folgt

L(i)11 x(i)1 = 0

L(i)11 ist regulär⇒ x(i)1 = 0.

Das heißt: Jede i-te Spalte von X = L−1 besitzt in den ersten i− 1 Einträgen lauterNullen. Folglich ist

L−1 = X =

∗ 0... . . .∗ · · · ∗

eine untere Dreiecksmatrix. Ist L zusätzlich normiert, das heißt lii = 1 ∀i = 1, . . . , n,so ergibt sich aus der ersten Gleichung von

L(i)22 x(i)2 = e(i)1 ,

dass das i-te Element von xi eine Eins sein muss:

xi =

(0

x(i)2

)=

01∗...∗

← i-te Zeile.

Der Beweis der Aussage (ii) folgt analog.

Satz 3.16 (Eindeutigkeit der LR-Zerlegung). Sei A ∈ Rn×n regulär mit

det(A[k]) 6= 0 ∀k = 1, . . . , n.

Dann existiert genau eine LR-Zerlegung von A, d.h.

A = LR

mit einer eindeutigen unteren normierten Dreiecksmatrix L ∈ Rn×n und einer eindeutigenoberen Dreiecksmatrix R ∈ Rn×n.

Beweis. Die Existenz wurde bereits gezeigt. Sei nun

A = L1R1 und A = L2R2

mit unteren normierten Dreiecksmatrizen L1, L2 ∈ Rn×n und oberen DreiecksmatrizenR1, R2 ∈ Rn×n. Dann gilt

L1R1 = L2R2

⇒ L1 = L2R2R−11

⇒ L−12 L1 = R2R−1

1 .

Die vorherigen Lemmata liefern:

37


• L−12 L1 ist eine untere normierte Dreiecksmatrix.

• R2R−11 ist eine obere Dreiecksmatrix.

Damit ist

L−12 L1 =

1 0. . .

∗ 1

=

∗ · · · ∗

. . . ...0 ∗

= R2R−1

1

und daraus ergibt sichL−1

2 L1 = I und R2R−11 = I.

Das ist gleichbedeutend mit

L−12 = L−1

1 , also L2 = L1 und R1 = R2.

Insgesamt folgt also die Behauptung.

Bemerkung 3.17. Der Beweis zeigt, dass die LR-Zerlegung immer eindeutig ist, soferndiese für eine reguläre Matrix existiert.

Bemerkung 3.18. Beachte, dass man in jeder Iteration des Algorithmus 3.7 eineFrobenius-Matrix Li abspeichern muss. Zur Vermeidung der Abspeicherung vonFrobenius-Matrizen betrachten wir das folgende Verfahren.

Algorithmus 3.8 Gauß-Elimination ohne Pivotisierung

1: for k = 1 : (n− 1) do2: for i = (k + 1) : n do3: aik := aik/akk4: for j = (k + 1) : n do5: aij := aij − aikakj6: end for7: end for8: end for

Ist Algorithmus 3.8 durchführbar, so erhalten wir am Ende des Prozesses die folgendeMatrix:

a11 · · · a1n... . . . ...

an1 · · · ann

−→

r11 r12 · · · rnn

l21. . . ...

... . . . ...ln1 · · · ln,m−1 rnn

.

Dann istA = LR

38

3.3 LR-Zerlegung mit Pivotisierung

mit

L =

1 0

l21. . .

... . . .ln1 · · · ln,m−1 1

und R =

r11 r12 · · · rnn. . . ...

. . . ...rnn

.

Zur Durchführung des Algorithmus 3.8 benötigt man ≈ 23 n3 Rechenoperationen.


Im vorherigen Abschnitt haben wir eine notwendige und hinreichende Bedingungzur Existenz und Eindeutigkeit der LR-Zerlegung für eine reguläre Matrix A ∈ Rn×n

hergeleitet. Diese lautet

det(A[k]) 6= 0 ∀k = 1, . . . , n.

Diese Bedingung ist aber zu stark und wird oft verletzt. Ein einfaches Beispiel dafür ist

A =

(0 11 1

).

Diese Matrix hat keine LR-Zerlegung, ist jedoch regulär. Also gibt es zu jedem Vektorb ∈ R2 genau eine Lösung x ∈ R2 von

Ax = b.

Unsere Idee besteht darin, die Zeilen von A durch eine Permutationsmatrix P zuvertauschen, so dass PA eine LR-Zerlegung besitzt, d.h.

PA = LR

mit einer unteren normierten Dreiecksmatrix L ∈ Rn×n und einer oberen Dreiecksma-trix R ∈ Rn×n. Für unser Beispiel ist

P =

(0 11 0

).

Dann ergibt sich

PA =

(1 10 1

)=

(1 00 1

)

︸︷︷︸=L

(1 10 1

)

︸︷︷︸=R

als LR-Zerlegung von PA.

39


Satz 3.19. Es sei A ∈ Rn×n regulär. Dann existiert eine Permutationsmatrix P ∈ Rn×n, sodass gilt

PA = LR,

mit einer unteren normierten Dreiecksmatrix L ∈ Rn×n und einer oberen DreiecksmatrixR ∈ Rn×n.

Beweis durch Induktion nach n.

Induktionsanfang n = 1: Wähle P = L = I =(1)

und R = A. Dann gilt

PA = LR.

Induktionsannahme: Die Behauptung gelte für alle regulären Matrizen der Dimension(n− 1)× (n− 1) für ein festes n ≥ 2.

Wir zeigen nun, dass die Aussage für alle regulären Matrizen der Dimensionn × n gilt. Sei also A ∈ Rn×n regulär. Die erste Spalte von A enthält mindes-tens ein von Null verschiedenes Element, da A regulär ist. Somit existiert einePermutationsmatrix Pn ∈ Rn×n, so dass

Pn A =

(α rT

l C

)

mit α ∈ R \ {0}, r, l ∈ Rn−1 und C ∈ R(n−1)×(n−1). Wir definieren nun

B = C− 1α

lrT ∈ R(n−1)×(n−1).

Dann ist

Pn A =

(α rT

l C

)=

(1 01α l I

)(α rT

0 B

).

Ferner giltdet(Pn A) = det(Pn)︸︷︷︸

6=0

det(A)︸︷︷︸6=0

6= 0.

Also ist

0 6= det(Pn A) = det(

1 01α l I

)det

(α rT

0 B

)= 1 det

(α rT

0 B

)

= α det(B).

Damit ist auchdet(B) 6= 0.

Also ist B ∈ R(n−1)×(n−1) regulär. Laut Induktionsannahme existieren eine Per-mutationsmatrix Pn−1 ∈ R(n−1)×(n−1), eine normierte untere Dreiecksmatrix

40


Ln−1 ∈ R(n−1)×(n−1) sowie eine obere Dreiecksmatrix Rn−1 ∈ R(n−1)×(n−1), sodass

Pn−1B = Ln−1Rn−1.

Auf Grund der Orthogonalität der Permutationsmatrix Pn−1 ist dies äquivalentzu

B = PTn−1Ln−1Rn−1.

Hieraus folgt

Pn A =

(1 01α l I

)(α rT

0 PTn−1Ln−1Rn−1

)=

(1 01α l I

)(1 00 PT

n−1Ln−1

)(α rT

0 Rn−1

)

=

(1 01α l PT

n−1Ln−1

)(α rT

0 Rn−1

)

=

(1 00 PT

n−1

)(1 0

1α Pn−1l Ln−1

)(α rT

0 Rn−1

).

Somit ist

Pn A =

(1 00 PT

n−1

)(1 0

1α Pn−1l Ln−1

)(α rT

0 Rn−1

),

und daher folgt

(1 00 PT

n−1

)T

Pn A =

(1 0

1α Pn−1l Ln−1

)(α rT

0 Rn−1

).

Setzen wir nun

P :=(

1 00 PT

n−1

)T

Pn (Permutationsmatrix)

L :=(

1 01α Pn−1l Ln−1

)(untere normierte Dreiecksmatrix)

R :=(

α rT

0 Rn−1

)(obere Dreiecksmatrix),

so erhalten wir die ZerlegungPA = LR.

Mit der Zerlegung PA = LR lässt sich das lineare Gleichungssystem

Ax = b

einfach lösen:Ax = b ⇔ PAx = Pb ⇔ LRx = Pb.

41


Algorithmus 3.9 Lösung von Ax = b durch LR-Zerlegung mit Pivotisierung

(S1) Bestimme eine Zerlegung PA = LR von A.(S2) Setze c := Pb (Algorithmus 2.5).(S3) Bestimme y als Lösung von Ly = c (Vorwärtseinsetzen).(S4) Bestimme x als Lösung von Rx = y (Rückwärtseinsetzen).

Definition 3.20. Für i, j ∈ {1, . . . , n} mit j ≥ i definieren wir die Matrix

Pij :=

1. . .

10 ··· ··· 0 1... 1 0... . . . ...0 1

...1 0 ··· ··· 0

1. . .

1

∈ Rn×n.

Die Matrix Pij entsteht aus der Einheitsmatrix I ∈ Rn×n durch Vertauschung der i-tenund j-ten Zeile. Das Produkt Pij A liefert die Matrix, die aus A durch Vertauschung deri-ten und j-ten Zeile hervorgeht.

Beispiel 3.21. Wir betrachten die Matrix

A =

0 2 −1 −22 −2 4 −11 1 1 1−2 1 −2 1

.

Das erste Element der ersten Zeile der Matrix A ist Null. Also müssen wir zuerstZeilen vertauschen, um eine LR-Zerlegung finden zu können. In unserem Fall tauschenwir nun die erste mit der zweiten Zeile:

P12A =

2 −2 4 −10 2 −1 −21 1 1 1−2 1 −2 1

.

Definiere nun eine Frobenius-Matrix

L1 =

10 1−1

2 11 1

,

42


dann ist

L1P12A =

2 −2 4 −12 −1 −22 −1 1.5−1 2 0

.

Das Pivot-Element ist nun 2, also insbesondere nicht 0. Wählen also als Permutations-matrix P22 = I und erhalten

P22L1P12A =

2 −2 4 −12 −1 −22 −1 1.5−1 2 0

.

Definieren wieder eine Frobenius-Matrix

L2 =

11−1 1

12 1

.

Es ergibt sich also

L2P22L1P12A =

2 −2 4 −12 −1 −2

0 3.51.5 −1

.

Wieder ist das Pivot-Element Null, vertausche also die dritte mit der vierten Zeile:

P34L2P22L1P12A =

2 −2 4 −12 −1 −2

1.5 −10 3.5

.

Abschließend wählen wir

L3 =

1. . .

1

und erhalten das Ergebnis

L3P34L2P22L1P12A =

2 −2 4 −12 −1 −2

1.5 −13.5

=: R.

Das lineare Gleichungssystem Ax = b für dieses Beispiel lösen wir wie folgt:

Ax = b ⇒ L3P34L2P22L1P12Ax = L3P34L2P22L1P12b =: c.

Daraus folgt Rx = c und wir können x bestimmen.

43


Algorithmus 3.10 Grundversion der LR-Zerlegung mit Pivotisierung

1: Setze A(1) := A.2: for i = 1 : (n− 1) do3: Wähle aus der i-ten Spalte von A(i) ein Element a(i)ji 6= 0 mit j ≥ i.

4: Setze A(i) = Pij A(i).5: Bestimme die Frobenius-Matrix Li ∈ Rn×n

Li =

1. . .

1

−li+1,i. . .

...−ln,i 1

mit lm,i = a(i)mi /a(i)ii für m = i + 1, . . . , n.6: Setze A(i+1) := Li A(i)

7: end for8: Setze R := A(n).

Bemerkung 3.22. Algorithmus 3.10 ist durchführbar, so lange A ∈ Rn×n regulär ist.Dieser Algorithmus liefert Frobenius-Matrizen

L1, . . . , Ln−1 ∈ Rn×n

sowie Permutationsmatrizen

P1, . . . , Pn−1 ∈ Rn×n, Pi = Pij(i) mit j(i) ≥ i,

so dassLn−1Pn−1Ln−2Pn−2 · . . . · L2P2L1P1 = R.

Wir wollen nun die DarstellungPA = LR

bestimmen, und zwar aus

Ln−1Pn−1Ln−2Pn−2 · . . . · L2P2L1P1 = R.

Für i = 1, . . . , n− 1 definieren wir

L′i = Pn−1 · . . . · Pi+1LiPTi+1 · . . . · PT

n−1 (wobei Pn = I).

Laut Konstruktion ist

Li =

1. . .

1

∗ . . ....∗ 1

= I − lieTi

44


mit einem geeigneten Vektor li ∈ Rn. Somit ist

L′i = Pn−1 · . . . · Pi+1(I − lieTi )PT

i+1 · . . . · PTn−1

= Pn−1 · . . . · Pi+1PTi+1 · . . . · PT

n−1︸︷︷︸=I

−Pn−1 · . . . · Pi+1 − lieTi PT

i+1 · . . . · PTn−1

= I − Pn−1 · . . . · Pi+1lieTi PT

i+1 · . . . · PTn−1

= I − l′i eTi

mit l′i = Pn−1 · . . . · Pi+1li.

Also ist L′i eine Frobenius-Matrix. Außerdem gilt nach Definition:

L′n−1L′n−2 · . . . · L′2L′1Pn−1 · . . . · P2P1 = Ln−1Pn−1Ln−2Pn−2 · . . . · P2L2P1L1.

Dies folgt aus

L′n−1L′n−2 = Pn︸︷︷︸=I

Ln−1 PTn︸︷︷︸

=I

Pn−1Ln−2PTn−1.


(L′n−1L′n−2 · . . . · L′2L′1)Pn−1 · . . . · P2P1A = R.

Daraus folgt

Pn−1 · . . . · P2P1A = (L′n−1L′n−2 · . . . · L′2L′1)−1R

und damit

PA = LR

mit

P = Pn−1 · . . . · P2P1 und L = (L′n−1L′n−2 · . . . · L′2L′1)−1.

In jeder Iteration des Algorithmus 3.10 werden sowohl Frobenius-Matrizen als aucheine Permutationsmatrix abgespeichert. Zur Vermeidung der Abspeicherung vondiesen Matrizen verwenden wir das Gauß-Verfahren mit Pivotisierung.

45


Algorithmus 3.11 Gauß-Elimination mit Pivotisierung (vgl. Algorithmus 3.8)

1: for k = 1 : (n− 1) do. Pivotsuche

2: z = 03: for i = k : n do4: if |aik| > z then5: z = |aik|6: s = i7: end if8: end for9: P(k) = s

10: if z = 0 then11: Warnung! (Matrix A ist singulär)12: end if

. Vertauschung der k-ten und s-ten Zeile13: if k < s then14: for j = 1 : n do15: z = akj16: akj = asj17: asj = z18: end for19: end if

. Berechnung der lik20: for i = (k + 1) : n do21: aik = aik/akk22: end for

. Spaltenweise Berechnung von aufdatierten A23: for j = (k + 1) : n do24: for i = (k + 1) : n do25: aij = aij − aikakj26: end for27: end for28: if ann = 0 then29: Warnung! (Matrix A ist singulär)30: end if31: end for

46


Algorithmus 3.12 Lösung von Ax = b nach erfolgter Gauß-Elimination mit Pivotisie-rung

IMPORT P(1), . . . , P(n− 1) und A aus Algorithmus 3.11.. Berechne b := Pb (vgl. Algorithmus 3.5)

1: for k = 1 : (n− 1) do2: if k < P(k) then3: z = bk4: bk = bP(k)5: bP(k) = z6: end if7: end for

. Berechne b := L−1b (vgl. Algorithmus 3.2)8: for k = 1 : (n− 1) do9: for i = (k + 1) : n do

10: bi = bi − aikbk11: end for12: end for

. Berechne b := R−1b (vgl. Algorithmus 3.3)13: for k = n : −1 : 1 do14: bk = bk/akk15: for i = 1 : (k− 1) do16: bi = bi − aikbk17: end for18: end for

47

3.4 Cholesky-Zerlegung


In diesem Abschnitt betrachten wir ein lineares Gleichungssystem

Ax = b

mit einer symmetrischen und positiv definiten (s.p.d.) Matrix A ∈ Rn×n. Insbesondereist A regulär und somit existiert die Zerlegung

PA = LR.

Nun zeigen wir, dass A auch die folgende Zerlegung

A = LLT

mit einer unteren Dreiecksmatrix L ∈ Rn×n besitzt.Bemerkung 3.23. Die Matrix L muss nicht unbedingt normiert sein.

Definition 3.24. Sei A ∈ Rn×n regulär. Eine Zerlegung der Gestalt

A = LLT

mit einer regulären unteren Dreiecksmatrix L ∈ Rn×n heißt Cholesky-Zerlegung.

Lemma 3.25 (Notwendige Bedingung für Cholesky-Zerlegung). Sei A ∈ Rn×n regulär.Es existiere eine Cholesky-Zerlegung für A, d.h.

A = LLT

mit einer regulären unteren Dreiecksmatrix L ∈ Rn×n. Dann ist A s.p.d.

Beweis.

• AT = (LLT)T = LLT = A.

• A = LLT ⇒ xT Ax = xT LLTx = (LTx)T LTx =∥∥LTx

∥∥22 > 0

∀x ∈ Rn×n \ {0}, denn LT ist regulär.

Nun wollen wir zeigen, dass die Bedingung

A ∈ Rn×n ist s.p.d.

hinreichend für die Cholesky-Zerlegung ist.

48


Lemma 3.26. Sei A ∈ Rn×n eine symmetrische und positiv definite Matrix. Dann gilt:

(i) aii > 0 für alle i = 1, . . . , n.

(ii) A[k] ist symmetrisch und positiv definit für alle k = 1, . . . , n.

Beweis.

Zu (i): aii = eTi Aei

A p.d.> 0 ∀i = 1, . . . , n.

Zu (ii): A = AT ⇔ A[k] = A[k]T und

yT A[k]y =

(y0

)T

A(

y0

)A p.d.> 0 ∀y ∈ Rk \ {0} ∀k = 1, . . . , n.

Satz 3.27 (Cholesky-Zerlegung). Sei A ∈ Rn×n symmetrisch und positiv definit. Dannexistiert genau eine untere Dreiecksmatrix L ∈ Rn×n mit positiven Diagonalelementen, so dass

A = LLT.

Beweis durch Induktion nach n ∈N.

Induktionsanfang n=1:

A = (a11) ⇒ a11 > 0 ⇒ L = +√

a11 ⇒ A = LLT.

Induktionsannahme: Die Aussage gelte für alle s.p.d. Matrizen der Dimension (n−1)× (n− 1).

Wir zeigen die Behauptung für alle s.p.d. Matrizen der Dimension n× n. Sei A ∈ Rn×n

s.p.d. Wir zerlegen A wie folgt:

A =

(An−1 b

bT ann

)

mit An−1 ∈ R(n−1)×(n−1), b ∈ Rn−1 und ann ∈ R. Da A s.p.d. ist, gilt

• ann > 0 und

• R(n−1)×(n−1) 3 An−1 = A[n− 1] s.p.d. (Lemma).

Nach Induktionsannahme existiert genau eine untere Dreiecksmatrix Ln−1 ∈ R(n−1)×(n−1)

mit positiven Diagonalelementen, so dass

An−1 = Ln−1LTn−1.

49


Wir machen den Ansatz

L :=(

L′n−1 0cT α

)

mit L′n−1 ∈ R(n−1)×(n−1), c ∈ Rn−1 und α ∈ R. Wir suchen L′n−1, c sowie α, so dass dieZerlegung

A = LLT

gilt. Diese ist äquivalent zu

A =

(An−1 b

bT ann

)=

(L′n−1 0

cT α

)((L′n−1)

T c0 α

),

was dem Gleichungssystem

An−1 = L′n−1(L′n−1)T

L′n−1c = b

cTc + α2 = ann

entspricht. Aus diesem ergibt sich nach Induktionsannahme

L′n−1 = Ln−1

sowiec = (Ln−1)

−1b

undα2 = ann − cTc.

Es bleibt zu zeigen, dass α > 0 ist. Dazu betrachte

0 6= det(A) = det(LLT) = det(L)2 Entwicklungssatz= (α det(Ln−1))

2 = α2 det(Ln−1)2.

Hieraus folgt0 < α2 = ann − cTc,

alsoα = +

√ann − cTc > 0.

Insgesamt giltA = LLT

mit

L =

(Ln−1 0

L′n−1b +√

ann − ((L′n−1)−1b)T(L′n−1)

−1b

).

50


Bemerkung 3.28. Die Eindeutigkeit der Cholesky-Zerlegung gilt nicht mehr, wennman untere Dreiecksmatrizen mit negativen Diagonalelementen zulässt.

Algorithmus 3.13 Lösung von Ax = b mit A ∈ Rn×n s.p.d.: Cholesky-Zerlegung

(S1) Bestimme eine Cholesky-Zerlegung A = LLT.(S2) Bestimme y als Lösung von Ly = b (Vorwärtseinsetzen).(S3) Bestimme x als Lösung von LTx = y (Rückwärtseinsetzen).

Die Matrix L = (lij) aus der Cholesky-Zerlegung bestimmt man wie folgt:

A = LLT

⇔

a11 a12 a13 · · · a1na21 a22 a23 · · · a2na31 a32 a33 · · · a3n...

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

an1 an2 an3 · · · ann

=

l11l21 l22l31 l32 l33...

...... . . .

ln1 ln2 ln3 · · · lnn

l11 l12 l13 · · · l1nl22 l23 · · · l2n

l33 · · · l3n. . . ...

lnn

.

Hieraus folgt:

1. Spalte von L:

a11 = l211 ⇒ l11 = +

√a11

a21 = l21l11 ⇒ l21 = a21/l11

...an1 = ln1l11 ⇒ ln1 = an1/l11

2. Spalte von L:

a22 = l221 + l2

22 ⇒ l22 = +√

a22 − l221

a32 = l31l21 + l32l22 ⇒ l32 = (a32 − l31l21)/l22

...an2 = ln1l21 + ln2l22 ⇒ ln2 = (an2 − ln1l21)/l22

n. Spalte von L:

ann = l2n1 + l2

n2 + . . . + l2nn ⇒ lnn =

√ann − l2

n1 − l2n2 − . . . l2

n,n−1

51


Algorithmus 3.14 Cholesky-Zerlegung für s.p.d. Matrizen A ∈ Rn×n

1: for j = 1 : n do

2: ljj :=√

ajj −∑j−1m=1 l2

jm (∑0m=1 l2

jm = 0 für j = 1)3: for i = (j + 1) : n do4: lij := (aij −∑

j−1m=1 limljm)/ljj

5: end for6: end for

Zusammenfassung 3.29. Sei A ∈ Rn×n regulär:

• A besitzt eine LR-Zerlegung genau dann, wenn det(A[k]) 6= 0 ∀k = 1, . . . , n.

• Die LR-Zerlegung ist eindeutig, sofern diese existiert.

• Es existiert immer eine Permutationsmatrix P ∈ Rn×n, so dass PA eine LR-Zerlegung besitzt.

• A besitzt eine Cholesky-Zerlegung (A = LLT mit einer regulären unteren Drei-ecksmatrix L ∈ Rn×n) genau dann, wenn A s.p.d. ist.

• Die Cholesky-Zerlegung ist eindeutig, falls man nur untere Dreiecksmatrizen Lmit positiven Diagonalelementen zulässt.

• Die Cholesky-Zerlegung (Algorithmus 3.14) benötigt ca. 13 n3 Rechenoperationen,

während die LR-Zerlegung ca. 23 n3 Rechenoperationen benötigt.

52

Kapitel 4Lineares Ausgleichsproblem

Seien A ∈ Rm×n und b ∈ Rm mit m ≥ n (oft m � n) gegeben. Wir betrachten dieAufgabe

Ax = b. (4.1)

Im Falle einer quadratischen Matrix A ∈ Rn×n (m = n) lässt sich (4.1) zum Beispielmit Hilfe der LR-Zerlegung mit Pivotisierung lösen, sofern A regulär ist.Ist m > n, so hat (4.1) mehr Gleichungen als Unbekannte und somit hat (4.1) imAllgemeinen keine Lösung mehr.Beispiel 4.1. Es seien

A =

(11

)∈ R2×1 und b =

(10

)

gegeben. Die Aufgabe (11

)x =

(10

)

besitzt dann keine Lösung.

Beim linearen Ausgleichsproblem suchen wir eine optimale Approximation x ∈ Rn, sodass

Ax ≈ b.

Mathematisch formulieren wir diese Problemstellung als ein Minimierungsproblem

(P) minx∈Rn

‖Ax− b‖,

wobei ‖·‖ : Rn → R eine Vektornorm bezeichnet. In diesem Kapitel setzen wir

‖·‖ = ‖·‖2 .

Beispiel 4.2 (Hasenpopulation). Zu verschiedenen Zeitpunkten

t1, . . . , tm

werden die Daten zur Hasenpopulation in Essen

b1, . . . , bm

53

Kapitel 4 Lineares Ausgleichsproblem

gesammelt. Gesucht ist eine Funktion p : R→ R, so dass

∥∥∥∥∥∥∥

p(t1)p(t2)

...p(tm)

−

b1b2...

bm

∥∥∥∥∥∥∥

minimiert wird. Die Funktion p : R→ R liefert eine kontinuierliche Approximationzur Hasenpopulation in Essen.

t1 t2 tm

b1

b2

bm

Für dieses Beispiel ist es sinnvoll (siehe Bild), eine quadratische Approximation

p(t) ≈ c0 + c1t + c2t2

zu wählen. Gesucht sind c0, c1, c2 ∈ R. Wir kommen auf das folgende Minimierungs-problem:

min

∥∥∥∥∥∥∥

c0+c1t1+c2t21

c0+c1t2+c2t22

...c0+c1tm+c2t2

m

−

b1b2...

bm

∥∥∥∥∥∥∥.

Dies entspricht gerade dem Problem

min ‖Ax− b‖

mit

A =

1 t1 t21

1 t2 t22

......

...1 tm t2

m

∈ Rm×3, b =

b1b2...

bm

, x =

c0c1c2

.

54


Definition 4.3. Seien A ∈ Rm×n und b ∈ Rm mit m ≥ n. Die Minimierungsaufgabe

minx∈Rn

‖Ax− b‖ (P)

heißt lineares Ausgleichsproblem. Ein Vektor x∗ ∈ Rn heißt (optimale) Lösung zu (P),falls gilt

‖Ax∗ − b‖ ≤ ‖Ax− b‖ ∀x ∈ Rn.

Bemerkung 4.4. Die Optimierungsaufgabe (P) ist endlich-dimensional und nicht linear,da die Normabbildung nicht linear ist.

Bei der Untersuchung des linearen Ausgleichsproblems wollen wir uns folgendeFragen stellen:

• Hat (P) eine Lösung?

• Ist die Lösung eindeutig?

• Gibt es eine Charakterisierung für die Lösung (Optimalitätsbedingung)?

• Numerische Approximation der Lösung?

Satz 4.5 (Notwendige und hinreichende Optimalitätsbedingung für (P)). Ein Vektorx∗ ∈ Rn ist genau dann eine Lösung zu (P), wenn x∗ den folgenden Normalgleichungengenügt:

AT Ax∗ = ATb.

Beweis.

„⇒“: Es sei x∗ ∈ Rn eine Lösung zu (P). Annahme:

r := AT(Ax∗ − b) 6= 0.

Mit xt := x∗ − tr ∈ Rn folgt:

‖Axt − b‖2 = ‖Axt − b + Ax∗ − Ax∗‖2

= ‖Ax∗ − b + A(xt − x∗)‖2

= ‖Ax∗ − b− tAr‖2

= 〈Ax∗ − b− tAr, Ax∗ − b− tAr〉= 〈Ax∗ − b, Ax∗ − b〉 − 2〈tAr, Ax∗ − b〉+ 〈tAr, tAr〉= ‖Ax∗ − b‖2 − 2t〈Ar, Ax∗ − b〉+ t2 ‖Ar‖2

= ‖Ax∗ − b‖2 − 2t〈r, AT(Ax∗ − b)︸︷︷︸=r

〉+ t2 ‖Ar‖2

= ‖Ax∗ − b‖2 − 2t ‖r‖2 + t2 ‖Ar‖2 .

55


Da r 6= 0 ist, gibt es ein t ∈ R+, so dass

−2t ‖r‖2 + t2 ‖Ar‖2 < 0 ∀t ∈ (0, t).

Somit gilt‖Axt − b‖2 < ‖Ax∗ − b‖2 ∀t ∈ (0, t).

Daraus ergibt sich der Widerspruch

‖Ax∗ − b‖2 ≤ ‖Axt − b‖2 < ‖Ax∗ − b‖2 ∀t ∈ (0, t)

und es folgt die Behauptung.

„⇐“: Es gelte AT Ax∗ = ATb für ein x∗ ∈ Rn. Zu zeigen: x∗ ist eine Lösung zu (P).Dazu sei x ∈ Rn beliebig aber fest. Dann gilt:

‖Ax− b‖2 = ‖A(x− x∗) + Ax∗ − b‖2

= ‖A(x− x∗)‖2 + 2〈A(x− x∗), Ax∗ − b〉+ ‖Ax∗ − b‖2

= ‖A(x− x∗)‖2 + 2〈x− x∗, AT(Ax∗ − b)︸︷︷︸=0

〉+ ‖Ax∗ − b‖2

= ‖A(x− x∗)‖2 + ‖Ax∗ − b‖2

≥ ‖Ax∗ − b‖2 .

Somit gilt‖Ax∗ − b‖ ≤ ‖Ax− b‖ ∀x ∈ Rn.

Es folgt also, dass x∗ optimal ist.

Bemerkung 4.6. Der Satz besagt, dass (P) und AT Ax = ATb äquivalent sind.

Satz 4.7 (Existenz einer optimalen Lösung zu (P)). Das lineare Ausgleichsproblem

(P) minx∈Rn

‖Ax− b‖

besitzt eine Lösung.

Beweis. Wir haben bereits gezeigt:

x∗ ∈ Rn ist eine Lösung zu (P) genau dann, wenn AT Ax = ATb gilt.

Also genügt es zu zeigen, dass AT Ax = ATb eine Lösung x ∈ Rn hat. Laut Definitionist ATb ∈ Bild(AT). In Satz 2.42 haben wir gezeigt, dass

Bild(AT) = Bild(AT A)

gilt. Also folgt insgesamt

ATb ∈ Bild(AT) = Bild(AT A) ⇒ ∃x ∈ Rn : AT Ax = ATb.

56


Bemerkung 4.8.

(i) Obwohl die Aufgabe Ax = b im Allgemeinen keine Lösung besitzt, hat dieAufgabe (P) stets eine Lösung.

(ii) Die Lösung von (P) ist im Allgemeinen nicht eindeutig, wie folgendes Beispielzeigt.

Beispiel 4.9. Es seien

A =

000

∈ R3×1 und b =

111

∈ R3×1

gegeben. Dann ist

AT A =(0 0 0

)

000

= 0 und ATb =

(0 0 0

)

111

= 0.

Es giltAT Ax = ATb ∀x ∈ Rn

und somit löst jedes x ∈ Rn die Aufgabe (P).

Satz 4.10 (Eindeutigkeit der Lösung). Das lineare Ausgleichsproblem hat genau dann eineeindeutige Lösung, wenn

Rang(A) = n.

Beweis.

„⇐“: Es sei Rang(A) = n.

Daraus folgt, dass AT A regulär ist (Übungsaufgabe). Aus der Regularität vonAT A ergibt sich, dass die Normalgleichungen AT Ax = ATb genau eine Lösungbesitzen. Dies ist gleichbedeutend mit der Eindeutigkeit des Problems (P).

„⇒“: (P) habe genau eine Lösung. Annahme: Rang(A) < n.

Daraus folgt, dass die Normalgleichungen AT Ax = ATb unendlich viele Lösun-gen besitzen. Dies ist äquivalent dazu, dass das Problem (P) unendlich vieleLösungen besitzt. Der Widerspruch zur Voraussetzung liefert die Behauptung.

Fazit 4.11. Das lineare Ausgleichsproblem

(P) min ‖Ax− b‖, A ∈ Rm×n, b ∈ Rm

hat stets eine Lösung und die Lösung ist genau dann eindeutig, wenn Rang(A) = n.Die Lösung ist gegeben durch

AT Ax = ATb.

57

4.1 QR-Zerlegung

4.1 QR-Zerlegung

Im Folgenden seien A ∈ Rm×n und b ∈ Rm mit m ≥ n und Rang(A) = n. Somit hat(P) genau eine Lösung x∗ ∈ Rn. Wir wollen eine numerische Lösung zu (P) mittelseiner QR-Zerlegung finden.Definition 4.12 (QR-Zerlegung).

(i) Die ZerlegungA = QR, Q ∈ Rm×n, R ∈ Rn×n

mit den Eigenschaften

– Q =

| |

q1 · · · qn| |

, qi ∈ Rm, 〈qi, qj〉 = δij (Spaltenvektoren sind ortho-

normiert)

– R =

r11 · · · r1n. . . ...

rnn

∈ Rn×n

heißt reduzierte QR-Zerlegung von A.

(ii) Die ZerlegungA = QR, Q ∈ Rm×m, R ∈ Rm×n


– QTQ = I ∈ Rm×m (Q ist orthogonal)

– R =

r11 · · · r1n. . . ...

rnn0

heißt (volle) QR-Zerlegung von A.

Bemerkung 4.13.

(i): Existiert eine QR-Zerlegung von A, so existiert auch eine reduzierte QR-Zerlegung.Denn laut Definition sind

Q =[Q Q

]

R =

[R0

]

und somit

A = QR =[Q Q

] [R0

]= QR.

58

4.1 QR-Zerlegung

(ii): Existiert eine reduzierte QR-Zerlegung von A, so existiert auch eine volle QR-Zerlegung von A, indem man die Spalten von Q zu einer Orthonormalbasis desRm erweitert.

Satz 4.14 (Eindeutigkeit der reduzierten QR-Zerlegung). Sei A ∈ Rm×n mit Rang(A) =n ≤ m. Ferner seien

A = Q1R1 und A = Q2R2

zwei reduzierte QR-Zerlegungen von A. Dann existiert eine Diagonalmatrix

Rn×n 3 D = diag(λ1, . . . , λn)

mit λi = ±1, so dassQ1 = Q2D und R2 = DR1.

Beweis. Es giltQ1R1 = Q2R2. (4.2)

Laut Definition sind die Spaltenvektoren von Q1 orthonormiert. Die gleiche Aussagegilt auch für Q2. Folglich ist

QT2 Q2 = QT

1 Q1 = I ∈ Rn×n.

Daraus folgt mit (4.2)

R1 = QT1 (Q1R1) = QT

1 Q2R2

R2 = QT2 (Q2R2) = QT

2 Q1R1

und

R1R−12 = QT

1 Q2R2R−1

1 = QT2 Q1,

}(4.3)

da R1 und R2 regulär sind (Rang(A) = n). Wir setzen

D := R2R−11 ∈ Rn×n.

Dann ist D eine obere Dreiecksmatrix, da R2 und R−11 obere Dreiecksmatrizen sind.

Andererseits gilt

DT Def.= (R2R−1

1 )T (4.3)= (QT

2 Q1)T = QT

1 Q2(4.3)= R1R−1

2 . (4.4)

Also ist DT eine obere Dreiecksmatrix. Damit ist D sowohl eine obere als auch eineuntere Dreiecksmatrix. Es gilt

D = diag λ1, . . . , λn

59

4.1 QR-Zerlegung

mit λi ∈ R. Es bleibt zu zeigen, dass λi = ±1 für alle i = 1, . . . , n. Hierzu verwendenwir (4.4):

DT = R1R−12 = (R2R−1

1 )−1 Def.= D−1.

Also gilt

I = DDT =

λ1. . .

λn

λ1. . .

λn

=

λ21

. . .λ2

n

und somit folgtλi = ±1 ∀i = 1, . . . , n.

Wir wollen nun das lineare Ausgleichsproblem

min ‖Ax− b‖, A ∈ Rn×n, b ∈ Rn (P)

mittels reduzierter QR-Zerlegung lösen. Auf Grund der Voraussetzung Rang(A) = nhat (P) genau eine Lösung und die Lösung erfüllt

AT Ax = ATb.

Existiert eine reduzierte QR-Zerlegung von A

A = QR, Q ∈ Rm×n, R ∈ Rn×n,

so folgtAT Ax = (QR)TQRx = RTQTQRx = RT Rx

undATb = (QR)Tb = RTQTb.

Folglich istAT Ax = b ⇔ RT Rx = RTQTb.

Da Rang(A) = n ist, ist R regulär. Somit ist auch RT regulär. Insgesamt gilt

AT Ax = ATb ⇔ Rx = QTb.

Algorithmus 4.1 Lösung von (P) mittels einer reduzierten QR-Zerlegung

(S1) Bestimme eine reduzierte QR-Zerlegung von A.(S2) Setze c := QTb.(S3) Löse Rx = c mittels Rückwärtseinsetzen.

Alternativ lässt sich (P) mittels voller QR-Zerlegung lösen. Sei A = QR mit Q ∈ Rm×m

und R ∈ Rm×n eine volle QR-Zerlegung. Das heißt:

60

4.1 QR-Zerlegung

(i) QTQ = I ∈ Rm×m

(ii) R =

(R0

)mit R =

r11 · · · r1n. . . ...

rnn

∈ Rn×n.

Wir setzen

c := QTb ∈ Rm

und zerlegen

c =(

cc

)mit c ∈ Rn und c ∈ Rm−n.

Folglich gilt

‖Ax− b‖ = ‖QRx− b‖ =∥∥∥QRx−QQTb

∥∥∥ =∥∥∥Q(Rx−QTb)

∥∥∥ .

Wegen

‖Qy‖2 = (Qy)TQy = yTQTQy = yTy = ‖y‖2 ∀y ∈ Rn,

folgt

‖Ax− b‖ =∥∥∥Rx−QTb

∥∥∥ =

∥∥∥∥(

R0

)x−

(cc

)∥∥∥∥ .

Insgesamt gilt

min ‖Ax− b‖ ⇔ min∥∥∥∥(

R0

)x−

(cc

)∥∥∥∥.

Da aber∥∥∥∥(

R0

)x−

(cc

)∥∥∥∥2 ‖·‖=‖·‖2=

∥∥Rx− c∥∥2

+ ‖c‖2

ist, kommen wir zur Folgerung:

x löst (P) ⇔ x löst min∥∥∥∥(

R0

)x−

(cc

)∥∥∥∥ ⇔ x löst min∥∥Rx− c

∥∥2 ⇔ Rx = c.

Das ist genau Schritt 3 im Algorithmus 4.1.

61

4.2 Gram-Schmidt-Orthogonalisierung


Im Folgenden sei A ∈ Rm×n mit Rang(A) = n ≤ m. Folglich hat (P) genau eineLösung, die wir mit Hilfe des Algorithmus 4.1 bestimmen wollen.In diesem Abschnitt zeigen wir die Existenz einer (reduzierten) QR-Zerlegung durchdie Gram-Schmidt-Orthogonalisierungsmethode.Aus A = QR mit Q ∈ Rm×n und R ∈ Rn×n erhalten wir (für eine bessere Lesbarkeitschreiben wir im Folgenden q(j) = q(j) für die Spaltenvektoren von Q)

A =

| |

a(1) · · · a(n)

| |

=

| |

q(1) · · · q(n)

| |

r11 · · · r1n. . . ...

rnn

.

Dies ist äquivalent zu dem System

a(1) = r11q(1)

a(2) = r12q(1) + r22q(2)

a(3) = r13q(1) + r23q(2) + r33q(3)

...

a(n) = r1nq(1) + . . . + rnnq(n),

also

a(j) =j

∑i=1

rijq(i) ∀j = 1, . . . , n.

Hieraus ergibt sich eine Formel für die orthonormalen Vektoren q(j) und die Zahlen rij.

Für j = 1:

a(1) = r11q(1) ⇒ q(1) =a(1)

r11und r11 = ±

∥∥∥a(1)∥∥∥ ,

damit∥∥∥q(1)

∥∥∥ = 1.

Für j = 2:

q(2) =1

r22(a(2) − r12q(1)) und r22 = ±

∥∥∥a(2) − r12q(1)∥∥∥ ,

damit∥∥∥q(2)

∥∥∥ = 1 und r12 bestimmen wir aus der Orthogonalitätsbedingung

0 = 〈q(1), q(2)〉 = 1r22

(〈q(1), a(2)〉 − r12 〈q(1), q(2)〉︸︷︷︸=1

) =1

r22(〈q(1), a(2)〉 − r12).

Also istr12 = 〈q(1), a(2)〉.

62


Für j ≥ 3:

q(j) =1rjj

(a(j) −

j−1

∑i=1

rijq(i))

und

rjj = ±∥∥∥∥∥a(j) −

j−1

∑i=1

rijq(i)∥∥∥∥∥ und

rij = 〈q(i), a(j)〉 ∀i = 1, . . . , j− 1.

Algorithmus 4.2 Gram-Schmidt-Verfahren

1: r11 :=∥∥∥a(1)

∥∥∥2: q(1) := a(1)/r113: for j = 2 : n do4: q(j) := a(j)

5: for i = 1 : (j− 1) do6: rij := 〈q(i), a(j)〉7: q(j) := q(j) − rijq(i)

8: end for9: rjj :=

∥∥∥q(j)∥∥∥

10: q(j) := q(j)/rjj11: end for

Satz 4.15 (Existenz der QR-Zerlegung). Es sei A ∈ Rm×n mit Rang(A) = n ≤ m. Dannexistiert eine QR-Zerlegung von A

A = QR

mit Q ∈ Rm×n und R ∈ Rn×n (siehe Definition). Somit existiert auch eine volle QR-Zerlegungvon A.

Beweis. Es bleibt zu zeigen, dass das Gram-Schmidt-Verfahren durchführbar ist. Wirmüssen also zeigen:

rjj 6= 0 ∀j = 1, . . . , n.

Annahme: Es gibt ein j ∈ {1, . . . , n}, so dass rii 6= 0 ∀i = 1, . . . , j− 1, aber rjj = 0.Folglich

0 = rjj =

∥∥∥∥∥a(j) −j−1

∑i=1

rijq(i)∥∥∥∥∥

⇔ a(j) =j−1

∑i=1

rijq(i)

⇔ a(j) ∈ Span{q(1), . . . , q(j−1)}.

63


Laut Konstruktion ist

Span{q(1), . . . , q(j−1)} = Span{a(1), . . . , a(j−1)}.

Somit gilta(j) ∈ Span{q(1), . . . , q(j−1)} = Span{a(1), . . . , a(j−1)}.

Mit anderen Worten lässt sich der j-te Spaltenvektor von A als Linearkombination vonden Spaltenvektoren a(1), . . . , a(j−1) darstellen. Dies ist ein Widerspruch zu Rang(A) =n. Also folgt die Behauptung.

Bemerkung 4.16. Das Gram-Schmidt-Verfahren ist numerisch instabil (vgl. Übung).Aufgrund der endlichen Rechengenauigkeit des Computers geht die Orthgonalität vonq(j) schnell verloren.

Wir wollen nun eine modifizierte (stabile) Variante des Gram-Schmidt-Verfahrensherleiten, indem wir weitere Hilfsvektoren p(j) einführen, um das Überschreiben derVektoren q(j) im Algorithmus 4.2 zu vermeiden.

1: p(1) = a(1)

2: r11 :=∥∥∥p(1)

∥∥∥3: q(1) := p(1)/r114: for j = 2 : n do5: for i = 1 : (j− 1) do6: rij := 〈q(i), a(j)〉7: end for8: p(j) := a(j) −∑

j−1i=1 rijq(i)

9: rjj :=∥∥∥p(j)

∥∥∥10: q(j) := p(j)/rjj11: end for

Die Berechnung für p(j) lautet:

p(j) = a(j) −j−1

∑i=1

rijq(i) = a(j) −j−1

∑i=1〈q(i), a(j)〉q(i) = a(j) −

j−1

∑i=1

q(i) 〈q(i), a(j)〉︸︷︷︸=(q(i))T a(j)

=

(I −

j−1

∑i=1

q(i)(q(i))T

)a(j).

Insgesamt gilt

p(j) =

(I −

j−1

∑i=1

q(i)(q(i))T

)a(j).

64


Wegen der Orthogonalitätsbedingung für q(i) gilt(

I −j−1

∑i=1

q(i)(q(i))T

)= (I − q(j−1)(q(j−1))T) · . . . · (I − q(2)(q(2))T)(I − q(1)(q(1))T)

und somit lässt sich p(j) auch wie folgt berechnen:

1: pj,1 := a(j)

2: for i = 1 : (j− 1) do3: rij := 〈q(i), pj,i〉4: pj,i+1 := pj,i − rijq(i)

5: end for6: p(j) = pj,j

Hier gilt:

pj,1 = a(j)

i = 1 : pj,2 = a(j) − 〈q(1), a(j)〉q(1)

= a(j) − q(1)(q(1))Ta(j)

= (I − q(1)(q(1))T)a(j)

i = 2 : pj,3 = (I − q(2)(q(2))T)(I − q(1)(q(1))T)a(j)

i = j− 1 : pj,j = (I − q(j−1)(q(j−1))T) · . . . · (I − q(2)(q(2))T)(I − q(1)(q(1))T)a(j).

Speichert man nun pj,i nicht explizit ab und überschreibt stattdessen einen Vektor q(j),so kommen wir auf das modifizierte Gram-Schmidt-Verfahren.

Algorithmus 4.3 Modifiziertes Gram-Schmidt-Verfahren

1: r11 :=∥∥∥a(1)

∥∥∥2: q(1) := a(1)/r113: for j = 2 : n do4: q(j) := a(j)

5: for i = 1 : (j− 1) do6: rij := 〈q(i), q(j)〉7: q(j) := q(j) − rijq(i)

8: end for9: rjj :=

∥∥∥q(j)∥∥∥

10: q(j) := q(j)/rjj11: end for

65

4.3 Householder-Spiegelungen

Bemerkung 4.17. Algorithmen 4.2 und 4.3 benötigen

4m(n− 1)n

2≈ 2mn2

Rechenoperationen.


In diesem Abschnitt befassen wir uns mit Householder-Spiegelungen zur Konstruktioneiner vollen QR-Zerlegung von A, das heißt

A = QR

mit einer Orthogonalmatrix Q ∈ Rm×m und einer Matrix R ∈ Rm×n der Gestalt

R =

r11 · · · r1n. . . ...

rnn0

.

Im Folgenden nehmen wir wieder an, dass

Rang(A) = n ≤ m.

Somit existiert eine volle QR-Zerlegung

A = QR ⇔ QT A = R.

Idee: Die Householder-Spiegelung liefert ein Produkt von einfachen Orthogonalmatri-zen Qi ∈ Rm×m mit

QT = Ql · . . . ·Q2Q1, l = min{n, m− 1}.

Jedes Qi sorgt dafür, dass in der i-ten Spalte von A geeignete Nullen erzeugt werden(vgl. Frobenius-Matrizen Li).

A =

∗ ∗ ∗∗ ∗ ∗∗ ∗ ∗∗ ∗ ∗

→ Q1A =

∗ ∗ ∗0 ∗ ∗0 ∗ ∗0 ∗ ∗

→ Q2Q1A =

∗ ∗ ∗0 ∗ ∗0 0 ∗0 0 ∗

→ Q3Q2Q1A =

∗ ∗ ∗0 ∗ ∗0 0 ∗0 0 0

=: R.

66


Definition 4.18 (Householder-Spiegelung). Eine Matrix H ∈ Rn×n heißt Householder-Spiegelung, falls gilt

H = I − 2uuT

für einen Vektor u ∈ Rn mit ‖u‖ = 1.

Lemma 4.19. Jede Householder-Spiegelung H ∈ Rn×n ist symmetrisch und orthogonal.

Beweis.

• HT = (I − 2uuT)T = I − 2uuT = H.

• HT H = HH = (I − 2uuT)(I − 2uuT) = I − 4uuT + 4uuTuuT = I − 4uuT +4uuT = I, da uTu = ‖u‖2 = 1 ist.

Bemerkung 4.20. Sei u ∈ Rn mit ‖u‖ = 1 und

H = I − 2uuT ∈ Rn×n (Householder-Spiegelung).

Geometrisch beschreibt die Abbildung

v 7→ Hv, Rn → Rn

eine Spiegelung des Vektors v an einer durch den Nullpunkt gehenden Ebene mit demNormalenvektor u.

u

v

Hv

Wir wollen zunächst untersuchen, wie man durch eine Householder-Spiegelung einenVektor v ∈ Rn \ {0} in ein Vielfaches von e1 ∈ Rn abbilden kann:

Hv = ρe1.

Sei alsoH = I − 2uuT ∈ Rn×n

67


mit u ∈ Rn, so dass ‖u‖ = 1. Weiter sei v ∈ Rn. Wir bestimmen nun u ∈ Rn undρ ∈ R, so dass

Hv = (I − 2uuT)v = v− 2uuTv != ρe1. (4.5)

Wegen‖v‖ = ‖Hv‖ = ‖ρe1‖ = |ρ|

istρ = ±‖v‖ .

Wir setzen nunγ = 2〈u, v〉,

und nehmen an, dass γ 6= 0 ist. Aus (4.5) erhalten wir

γu Def.= 2u〈u, v〉 = 2uuTv = v− ρe1

und daher

u =1γ(v− ρe1)

sowieγ = ±‖v− ρe1‖ , da ‖u‖ = 1.

Im Folgenden wählen wir γ = + ‖v− ρe1‖. Insgesamt haben wir das folgende Resultat.

Lemma 4.21. Sei v ∈ Rn \ {0} und

ρ := ±‖v‖ , γ = + ‖v− ρe1‖ , u =1γ(v− ρe1).

Ist γ 6= 0, so gilt für die Householder-Matrix H = I − 2uuT ∈ Rn×n:

Hv = ρe1.

Bemerkung 4.22. Aus numerischen Gründen ist die Wahl

ρ =

{−‖v‖ falls v1 ≥ 0+ ‖v‖ falls v1 < 0

besser geeignet, um den Effekt der Auslöschung bei der Berechnung von

v− ρe1

zu vermeiden. Bei dieser Wahl ist die Bedingung γ 6= 0 stets gesichert.

Das obige Lemma ist eine wichtige Grundlage für die volle QR-Zerlegung mittelsHouseholder-Spiegelungen. Dazu betrachten wir ein Beispiel.

68


Beispiel 4.23. Es sei

A =

0 1 00 0 −

√2

2 1 00 1 2

.

1. Konstruktion von Q1 (vgl. Lemma 4.21):

v =

0020

.

Damit ergibt sich

‖v‖ = 2ρ = −‖v‖ = −2 (siehe Bemerkung 4.22)

γ = ‖v− ρe1‖ =∥∥∥∥( 0

020

)+

( 2000

)∥∥∥∥ =

∥∥∥∥( 2

020

)∥∥∥∥ =√

8

Es folgt

u =1γ(v− ρe1) =

1√8

( 2020

)

und

H1 = I − 2uuT = I − 218

2020

(2 0 2 0

)

=

11

11

−

1 0 1 00 0 0 01 0 1 00 0 0 0

.


H1 =

0 0 −1 00 1 0 0−1 0 0 0

0 0 0 1

.

Setze Q1 = H1 und erhalte

Q1A =

0 0 −1 00 1 0 0−1 0 0 0

0 0 0 1

0 1 00 0 −

√2

2 1 00 1 2

=

−2 −1 00 0 −

√2

0 −1 00 1 2

.

69


Kontrolle:

Hv = ρe1 =

−2000

.

2. Konstruktion von Q2:

v =

0−1

1

.

Damit ergibt sich

‖v‖ =√

2ρ = −‖v‖ = −

√2

γ = ‖v− ρe1‖ =∥∥∥∥( 0−11

)+

(√2

00

)∥∥∥∥ =

∥∥∥∥(√

2−11

)∥∥∥∥ = 2

Es folgt

u =1γ(v− ρe1) =

12

√

2−11

und

H2 = I − 2uuT = I − 214

√

2−1

1

(√2 −1 1

)

=

11

1

− 1

2

2 −√

2√

2−√

2 1 −1√2 −1 1

=

0 12

√2 −1

2

√2

12

√2 1

2 −12

−12

√2 1

212

.

Setze

Q2 =

1 0 0 000 H20

=

1 0 0 00 0 1

2

√2 −1

2

√2

0 12

√2 1

212

0 −12

√2 1

212

und

Q2Q1A =

1 0 0 00 0 1

2

√2 −1

2

√2

0 12

√2 1

212

0 −12

√2 1

212

2 −1 00 0 −

√2

0 −1 00 1 2

=

−2 −1 00 −

√2 −

√2

0 0 00 0 2

.

70


Kontrolle:

Hv = ρe1 =

−√

200

.

3. Konstruktion von Q3:

v =

(02

).

Damit ergibt sich

‖v‖ = 2ρ = −‖v‖ = −2

γ = ‖v− ρe1‖ =∥∥∥(

22

)∥∥∥ =√

8

Es folgt

u =1γ(v− ρe1) =

1√8

(22

)

und

H3 = I − 2uuT = I − 218

(22

) (2 2

)

=

(1

1

)−(

1 11 1

)

=

(0 −1−1 0

).

Setze

1 0 0 00 1 0 00 0

H30 0

=

1 0 0 00 1 0 00 0 0 −10 0 −1 0

.

Somit gilt:

Q3Q2Q1A =

1 0 0 00 1 0 00 0 0 −10 0 −1 0

−2 −1 00 −

√2 −

√2

0 0 00 0 2

=

−2 −1 00 −

√2 −

√2

0 0 −20 0 0

=

(R0

)=: R.

71


Kontrolle:

Hv = ρe1 =

(−20

).

Insgesamt haben wir eine QR-Zerlegung gefunden:

QT : = Q3Q2Q1

A = QR ⇔ QT A = R.

Die Orthogonalmatrix lautet

Q =

0 −12

√2 1

212

0 0 12

√2 −1

2

√2

−1 0 0 00 −1

2

√2 −1

2 −12

.

Kontrolle:

QR =

0 −12

√2 1

212

0 0 12

√2 −1

2

√2

−1 0 0 00 −1

2

√2 −1

2 −12

−2 −1 00 −

√2 −

√2

0 0 00 0 2

=

0 1 00 0 −

√2

2 1 00 1 2

= A.

Das Verfahren verläuft wie folgt:Seien Orthogonalmatrizen Qi, . . . , Q1 ∈ Rm×n bereits berechnet mit

QiQi−1 · . . . ·Q1A =

(Ri Ci0 Mi

)

für eine obere Dreiecksmatrix Ri ∈ Ri×i und geeignete Matrizen Ci ∈ Ri×(n−i) undMi ∈ R(m−i)×(m−i).Mit v ∈ Rm−i bezeichnen wir den ersten Spaltenvektor von Mi. Hieraus bestimmenwir Hi+1 ∈ R(m−i)×(m−i) und setze

Qi+1 =

(Ii 00 Hi+1

)∈ Rm×m

mit Ii ∈ Ri×i. Dann ist

Qi+1Qi · . . . ·Q1A =

(Ii 00 Hi+1

)(Ri Ci0 Mi

)=

(Ri+1 Ci+1

0 Mi+1

)

mit einer neuen oberen Dreiecksmatrix Ri+1 ∈ R(i+1)×(i+1) und geeigneten MatrizenCi+1 und Mi+1.Nach spätestens l = min{m− 1, n} Schritten ist dann

QlQl−1 · . . . ·Q1A =

(R0

)= R

72


mit R ∈ Rn×n obere Dreicksmatrix.

Algorithmus 4.4 Householder-Spiegelung

1: Setze l := min{m− 1, n}.2: for i = 1 : l do3: v := A(i : m, i)4: if v1 ≥ 0 then5: ρ = −‖v‖6: else7: ρ = + ‖v‖8: end if9: u(i) := v− ρe1

10: u(i) := u(i)/∥∥∥u(i)

∥∥∥11: A(i : m, i : m) := A(i : m, i : m)− 2u(i)(u(i))T A(i : m, i : m)12: end for

Bemerkung 4.24. Algorithmus 4.4 überschreibt die Matrix A mit der oberen Dreiecks-matrix R und speichert u(i) zur Berechnung von Q.

73

Kapitel 5CG-Verfahren

Dieses Kapitel befasst sich mit dem bekannten CG-Verfahren (CG = Conjugate Gradi-ent) zur Lösung von

Ax = b

mit einer symmetrischen und positiv definiten Matrix A ∈ Rn×n.

5.1 Herleitung des CG-Verfahrens

Das CG-Verfahren ist eine iterative Methode zur Lösung von Ax = b. Die Herlei-tung des CG-Verfahrens erfolgt über ein quadratisches Optimierungsproblem. Dazubetrachten wir ein quadratisches Funktional

f : Rn → R, f (x) :=12

xT Ax− bTx

und das Optimierungsproblemminx∈Rn

f (x).

Lemma 5.1. Es sei A ∈ Rn×n symmetrisch und positiv definit, b ∈ Rn, sowie

f : Rn → R, f (x) :=12

xT Ax− bTx.

Dann ist x∗ ∈ Rn genau dann eine Lösung des linearen Gleichungssystems

Ax = b,

wenn x∗ das Minimierungsproblemminx∈Rn

f (x)

löst.

Beweis. Da A ∈ Rn×n positiv definit ist, ist die Inverse A−1 ebenso positiv definit. Wirdefinieren ein Hilfsfunktional

g : Rn → R, g(x) := f (x) +12

bT A−1b.

75


Die Funktion g unterscheidet sich von f nur um einen konstanten Term 12 bT A−1b und

somit giltx∗ löst min

x∈Rnf (x) ⇔ x∗ löst min

x∈Rng(x).

Ferner ist

12(Ax− b)T A−1(Ax− b) =

12(Ax− b)T(x− A−1b)

=12(xT Ax− xT AA−1b− bTx + bT A−1b)

=12(xT Ax− 2bTx + bT A−1b)

=12

xT Ax− bTx +12

bT A−1b

= f (x) +12

bT A−1b = g(x).

Folglich ist

x∗ löst minx∈Rn

f (x) ⇔ x∗ löst minx∈Rn

12(Ax− b)T A−1(Ax− b).

Da A−1 positiv definit ist, gilt

12(Ax− b)T A−1(Ax− b) ≥ 0 ∀x ∈ Rn

und12(Ax− b)T A−1(Ax− b) = 0 ⇔ Ax = b.

Insgesamt gilt alsoAx∗ = b ⇔ x∗ löst min

x∈Rnf (x).

Bemerkung 5.2. Die Aussage

Ax∗ = b ⇔ x∗ löst minx∈Rn

f (x)

folgt auch unmittelbar aus den klassischen notwendigen und hinreichenden Optimali-tätsbedingungen für differenzierbare Minimierungsaufgaben.

Definition 5.3. Es sei A ∈ Rn×n symmetrisch und positiv definit.

(i): Zwei Vektoren x, y ∈ Rn \ {0} heißen A-orthogonal (A-konjugiert), wenn gilt

xT Ay = 0.

76


(ii): Es seien d0, d1, . . . , dn−1 ∈ Rn \ {0} paarweise A-orthogonale Vektoren, d.h.

(di)T Adj = 0 ∀i, j ∈ {0, . . . , n− 1} : i 6= j.

Das Verfahren der sukzessiven 1D-Minimierung entlang d0, d1, . . . , dn−1 ist wiefolgt definiert:

xk+1 = xk + tkdk

f (xk + tkdk) = mint∈R f (xk + tdk)

k = 0, . . . , n− 1

mit x0 ∈ Rn Startvektor.

Bemerkung 5.4. Die paarweise A-orthogonalen Vektoren d0, d1, . . . , dn−1 ∈ Rn \ {0}lassen sich zum Beispiel mittels Gram-Schmidt bestimmen.

Lemma 5.5. Es sei A ∈ Rn×n symmetrisch und positiv definit. Ferner seien d0, d1, . . . , dn−1 ∈Rn \ {0} paarweise A-orthogonale Vektoren. Dann sind d0, d1, . . . , dn−1 linear unabhängig.

Beweis. Annahme: d0, d1, . . . , dn−1 ∈ Rn \ {0} wären nicht linear unabhängig. Dannexistiert ein Index m ∈ {0, . . . , n− 1}, so dass

dm =n−1

∑j=0j 6=m

λjdj

mit λj ∈ R, j ∈ {0, . . . , n− 1} und es gibt (mindestens) einen Index i 6= m mit λi 6= 0.Folglich:

0 = (di)T Adm = (di)T An−1

∑j=0j 6=m

λjdj =n−1

∑j=0j 6=m

λj (di)T Adj︸︷︷︸=0 für i 6=j

= λi(di)T Adi,

also0 = (di)T Adi ⇔ di = 0.

Dies ist ein Widerspruch zur Annahme, also folgt die Behauptung.

Satz 5.6. Es seien A ∈ Rn×n symmetrisch und positiv definit, b ∈ Rn, und d0, d1, . . . , dn−1 ∈Rn \ {0} paarweise A-orthogonale Vektoren. Dann konvergiert das Verfahren der sukzessiven1D-Minimierung entlang d0, d1, . . . , dn−1 nach spätestens n Schritten gegen die Lösungx∗ ∈ Rn von

Ax = b.

Für jedes k = 0, . . . , n− 1 gilt mit gk = Axk − b

tk = −(gk)Tdk

(dk)T Adk

77


und(gk+1)Tdj = 0 ∀j = 0, . . . , k.

Beweis. Laut Konstruktion löst jedes tk ∈ R die Minimierungsaufgabe

mint∈R

ϕ(t) := f (xk + tdk).

Aus der notwendigen und hinreichenden Optimalitätsbedingung

ϕ′(tk) = 0

folgt

tk = −(gk)Tdk

(dk)T Adk . (5.1)

Beachte, dass (dk)T Adk 6= 0 ist, da dk 6= 0 und A positiv definit ist. Für allek = 0, . . . , n− 1 gilt:

(gk+1)Tdk = (Axk+1 − b)Tdk = (Axk + tk Adk − b)Tdk

= (Axk − b + tk Adk)Tdk

= (gk)Tdk + tk(dk)T Adk.

Also folgt

(gk+1)Tdk (5.1)= 0. (5.2)

Laut Voraussetzung ist (dj)T Adi = 0 für alle i 6= j und somit

(gi+1 − gi)Tdj = (Axi+1 − b− (Axi − b))Tdj

= (Axi+1 − Axi)Tdj

= (Axi − Axi + ti Adi)Tdj

= ti(di)T Adj = 0 ∀i 6= j. (5.3)

Aus (5.2) und (5.3) erhalten wir

(gk+1)Tdj = (gj+1)Tdj︸︷︷︸=0 (5.2)

+k

∑i=j+1

(gi+1 − gi)Tdj︸︷︷︸

=0 (5.3)

= 0 ∀j = 0, . . . , k.

Es bleibt zu zeigen, dass das Verfahren nach spätestens n Schritten die Lösung vonAx = b liefert. Wir wissen aber

(gn)Tdj = 0 ∀j = 0, . . . , n− 1.

Daher istgn = 0,

78


da d0, d1, . . . , dn−1 linear unabhängig sind. Es folgt

0 = gn = Axn − b

und insgesamtAxn = b.

Bemerkung 5.7. Das Verfahren der sukzessiven 1D-Minimierung entlang d0, d1, . . . , dn−1

konvergiert stets gegen die Lösung von Ax = b. Der Nachteil dieses Verfahrens bestehtdarin, dass man vorab die paarweisen A-orthogonalen Vektoren d0, d1, . . . , dn−1 ∈Rn \ {0} bestimmen muss. Das kann numerisch teuer sein.

Beim CG-Verfahren bestimmt man die Vektoren d0, d1, . . . , dn−1 ∈ Rn \ {0} nicht vorab.Diese werden sukzessiv im Verfahren mit der Bedingung

∇ f (xk)Tdk < 0

erzeugt, d.h. dk ∈ Rn \ {0} ist eine Abstiegsrichtung für f in xk. Beachte:

∇ f (xk) = Axk − b = g(k).

Wir starten mitd0 := −∇ f (x0) = −(Ax0 − b) = −g0.

Nun gehen wir davon aus, dass d0, d1, . . . , dl (l ∈ {0, . . . , n− 2}) mit

(di)T Adj = 0 ∀i 6= j

vorliegen. Analog zum Verfahren der sukzessiven 1D-Minimierung setzen wir

xl+1 = xl + tldl,

tl = −(gl)Tdl

(dl)T Adl .

Istgl+1 Def.

= Axl+1 − b = 0,

so sind wir fertig. Ist gl+1 6= 0, dann betrachten wir den Ansatz

dl+1 := −gl+1 +l

∑i=0

βlid

i (A)

mit geeigneten Koeffizienten βli ∈ R, so dass

(dl+1)T Adj = 0 ∀j = 0, . . . , l.

Aus(di)T Adj = 0 ∀i 6= j = 0, . . . , l

folgt unmittelbar

βlj =

(gl+1)T Adj

(dj)T Adj ∀j = 0, . . . , l.

79


Bemerkung 5.8. Wir sehen leicht, dass dl+1 aus dem Ansatz (A) eine Abstiegsrichtungfür f in xl+1 ist, denn:

∇ f (xl+1)Tdl+1 = (Axl+1 − b)dl+1 Def.= (gl+1)Tdl+1

(A)= (gl+1)T(−gl+1 +

l

∑i=0

βlid

i)

= −‖gl+1‖22 +

l

∑i=0

βli (gl+1)Tdi︸︷︷︸

=0

= −‖gl+1‖22

Die gleiche Aussage gilt auch für die alten Iterationen:

(gj)T︸︷︷︸=∇ f (xj)

dj = −‖gj‖22 ∀j = 0, . . . , l + 1.

Lemma 5.9. Es gilt:

βlj = 0 ∀j = 0, . . . , l − 1 und

βll =‖gl+1‖2

2‖gl‖2

2.

Beweis. Für j = 0, . . . , l ist

(gl+1)Tgj (A)= (gl+1)T(

j−1

∑i=0

βj−1i di − dj) =

j−1

∑i=0

βj−1i (gl+1)Tdi − (gl+1)Tdj = 0. (5.4)

Weiter ist

βlj =

(gl+1)T Adj

(dj)T Adj =(gl+1)T(gj+1 − gj)

tj(dj)T Adj , (5.5)

weil:(gj+1 − gj)

Def.= (Axj+1 − b− Axj + b) = Axj + tj Adj − Axj = tj Adj.

Aus (5.4) und (5.5) folgtβl

j = 0 ∀j = 0, . . . , l − 1

und

βll =

(gl+1)T Adl

(dl)T Adl =

∥∥gl+1∥∥2

2∥∥gl∥∥2

2

,

denn:

tl(dl)T Adl = (gl+1 − gl)Tdl = (gl+1)Tdl︸︷︷︸

=0

−(gl)Tdl = −(gl)Tdl Bem.= ‖gl‖2

2.

80


Algorithmus 5.1 CG-Verfahren zur Lösung von Ax = b mit A s.p.d.

1: Wähle x0 ∈ Rn.2: Setze g0 = Ax0 − b und d0 = −g0.3: for k = 0, 1, 2, . . . do

4: tk = − (gk)Tdk

(dk)T Adk =‖gk‖2

2(dk)T Adk

5: xk+1 = xk + tkdk

6: gk+1 = Axk+1 − b = gk + tk Adk

7: βk =‖gk+1‖2

2

‖gk‖22

8: dk+1 = −gk+1 + βkdk

9: if∥∥gk+1

∥∥2 = 0 then

10: STOP11: end if12: end for

Bemerkung 5.10. In der Praxis ersetzt man das Abbruchkriterium

‖gk+1‖2 = 0

durch‖gk+1‖2

‖g0‖2≤ ε

mit ε > 0 „klein“.

Satz 5.11. Es sei A ∈ Rn×n symmetrisch und positiv definit, b ∈ Rn, und

f : Rn → R, f (x) :=12

xT Ax− bTx.

Dann liefert der Algorithmus 5.1 (CG-Verfahren) nach spätestens n Schritten die Lösungx∗ ∈ Rn von

Ax = b ⇔ minx∈Rn

f (x).

Ist m ∈ {0, . . . , n} die kleinste Zahl mit

xm = x∗,

dann gilt

(dk)T Adj = 0 ∀k = 1, . . . , m ∀j = 0, . . . , k− 1 (A-Orthogonalität)

(gk)Tgj = 0 ∀k = 1, . . . , m ∀j = 0, . . . , k− 1

(gk)Tdj = 0 ∀k = 1, . . . , m ∀j = 0, . . . , k− 1

(gk)Tdk = −‖gk‖22 ∀k = 1, . . . , m (Abstiegsrichtung, siehe Bemerkung).

81

5.2 Charakterisierung des CG-Verfahrens

Beweis. Es bleibt nur noch zu zeigen, dass Algorithmus 4.1 nach spätestens n Schrittenkonvergiert. Ist gm = 0 für m < n, dann sind wir fertig. Andernfalls erzeugt der Algo-rithmus paarweise A-orthogonale Vektoren d0, d1, . . . , dn−1 ∈ Rn \ {0}, und somit istder Algorithmus 5.1 äquivalent zum Verfahren der sukzessiven 1D-Minimierungsaufgabeentlang d0, d1, . . . , dn−1. Folglich ist xn = x∗ wegen der Konvergenz der sukzessiven1D-Minimierung entlang d0, d1, . . . , dn−1.


Im Folgenden sei A ∈ Rn×n symmetrisch und positiv definit und b ∈ Rn. Wir untersu-chen weiter das CG-Verfahren.Definition 5.12. Sei A ∈ Rn×n symmetrisch und positiv definit. Dann definieren wirdie A-gewichtete Norm auf Rn wie folgt:

‖·‖RnA

: Rn → R, ‖x‖RnA

:=√

xT Ax.

Bemerkung 5.13. Da A symmetrisch und positiv definit ist, definiert

〈·, ·〉RnA

: Rn ×Rn → Rn, 〈x, y〉RnA

:= xT Ay

ein Skalarprodukt auf Rn und

‖x‖RnA=√〈x, x〉Rn

A

ist eine Norm auf Rn.

Wir verwenden folgende Notation:

gk = Axk − b (Gradient)

rk = b− Axk = −gk (Residuum).

Definition 5.14 (Krylov-Raum). Der folgende Raum

Kk := Kk(r0, A) := Span{r0, Ar0, . . . , Ak−1r0}

heißt der von dem Anfangsresiduum r0 und der Matrix A aufgespannte k-te Krylov-Raum.

Lemma 5.15. Bricht der Algorithmus 5.1 nicht vorzeitig ab, so gelten die folgenden Identitäten:

Span{d0, d1, . . . , dk−1} = Span{g0, g1, . . . , gk−1} = Span{r0, r1, . . . , rk−1}= Span{r0, Ar0, . . . , Ak−1r0}= Kk(r0, A) = Kk,

für k = 1, 2, . . ..

82


Aufgrund der Definition ist

xk = xk−1 + tk−1dk−1 = (xk−2 + tk−2dk−2) + tk−1dk−1

= . . . = x0 +k−1

∑i=0

tidi.

Daher ist

xk = x0 +k−1

∑i=0

tidi ∈ x0 + Span{d0, . . . , dk−1}

Lemma 5.15= x0 +Kk.

Mit anderen Worten ist die k-te Iterierte xk ∈ Rn des CG-Verfahrens stets ein Elementdes affinen Raums x0 +Kk.Lemma 5.16. Die k-te Iterierte xk ∈ Rn des CG-Verfahrens löst das folgende Minimierungs-problem

minx∈x0+Kk

f (x)

mit dem Zielfunktional f (x) = 12 xT Ax− bTx.

Satz 5.17. Es sei A ∈ Rn×n symmetrisch und positiv definit, b ∈ Rn, und x∗ := A−1b.Dann sind die folgenden Aussagen äquivalent:

(i) xk löst minx∈x0+Kk

f (x),

(ii) xk löst minx∈x0+Kk

‖x− x∗‖2Rn

A(xk minimert den Fehler ‖x− x∗‖2

RnA

auf x0 +Kk.),

(iii) xk löst minx∈x0+Kk

‖Ax− b‖2Rn

A−1(xk minimert das Residuum ‖Ax− b‖2

RnA−1

auf x0 +

Kk.),

wobei ‖y‖RnA=√

yT Ay und ‖y‖RnA−1

=√

yT A−1y.

Beweis.

(i)⇔ (ii): Laut Definition gilt

‖x− x∗‖2Rn

A= (x− x∗)T A(x− x∗) = xT Ax− 2xT Ax∗ + (x∗)T Ax∗

Ax∗=b= xT Ax− 2bTx + bTx∗

De f .= 2 f (x) + bTx∗.

Somit gilt:

xk löst minx∈x0+Kk

f (x) ⇔ xk löst minx∈x0+Kk

‖xk − x∗‖2Rn

A.

83

5.3 Konvergenzanalyse des CG-Verfahrens

(ii)⇔ (iii): Laut Definition gilt

‖x− x∗‖2Rn

A= (x− x∗)T A(x− x∗) = (x− x∗)T AA−1A(x− x∗)

A=AT= (A(x− x∗))T A−1(A(x− x∗))

Ax∗=b= (Ax− b)T A−1(Ax− b)

De f .= ‖Ax− b‖2

RnA−1

.

Somit gilt:

xk löst minx∈x0+Kk

‖xk − x∗‖2Rn

A⇔ xk löst min

x∈x0+Kk

‖Ax− b‖2Rn

A−1.

Folgerung 5.18. Jede k-te Iterierte xk des CG-Verfahrens minimiert sowohl den Fehler‖x− x∗‖Rn

Aals auch das Residuum ‖Ax− b‖Rn

A−1auf x0 +Kk.


Ist A ∈ Rn×n symmetrisch und positiv definit, so wissen wir, dass das CG-Verfahrennach spätestens n Schritten die exakte Lösung x∗ ∈ Rn von

Ax = b

liefert. In diesem Abschnitt wollen wir den Fehler

‖x− x∗‖RnA

für jede CG-Iterierte xk ∈ Rn analysieren. Dazu betrachten wir zuerst die Tschebyscheff-Polynome.Definition 5.19 (Tschebyscheff-Polynom). Es sei n ∈N∪ {0}. Dann heißt

Tn : [−1, 1]→ R, Tn(x) := cos(n arccos(x))

n-tes Tschebyscheff-Polynom auf dem Intervall [−1, 1].

Satz 5.20 (Eigenschaften der Tschebyscheff-Polynome). Es gelten die folgenden Aussagen:

(i) Tn+1(x) = 2xTn(x)− Tn−1(x) für n = 1, 2, . . ., wobei T0 ≡ 1 und T1(x) = x.

(ii) Tn ist ein Polynom vom Grad n.

(iii) Tn+1 besitzt n+1 Nullstellen bei

xi = cos(

2i + 12n + 2

π

), i = 0, 1, . . . , n.

84


(iv) Tn+1 besitzt die folgende Darstellung:

Tn+1(x) = 2nn

∏i=0

(x− xi), n = 0, 1, . . . .

(v) Es gilt |Tn+1(x)| ≤ 1 für alle x ∈ [−1, 1] und

Tn+1(xi) = (−1)i bei xi = cos(

iπn + 1

)für i = 0, 1, . . . , n + 1.

Beweis.

zu (i): T0(x)De f .= cos(0 arccos(x)) = cos(0) = 1 ∀x ∈ [−1, 1]. Außerdem gilt

T1(x)De f .= cos(1 arccos(x)) = x ∀x ∈ [−1, 1].

Wir benutzen das Additionstheorem

cos(a + b) = cos(a) cos(b)− sin(a) sin(b)

mit a := n arccos(x) und b := arccos(x), und somit erhalten wir

Tn+1(x)De f .= cos((n + 1) arccos(x))= cos(n arccos(x)) cos(arccos(x))− sin(n arccos(x)) sin(arccos(x))= Tn(x)x− sin(n arccos(x)) sin(arccos(x)).

Nun benutzen wir ein anderes Additionstheorem

sin(

a + b2

)sin(

a− b2

)=

12(cos(b)− cos(a))

mit a := (n + 1) arccos(x) und b := (n− 1) arccos(x), und somit:

sin(n arccos(x)) sin(arccos(x)) =12(cos((n− 1) arccos(x))

− cos((n + 1) arccos(x)))

=12(Tn−1(x)− Tn+1(x)) .

Somit erhalten wir

Tn+1(x) = xTn(x)− 12(Tn−1(x)− Tn+1(x))


Tn+1(x) = 2xTn(x)− Tn−1(x).

85


zu (ii): Dies folgt unmittelbar aus (i).

zu (iii): Dies folgt unmittelbar aus der Definition

Tn+1(x) = cos((n + 1) arccos(x)).

zu (iv): Dies folgt aus (iii) und (i).

zu (v): Dies folgt unmittelbar aus der Definition.

Aus der Rekursivformel{

Tn+1(x) = 2xTn(x)− Tn−1(x), x ∈ [−1, 1], n = 1, 2, . . .T0 ≡ 1, T1(x) = x

ergibt sich

T2(x) = 2x2 − 1

T3(x) = 4x3 − 3x

T4(x) = 8x4 − 8x2 + 1.

−1 −0.8−0.6−0.4−0.2 0 0.2 0.4 0.6 0.8 1−1

−0.8

−0.6

−0.4

−0.2

0

0.2

0.4

0.6

0.8

1x2

T2T3T4

Definition 5.21. Das Tschebyscheff-Polynom auf R ist definiert durch die Rekursivfor-mel

Tn+1(x) = 2xTn(x)− Tn−1(x)T0 ≡ 1, T1(x) = xfür x ∈ R, n = 1, 2, . . .

86


Lemma 5.22. Das Tschebyscheff-Polynom auf R \ [−1, 1] besitzt die folgende Darstellung

Tn(x) =12

[(x +

√x2 − 1)n + (x−

√x2 − 1)n

]

für alle x ∈ R mit |x| > 1.

Beweis. Wir verwenden

cos(nθ) + i sin(nθ) = (cos(θ) + i sin(θ))n ,cos(nθ)− i sin(nθ) = (cos(θ)− i sin(θ))n .

Addition beider Gleichungen liefert:

cos(nθ) =12[(cos(θ) + i sin(θ))n + (cos(θ)− i sin(θ))n] .

Wir setzen θ := arccos(x) für x ∈ [−1, 1] und erhalten mit sin2(θ) + cos2(θ) = 1:

cos(n arccos(x)) =12[(cos(θ) + i sin(θ))n + (cos(θ)− i sin(θ))n]

=12

[(x + i

√1− x2)n + (x− i

√1− x2)n

].

Daher ist

Tn(x) =12

[(x + i

√1− x2)n + (x− i

√1− x2)n

]∀x ∈ [−1, 1].

Da sowohl Tn als auch die rechte Seite ein Polynom vom Grad n ist und diese beidenPolynome für alle x ∈ [−1, 1] übereinstimmen, so gilt diese Gleichung auch für allex ∈ R. Mit

i√

1− x2 =√

x2 − 1 ∀x ∈ R mit |x| > 1

folgt

Tn(x) =12

[(x +

√x2 − 1)n + (x−

√x2 − 1)n

]∀x ∈ [−1, 1].

Folgerung 5.23. Das Tschebyscheff-Polynom genügt der folgenden Ungleichung

∣∣∣∣Tn

(α + 1α− 1

)∣∣∣∣ ≥12

(√α + 1√α− 1

)n

für alle α > 1.

87


Beweis. Sei α > 1. Dann gilt für x := α+1α−1 > 1

x±√

x2 − 1 =

√α± 1√α∓ 1

und somit folgt

Tn

(α + 1α− 1

)Lemma=

12

[(√α + 1√α− 1

)n

+

(√α− 1√α + 1

)n]

≥ 12

(√α + 1√α− 1

)n

.

Folgerung 5.24. Es gilt

Tn(−x) = (−1)nTn(x) für alle n ∈N.

Beweis mittels vollständiger Induktion:

Induktionsanfang n = 1:

T1(−x) Def.= −x Def.

= (−1)1T1(x).

Induktionsannahme: Die Aussage gelte für alle Indizes k = 1, . . . , n mit n ∈N.Nun zeigen wir, dass die Aussage auch für n + 1 gilt. Lauf Definition ist aber

Tn+1(−x) Def.= (−2x)Tn(−x)− Tn−1(−x)

IA= (−2x)(−1)nTn(x)− (−1)n−1Tn−1(x)

= (−1)n+12xTn(x)− (−1)n−1(−1)2Tn−1(x)

= (−1)n+1 [2xTn(x)− Tn−1(x)]

= (−1)n+1Tn+1(x).

Satz 5.25 (Fehlerabschätzung für das CG-Verfahren). Es sei A ∈ Rn×n symmetrischund positiv definit, b ∈ Rn, sowie x∗ = A−1b. Dann erfüllt die durch Algorithmus 5.1(CG-Verfahren) erzeugte k-te Iterierte die folgende Fehlerabschätzung:

‖xk − x∗‖RnA≤ 2

(√α− 1√α + 1

)k

‖x0 − x∗‖RnA

für alle k = 1, 2, . . . ,

mit√

α− 1√α + 1

∈ [0, 1), und α := cond2(A) ∈ [1, ∞) der Spektralkondition von A.

88


Bemerkung 5.26. Es ist

cond2(A) = ‖A‖2 ‖A−1‖2

(=

λmax(A)

λmin(A), falls A symmetrisch und positiv definit

).

Beweis. Wir beweisen die Aussage für α > 1. Die Aussage für α = 1 erhalten wir späterals Folgerung.Wir setzen

ek := xk − x∗ (e : error)

für den Fehler in der k-ten Iteration. Aus

r0 De f .= b− Ax0 b=Ax∗

= A(x∗ − x0) = −Ae0

folgtAjr0 = −Aj+1e0.

AlsoKk

De f .= Span{r0, Ar0, . . . , Ak−1r0} = Span{Ae0, . . . , Ake0}.

Somit gilt:

‖ek‖RnA= ‖xk − x∗‖Rn

A

Satz= min

x∈x0+Kk

‖x− x∗‖RnA= min

z∈Kk‖x0 − x∗ + z‖Rn

A

= minµ1,...,µk∈R

‖x0 − x∗ +k

∑j=1

µj Aje0‖RnA

= minµ1,...,µk∈R

‖e0 +k

∑j=1

µj Aje0‖RnA

= minp∈∏k

p(0)=1

‖p(A)e0‖RnA

, (5.6)

wobei ∏k den Raum aller Polynome vom Grad k bezeichnet. Da A symmetrisch ist,existiert eine Orthonormalbasis {v1, . . . , vn} aus Eigenvektoren von A zu Eigenwertenλi ∈ R+, i = 1, . . . , n. Beachte, dass λi > 0 für alle i = 1, . . . , n gilt, da A symmetrischund positiv definit ist. Somit hat e0 die Darstellung

e0 =n

∑i=1

γivi

und folglich

‖e0‖2Rn

A

De f .= 〈e0, Ae0〉 =

⟨n

∑i=1

γivi, A( n

∑j=1

γjvj)⟩

Avj=λjvj

=

⟨n

∑i=1

γivi,n

∑j=1

γjλjvj

⟩

〈vi,vj〉=δij=

n

∑i=1

γ2i λi. (5.7)

89


Analog gilt:

‖p(A)e0‖2Rn

A

De f .=⟨

p(A)e0, A(p(A)e0)⟩=

⟨n

∑i=1

γi p(A)vi, A( n

∑j=1

γj p(A)vj)⟩

p(A)vj=p(λj)vj

=

⟨n

∑i=1

γi p(λi)vi, A( n

∑j=1

γj p(λj)vj)⟩

Avj=λjvj

=

⟨n

∑i=1

γi p(λi)vi,n

∑j=1

γj p(λj)λjvj

⟩

〈vi,vj〉=δij=

n

∑i=1

γ2i p(λi)

2λi. (5.8)

Zusammen ergibt sich

‖ek‖RnA

(5.6)= min

p∈∏kp(0)=1

‖p(A)e0‖RnA

(5.8)= min

p∈∏kp(0)=1

(n

∑i=1

γ2i p(λi)

2λi

) 12

≤ minp∈∏k

p(0)=1

max1≤j≤n

|p(λj)|(

n

∑i=1

γ2i λi

) 12

(5.7)= min

p∈∏kp(0)=1

max1≤j≤n

|p(λj)|‖e0‖Rn A.

Es bleibt nur noch zu zeigen, dass

minp∈∏k

p(0)=1

max1≤j≤n

|p(λj)| ≤ 2

(√2− 1√2 + 1

)k

gilt. Dazu betrachten wir ein spezielles Polynom:

pk(λ) :=Tk(F(λ))Tk(F(0))

,

mit

F(λ) =2λ− (λmax + λmin)

λmax − λmin.

Da Tk ein Polynom vom Grad k ist, gilt pk ∈ ∏k und auf Grund der Konstruktion gilt

90


p(0) = 1. Somit gilt

minp∈∏k

p(0)=1

max1≤j≤n

|p(λj)| ≤ max1≤j≤n

|pk(λj)| = max1≤j≤n

|Tk(F(λj))||Tk(F(0))|

≤ 1 1|Tk(F(0))|

=

∣∣∣∣Tk

(−λmax + λmin

λmax − λmin

) ∣∣∣∣−1

Folgerung=

∣∣∣∣Tk

(λmax + λmin

λmax − λmin

) ∣∣∣∣−1

α= λmaxλmin=

∣∣∣∣Tk

(α + 1α− 1

) ∣∣∣∣−1

Folgerung≤

(12

(√α + 1√α− 1

)k)−1

= 2(√

α− 1√α + 1

)k

.

Insgesamt gilt also

‖ek‖Rn A ≤ minp∈∏k

p(0)=1

max1≤j≤n

|p(λj)|‖e0‖Rn A ≤ max1≤j≤n

|pk(λj)|‖e0‖Rn A ≤ 2(√

α− 1√α + 1

)k

‖e0‖Rn A.

Bemerkung 5.27. Ist A schlecht konditioniert (cond2(A)� 1), so ist

(√α− 1√α + 1

)≈ 1− ε

mit ε > 0 „klein“. In diesem Fall kann das CG-Verfahren sehr langsam konvergieren.Ist umgekehrt

cond2(A) ≈ 1 + ε,

so ist(√

α− 1√α + 1

)nahe Null und somit konvergiert das CG-Verfahren in diesem Fall

recht schnell.

Folgerung 5.28. Besitzt A ∈ Rn×n insgesamt m ≤ n verschiedene Eigenwerte, so brichtdas CG-Verfahren nach spätestens m Schritten mit der Lösung von Ax = b ab.

1|Tk(x)| ≤ 1 ∀x ∈ [−1, 1]

91

5.4 Das präkonditionierte CG-Verfahren

Beweis. Im Beweis vom obigen Satz haben wir bereits gezeigt, dass


p(0)=1

max1≤j≤n

|p(λj)|‖e0‖Rn A

mit ek = xk − x∗ und λ1, . . . , λn Eigenwerte von A gilt. Hat A insgesamt nur m ≤ nverschiedene Eigenwerte, so definieren wir

pm(x) := κm

∏i=1

(x− λi)

mit κ ∈ R, so dass pm(0) = 1 gilt. Folglich ist


p(0)=1

max1≤j≤n

|p(λj)|‖e0‖Rn Ap=pm≤ 0.

Dies bedeutet‖xm − x∗‖ = 0,

alsoxm = x∗.

Fazit 5.29. Die Konditionszahl der Matrix A ist wichtig für die Konvergenzgeschwin-digkeit des CG-Verfahrens.


Im Folgenden sei A ∈ Rn×n symmetrisch und positiv definit, und b ∈ Rn. UnsereKonvergenzanalyse zeigt, dass das CG-Verfahren „schneller“ konvergiert, wenn dieMatrix A gut konditioniert ist. Aus diesem Grund ist es naheliegend, das lineareGleichungssystem Ax = b umzuschreiben in der folgenden Gestalt:

C−1Ax = C−1b (5.9)

mit einer regulären Matrix C ∈ Rn×n, so dass

cond2(C−1A)De f .= ‖C−1A‖2‖A−1C‖2 < ‖A‖2 ‖A−1‖2 = cond2(A).

Beachte, dass C−1A im Allgemeinen nicht symmetrisch ist, so dass das CG-Verfahrennicht direkt auf (5.9) angewendet werden kann. Deshalb benutzen wir die Umformu-lierung

C−1AC−TCTx = C−1b, mit C−T := (C−1)T = (CT)−1

⇔ Ax = b, (5.10)

92


mit

A := C−1AC−T symmetrisch und positiv definit,x := CTx,b := C−1b.

(5.11)

Nun ist A ∈ Rn×n symmetrisch und positiv definit, und das auf (5.10) angewandteCG-Verfahren lautet dann wie folgt (siehe Algorithmus 5.1):

1: Wähle x0 ∈ Rn.2: Setze g0 = Ax0 − b und d0 = −g0.3: for k = 0, 1, 2, . . . do

4: tk =‖gk‖2

2(dk)T Adk

5: xk+1 = xk + tkdk

6: gk+1 = Axk+1 − b = gk + tk Adk

7: βk =‖gk+1‖2

2

‖gk‖22

8: dk+1 = −gk+1 + βkdk

9: end for

Setzen wir

xk := C−T xk und

dk := C−T dk,

so ergibt sich aus (5.11):

(i): gk = Axk − b = C−1AC−TCTxk − C−1b = C−1(Axk − b) = C−1gk,

(ii): tk(i)= (gk)T gk

(dk)T Adk = (gk)T(C−TC−1gk)(CTdk)TC−1 AC−TCTdk = (gk)T(C−TC−1gk)

(dk)T Adk ,

(iii): CTxk+1 = xk+1 = xk + tkdk = CTxk + tkCTdk,

(iv): C−1gk+1 (i)= gk+1 = gk + tk Adk (i)

= C−1gk + tkC−1AC−TCTdk = C−1(gk + tk Adk),

(v): βk =‖gk+1‖2

2

‖gk‖22

(i)= (gk+1)T(C−TC−1gk+1)

(gk)T(C−TC−1gk),

(vi): CTdk+1 = dk+1 = −gk+1 + βkdk = −C−1gk+1 + βkCTdk.

Multiplizieren wir die Aufdatierungsformeln (iii) und (vi) von links mit C−T und dieFormel (iv) von links mit C, so kommen wir auf das folgende CG-Verfahren für (5.10):

1: Wähle x0 ∈ Rn.

93


2: Setze g0 = Ax0 − b und d0 = −C−TC−1g0.3: for k = 0, 1, 2, . . . do4: tk =

(gk)T(C−TC−1gk)(dk)T Adk . (vgl. (ii))

5: xk+1 = xk + tkdk . (vgl. (iii))6: gk+1 = gk + tk Adk . (vgl. (iv))

7: βk =(gk+1)T(C−TC−1gk+1)(gk)T(C−TC−1gk)

. (vgl. (v))

8: dk+1 = −C−TC−1gk+1 + βkdk . (vgl. (vi))9: end for

Setzen wir P := CCT, so erhalten wir das präkonditionierte CG-Verfahren.

Algorithmus 5.2 Präkonditioniertes CG-Verfahren

(S1) Wähle x0 ∈ Rn und P ∈ Rn×n symmetrisch und positiv definit.(S2) Setze g0 = Ax0 − b und bestimme d0 aus Pd0 = −g0.(S3) Setze z0 = d0.

1: for k = 0, 1, 2, . . . do2: tk =

(gk)Tzk

(dk)T Adk

3: xk+1 = xk + tkdk

4: gk+1 = gk + tk Adk

5: Bestimme zk+1 aus Pzk+1 = −gk+1.

6: βk =(gk+1)Tzk+1

(gk)Tzk

7: dk+1 = zk + βkdk

8: if∥∥gk+1

∥∥2 = 0 then

9: STOP10: end if11: end for

Bemerkung 5.30. Die Matrix P wird als Präkonditionierer bezeichnet. VerschiedeneStrategien zur effizienten Konstruktion von P findet man in der Literatur. Beispieledafür sind

• die unvollständige Cholesky-Zerlegung,

• die unvollständige LR-Zerlegung,

• die splittingbasierte Zerlegung.

94

Kapitel 6Nichtlineares Gleichungssystem

In diesem Kapitel beschäftigen wir uns mit einem nichtlinearen Gleichungssystem derGestalt

F(x) = 0

mit einer nichtlinearen Funktion F : Rn → Rn.

6.1 Differentialrechnung

Definition 6.1. Es sei U ⊂ Rn offen und F : Rn ⊃ U → Rm.

(i) F heißt in x ∈ U differenzierbar, wenn eine lineare Abbildung F′(x) : Rn → Rm

existiert, so dass für alle y ∈ U gilt:

F(y) = F(x) + F′(x)(y− x) + r(x, y),

und das Restglied r(x, y) genügt der Bedingung:

limy→x

r(x, y)‖y− x‖2

= 0.

(ii) Die lineare Abbildung F′(x) : Rn → Rm heißt Ableitung von F in x ∈ U.

(iii) F heißt (auf U) differenzierbar, falls F in allen Punkten x ∈ U differenzierbar ist.

Bemerkung 6.2.

(i) Ist F : Rn ⊃ U → Rm in x ∈ U differenzierbar, so ist die AbleitungF′(x) : Rn → Rm eindeutig bestimmt.

(ii) Die Ableitung F′(x) : Rn → Rm ist laut Definition linear und somit stetig, dajede lineare Abbildung auf endlichdimensionalen Vektorräumen stetig ist.

(iii) Das Restglied ist wiederum eine Funktion r(x, ·) : Rn ⊃ U → Rm. Diese hat dieexplizite Darstellung

r(x, y) = F(y)− F(x)− F′(x)(y− x) ∀y ∈ U,

95


und

limy→x

r(x, y)‖y− x‖2

= 0 ⇔ limy→x

‖r(x, y)‖2‖y− x‖2

= 0 ⇔ limy→x

ri(x, y)‖y− x‖2

= 0

für alle i = 1, . . . , m mit r(x, y) =

r1(x, y)...

rm(x, y)

.

Geometrische Interpretation: Im Falle F : R2 → R ist der Graph von

y 7→ F(x) + F′(x)(y− x)

die Tangentialebene im Punkt F(x) an den Graphen von F.Bemerkung 6.3 (Komponentenweise Differentiation).

F =

F1...

Fm

: Rn ⊃ U → Rm

mit Komponentenfunktionen Fi : Rn ⊃ U → R. Dann ist F genau dann in x ∈ Udifferenzierbar, wenn alle Fi in x differenzierbar sind. In diesem Fall gilt

F′(x)h =

F′1(x)h...

F′m(x)h

∀h ∈ Rn.

Definition 6.4 (Richtungsableitung). Es sei U ⊂ Rn offen und F : Rn ⊃ U → Rm.Dann heißt F in x ∈ U in Richtung h ∈ Rn richtungsdifferenzierbar, wenn der Limes

∂F∂h

(x) := limt→0

F(x + th)− F(x)t

∈ Rm

existiert. Der Limes ∂F∂h (x) ∈ Rm heißt Richtungsableitung von F in x in Richtung h.

Existiert ∂F∂h (x) ∈ Rm für alle Richtungen h ∈ Rn, so heißt F in x richtungsdifferenzier-

bar.

Bemerkung 6.5. Im Falle h = ei heißt

∂F∂ei

(x)

die i-te partielle Ableitung von F in x. Andere Notationen:

∂F∂ei

(x) =∂F∂xi

(x) = ∂iF(x).

96


Lemma 6.6. Es sei U ⊂ Rn offen und F : Rn ⊃ U → Rm. Ist F in x ∈ U differenzierbar, soist F in x richtungsdifferenzierbar mit

∂F∂h

(x) = F′(x)h ∀h ∈ Rn.

Beweis. Sei h ∈ Rn beliebig aber fest. Dann ist

x + th ∈ U

für hinreichend kleines |t|, weil U offen ist. Somit gilt für hinreichend kleines |t|:F(x + th)− F(x)

t=

F(x) + F′(x)(th) + r(x, x + th)− F(x)t

= F′(x)h +r(x, x + th)

t.

Aberr(x, x + th)

t=

r(x, x + th)‖x + th− x‖2︸︷︷︸→0 für t→0

‖th‖2t︸︷︷︸

=±‖h‖2

→ 0 für t→ 0.

Und somit

limt→0

F(x + th)− F(x)t

= F′(x)h.

Lemma 6.7. Ist F : Rn ⊃ U → Rm (U ⊂ Rn offen) in x ∈ U differenzierbar, so ist F in xstetig.

Beweis. Sei {xk}∞k=1 ⊂ Rn mit

limk→∞‖xk − x‖2 = 0.

Dann gilt:

‖F(xk)− F(x)‖2 =∥∥F′(x)(xk − x) + r(x, xk)

∥∥2 → 0 für k→ ∞,

denn F′(x) : Rn → Rm ist stetig und limk→∞‖r(x, xk)‖2 = 0.

Fazit 6.8. Ist F : Rn ⊃ U → Rm in x ∈ U differenzierbar, so ist

• F in x stetig und

• F in x richtungsdifferenzierbar mit ∂F∂h (x) = F′(x)h ∀h ∈ Rn.

Bemerkung 6.9. Ist F : Rn ⊃ U → Rm in x ∈ U richtungsdifferenzierbar, so muss Fin x nicht notwendigerweise differenzierbar sein (F muss auch nicht unbedingt in xstetig sein). Ein Beispiel dazu ist:

F : R2 → R, F(x1, x2) =

{1 falls x1 = x2

2 und x2 6= 00 sonst

.

Diese Abbildung ist im Nullpunkt richtungsdifferenzierbar, ist dort aber nicht stetig.

97


Beispiel 6.10. Jede lineare Abbildung F : Rn → Rm ist differenzierbar mit der Ablei-tungF′(x) = F ∀x ∈ Rn, denn

F(y) = F(x) + F(y− x)︸︷︷︸=F′(x)(y−x)

+ 0︸︷︷︸=r(x,y)

∀x, y ∈ Rn.


F : Rn → Rn−1, F(x) = x1

x2...

xn

.

Diese Abbildung ist auf Rn differenzierbar und für jedes x ∈ Rn gilt

F′(x)y = x1

y2...

yn

+ y1

x2...

xn

∀y ∈ Rn.

Beweis. Offenbar ist F′(x) : Rn → Rn−1, wie oben definiert, linear. Es bleibt also zuzeigen, dass das Restglied

r(x, y) = F(y)− F(x)− F′(x)(y− x)

die Bedingung

limy→x

r(x, y)‖y− x‖2

= 0

erfüllt (siehe Definition). Laut Definition gilt aber:

r(x, y)‖y− x‖2

=F(y)− F(x)− F′(x)(y− x)

‖y− x‖2

=

y1

( y2...

yn

)− x1

( x2...

xn

)− x1

(y2−x2

...yn−xn

)− (y1 − x1)

( x2...

xn

)

‖y− x‖2

=

y1

(y2−x2

...yn−xn

)− x1

(y2−x2

...yn−xn

)

‖y− x‖2

=

(y1 − x1)

(y2−x2

...yn−xn

)

‖y− x‖2.

98


Daraus folgt

‖r(x, y)‖2‖y− x‖2

=

∥∥∥∥∥(y1 − x1)

(y2−x2

...yn−xn

)∥∥∥∥∥2

‖y− x‖2≤ ‖y− x‖2 ‖y− x‖2

‖y− x‖2= ‖y− x‖2


limy→x

r(x, y)‖y− x‖2

= 0.

Bemerkung 6.12. Ist F : Rn ⊃ U → Rm in x ∈ U differenzierbar, so ist laut DefinitionF′(x) : Rn → Rm linear. Nach linearer Algebra kann somit F′(x) : Rn → Rm eindeutigdurch eine Darstellungsmatrix dargestellt werden:

F′(x)y =[F′(x)e1 · · · F′(x)en

]y

Lemma=

[∂F∂e1

(x) · · · ∂F∂en

(x)]

y

=

∂F1∂e1

(x) · · · ∂F1∂en

(x)... . . . ...

∂Fm∂e1

(x) · · · ∂Fm∂en

(x)

y1...

yn

∀y ∈ Rn.

Definition 6.13. Ist F : Rn ⊃ U → Rm in x ∈ U differenzierbar, so heißt die Darstel-lungsmatrix

Rm×n 3(

∂Fi

∂ej(x)

)Andere Notation

=

(∂Fi

∂xj(x)

)=(∂jFi(x)

)

Jacobi-Matrix bzw. Funktionalmatrix von F in x.

Beispiel 6.14. Ist F : Rn ⊃ U → R (m = 1) in x ∈ U differenzierbar, so ist

F′(x)y = ∇F(x)Ty = (∂1F(x), . . . , ∂nF(x))

y1...

yn

.

Mit anderen Worten: Im Falle m = 1 stimmt die Jacobi-Matrix mit ∇F(x)T überein.

Der folgende Satz liefert ein überaus wichtiges Kriterium für die Differenzierbarkeiteiner Funktion F : Rn ⊃ U → Rm.

99


Lemma 6.15 (Differenzierbarkeit und partielle Ableitungen). Es sei U ⊂ Rn offen undF : Rn ⊃ U → Rm. Existieren alle partiellen Ableitungen von F auf U und sind dieAbbildungen

∂iF : Rn ⊃ U → Rm, x 7→ ∂iF(x) ∀i = 1, . . . , n

stetig, so ist F auf U differenzierbar, und die Jacobi-Matrix von F in jedem x ∈ U ist gegebendurch (∂jFi(x)) ∈ Rm×n.

Beweis. Wir zeigen die Aussage für m = 1, d.h. F : Rn ⊃ U → R. Die Aussage fürm > 1 folgt danach mittels komponentenweiser Differentiation. Sei x ∈ U beliebig aberfest. Wir definieren:

F′(x) : Rn → Rm, F′(x)y :=n

∑i=1

yi ∂iF(x)︸︷︷︸∈R

.

Da U ⊂ Rn offen ist, existiert ein ε > 0 mit

Bε(x) := {y ∈ Rn | ‖y− x‖2 < ε} ⊂ U.

Sei y ∈ Bε(x) mit y 6= x. Dann gilt:

F(y)− F(x) = F(y1, . . . , yn)− F(x1, . . . , xn)

= F(y1, . . . , yn)− F(x1, y2, . . . , yn)

+ F(x1, y2, . . . , yn)− F(x1, x2, y3 . . . , yn)

+ F(x1, x2, y3 . . . , yn)− F(x1, x2, x3, y4 . . . , yn)

± . . .+ F(x1, x2, . . . , xn−1, yn)− F(x1, . . . , xn).

In jeder Zeile verwenden wir nun den Mittelwertsatz aus der Analysis I:

F(y)− F(x) = ∂1F(ξ1, y2, . . . , yn)(y1 − x1)

+ ∂2F(x1, ξ2, y3, . . . , yn)(y2 − x2)

+ . . .+ ∂nF(x1, x2, . . . , xn−1, ξn)(yn − xn)

mit ξi zwischen yi und xi. Folglich:

r(x, y)‖y− x‖2

=F(y)− F(x)− F′(x)(y− x)

‖y− x‖2

=∑n

i=1 (yi − xi) [∂iF(x1, . . . , xi−1, ξi, yi+1, . . . , yn)− ∂iF(x)]‖y− x‖2

.

100


Daraus folgt

r(x, y)‖y− x‖2

=n

∑i=1

yi − xi

‖y− x‖2[∂iF(x1, . . . , xi−1, ξi, yi+1, . . . , yn − ∂iF(x1, . . . , xn)]︸︷︷︸

→0 für y→x

→ 0 für y→ x,

da die Abbildungen ∂iF : Rn → R laut Voraussetzung stetig sind. Beachte dabei, dassaus der Konvergenz von y gegen x die Konvergenz von ξ gegen x folgt. Insgesamtergibt sich damit die Behauptung.


F : R2 → R2, F(x) =(

sin(x1) sin(x2)cos(x2)

).

Alle partiellen Ableitungen existieren:

∂1F1(x) = cos(x1) sin(x2),∂1F2(x) = 0,∂2F1(x) = sin(x1) cos(x2),∂2F2(x) = − sin(x2).

Die Abbildungenx 7→ ∂1F(x), x 7→ ∂2F(x)

sind stetig. Also ist F : R2 → R2 differenzierbar mit

F′(x)y =

(cos(x1) sin(x2) sin(x1) cos(x2)

0 − sin(x2)

)(y1y2

)∀y ∈ R2.

Definition 6.17. Mit L (Rn, Rm) bezeichnen wir den Raum aller linearen Funktionenzwischen Rn und Rm, d.h.

L (Rn, Rm) = { f : Rn → Rm linear}.

Eine differenzierbare Abbildung F : Rn ⊃ U → Rm heißt (auf U) stetig differenzierbar,falls die Abbildung

F′ : Rn ⊃ U → L (Rn, Rm), x 7→ F′(x)

stetig ist. Dies ist äquivalent dazu, dass die Abbildung

(∂jFi) : Rn ⊃ U → Rm×n, x 7→ (∂jFi(x))

stetig ist.

101

6.2 Newton-Verfahren

Korollar 6.18. Es sei U offen und F : Rn ⊃ U → Rm. Existieren alle partiellen Ableitungenvon F auf U und sind die Ableitungen

∂iF : Rn ⊃ U → Rm, x 7→ ∂iF(x) ∀i = 1, . . . , n

stetig, so ist F stetig differenzierbar.

Beweis. Die Aussage folgt aus dem vorherigen Satz und der Definition der stetigenDiffernzierbarkeit.


Im Folgenden untersuchen wir das nichtlineare Gleichungssystem

F(x) = 0

mit einer stetig differenzierbaren Funktion F : Rn ⊃ U → Rn (m = n, U = Rn).Idee des Newton-Verfahrens: Angenommen es sei xk ∈ Rn bereits berechnet. Danngilt

F(x) = F(xk) + F′(xk)(x− xk) + r(xk, x).

Ist xk ≈ x, dann ist r(xk, x) sehr klein, so dass

F(x) ≈ F(xk) + F′(xk)(x− xk)

gilt. Es ist also naheliegend, das lineare Gleichungssystem

F(xk) + F′(xk)(x− xk) = 0

zu lösen. Diese Aufgabe ist äquivalent zu:Finde xk+1 ∈ Rn als Lösung von

F′(xk)(x− xk) = −F(xk)

⇔ F′(xk)dk = −F(xk) mit

xk+1 = xk + dk.

Algorithmus 6.1 Newton-Verfahren

(S0) Wähle einen Startwert x0 ∈ Rn und setze k = 0.(S1) if F(xK) = 0 then STOP end if(S2) Bestimme dk ∈ Rn als Lösung von

F′(xk)dk = −F(xk).

(S3) Setze xk+1 = xk + dk, k = k + 1 und gehe zu (S1).

102


Bemerkung 6.19. In der Praxis ersetzen wir das Abbruchkriterium F(xk) = 0 durch

‖F(xk)‖ < ε mit ε „klein“.

Illustration des Newton-Verfahrens: (F(x) = x2)

x2 x1 x0

F F(x0)

F(x0) + F′(x0)(x− x0)

Bemerkung 6.20. Die Konvergenz des Newton-Verfahrens hängt im wesentlichen vondem Startwert ab! Ist der Startwert ungünstig gewählt, so kann das Newton-Verfahrenunter Umständen nicht konvergieren! Dazu betrachten wir das Beispiel

F : R→ R, F(x) = arctan(x).

x1 x0

− π2

π2

Unser Ziel ist es, die lokale Konvergenz des Newton-Verfahrens zu analysieren.

103

6.3 Lokale Konvergenz


Definition 6.21. Es sei {xk}∞k=1 eine gegen x ∈ Rn konvergente Folge.

(i) Die Folge {xk}∞k=1 konvergiert linear mit Rate α ∈ (0, 1) gegen x, falls ein n ∈N

existiert, so dass‖xk+1 − x‖2 ≤ α‖xk − x‖2 ∀k ≥ n.

(ii) Die Folge {xk}∞k=1 konvergiert superlinear gegen x, falls

limk→∞

∥∥xk+1 − x∥∥

2∥∥xk − x∥∥

2= 0.

Andere Notation: ‖xk+1 − x‖2 = O(‖xk − x‖2

)für k→ ∞.

(iii) Die Folge {xk}∞k=1 konvergiert quadratisch gegen x, falls es eine Konstante c > 0

gibt, so dass‖xk+1 − x‖2 ≤ c‖xk − x‖2

2 ∀k ∈N.

Andere Notation: ‖xk+1 − x‖2 = O(‖xk − x‖2

2)∀k ∈N.

Lemma 6.22. Es sei F : Rn → Rn stetig differenzierbar und x ∈ Rn eine Nullstelle von F.Ferner sei die Jacobi-Matrix von F in x

(∂jFi(x)) ∈ Rn×n

regulär. Dann gibt es ein ε > 0 und ein β > 0 mit

β ‖x− x‖2 ≤ ‖F(x)‖2 ∀x ∈ Bε(x).

Insbesondere ist x ∈ Rn die einzige Nullstelle von F in Bε(x).

Beweis. Laut Voraussetzung gilt

‖x− x‖2 = ‖(∂jFi(x))−1(∂jFi(x))(x− x)‖2 ≤ ‖(∂jFi(x))−1‖2‖(∂jFi(x))(x− x)‖2

= ‖(∂jFi(x))−1‖2‖F′(x)(x− x)‖2.

Wir setzenβ :=

12∥∥(∂jFi(x))−1

∥∥2

∈ R+.

Dann gilt2β ‖x− x‖2 ≤

∥∥F′(x)(x− x)∥∥

2 . (6.1)

Laut Definition der Differenzierbarkeit gilt

‖r(x, x)‖2‖x− x‖2

=‖F(x)− F(x)− F′(x)(x− x)‖2

‖x− x‖2→ 0 für x → x.

104


Folglich gibt es ein ε > 0, so dass

‖F(x)− F(x)− F′(x)(x− x)‖2‖x− x‖2

≤ β ∀x ∈ Bε(x). (6.2)

Mit (6.1), (6.2) und F(x) = 0 folgt

2β ‖x− x‖2 ≤∥∥F′(x)(x− x)

∥∥2 =

∥∥F(x)− (F(x)− F(x)− F′(x)(x− x))∥∥

2 ‖x− x‖2

≤ ‖F(x)‖2 +∥∥F(x)− F(x)− F′(x)(x− x)

∥∥2

≤ ‖F(x)‖2 + β ‖x− x‖2 ∀x ∈ Bε(x).


β ‖x− x‖2 ≤ ‖F(x)‖2 ∀x ∈ Bε(x).

Lemma 6.23. Es seien A, B ∈ Rn×n mit ‖I − BA‖2 < 1. Dann sind A und B regulär undes gilt

‖A−1‖2 ≤‖B‖2

1− ‖I − BA‖2.

Lemma 6.24. Es sei F : Rn → Rn stetig differenzierbar und x ∈ Rn eine Nullstelle von F.Ferner sei die Jacobi-Matrix von F in x


regulär. Dann gibt es wieder ein ε > 0, so dass die Jacobi-Matrizen

(∂jFi(x)) ∈ Rn×n ∀x ∈ Bε(x)

regulär sind, und es gilt

‖(∂jFi(x))−1‖2 ≤ 2‖(∂jFi(x))−1‖2 ∀x ∈ Bε(x).

Beweis. Da F : Rn → Rn stetig differenzierbar ist, gibt es ein ε > 0, so dass

∥∥(∂jFi(x))− (∂jFi(x))∥∥

2 ≤1

2∥∥(∂jFi(x))−1

∥∥2

∀x ∈ Bε(x).

Also ist

‖I − (∂jFi(x))−1(∂jFi(x))‖2 = ‖(∂jFi(x))−1 [(∂jFi(x))− (∂jFi(x))]‖2

≤ ‖(∂jFi(x))−1‖2‖(∂jFi(x))− (∂jFi(x))‖2

≤ 12∀x ∈ Bε(x).

105


Das vorherige Lemma mit

B := (∂jFi(x))−1 und A := (∂jFi(x))

impliziert, dass (∂jFi(x)) ∈ Rn×n für alle x ∈ Bε(x) regulär ist und es gilt

‖(∂jFi(x))−1‖2 ≤∥∥(∂jFi(x))−1

∥∥2

1−∥∥I − (∂jFi(x))−1(∂jFi(x))

∥∥2

≤ 2‖(∂jFi(x))−1‖2 ∀x ∈ Bε(x).

Satz 6.25 (Lokale superlineare Konvergenz des Newton-Verfahrens). Es sei F : Rn →Rn stetig differenzierbar und x ∈ Rn eine Nullstelle von F. Ferner sei die Jacobi-Matrix


regulär. Dann gibt es ein ε > 0, so dass für jeden Startwert x0 ∈ Bε(x) das Newton-Verfahrensuperlinear gegen x konvergiert.

Beweis. Laut Definition ist

r(x, y) = F(y)− F(x)− F′(x)(y− x) ∀x, y ∈ Rn.

Wir zeigen zunächst, dass das Restglied r auch die folgende Integraldarstellung besitzt:

r(x, y) =∫ 1

0F′(x + t(y− x))(y− x)− F′(x)(y− x)dt ∀x, y ∈ Rn.

Seien x, y ∈ Rn beliebig aber fest. Wir definieren die Hilfsfunktion

G : R→ Rn, G(t) = F(x + t(y− x)).

Nach Kettenregel istG′(t) = F′(x + t(y− x))(y− x).

Somit gilt∫ 1

0F′(x + t(y− x))(y− x)dt =

∫ 1

0G′(t)dt = G(1)− G(0) = F(y)− F(x).

Daraus folgt ∫ 1

0F′(x + t(y− x))(y− x)dt = F(y)− F(x)

und wir erhalten

r(x, y) =∫ 1

0F′(x + t(y− x))(y− x)dt− F′(x)(y− x)

=∫ 1

0F′(x + t(y− x))(y− x)− F′(x)(y− x)dt.

106


Aus der Integraldarstellung folgt

‖r(x, x)‖2 =

∥∥∥∥∫ 1

0F′(x + t(x− x))(x− x)− F′(x)(x− x)dt

∥∥∥∥2

≤∫ 1

0

∥∥F′(x + t(x− x))(x− x)− F′(x)(x− x)∥∥

2 dt

=∫ 1

0

∥∥[(∂jFi(x + t(x− x)))− (∂jFi(x))](x− x)

∥∥2 dt

≤∫ 1

0

∥∥(∂jFi(x + t(x− x)))− (∂jFi(x))∥∥

2 ‖x− x‖2 dt

=∫ 1

0

∥∥(∂jFi(x + t(x− x)))− (∂jFi(x))∥∥

2 dt ‖x− x‖2 (6.3)

Sei α ∈ (0, 1) beliebig aber fest. Da F : Rn → Rn stetig differenzierbar ist, existiert einε > 0, so dass

∫ 1

0

∥∥(∂jFi(x + t(x− x)))− (∂jFi(x))∥∥

2 dt ≤ α

2∥∥(∂jFi(x))−1

∥∥2

∀x ∈ Bε(x).

Folglich

‖r(x, x)‖2 =α

2∥∥(∂jFi(x))−1

∥∥2

‖x− x‖2 ∀x ∈ Bε(x). (6.4)

Wir können nun ε > 0 verkleinern, so dass gilt:

‖(∂jFi(x))−1‖2 ≤ 2‖(∂jFi(x))−1‖2 ∀x ∈ Bε(x). (6.5)

Sei nun xk ∈ Bε(x) und xk+1 ∈ Rn die durch das Newton-Verfahren erzeugte neueIterierte. Dann gilt:

xk+1 − x = xk+1 − xk + xk − x

= (∂jFi(xk))−1(∂jFi(xk))(xk+1 − xk) + xk − x

= (∂jFi(xk))−1 F′(xk)(xk+1 − xk)︸︷︷︸=−F(xk)

+xk − x

= −(∂jFi(xk))−1F(xk) + xk − x

= (∂jFi(xk))−1(−F(xk) + (∂jFi(xk))(xk − x))F(x)=0= (∂jFi(xk))−1(F(x)− F(xk) + F′(xk)(xk − x))

= (∂jFi(xk))−1(F(x)− F(xk)− F′(xk)(x− xk))

= (∂jFi(xk))−1r(xk, x).

107


Aus (6.4) und (6.5) folgt

‖xk+1 − x‖2 = ‖(∂jFi(xk))−1r(xk, x)‖2

≤ ‖(∂jFi(xk))−1‖2‖r(xk, x)‖2

(6.4)≤ ‖(∂jFi(xk))−1‖2

α

2‖(∂jFi(x))−1‖2‖x− x‖2

(6.5)≤ α‖xk − x‖2.

Ist x0 ∈ Bε(x), so folgt aus der obigen Abschätzung:

xk ∈ Bε(x) ∀k ∈N (denn ‖xk − x‖2 ≤ αk‖x0 − x‖2 < ε)

und‖xk+1 − x‖2 ≤ α‖xk − x‖2 ∀k ∈N.

Mit anderen Worten: Ist x0 ∈ Bε(x), so konvergiert die durch das Newton-Verfahrenerzeugte Folge {xk}∞

k=1 linear (mit Rate α) gegen x. Nun zeigen wir die superlineareKonvergenz:

‖xk+1 − x‖2 = ‖(∂jFi(xk))−1r(xk, x)‖2

(6.5)≤ 2‖(∂jFi(xk))−1‖2‖r(xk, x)‖2

(6.3)≤ 2‖(∂jFi(xk))−1‖2

∫ 1

0

∥∥∥(∂jFi(xk + t(x− xk)))− (∂jFi(xk))∥∥∥

2dt‖xk − x‖2

mit ∫ 1

0

∥∥∥(∂jFi(xk + t(x− xk)))− (∂jFi(xk))∥∥∥

2dt→ 0 für k→ ∞,

da F : Rn → Rn stetig differenzierbar ist und {xk}∞k=1 gegen x konvergiert. Insgesamt

erhalten wir: ∥∥xk+1 − x∥∥

2∥∥xk − x∥∥

2→ 0 für k→ ∞.

Satz 6.26 (Lokale quadratische Konvergenz des Newton-Verfahrens). Es sei F : Rn →Rn stetig differenzierbar und x ∈ Rn eine Nullstelle von F. Ferner sei die Jacobi-Matrix


regulär. Es existiere ein L > 0, so dass∥∥(∂jFi(x))− (∂jFi(y))

∥∥2 ≤ L ‖x− y‖2 ∀x, y ∈ Rn,

d.h. F′ : Rn → L (Rn, Rn) ist Lipschitz-stetig. Dann existiert ein ε > 0, so dass für jedenStartwert x0 ∈ Bε(x) das Newton-Verfahren quadratisch gegen x konvergiert.

108

6.4 Das vereinfachte Newton-Verfahren

Beweis. Im Beweis des vorherigen Satzes haben wir gezeigt, dass ein ε > 0 existiert, sodass für jeden Startwert x0 ∈ Bε(x) gilt:

(i): xk ∈ Bε(x) ∀k ∈N.

(ii): {xk}∞k=1 konvergiert superlinear gegen x.

Außerdem gilt

‖xk+1 − x‖2 ≤ 2∥∥∥(∂jFi(xk))−1

∥∥∥2

∫ 1

0

∥∥∥(∂jFi(xk + t(x− xk)))− (∂jFi(xk))∥∥∥

2dt‖xk − x‖2

≤ 2∥∥∥(∂jFi(xk))−1

∥∥∥2

∫ 1

0L∥∥∥xk + t(x− xk)− xk

∥∥∥2

dt‖xk − x‖2

≤ 2∥∥∥(∂jFi(xk))−1

∥∥∥2

L∫ 1

0t dt‖xk − x‖2

2

= c‖xk − x‖22

mit c := 2∥∥(∂jFi(xk))−1

∥∥2 L 1

2 = L∥∥(∂jFi(xk))−1

∥∥2. Folglich gilt

‖xk+1 − x‖2 ≤ c‖xk − x‖22 ∀k ∈N.

Bemerkung 6.27. Die Aussage bleibt auch richtig, falls die Abbildung F′ : Rn →L (Rn, Rn) nur lokal Lipschitz-stetig ist.


Es sei F : Rn → Rn stetig differenzierbar. Beim Newton-Verfahren zur Lösung vonF(x) = 0 lösen wir in jeder Iteration das lineare Gleichungssystem

F′(xk)dk = −F(xk).

Für jedes k ∈N müssen wir also:

(i) Die Jacobi-Matrix von F in xk berechnen.

(ii) Das lineare Gleichungssystem

(∂jFi(xk))dk = −F(xk)

zum Beispiel mittels LR-Zerlegung lösen.

109


Die Berechnungen von (i) und (ii) müssen in jedem Schritt des Newton-Verfahrensdurchgeführt werden. Das kann aber manchmal sehr teuer sein.Beim vereinfachten Newton-Verfahren bestimmen wir die dk ∈ Rn als Lösung von

F′(x0)dk = −F(xk) (Newton-Verfahren: F′(xk)dk = −F(xk)).

Diese Strategie hat den Vorteil, dass man im gesamten Algorithmus nur einmal dieJacobi-Matrix von F in x0 und deren LR-Zerlegung bestimmen muss. Der Rechenauf-wand beim vereinfachten Newton-Verfahren ist somit meist überaus geringer als beimNewton-Verfahren.

Algorithmus 6.2 Vereinfachtes Newton-Verfahren

(S0) Wähle einen Startwert x0 ∈ Rn und setze k = 0.(S1) if F(xK) = 0 then STOP end if(S2) Bestimme dk ∈ Rn als Lösung von

F′(x0)dk = −F(xk).

(S3) Setze xk+1 = xk + dk, k = k + 1 und gehe zu (S1).

Der Nachteil des Verfahrens besteht darin, dass dieses nur linear lokal konvergiert(also insbesondere nicht superlinear oder quadratisch).

Satz 6.28 (Lokale lineare Konvergenz des vereinfachten Newton-Verfahrens). Es seiF : Rn → Rn stetig differenzierbar und x ∈ Rn eine Nullstelle von F. Ferner sei dieJacobi-Matrix von F in x


regulär. Dann gibt es ein ε > 0, so dass für jeden Startwert x0 ∈ Bε(x) das vereinfachteNewton-Verfahren linear gegen x konvergiert.

Beweis. Da die Jacobi-Matrix (∂jFi(x)) ∈ Rn×n regulär ist, gibt es ein ε1 > 0, so dassgilt

‖(∂jFi(x))−1‖2 ≤ c ∀x ∈ Bε1(x) (6.6)

mit c := 2∥∥(∂jFi(x))−1

∥∥2 (siehe Lemma). Weiter existiert ein ε2 > 0, so dass gilt

‖r(x, x)‖2 ≤14

c ‖x− x‖2 ∀x ∈ Bε2(x), (6.7)

siehe Beweis der lokalen superlinearen Konvergenz des Newton-Verfahrens (Integ-raldarstellung von r). Außerdem gibt es ein ε3 > 0, so dass

∥∥(∂jFi(x))− (∂jFi(y))∥∥

2 ≤14c∀x, y ∈ Bε3(x), (6.8)

110


da F : Rn → Rn stetig differenzierbar ist. Setze nun ε = min{ε1, ε2, ε3} und wählex0 ∈ Bε(x). Dann gilt für jedes xk ∈ Bε(x) und die durch das vereinfachte Newton-Verfahren erzeugte neue Iterierte xk+1 ∈ Rn:

‖xk+1 − x‖2 = ‖xk+1 − xk + xk − x‖2

= ‖F′(x0)−1F′(x0)(xk+1 − xk) + xk − x‖2

= ‖−F′(x0)−1F(xk) + xk − x‖2

= ‖(∂jFi(x0))−1(−F(xk) + (∂jFi(x0))(xk − x))‖2

≤ ‖(∂jFi(x0))−1‖2‖−F(xk) + F′(x0)(xk − x)‖2

(6.6)≤ c‖−F(xk) + F′(x0)(xk − x)‖2

= c‖−F(xk) + F(x)− F′(x0)(x− xk)‖2

≤ c‖F(x)− F(xk)− F′(xk)(x− xk)‖2 + c‖F′(xk)(x− xk)− F′(x0)(x− xk)‖2

(6.7)≤ 1

4‖x− xk‖2 + c‖

[(∂jFi(xk))− (∂jFi(x0))

](x− xk)‖2

≤ 14‖x− xk‖2 + c‖(∂jFi(xk))− (∂jFi(x0))‖2‖(x− xk)‖2

(6.8)≤ 1

4‖x− xk‖2 +

14‖x− xk‖2

=12‖x− xk‖2.

Folglich gilt für jeden Startwert x0 ∈ Bε(x):{

xk ∈ Bε(x) ∀k ∈N (denn∥∥xk − x

∥∥2 ≤ (1

2)k∥∥x0 − x

∥∥2 < ε),∥∥xk+1 − x

∥∥2 ≤ 1

2

∥∥xk − x∥∥

2 ∀k ∈N.

Daraus ergibt sich die Behauptung.

111

Kapitel 7Interpolation

Die Interpolationsaufgabe besteht darin, zu gegebenen Punktepaaren (xk, yk) ∈ R2,k = 0, . . . , n eine Funktion p : R→ R zu finden, so dass

p(xk) = yk ∀k = 0, . . . , n

gilt. Mögliche Ansätze für p : R→ R sind:

(i) Polynome p(x) = a0 + a1x + . . . + anxn, also z.B.

x0 x1 x2 x3

y0

y1

y2 y3

113

7.1 Polynominterpolation

(ii) Splines (stückweise Polynome), also z.B.

x0 x1 x2 x3

y0

y1

y2 y3

(iii) Trigonometrische Polynome.


Definition 7.1. Mit Πn bezeichnen wir die Menge aller reeller Polynome vom Grad ≤n, das heißt

Πn := {p : R→ R | ∃a0, . . . , an ∈ R : p(x) = a0 + a1x + . . . + anxn ∀x ∈ R}.

Satz 7.2. Es sei n ∈ N und (xk, yk) ∈ R2, k = 0, . . . , n mit xi 6= xj ∀i 6= j. Dann gibt esgenau ein Polynom p ∈ Πn mit

p(xk) = yk ∀k = 0, . . . , n.

Beweis. Wir machen den Ansatz

p(x) = a0 + a1x + . . . + anxn.

Gesucht sind a0, . . . , an ∈ R, so dass gilt

p(xk) = yk ∀k = 0, . . . , n.

Dies ist äquivalent zu:

a0 + a1x0 + a2x20 + . . . + anxn

0 = y0

a0 + a1x1 + a2x21 + . . . + anxn

1 = y1

...

a0 + a1xn + a2x2n + . . . + anxn

n = yn.

114


Dieses Gleichungssystem entspricht

1 x0 x20 · · · xn

01 x1 x2

1 · · · xn1

......

......

1 xn x2n · · · xn

n

a0a1...

an

=

y0y1...

yn

.

Da die xi paarweise verschieden sind, ist die obige Matrix regulär. Somit gibt es genaua0, . . . , an ∈ R mit

p(xk) = yk ∀k = 0, . . . , n.

x0 x1

y0

y1

n = 1

x0 x1 x2

y0 y1

y2

n = 2

Definition 7.3 (Interpolationspolynom). Das nach dem Satz eindeutig bestimmte Poly-nom p ∈ Πn mit der Eigenschaft

p(xk) = yk ∀k = 0, . . . , n

für die vorgegebenen Punktepaare

(xk, yk) ∈ R2, k = 0, . . . , n mit xi 6= xj ∀i 6= j

heißt Interpolationspolynom.

Eine kanonische Basis für Πn ist gegeben durch

{1, x, x2, . . . , xn} (Monombasis).

Wir wollen nun weitere Basis-Ansätze für Πn untersuchen.

7.1.0.1 Lagrangesche Interpolation

Im Folgenden seien Punktepaare

(xk, yk) ∈ R2, k = 0, . . . , n

mit paarweise verschiedenen xk gegeben. Die Interpolationsaufgabe lautet:

Finde p ∈ Πn, so dass p(xk) = yk ∀k = 0, . . . , n.

115


Definition 7.4 (Lagrange-Basispolynome). Die Polynome

Lk(x) =(x− x0) · . . . · (x− xk−1)(x− xk+1) · . . . · (x− xn)

(xk − x0) · . . . · (xk − xk−1)(xk − xk+1) · . . . · (xk − xn)

heißen Lagrange-Basispolynome.

Beispiel 7.5 (n=2).

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

L0(x0) = 1, L0(x1) = 0, L0(x2) = 0,L1(x0) = 0, L1(x1) = 1, L1(x2) = 0,L2(x0) = 0, L2(x1) = 0, L2(x2) = 1.

Lemma 7.6. Es giltLk(xj) = δkj ∀k, j = 0, . . . , n.

Beweis. Klar nach Definition der Lagrange-Basispolynome.

Satz 7.7 (Lagrangesche Interpolationsformel). Gegeben seien Punktepaare (xk, yk) ∈ R2,k = 0, . . . , n mit paarweise verschiedenen xk. Dann ist das Lagrangesche Interpolationspolynom

pL(x) :=n

∑k=0

ykLk(x)

die eindeutige Lösung der Interpolationsaufgabe:

Finde p ∈ Πn, so dass p(xk) = yk ∀k = 0, . . . , n.

Beweis. Da Lj(xi) = δij ∀i, j = 0, . . . , n gilt, so ist

pL(xk) =n

∑j=0

yjLj(xk) =n

∑j=0

yjδjk = yk

für alle k = 0, . . . , n.

116


7.1.0.2 Newtonsche Interpolation

Definition 7.8 (Newton-Basispolynome). Die Polynome

Nk(x) := (x− x0) · . . . · (x− xk−1) ∀k = 1, . . . , nN0(x) ≡ 1

heißen Newton-Basispolynome.

Wir wollen nun die Interpolationsaufgabe mittels Newton-Basispolynomen lösen.Gesucht sind

a0, . . . , an ∈ R,

so dass gilt:

yk =n

∑j=0

ajNj(xk) ∀k = 0, . . . , n.

Lemma 7.9 (Aitken). Es seien Punktepaare (xk, yk) ∈ R2 mit paarweise verschiedenen xk.Ferner seien P[0], P[n] ∈ Πn−1 gegeben mit

P[0](xk) = yk ∀k = 0, . . . , n− 1,

P[n](xk) = yk ∀k = 1, . . . , n.

Dann löst das Polynom

p(x) :=P[0](x)(x− xn)− P[n](x)(x− x0)

x0 − xn

die Interpolationsaufgabe

Finde p ∈ Πn mit p(xk) = yk ∀k = 0, . . . , n.

Beweis. Laut Definition ist

p(x) =P[0](x)(x− xn)− P[n](x)(x− x0)

x0 − xn

ein Polynom vom Grad ≤ n, denn

P[0], P[n] ∈ Πn.

Für k ∈ {1, . . . , n− 1} gilt

p(xk) =P[0](xk)(x− xn)− P[n](xk)(x− x0)

x0 − xn=

yk(x− xn)− yk(x− x0)

x0 − xn

=yk(x0 − xn)

x0 − xn

= yk.

117


Für k ∈ {0, n} gilt

p(x0) =P[0](x0)(x− xn)− P[n](x0)(x− x0)

x0 − xn=

y0(x0 − xn)

x0 − xn= y0,

p(xn) =P[0](xn)(x− xn)− P[n](xn)(x− x0)

x0 − xn=

yn(x0 − xn)

x0 − xn= yn,

und somit folgt die Behauptung.

Wir setzen

P00(x) ≡ y0,P11(x) ≡ y1,

...Pnn(x) ≡ yn.

Wir definieren Pi j ∈ Πj−i, 0 ≤ i < j ≤ n rekursiv wie folgt:

Pij(x) :=Pi,j−1(x)(x− xj)− Pi+1,j(x)(x− xi)

xi − xj, 0 ≤ i < j ≤ n.

Aus dem Lemma von Aitken wird die Interpolationsaufgabe

Finde p ∈ Πn mit p(xk) = yk ∀k = 0, . . . , n

eindeutig durch P0n gelöst. Das Interpolationspolynom an der Stelle x lässt sich durchdas sogenannte Neville-Aitken-Schema berechnen:

x0 y0 = P00(x)↘

x1 y1 = P11(x) → P01(x) ↘↘

x2 y2 = P22(x) → P12(x) → P02(x)...

...... . . .

xn yn = Pnn(x) → Pn−1,n(x) · · · · · · · · · P0n(x)

Satz 7.10 (Newtonsche Interpolationsformel). Gegeben seien Punktepaare

(xk, yk) ∈ R2, k = 0, . . . , n

mit paarweise verschiedenen xk. Dann ist das Newtonsche Interpolationspolynom

pN(x) :=n

∑j=0

[x0, . . . , xn]Nj(x),

118


wobei[x0, . . . , xn] ∈ R

rekursiv definiert sind durch

[xj] := yj ∀j = 0, . . . , n

[xk, . . . , xj] :=[xk+1, . . . , xj]− [xk, . . . , xj−1]

xj − xk∀j > k ≥ 0,

die eindeutige Lösung der Interpolationsaufgabe.

Beweis. Wir haben oben bereits gezeigt, dass P0n ∈ Πn die eindeutige Lösung derInterpolationsaufgabe ist, und P0n ist definiert durch

Pij(x) =Pi,j−1(x)(x− xj)− Pi+1,j(x)(x− xi)

xi − xj

= Pi+1,j(x) +(

Pi+1,j(x)− Pi,j−1(x)) x− xj

xj − xi, 0 ≤ i < j ≤ n

undPkk(x) ≡ yk ∀k = 0, . . . , n.

Durch Koeffizientenvergleich folgt

P0n(x) = pN(x).

7.1.0.3 Interpolationsfehler

Definition 7.11. Es seien a, b ∈ R mit a < b und n ∈N. Wir definieren

C[a, b] := { f : [a, b]→ R stetig}Cn[a, b] :=

{f : (a, b)→ R n-mal differenzierbar | f (k) ∈ C[a, b] ∀k = 0, . . . , n

},

wobei f (k) die k-te Ableitung von f bezeichnet, und f (0) = f .

Lemma 7.12 (Rolle). Es seien a, b ∈ R mit a < b und f ∈ C1[a, b]. Ferner seien x1, x2 ∈[a, b] mit x1 6= x2 und f (x1) = f (x2). Dann gibt es ein ξ ∈ (x1, x2) mit f ′(ξ) = 0.

Satz 7.13 (Interpolationsfehler). Seien a, b ∈ R mit a < b und

a = x0 < x1 < . . . < xn = b.

Ferner seien f ∈ Cn+1[a, b], n ∈N und p ∈ Πn mit

p(xk) = f (xk) ∀k = 0, . . . , n.

119


Dann gibt es zu jedem x ∈ [a, b] ein ξ ∈ (a, b) mit

f (x)− p(x) =w(x)

(n + 1)!f (n+1)(ξ)

mit w(x) = (x− x0) · . . . · (x− xn).

Beweis. Die Aussage ist bereits richtig für x ∈ {x0, . . . , xn}, denn

f (xk) = p(xk) ∀k = 0, . . . , n

undw(xk) = 0 ∀k = 0, . . . , n.

Sei nun x ∈ [a, b] \ {x0, . . . , xn}. Dann gilt

w(x) Def.= (x− x0) · . . . · (x− xn) 6= 0

und wir setzen

α :=f (x)− p(x)

w(x)∈ R.

Dann istf (x)− p(x)− αw(x) = 0. (7.1)

Nun sei F : [a, b]→ R, F(x) := f (x)− p(x)− αw(x) mit w(x) = (x− x0) · . . . · (x−xn). Dann ist F ∈ Cn+1[a, b], da f ∈ Cn+1[a, b], p ∈ Πn und w ∈ Πn+1. Außerdem hatF gerade n + 2 Nullstellen

x, x0, x1, . . . , xn.

Nach dem Lemma von Rolle hat die Ableitung F′ also n + 1 Nullstellen, F(2) n Null-stelle, . . . , F(n+1) hat eine Nullstelle ξ ∈ (a, b). Also

0 = F(n+1)(ξ) = f (n+1)(ξ)− p(n+1)(ξ)− αw(n+1)(ξ) = f (n+1)(ξ)− 0− α(n + 1)!.

Daraus ergibt sich

α =1

(n + 1)!f (n+1)(ξ).

Einsetzen in (7.1) liefert

f (x)− p(x) =w(x)

(n + 1)!f (n+1)(x).

120

7.2 Hermite-Interpolation

7.2 Hermite-Interpolation

Im Folgenden sei f ∈ C1[a, b] mit −∞ < a < b < ∞. Ferner seien x0, . . . , xn ∈ [a, b]paarweise verschieden. Die Hermitesche Interpolationsaufgabe lautet:Finde p ∈ Π2n+1 mit

p(xk) = f (xk) ∀k = 0, . . . , n,p′(xk) = f ′(xk) ∀k = 0, . . . , n.

Wir suchen Basis-Funktionen

Φk, Ψk ∈ Π2n+1 ∀k = 0, . . . , n


Φk(xj) = δkj ∀k, j = 0, . . . , n,

Φ′k(xj) = 0 ∀k, j = 0, . . . , n,

Ψk(xj) = 0 ∀k, j = 0, . . . , n,

Ψk(xj) = δkj ∀k, j = 0, . . . , n.

Wir wissen bereits, dass die Lagrange-Basisfunktion Lk ∈ Πn mit

Lk(x) =(x− x0) · . . . · (x− xk−1)(x− xk+1) · . . . · (x− xn)

(xk − x0) · . . . · (xk − xk−1)(xk − xk+1) · . . . · (xk − xn)

die EigenschaftLk(xj) = δkj ∀k, j = 0, . . . , n

erfüllt. Nun definieren wir

Φk(x) := (1− 2L′k(xk)(x− xk))L2k(x)

Ψk(x) := (x− xk)L2k(x)

für alle k = 0, . . . , n. Nun gilt

(i) Φk(xj) = (1− 2L′k(xk)(xj − xk)) L2k(xj)︸︷︷︸=δkj

= δkj ∀k, j = 0, . . . , n.

(ii) Φ′k(xj) = −2L′k(xk)L2k(xj) + (1− 2L′k(xk)(xj − xk))2Lk(xj)L′k(xj)

= −2L′k(xk)δkj + (1− 2L′k(xk)(xj − xk))2δkjL′k(xj)

= −2L′k(xk)δkj + 2δkjL′k(xj)

= 0 ∀k, j = 0, . . . , n.

(iii) Ψk(xj) = (xj − xk)L2k(xj) = (xj − xk)δkj = 0 ∀k, j = 0, . . . , n.

121

7.3 Splines

(iv) Ψ′k(xj) = L2k(xj)︸︷︷︸=δkj

+(xj − xk)2= Lk(xj)︸︷︷︸δkj

L′k(xj)

= δkj + 2(xj − xk)δkjL′k(xj)

= δkj ∀k, j = 0, . . . , n.

Die Lösung der Hermiteschen Interpolationsaufgabe ist dann gegeben durch

pH(x) :=n

∑j=0

f (xj)Φj(x) +n

∑j=0

f ′(xj)Ψj(x),

denn laut Konstruktion sind

Φj, Ψj ∈ Π2n+1 ∀j = 0, . . . , n

und

pH(xk) =n

∑j=0

f (xj)Φj(xk)︸︷︷︸=δkj

+n

∑j=0

f ′(xj)Ψj(xk)︸︷︷︸=0

= f (xk),

p′H(xk) =n

∑j=0

f (xj)Φ′j(xk)︸︷︷︸=0

+n

∑j=0

f ′(xj)Ψ′j(xk)︸︷︷︸=δkj

= f ′(xk).

7.3 Splines

Der Nachteil der Polynominterpolation besteht darin, dass man bei der Verfeinerungder Zerlegung keine gleichmäßige Konvergenz erwarten kann. Deshalb untersuchenwir nun ein anderes Verfahren mit einer besseren Konvergenzeigenschaft. Im Folgen-den betrachten wir die Zerlegung

∆ = {a = x0 < x1 < . . . < xn = b}.

Definition 7.14. Seien p, q ∈ N ∪ {0} mit der Eigenschaft 0 ≤ q < p. Wir definierenden Raum

S(∆, p, q) :={

s ∈ Cq[a, b] | s|[xk−1,xk]∈ Πp ∀k = 1, . . . , n

}.

Eine Funktion s ∈ S(∆, p, q) wird als Spline vom Grad p der Diffenzierbarkeitsklasse qzur Zerlegung ∆ bezeichnet.

Die Spline-Interpolationsaufgabe lautet wie folgt:Zu gegebener Funktion f : [a, b]→ R suchen wir ein s ∈ S(∆, p, q) mit

s(xk) = f (xk) ∀k = 0, . . . , n.

122

7.3 Splines

Nachdem wir diese Aufgabe gelöst haben, wollen wir den Fehler

‖ f − s‖C[a,b]

untersuchen.

Definition 7.15 (Stückweise lineare Splines). Mit

S1 := S(∆, 1, 0) ={

s ∈ C[a, b] | s|[xk−1,xk]∈ Π1 ∀k = 1, . . . , n

}

bezeichnen wir den Raum aller stückweise linearer Splines.

Beispiel 7.16.

a x1 x2 x3 x4 b

s ∈ S1

Bemerkung 7.17. Ist s ∈ S1, so gilt für jedes k ∈ {1, . . . , n}

s|[xk−1,xk]∈ Π1,

also

∃ak, bk ∈ R : s(x) = ak + bkx ∀x ∈ [xk−1, xk].

Definition 7.18 (Hutsche Basisfunktion). Wir definieren die Hutsche Basisfunktion Φj ∈S1 für j = 0, . . . , n wie folgt:

Φj(x) :=

x−xj−1xj−xj−1

falls x ∈ [xj−1, xj] und j > 0xj+1−xxj+1−xj

falls x ∈ [xj, xj+1] und j < n

0 sonst

.

123

7.3 Splines

x0 x1 x2

1

Φ0

j = 0

x0 x1 x2

1Φ1

j = 1

xn−1 xn

1Φn

j = n

Offenbar gilt für die Hutsche Basisfunktion

Φj(xi) = δji ∀i, j = 0, . . . , n.

Definition 7.19. Die Abbildung

I1 : C[a, b]→ S1, (I1 f )(x) :=n

∑j=0

f (xj)Φj(x)

heißt S1-Interpolationsoperator.

Sei f ∈ C[a, b]. Dann gilt per Definition

I1 f ∈ S1

und

(I1 f )(xk) =n

∑j=0

f (xj)Φj(xk)︸︷︷︸=δjk

= f (xk) ∀k = 0, . . . , n.

Somit löst I1 f ∈ S1 die Spline-Interpolationsaufgabe auf S1.

Satz 7.20 (Fehlerabschätzung für den S1-Interpolationsoperator). Es sei f ∈ C2[a, b]und hmax := max1≤k≤n |xk−1 − xk|. Dann gilt

‖ f − I1 f ‖C[a,b] ≤18

h2max

∥∥∥ f (2)∥∥∥

C[a,b],

wobei‖·‖C[a,b] : C[a, b]→ R, ‖v‖C[a,b] = max

t∈[a,b]|v(t)|.

Beweis. Nach Definition gilt

(I1 f )(xk) = f (xk) ∀k = 0, . . . , n

124

7.3 Splines

und für jedes k ∈ {1, . . . , n} ist I1 f |[xk−1,xk]∈ Π1. Sei nun k ∈ {1, . . . , n} beliebig aber

fest. Das Polynom I1 f |[xk−1,xk]∈ Π1 löst die Aufgabe:

Finde p ∈ Π1 mit

p(xk−1) = f (xk−1)

p(xk) = f (xk).

Nach unserem Satz über den Interpolationsfehler (Kapitel 7.1.3) gilt:Zu jedem x ∈ [xk−1, xk] gibt es ein ξ ∈ (xk−1, xk) mit

f (x)− (I1 f )(x) =w(x)

2!f (2)(ξ)

mit w(x) = (x− xk−1)(x− xk). Es folgt

| f (x)− (I1 f )(x)| =∣∣∣∣(x− xk−1)(x− xk)

2

∣∣∣∣ | f (2)(ξ)|x∈[xk−1,xk]

≤ |xk − xk−1|28

| f (2)(ξ)|

≤ 18

h2max| f (2)(ξ)|

≤ 18

h2max max

x∈[a,b]| f (2)(x)|

=18

h2max

∥∥∥ f (2)(x)∥∥∥

C[a,b].

Daraus folgt

maxx∈[xk−1,xk]

| f (x)− (I1 f )(x)| ≤ 18

h2max

∥∥∥ f (2)(x)∥∥∥

C[a,b].

Da k ∈ {1, . . . , n} beliebig aber fest gewählt wurde, folgt daraus

maxx∈[a,b]

| f (x)− (I1 f )(x)| ≤ 18

h2max

∥∥∥ f (2)(x)∥∥∥

C[a,b].

Mit anderen Worten:

‖ f − I1 f ‖C[a,b] ≤18

h2max

∥∥∥ f (2)(x)∥∥∥

C[a,b].

Definition 7.21 (Kubische Splines). Mit

S3 := S(∆, 3, 1) ={

s ∈ C1[a, b] | s|[xk−1,xk]∈ Π3 ∀k = 1, . . . , n

}

bezeichnen wir den Raum der kubischen Hermite-Splines.

125

7.3 Splines

Beispiel 7.22.

x0 x1 x2

f

s1 ∈ S1

x0 x1 x2

f

s3 ∈ S3

Definition 7.23 (Basisfunktion für S3). Wir definieren ηj, Ψj ∈ S3, j = 0, . . . , n wie folgt

ηj(x) :=

(x−xj−1)2(3xj−xj−1−2x)

(xj−xj−1)3 falls x ∈ [xj−1, xj] und j > 0(xj+1−x)2(xj+1−3xj+2x)

(xj+1−xj)3 falls x ∈ [xj, xj+1] und j < n

0 sonst

.

Ψj(x) :=

(x−xj−1)2(x−xj)

(xj−xj−1)2 falls x ∈ [xj−1, xj] und j > 0(xj+1−x)2(x−xj)

(xj+1−xj)2 falls x ∈ [xj, xj+1] und j < n

0 sonst

.

Die Basisfunktionen ηj und Ψj erfüllen

ηj(xi) = δij, η′j(xi) = 0 ∀i, j = 0, . . . , n,

Ψj(xi) = 0, Ψ′j(xi) = δij ∀i, j = 0, . . . , n.

Definition 7.24 (S3-Interpolationsoperator). Die Abbildung

I3 : C1[a, b]→ S3, (I3 f )(x) :=n

∑j=0

f (xj)ηj(x) + f ′(xj)Ψj(x)

heißt S3-Interpolationsoperator.

Ist f ∈ C1[a, b], so erfüllt I3 f

(I3 f )(xk) =n

∑j=0

f (xj) ηj(xk)︸︷︷︸=δjk

+ f ′(xj)Ψj(xk)︸︷︷︸=0

= f (xk),

(I3 f )′(xk) =n

∑j=0

f (xj) η′j(xk)︸︷︷︸=0

+ f ′(xj)Ψ′j(xk)︸︷︷︸=δjk

= f ′(xk).

126

7.3 Splines

Insbesondere löst I3 f die Spline-Interpolationsaufgabe auf S3. Nun wollen wir denFehler

‖ f − I3 f ‖C[a,b]

untersuchen.

Lemma 7.25 (Hermite-Interpolationsfehler). Es sei n ∈ N und f ∈ C2n+2[a, b]. Fernersei p ∈ Π2n+1 eine Lösung der Hermite-Interpolationsaufgabe, d.h.

p(xk) = f (xk) ∀k = 0, . . . , n,p′(xk) = f ′(xk) ∀k = 0, . . . , n.

Dann gibt es zu jedem x ∈ [a, b] ein ξ ∈ (a, b), so dass

f (x)− p(x) =w2(x)

(2n + 2)!f (2n+2)(ξ)

mit w(x) = (x− x0) · . . . · (x− xn).

Beweis. Lemma von Rolle, siehe Kapitel 7.1.3.

Satz 7.26 (Fehlerabschätzung für den S3-Interpolationsoperator). Es sei f ∈ C4[a, b]und hmax := max1≤k≤n |xk−1 − xk|. Dann gilt

‖ f − I3 f ‖C[a,b] ≤1

384h4

max

∥∥∥ f (4)∥∥∥

C[a,b].

Beweis. Nach Definition ist

(I3 f )(xk) = f (xk) ∀k = 0, . . . , n,(I3 f )′(xk) = f ′(xk) ∀k = 0, . . . , n

(7.2)

undI3 f |[xk−1,xk]

∈ Π3 ∀k = 0, . . . , n. (7.3)

Sei nun k ∈ {1, . . . , n} beliebig aber fest. Laut (7.2) und (7.3) löst I3 f |[xk−1,xk]∈ Π2n+1

(mit n = 1) die Hermite-Interpolationsaufgabe:Finde p ∈ Π3, so dass gilt

p(xk−1) = f (xk−1), p(xk) = f (xk)

p′(xk−1) = f ′(xk−1), p′(xk) = f ′(xk).

Aus dem obigen Lemma (mit n = 1, a = xk−1, b = xk) folgt:Zu jedem x ∈ [xk−1, xk] gibt es ein ξ ∈ (xk−1, xk), so dass

f (x)− p(x) =w2(x)

4!f (4)(ξ).

127

7.3 Splines

Folglich gilt für jedes x ∈ [xk−1, xk]

| f (x)− (I3 f )(x)| = 14!|x− xk−1|2|x− xk|2| f (4)(ξ)|

≤ 14!|x− xk−1|2|x− xk|2

∥∥∥ f (4)(x)∥∥∥

C[a,b]

x∈[xk−1,xk]

≤ 14!

14|xk − xk−1|2

14|xk − xk|2

∥∥∥ f (4)(x)∥∥∥

C[a,b]

=1

384|xk − xk−1|4

∥∥∥ f (4)(x)∥∥∥

C[a,b]

≤ 1384

h4max

∥∥∥ f (4)(x)∥∥∥

C[a,b].

Daraus folgt

maxx∈[xk−1,xk]

| f (x)− (I3 f )(x)| ≤ 1384

h4max

∥∥∥ f (4)(x)∥∥∥

C[a,b].

Da k ∈ {1, . . . , n} beliebig gewählt wurde, folgt

‖ f − I3 f ‖C[a,b] ≤1

384h4

max

∥∥∥ f (4)(x)∥∥∥

C[a,b].

Fazit 7.27. Die Konvergenzrate beim S3-Interpolationsoperator ist besser (2 Ordnungermehr) als beim S1- Interpolationsoperator. Der Nachteil ist, dass f ∈ C4[a, b] sein muss(bei I1 nur f ∈ C2[a, b]).

128

Kapitel 8Numerische Integration

Wir untersuchen in diesem Kapitel numerische Verfahren zur Approximation von

∫ b

af (x)dx

für eine gegebene Funktion f ∈ C[a, b].

8.1 Newton-Cotes-Formel

Idee:

(i) Wir betrachten eine äquidistante Zerlegung von [a, b]

xk := a + kh, k = 0, . . . , n, n ∈N

mit

h :=b− a

n.

(ii) Wir lösen das Interpolationspolynom für die Punktepaare (xk, f (xk)), d.h. wirsuchen ein Polynom p ∈ Πn mit

p(xk) = f (xk) ∀k = 0, . . . , n.

(iii) Wir betrachten die Approximation

∫ b

af (x)dx ≈

∫ b

ap(x)dx.

129


Illustration:

a b

s ∈ S1

b∫a

p1(x)dx

Wir wissen bereits, dass die Interpolationsaufgabe in (ii) genau eine Lösung besitzt:

p(x) =n

∑k=0

f (xk)Lk(x),

wobei Lk ∈ Πn die Lagrange-Basispolynome bezeichnen. Folglich

∫ b

ap(x)dx =

n

∑k=0

f (xk)∫ b

aLk(x)dx

=n

∑k=0

f (xk)∫ b

a

n

∏j=0j 6=k

x− xj

xk − xjdx

=n

∑k=0

f (xk)∫ b

a

n

∏j=0j 6=k

x− (a + jh)(a + kh)− (a + jh)

dx

=n

∑k=0

f (xk)∫ b

a

n

∏j=0j 6=k

x−ah − jk− j

dx.

Mit der Substitution

z =x− a

h(also dz =

1h

dx)

folgt

∫ b

ap(x)dx =

n

∑k=0

f (xk)h∫ n

0

n

∏j=0j 6=k

z− jk− j

dz

= (b− a)n

∑k=0

f (xk)1n

∫ n

0

n

∏j=0j 6=k

z− jk− j

dz.

130


Definition 8.1 (Newton-Cotes-Quadraturformel). Es sei f ∈ C[a, b] und n ∈ N. DieNäherungsformel

Qn( f ) := (b− a)n

∑k=0

f (xk)σk

mit den Gewichten

σk :=1n

∫ n

0

n

∏j=0j 6=k

z− jk− j

dz, k = 0, . . . , n

heißt Newton-Cotes-Quadraturformel.

Beispiel 8.2.

(i) n = 1: Trapezregel

σ0 =11

∫ 1

0

z− 10− 1

dz = −12(z− 1)2

∣∣∣∣1

0=

12

,

σ1 =11

∫ 1

0

z− 01− 0

dz = −12(z)2

∣∣∣∣1

0=

12

.

Daraus folgt

∫ b

af (x)dx ≈ Q1( f ) = (b− a)

1

∑k=0

f (xk)σk =(b− a)

2( f (a) + f (b)) ,

denn x0 = a, x1 = b.

(ii) n = 2: Simpson-Regel

σ0 =12

∫ 2

0

z− 10− 1

z− 20− 2

dz =14

∫ 2

0z2 − 3z + 3 dz =

16

,

σ1 =12

∫ 2

0

z− 01− 0

z− 21− 2

dz =46

,

σ2 =16

.

(iii) n = 3: Newtonsche 38 -Regel

Q3( f ) =b− a

8

(f (a) + 3 f (

2a + b3

) + 3 f (a + 2b

3) + f (b)

).

Beachte, dass für die Newton-Cotes-Formel für n ≥ 8 negative Gewichte σk auftreten.Deshalb ist die Formel nur für kleines n ∈ N gut geeignet. Um die Genauigkeit zuerhöhen, verwendet man die Newton-Cotes-Formel für kleines n (n ∈ {1, 2, 3}) aufTeilintervallen von [a, b] und summiert auf. Dies führt auf die summierte Newton-Cotes-Formel.

131


n σk Name

112

12

Trapezregel

216

46

16

Simpson-Regel

318

38

38

18

Newtonsche38

-Regel

Beispiel 8.3 (summierte Trapezregel).

(i) Wir betrachten eine äquidistante Zerlegung von [a, b]

xk = a + kh, k = 0, . . . , m, k =b− a

m.

(ii) Das Integral von f auf jedem Teilintervall [xk−1, xk] approximieren wir durch dieTrapezregel ∫ xk

xk−1

f (x)dx ≈ xk − xk−1

2( f (xk−1) + f (xk)) .

(iii) Wir summieren die obigen Formeln auf∫ b

af (x)dx =

∫ x1

x0=af (x)dx + . . . +

∫ xm=b

xm−1

f (x)dx

≈ h2( f (a) + f (x1)) + . . . +

h2( f (xm−1) + f (b))

= h

(12

f (a) +m−1

∑j=1

f (xj) +12

f (b)

)

=: Tm(F).

Illustration:

a x1 x2 b

132

8.2 Fehler von allgemeinen Quadraturformeln


Im Folgenden betrachten wir das Integral

I( f ) =∫ 1

0f (x)dx

und eine allgemeine Quadraturformel der Gestalt

Q( f ) :=n

∑j=0

ωj f (xj) (8.1)

mit xj ∈ [0, 1] und ωj ∈ R für alle j = 0, . . . , n.

Definition 8.4. Eine Quadraturformel der Gestalt (8.1) hat die Fehlerordnung m ∈ N

(den Genauigkeitsgrad m), falls gilt:

Q(p) =∫ 1

0p(x)dx ∀p ∈ Πm−1

Q( p) 6=∫ 1

0p(x)dx

für ein Polynom p ∈ Πm.

Beispiel 8.5. Die Trapezregel ist eine Quadraturformel der Ordnung 2, denn

T( f ) =12( f (0) + f (1)) .

Also

T(p) =∫ 1

0p(x)dx ∀p ∈ Π1 und T(x2) =

12(02 + 12) =

126=∫ 1

0x2 dx =

13

.

Satz 8.6. Es sei Q eine Quadraturformel der Ordnung m ∈N. Dann gibt es eine Funktionk : [0, 1]→ R, so dass gilt

∫ 1

0f (x)dx−Q( f ) =

∫ 1

0k(x) f (m)(x)dx ∀ f ∈ Cm[0, 1].

Die Funktion k heißt Peano-Kern der Quadraturformel.

Beweis. Es sei f ∈ Cm[0, 1]. Der Satz von Taylor liefert

f (x) =m−1

∑k=0

f (k)(0)k!

xk

︸︷︷︸=:p(x)

+1

(m− 1)!

∫ x

0f (m)(t)(x− t)m−1 dt

︸︷︷︸=:r(x)

.

133


Das heißt,f (x) = p(x) + r(x).

Da p ∈ Πm−1 ist und Q die Fehlerordnung m hat, gilt

Q(p) =∫ 1

0p(x)dx.

Folglich gilt

∫ 1

0f (x)dx−Q( f ) =

∫ 1

0p(x) + r(x)dx−Q(p + r) =

∫ 1

0r(x)dx−Q(r) =: I.

Für I gilt

I =1

(m− 1)!

(∫ 1

0

∫ x

0f (m)(t)(x− t)m−1 dt dx−

n

∑j=0

ωj

∫ xj

0f (m)(t)(xj − t)m−1 dt

)

=1

(m− 1)!

(∫ 1

0

∫ 1

tf (m)(t)(x− t)m−1 dt dx−

n

∑j=0

ωj

∫ 1

0f (m)(t)(xj − t)m−1

+ dt

),

wobei (c)+ = max{c, 0}. Insgesamt gilt

I =1

(m− 1)!

(∫ 1

0f (m)(t)

1m(1− t)m−1 dt−

n

∑j=0

ωj

∫ 1

0f (m)(t)(xj − t)m−1

+ dt

)

=1

(m− 1)!

∫ 1

0

(1m(1− t)m −

n

∑j=0

ωj(xj − t)m−1+

)

︸︷︷︸=:k(t)

f (m)(t)dt.

Es folgt also die Behauptung

∫ 1

0f (x)dx−Q(t) =

∫ 1

0k(x) f (m)(x)dx ∀ f ∈ Cm[0, 1]

mit

k(x) =1

(m− 1)!

(1m(1− t)m −

n

∑j=0

ωj(xj − t)m−1+

).

Beispiel 8.7. Die Trapezregel

T( f ) =12( f (0) + f (1))

134


ist eine Quadraturformel der Ordnung 2. Der Satz ist also anwendbar, und der Peano-Kern k ist gegeben durch

kT(x) =11!

(12(1− x)2 − 1

2(−x)+ −

12(1− x)+

)

=12(1− x)2 − 1

2(1− x) ∀x ∈ [0, 1]

(= −1

2x(1− x)

).

Es gilt für jede Funktion f ∈ C2[0, 1]:

∫ 1

0f (x)dx− T( f ) = −1

2

∫ 1

0x(1− x) f (2)(x)dx.

135

Numerische Mathematik 1 - uni-due.de · Kapitel 1 Einleitung Die Aufgabe der Numerik ist die...

Documents

Transcript of Numerische Mathematik 1 - uni-due.de · Kapitel 1 Einleitung Die Aufgabe der Numerik ist die...