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
2.1 Zahlendarstellung und Rundungsfehler
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
2.1 Zahlendarstellung und Rundungsfehler
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.
2.2 Vektor- und Matrixnormen
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
2.2 Vektor- und Matrixnormen
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
2.2 Vektor- und Matrixnormen
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
2.2 Vektor- und Matrixnormen
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
2.2 Vektor- und Matrixnormen
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
2.2 Vektor- und Matrixnormen
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
2.2 Vektor- und Matrixnormen
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
2.2 Vektor- und Matrixnormen
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).
2.3 Kondition und Fehlerabschätzung
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
2.3 Kondition und Fehlerabschätzung
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
2.3 Kondition und Fehlerabschätzung
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
2.3 Kondition und Fehlerabschätzung
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
2.3 Kondition und Fehlerabschätzung
(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
2.3 Kondition und Fehlerabschätzung
(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
2.3 Kondition und Fehlerabschätzung
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
2.3 Kondition und Fehlerabschätzung
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
2.3 Kondition und Fehlerabschätzung
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
2.3 Kondition und Fehlerabschätzung
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
3.1 Lineare Gleichungssysteme einfacher Strukturen
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
3.1 Lineare Gleichungssysteme einfacher Strukturen
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
3.2 LR-Zerlegung ohne Pivotisierung
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
3.2 LR-Zerlegung ohne Pivotisierung
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
3.2 LR-Zerlegung ohne Pivotisierung
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
3.2 LR-Zerlegung ohne Pivotisierung
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
3.2 LR-Zerlegung ohne Pivotisierung
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
3.2 LR-Zerlegung ohne Pivotisierung
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
3.2 LR-Zerlegung ohne Pivotisierung
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
3.2 LR-Zerlegung ohne Pivotisierung
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
3.2 LR-Zerlegung ohne Pivotisierung
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
3.2 LR-Zerlegung ohne Pivotisierung
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
3.2 LR-Zerlegung ohne Pivotisierung
• 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.
3.3 LR-Zerlegung mit Pivotisierung
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
3.3 LR-Zerlegung mit Pivotisierung
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
3.3 LR-Zerlegung mit Pivotisierung
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
3.3 LR-Zerlegung mit Pivotisierung
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
3.3 LR-Zerlegung mit Pivotisierung
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
3.3 LR-Zerlegung mit Pivotisierung
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
3.3 LR-Zerlegung mit Pivotisierung
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.
Insgesamt ergibt sich also
(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
3.3 LR-Zerlegung mit Pivotisierung
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
3.3 LR-Zerlegung mit Pivotisierung
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
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
3.4 Cholesky-Zerlegung
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
3.4 Cholesky-Zerlegung
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
3.4 Cholesky-Zerlegung
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
3.4 Cholesky-Zerlegung
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
Kapitel 4 Lineares Ausgleichsproblem
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
Kapitel 4 Lineares Ausgleichsproblem
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
Kapitel 4 Lineares Ausgleichsproblem
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
mit den Eigenschaften
– 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
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
4.2 Gram-Schmidt-Orthogonalisierung
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
4.2 Gram-Schmidt-Orthogonalisierung
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
4.2 Gram-Schmidt-Orthogonalisierung
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.
4.3 Householder-Spiegelungen
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
4.3 Householder-Spiegelungen
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
4.3 Householder-Spiegelungen
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
4.3 Householder-Spiegelungen
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
.
Insgesamt ergibt sich also
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
4.3 Householder-Spiegelungen
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
4.3 Householder-Spiegelungen
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
4.3 Householder-Spiegelungen
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
4.3 Householder-Spiegelungen
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
5.1 Herleitung des CG-Verfahrens
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
5.1 Herleitung des CG-Verfahrens
(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
5.1 Herleitung des CG-Verfahrens
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
5.1 Herleitung des CG-Verfahrens
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
5.1 Herleitung des CG-Verfahrens
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
5.1 Herleitung des CG-Verfahrens
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.
5.2 Charakterisierung des CG-Verfahrens
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
5.2 Charakterisierung des CG-Verfahrens
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.
5.3 Konvergenzanalyse des CG-Verfahrens
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
5.3 Konvergenzanalyse des CG-Verfahrens
(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))
und insgesamt ergibt sich
Tn+1(x) = 2xTn(x)− Tn−1(x).
85
5.3 Konvergenzanalyse des CG-Verfahrens
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
5.3 Konvergenzanalyse des CG-Verfahrens
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
5.3 Konvergenzanalyse des CG-Verfahrens
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
5.3 Konvergenzanalyse des CG-Verfahrens
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
5.3 Konvergenzanalyse des CG-Verfahrens
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
5.3 Konvergenzanalyse des CG-Verfahrens
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
‖ek‖Rn A ≤ minp∈∏k
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
‖ek‖Rn A ≤ minp∈∏k
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.
5.4 Das präkonditionierte CG-Verfahren
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
5.4 Das präkonditionierte CG-Verfahren
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
5.4 Das präkonditionierte CG-Verfahren
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
6.1 Differentialrechnung
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
6.1 Differentialrechnung
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
6.1 Differentialrechnung
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.
Beispiel 6.11. Betrachte
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
6.1 Differentialrechnung
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
und insgesamt ergibt sich
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
6.1 Differentialrechnung
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
6.1 Differentialrechnung
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.
Beispiel 6.16. Betrachte
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.
6.2 Newton-Verfahren
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
6.2 Newton-Verfahren
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
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
6.3 Lokale Konvergenz
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).
Insgesamt ergibt sich also
β ‖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
(∂jFi(x)) ∈ Rn×n
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
6.3 Lokale Konvergenz
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
(∂jFi(x)) ∈ Rn×n
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
6.3 Lokale Konvergenz
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
6.3 Lokale Konvergenz
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
(∂jFi(x)) ∈ Rn×n
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.
6.4 Das vereinfachte Newton-Verfahren
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
6.4 Das vereinfachte Newton-Verfahren
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
(∂jFi(x)) ∈ Rn×n
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
6.4 Das vereinfachte Newton-Verfahren
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.
7.1 Polynominterpolation
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
7.1 Polynominterpolation
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
7.1 Polynominterpolation
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 Polynominterpolation
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
7.1 Polynominterpolation
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
7.1 Polynominterpolation
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
7.1 Polynominterpolation
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
mit den Eigenschaften
Φ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
8.1 Newton-Cotes-Formel
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
8.1 Newton-Cotes-Formel
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
8.1 Newton-Cotes-Formel
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
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
8.2 Fehler von allgemeinen Quadraturformeln
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
8.2 Fehler von allgemeinen Quadraturformeln
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