Plaintext
Einführung in die Numerik
Dr. Ole Klein
Interdisziplinäres Zentrum für Wissenschaftliches Rechnen (IWR)
Universität Heidelberg
email: ole.klein@iwr.uni-heidelberg.de
Basierend auf den Vorlesungsunterlagen
von Prof. Dr. Peter Bastian
Wintersemester 2021/22
Lizenz und Copyright
Dieses Werk ist unter der Creative Commons-Lizenz 4.0 (CC BY-SA)
veröffentlicht.
Bis auf einzelne Ausnahmen, die jeweils entsprechend gekennzeichnet sind, gilt für den Inhalt
dieses Skripts:
Text, Illustrationen © 2021 Ole Klein, Peter Bastian
Falls Sie in Ihren wissenschaftlichen Arbeiten Bezug auf den Inhalt dieses Skripts nehmen
möchten, referenzieren Sie es bitte via “Klein, O., and Bastian, P., 2021: Einführung in die
Numerik, Lecture Notes, Heidelberg University” unter Angabe des Links, unter dem dieses
Dokument online zu finden ist:
https://conan.iwr.uni-heidelberg.de/people/oklein/
Inhaltsverzeichnis
Inhaltsverzeichnis
1 Motivation 3
1.1 Die Wissenschaftliche Methode . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.2 Ein Beispiel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.3 Lösungsansätze . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.4 Numerische Resultate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.5 Weitere Anwendungen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.6 Ausblick . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2 Fließkommazahlen 15
2.1 Zahlendarstellung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.2 Runden und Rundungsfehler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.3 Fließkommaarithmetik . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.4 Fehleranalyse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
2.5 Die quadratische Gleichung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
2.6 Auslöschung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.7 Die Exponentialfunktion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
2.8 Rekursionsformel für Integrale . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
3 Motivation linearer Gleichungssysteme 37
3.1 Strömung in Rohrleitungsnetzen . . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.2 Die Kirchhoffschen Gesetze . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
3.3 Das Knotenpotentialverfahren . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.4 Partielle Differentialgleichungen . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
3.5 Radiosity-Methode in der Computergraphik . . . . . . . . . . . . . . . . . . . . 45
4 Kondition linearer Gleichungssysteme 50
4.1 Lösbarkeit linearer Gleichungssysteme . . . . . . . . . . . . . . . . . . . . . . . 50
4.2 Vektornormen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
4.3 Matrixnormen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
4.4 Eigenwerte und Eigenvektoren . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
4.5 Die Spektralnorm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
4.6 Positiv definite Matrizen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
4.7 Störungstheorie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
5 Eliminationsverfahren zur Lösung linearer Gleichungssysteme 61
5.1 Dreieckssysteme (gestaffelte Gleichungssysteme) . . . . . . . . . . . . . . . . . . 61
5.2 Gauß-Elimination . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
5.3 LR-Zerlegung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
5.4 Rundungsfehleranalyse der LR-Zerlegung . . . . . . . . . . . . . . . . . . . . . . 71
5.5 Pivotisierung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
5.6 Spezielle Systeme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
5.7 Nichtreguläre Systeme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
1
Inhaltsverzeichnis
6 Interpolation und Approximation 92
6.1 Einführung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
6.2 Polynominterpolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
6.3 Anwendungen der Polynominterpolation . . . . . . . . . . . . . . . . . . . . . . 100
6.4 Bernstein-Polynome und Kurvendarstellung . . . . . . . . . . . . . . . . . . . . 104
6.5 Splines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
6.6 Trigonometrische Interpolation . . . . . . . . . . . . . . . . . . . . . . . . . . . 114
6.7 Gauß-Approximation von Funktionen . . . . . . . . . . . . . . . . . . . . . . . . 120
6.8 Tschebyscheff-Approximation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128
7 Numerische Integration 131
7.1 Einführung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131
7.2 Newton-Cotes-Formeln . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131
7.3 Summierte Quadraturformeln . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135
7.4 Quadraturen höherer Ordnung . . . . . . . . . . . . . . . . . . . . . . . . . . . 137
7.5 Ausblick . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139
8 Iterative Lösung von Gleichungssystemen 141
8.1 Einführung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141
8.2 Newton-Verfahren . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142
8.3 Sukzessive Approximation (Fixpunktiteration) . . . . . . . . . . . . . . . . . . . 145
8.4 Iterationsverfahren zur Lösung linearer Gleichungssysteme . . . . . . . . . . . . 146
9 Schlussbemerkungen 156
2
1 Motivation
1.1 Die Wissenschaftliche Methode
• Experiment: Beobachte (und messe) was passiert.
• Theorie: Versuche die Beobachtung mit Hilfe von Modellen zu erklären.
• Theorie und Experiment werden sukzessive verfeinert und verglichen, bis eine akzeptable
Übereinstimmung vorliegt.
• In Naturwissenschaft und Technik liegen Modelle oft in Form mathematischer Gleichun-
gen vor, z. B. gewöhnliche oder partielle Differentialgleichungen.
• Oft können die Modellgleichungen nicht geschlossen (mit Papier und Bleistift, oder Ma-
thematica . . . ) gelöst werden.
⇒ Numerische Simulation und Optimierung
Simulation
Simulation: Gleichungen numerisch lösen.
• Undurchführbare Experimente ermöglichen (z.B. Galaxienkollisionen).
• Einsparung teurer Experimente (z.B. Windkanal).
• Parameterstudien schneller durchführbar.
• (Automatische) Optimierung von Prozessen.
• Identifikation von Modellparametern aus Messungen.
• Abschätzung von Unsicherheiten.
• Vielfältiger Einsatz in Naturwissenschaft, Technik und Industrie: Strömungsberechnung
(Wetter, Klima), Festigkeit von Bauwerken . . .
• Grundlage für alle diese Anwendungen sind numerische Algorithmen!
• Auch Informatiker sind beteiligt: Software-Engineering, High-Performance (Parallel) Com-
puting, . . .
Die prinzipielle Herangehensweise im Wissenschaftlichen Rechnen zeigt Abbildung 1. Die er-
folgreiche Durchführung einer Simulation erfordert die interdisziplinäre Zusammenarbeit von
Naturwissenschaftlern oder Ingenieuren mit Mathematikern und Informatikern.
3
1 Motivation
konzeptionelles Modell mathematisches Modell
wesentliche Prozesse algebraische Gleichungen
Differentialgleichungen
Realität Transport von Materie Wahrscheinlichkeiten
Wellenausbreitung
Reaktion Objekte: reelle Zahlen,
Phasenübergänge Funktionen, ...
...
?
Computerprogramm numerisches Modell
Komplexe SW Näherungsverfahren
SW−Engineering, Qualität zur Lösung oben
Simulation
Effizienz ("Teraflop") genannter Gleichungen
High Performance Comp.
Visualisierung
Abbildung 1: Prinzipielles Vorgehen im Wissenschaftlichen Rechnen.
4
1.2 Ein Beispiel
Fehlerquellen
Unterschiede zwischen Experiment und Simulation haben verschiedene Gründe:
• Modellfehler: Ein relevanter Prozess wurde nicht oder ungenau modelliert (Temp. kon-
stant, Luftwiderstand vernachlässigt, . . . )
• Datenfehler: Messungen von Anfangsbedingungen, Randbedingungen, Werte für Parame-
ter sind fehlerbehaftet.
• Abschneidefehler: Abbruch von Reihen oder Iterationsverfahren, Approximation von Funk-
tionen (z.B. stückweise Polynome).
• Rundungsfehler: Reelle Zahlen werden im Rechner genähert dargestellt.
Untersuchung von Rundungsfehlern und Abschneidefehlern ist ein zentraler Aspekt der Vorle-
sung!
1.2 Ein Beispiel
Pisa, 1582
Der Student Galileo Galilei sitzt in der Kirche und ihm ist langweilig. Er beobachtet den
langsam über ihm pendelnden Kerzenleuchter und denkt: „Wie kann ich nur die Bewegung
dieses Leuchters beschreiben?“.
Konzeptionelles Modell
Welche Eigenschaften (physikalischen Prozesse) sind für die gestellte Frage relevant?
• Leuchter ist ein Massenpunkt mit der Masse m.
• Der Faden der Länge ` wird als rigide und masselos angenommen.
• Der Luftwiderstand wird vernachlässigt.
Nun entwickle mathematisches Modell.
5
1 Motivation
(0, 0)
φ
`
m
Kräfte
• Annahme: Pendel läuft auf Kreisbahn: Nur Tangentialkraft ist relevant.
• Tangentialkraft bei Auslenkung φ:
cos(φ)
F~T (φ) = −m g sin(φ) .
sin(φ)
• Also etwa:
0
F~T (0) = −mg ,
0
0
F~T (π/2) = −mg .
1
• Vorzeichen kodiert Richtung.
(0, 0)
φ
`
F~T m
F~N
F~
6
1.3 Lösungsansätze
Weg, Geschwindigkeit, Beschleunigung
• Weg s(t), Geschwindigkeit v(t), Beschleunigung a(t) erfüllen:
ds(t) dv(t)
v(t) = , a(t) = .
dt dt
• Für den zurückgelegten Weg (mit Vorzeichen!) gilt s(t) = `φ(t).
• Also für die Geschwindigkeit
d s(φ(t)) d `φ(t) dφ(t)
v(t) = = =`
dt dt dt
• und die Beschleunigung
d v(φ(t)) d2 φ(t)
a(t) = =` .
dt dt2
Bewegungsgleichung
• Einsetzen in das 2. Newton’sche Gesetz m a(t) = F (t) liefert nun:
d2 φ(t)
m` = −mg sin(φ(t)) ∀t > t0 .
dt2
• Die Kraft ist hier skalar (vorzeichenbehafteter Betrag der Tangentialkraft), da wir nur
den zurückgelegten Weg betrachten.
• Ergibt gewöhnliche Differentialgleichung 2. Ordnung für die Auslenkung φ(t):
d2 φ(t) g
2
= − sin(φ(t)) ∀t > t0 . (1.1)
dt `
• Eindeutige Lösung erfordert zwei Anfangsbedingungen (t0 = 0):
dφ
φ(0) = φ0 , (0) = u0 . (1.2)
dt
1.3 Lösungsansätze
Lösung bei kleiner Auslenkung
• Allgemeine Gleichung für das Pendel ist schwer „analytisch“ zu lösen.
• Für kleine Winkel φ gilt
sin(φ) ≈ φ,
z.B. sin(0, 1) = 0, 099833417.
7
1 Motivation
• Diese Näherung reduziert die Gleichung auf
d2 φ(t) g
= − φ(t).
dt2 `
• Ansatz φ(t) = A cos(ωt) liefert mit φ(0) = φ0 , dφ
dt (0) = 0 dann die aus der Schule bekannte
Formel r
g
φ(t) = φ0 cos t (1.3)
`
Volles Modell; Verfahren 1
• Löse das volle Modell mit zwei numerischen Verfahren.
• Ersetze Gleichung zweiter Ordnung durch zwei Gleichungen erster Ordnung:
dφ(t) d2 φ(t) du(t) g
= u(t), = = − sin(φ(t)).
dt dt2 dt `
• Ersetze Ableitungen durch Differenzenquotienten:
φ(t + ∆t) − φ(t) dφ(t)
≈ = u(t),
∆t dt
u(t + ∆t) − u(t) du(t) g
≈ = − sin(φ(t)).
∆t dt `
• Mit φn = φ(n∆t), un = u(n∆t) erhält man Rekursion (Euler):
φn+1 = φn + ∆t un φ0 = φ0 (1.4)
un+1 = un − ∆t (g/`) sin(φn ) u0 = u0 (1.5)
Volles Modell; Verfahren 2
• Nutze Näherungsformel für die zweite Ableitung, sogenannter zentraler Differenzenquo-
tient:
φ(t + ∆t) − 2φ(t) + φ(t − ∆t) d2 φ(t) g
2
≈ 2
= − sin(φ(t)).
∆t dt `
• Auflösen nach φ(t + ∆t) ergibt Rekursionsformel (n ≥ 2):
φn+1 = 2φn − φn−1 − ∆t2 (g/`) sin(φn ) (1.6)
mit der Anfangsbedingung
φ0 = φ0 , φ1 = φ0 + ∆t u0 . (1.7)
(Die zweite Bedingung kommt aus dem Eulerverfahren oben).
8
1.4 Numerische Resultate
Konvergenz Differenzenquotient (Euler), phi=0.1
0.2
vereinfachtes Modell
dt=0.2
0.15 dt=0.1
dt=0.01
dt=0.001
0.1
0.05
Auslenkung
0
-0.05
-0.1
-0.15
-0.2
0 0.5 1 1.5 2 2.5 3 3.5 4
Zeit
Abbildung 2: Simulation des Fadenpendels (volles Modell) bei φ0 = 0.1 ≈ 5.7◦ mit dem Euler-
verfahren.
1.4 Numerische Resultate
Abbildung 2 zeigt das Eulerverfahren in Aktion. Für festen Zeitpunkt t und ∆t → 0 konvergiert
das Verfahren. Für festes ∆t und t → ∞ nimmt das Verfahren immer größere Werte an.
Abbildung 3 zeigt zum Vergleich das zentrale Verfahren für die gleiche Anfangsbedingung. Im
Unterschied zum expliziten Euler scheint der Fehler bei festem ∆t und t → ∞ nicht unbe-
schränkt zu wachsen.
Nun können wir das volle Modell mit dem vereinfachten Modell vergleichen und sehen welche
Auswirkungen die Annahme sin φ ≈ φ auf das Ergebnis hat. Abbildung 4 zeigt die numerische
Simulation. Selbst bei 28.6◦ ist die Übereinstimmung noch einigermaßen passabel. Für größere
Auslenkungen ist das vereinfachte Modell völlig unbrauchbar, die Form der Schwingung ist kein
Kosinus mehr.
1.5 Weitere Anwendungen
Eine Geothermieanlage
Wir betrachten die Modellierung und Simulation einer Geothermieanlage. Den schematischen
Aufbau zeigt die Abbildung 5. Kaltes Wasser fließt in einer Bohrung nach unten, wird erwärmt
und in einem isolierten Innenrohr wieder nach oben geführt.
• Grundwasserströmung gekoppelt mit Wärmetransport.
9
1 Motivation
Konvergenz zentriertes Verfahren, phi=0.1
0.15
vereinfachtes Modell
dt=0.2
dt=0.1
0.1 dt=0.01
dt=0.001
0.05
Auslenkung
0
-0.05
-0.1
-0.15
0 0.5 1 1.5 2 2.5 3 3.5 4
Zeit
Abbildung 3: Simulation des Fadenpendels (volles Modell) bei φ0 = 0.1 ≈ 5.7◦ mit dem zen-
tralen Verfahren.
• Welche Leistung erzielt so eine Anlage?
Modell für eine Geothermieanlage
• Strömung des Wassers in und um das Bohrloch
∇ · u = f,
K
u=− (∇p − %w G)
µ
• Transport der Wärme durch Konvektion und Wärmeleitung
∂(ce %e T )
+ ∇ · q + cw %w f T = g,
∂t
q = cw %w uT − λ∇T
in Abhängigkeit diverser Parameter: Bodendurchlässigkeit, Wärmekapazität, Dichte, Wär-
meleitfähigkeit, Pumprate sowie Rand- und Anfangsbedingungen.
Abbildung 5 zeigt die Entzugsleistung so einer Anlage für die ersten hundert Tage Betriebszeit.
Dabei wurden plausible Werte für die Modellparameter gewählt.
10
1.5 Weitere Anwendungen
Zentriertes Verfahren, phi=0.5
0.6
vereinfachtes Modell
dt=0.01
dt=0.0001
0.4
0.2
Auslenkung
0
-0.2
-0.4
-0.6
0 0.5 1 1.5 2 2.5 3 3.5 4
Zeit
Zentriertes Verfahren, phi=3.0
4
vereinfachtes Modell
dt=0.0001
3
2
1
Auslenkung
0
-1
-2
-3
-4
0 0.5 1 1.5 2 2.5 3 3.5 4
Zeit
Zentriertes Verfahren, phi=3.14
4
vereinfachtes Modell
dt=0.0001
3
2
1
Auslenkung
0
-1
-2
-3
-4
0 2 4 6 8 10
Zeit
Abbildung 4: Vergleich von vollem und vereinfachtem Modell (jeweils in rot) bei den Winkeln
φ = 0.5, 3.0, 3.14 gerechnet mit dem zentralen Verfahren.
11
1 Motivation
Ruecklauf Vorlauf
Wasserstrom
Waermestrom
Abbildung 5: Schematischer Aufbau einer Geothermieanlage.
Entzugsleistung im kalibrierten System
1.3e+06
1.2e+06
1.1e+06
Entzugsleistung [Watt]
1e+06
900000
800000
700000
600000
500000
400000
300000
0 10 20 30 40 50 60 70 80 90 100
Zeit [Tage]
Abbildung 6: Entzugsleistung der Geothermieanlage über die Zeit für bestimmte Systempara-
meter.
12
1.6 Ausblick
Abbildung 7: Ausbreitung eines mit Wasser nicht mischbaren Schadstoffes im Boden.
Schadstoffausbreitung
Als weiteres Beispiel nennen wir die Modellierung und Simulation einer Schadstoffausbreitung
im Boden, siehe Abbildung 7. Der Schadstoff wird hier als nicht mit Wasser mischbar und
schwerer als Wasser angenommen (Bestimmte Reinigungsmittel haben diese Eigenschaften).
• Wo erreicht der Schadstoff welche Konzentrationen?
• Wie bekommt man den Schadstoff wieder weg?
• Wohin bewegt sich gelöster Schadstoff?
1.6 Ausblick
Wir werden in dieser Vorlesung die folgenden Themengebiete behandeln:
• Grundbegriffe, Gleitpunktzahlen, Gleitpunktarithmetik
• Direkte Methoden zur Lösung linearer Gleichungssysteme
• Interpolation und Approximation
• Numerische Integration
• Iterationsverfahren zur Lösung linearer Gleichungssysteme
13
1 Motivation
• Iterationsverfahren zur Lösung nichtlinearer Gleichungssysteme
• Eigenwerte und Eigenvektoren
14
2 Fließkommazahlen
2.1 Zahlendarstellung
Stellenwertsystem
Definition 2.1 (Stellenwertsystem / Positionssystem). System, das Zahlen x ∈ R wie folgt
darstellt:
x = ± . . . m2 β 2 + m1 β + m0 + m−1 β −1 + m−2 β −2 . . .
X
= mi β i
i∈Z
β ∈ N, β ≥ 2 heißt Basis, mi ∈ {0, 1, 2, . . . , β − 1} heißen Ziffern
Geschichte: [Knuth, Band 2, S. 194]
• Babylonier (≈ −1750), β = 60
• Basis 10 ab ca. 1580
• Pascal: alle β ≥ 2 möglich
Festkommazahlen
Festkommazahlen: brich Reihe nach endlich vielen Einträgen ab
n
X
x=± mi β i
i=−k
Problem: wissenschaftl. Anwendungen brauchen Zahlen sehr unterschiedlicher Größe
Plancksches Wirkungsquantum h: 6.626093 · 10−34 Js
Avogadro-Konstante NA : 1
6.021415 · 1023 mol
Ruhemasse Elektron me : 9.109384 · 10−31 kg
Lichtgeschwindigkeit c: 2.997925 · 108 ms
Fließkommazahlen können unterschiedlich große Zahlen hinreichend genau und effizient dar-
stellen
15
2 Fließkommazahlen
Fließkommazahlen
Definition 2.2 (Fließkommazahlen). Sei β, r, s ∈ N und β ≥ 2. Das Fließkommagitter F(β, r, s) ⊂
R besteht aus den Zahlen mit folgenden Eigenschaften:
1. ∀x ∈ F(β, r, s) : x = m(x) · β e(x) mit
r
X s−1
X
−i
m(x) = ± mi β , e(x) = ± ej β j
i=1 j=0
mit Ziffern mi und ej . m heißt Mantisse, e heißt Exponent.
2. ∀x ∈ F(β, r, s) : x = 0 ∨ m1 6= 0. m1 6= 0 heißt Normierung und macht die Darstellung
eindeutig.
Eigenschaften
Ist x ∈ F(β, r, s) und x 6= 0, dann gilt wegen Punkt 2:
β −1 ≤ |m(x)| < 1 und damit β e(x)−1 ≤ |x| < β e(x) (2.1)
Größte / kleinste darstellbare Zahl in F(β, r, s):
s −1
X+/− = ±(1 − β −r ) · β β
Kleinste positive / größte negative Zahl:
s
x+/− = ±β −β
Beweis: Tafel
Beispiel
Beispiel 2.3. 1. F(10, 3, 1) besteht aus Zahlen der Form:
x = ±(m1 · 0.1 + m2 · 0.01 + m3 · 0.001) · 10±e0
mit m1 6= 0 ∨ (m1 = m2 = m3 = 0), z.B. 0, 0.999 · 104 , und 0.123 · 10−1 , aber nicht
0.140 · 10−10 , da Exponent zu klein
2. F(2, 2, 1) besteht aus Zahlen der Form:
x = ±(m1 · 0.5 + m2 · 0.25) · 2±e0
3 3 1 3 1 1 3 1 3 3
=⇒ F(2, 2, 1) = − , −1, − , − , − , − , 0, , , , , 1,
2 4 2 8 4 4 8 2 4 2
16
2.1 Zahlendarstellung
Abstände zwischen Fließkommazahlen
Sei x0 ∈ F die nächstliegende Fließkommazahl zu x ∈ F, x > 0.
Dann gilt:
|x|
1. Für m(x) 6= β −1 : |x − x0 | = |m(x)| β
−r
2. Für m(x) = β −1 : |x − x0 | = |x|β −r
Für den relativen Abstand gilt:
|x − x0 |
β −r ≤ ≤ β 1−r
|x|
Schwankt um Faktor β je nach Größe von |x|. Nennt man wobble. Daher: kleine β besser!
Beweis: Tafel
Standard: IEEE 754 / IEC 559
Ziel: Portabilität von Programmen mit Fließkommaarithmetik
Standard IEEE 754 / IEC 559, verabschiedet 1985
β = 2, mit vier Genauigkeitsstufen und normierter Darstellung:
single single-ext double double-ext
emax 127 ≥ 1024 1023 ≥ 16384
emin -126 ≤ -1021 -1022 ≤ -16381
Bits Expon. 8 ≤ 11 11 ≥ 15
Bits total 32 ≥ 43 64 ≥ 79
Der Standard definiert vier Rundungsarten:
nach −∞, nach +∞, nach 0, nearest
Seit 2008: zusätzlich half precision und quadruple precision
17
2 Fließkommazahlen
Double Precision
Betrachte double precision genauer:
• Insgesamt 64 bit
• 11 bit für Exponent, der vorzeichenlos als Zahl c ∈ [1, 2046] gespeichert wird
• Setze e := c − 1023 =⇒ e ∈ [−1022, 1023], kein Vorzeichen nötig
• Die Werte c ∈ {0, 2047} werden anderweitig genutzt:
– c = 0 ∧ m = 0 kodiert die Null
– c = 0 ∧ m 6= 0 kodiert denormalisierte Darstellung
– c = 2047 ∧ m = 0 kodiert ∞ (Überlauf)
– c = 2047 ∧ m 6= 0 kodiert NaN = “not a number”, z.B. bei Division durch Null
• 64 − 11 = 53 bit für Mantisse, eines für Vorzeichen, bleiben 52 Nachkommastellen
• Wegen β = 2 gilt immer m1 = 1
• Diese Ziffer heißt hidden bit und wird nicht gespeichert
• Somit r = 53 im Sinne der Definition von Fließkommazahlen
Bis auf die Spezialcodes entspricht double precision somit F(2, 53, 10).
2.2 Runden und Rundungsfehler
Rundungsfunktion
Um x ∈ R in F(β, r, s) zu approximieren, brauchen wir
rd : D(β, r, s) → F(β, r, s), (2.2)
wobei D(β, r, s) ⊂ R der Bereich ist, der F(β, r, s) enthält:
D := [X− , x− ] ∪ {0} ∪ [x+ , X+ ]
Achtung: das setzt voraus, dass x im darstellbaren Bereich liegt! Bei Über-/Unterlauf sind r,s
zu ändern.
Sinnvollerweise soll für rd gelten:
∀x ∈ D : |x − rd(x)| = min |x − y| (Bestapproximation)
y∈F
18
2.2 Runden und Rundungsfehler
Mit l(x) := max{y ∈ F|y ≤ x} und r(x) := min{y ∈ F|y ≥ x} gilt dann:
x l(x) = r(x), x ∈ F
l(x) |x − l(x)| < |x − r(x)|
rd(x) =
r(x) |x − l(x)| > |x − r(x)|
?
|x − l(x)| = |x − r(x)|
Im letzten Fall ist eine Rundung nötig. Dafür gibt es verschiedene Möglichkeiten.
Natürliche Rundung
P∞
Definition 2.4 (Natürliche Rundung). Sei x = sign(x) · −i β e die normierte Dar-
i=1 mi β
stellung von x ∈ D. Setze
( Pr −i β e
l(x) = sign(x) · i=1 mi β falls 0 ≤ mr+1 < β/2
rd(x) :=
r(x) = l(x) + β e−r (letzte Stelle) falls β/2 ≤ mr+1 < β
Das ist die übliche Rundung, die jeder kennt.
Gerade Rundung
Definition 2.5 (Gerade Rundung). Setze (mit Notation wie bei natürlicher Rundung)
l(x)
falls |x − l(x)| < |x − r(x)|
rd(x) := l(x) falls |x − l(x)| = |x − r(x)| ∧ mr gerade
r(x) sonst
Damit ist mr in rd(x) immer gerade, wenn gerundet wurde.
• Für rd(x) = l(x) ist das nach Definition so.
• Sonst rd(x) = r(x) = l(x) + β e−r , mr in l(x) ungerade, und Addition von β e−r ändert
die letzte Stelle um Eins.
Diese Wahl von Rundung vermeidet eine mögliche Drift beim Aufrunden. Entspricht “round to
nearest” im Standard.
19
2 Fließkommazahlen
Beispiel
Beispiel 2.6. Stelle die Zahl x = (106)10 in F(4, 3, 2) dar, d.h. konvertiere ins Vierersystem
(β = 4) und runde.
• log4 (106) = 3.36 Exponent 4, damit Mantisse < 1 ist
• Berechne r + 1 = 4 Ziffern, um runden zu können
m1 = 1, m2 = 2, m3 = 2, m4 = 2
• =⇒ (106)10 = 0.1222 · (10)10
4 (exakt, Ziffern zur Basis 4)
• Natürliche Rundung: (x0 )4 = 0.123 · 410 , da mr+1 = 2 = 4/2
• Gerade Rundung: (x0 )4 = 0.122 · 410 , da mr = 2 gerade
Absoluter und relativer Fehler
Definition 2.7 (absoluter und relativer Fehler). Sei x0 ∈ R eine Näherung von x ∈ R. Dann
heißt
∆x := x0 − x absoluter Fehler
und für x 6= 0
∆x
x0 := relativer Fehler
x
Umformen liefert:
0 ∆x
x = x + ∆x = x · 1 + = x · (1 + x0 )
x
Motivation
Motivation:
Es sei ∆x = x0 − x = 100 km.
Für x = Entfernung Erde—Sonne ≈ 1.5 · 108 km ist
102 km
x0 = ≈ 6.6 · 10−7
1.5 · 108 km
relativ klein.
Für x = Entfernung Heidelberg—Paris ≈ 460 km ist
102 km
x0 = ≈ 0.22 (22%)
4.6 · 102 km
dagegen relativ groß.
20
2.3 Fließkommaarithmetik
Fehlerabschätzung
Lemma 2.8 (Rundungsfehler). Bei der Rundung in F(β, r, s) gilt für den absoluten Fehler
1
|x − rd(x)| ≤ β e(x)−r (2.3)
2
und für den relativen Fehler (x 6= 0)
|x − rd(x)| 1
≤ β 1−r .
|x| 2
Diese Abschätzung ist scharf (d.h. der Fall “=” existiert).
Die Zahl eps := 12 β 1−r heißt Maschinengenauigkeit.
Beweis: Tafel
2.3 Fließkommaarithmetik
Wir benötigen eine Arithmetik auf F:
~ : F × F → F mit ~ ∈ {⊕, , , }
die den bekannten Operationen ∗ ∈ {+, −, ·, /} auf R entsprechen.
Problem: in der Regel x, y ∈ F =⇒
6 x∗y ∈F
Daher muss das Ergebnis wieder gerundet werden. Wir definieren
∀x, y ∈ F : x ~ y := rd(x ∗ y) (2.4)
Man sagt: ~ ist “exakt gerundet”. Das ist nicht trivial!
Guard digit
Beispiel 2.9 (Guard digit). Sei F = F(10, 3, 1), x = 0.215·108 , y = 0.125·10−5 . Wir betrachten
die Subtraktion x y = rd(x − y).
1. Subtraktion und anschließendes Runden benötigt extrem hohe Stellenzahl O(β s )!
2. Runden vor der Subtraktion liefert das selbe Ergebnis. Gute Idee?
3. Aber: betrachte z.B. x = 0.101 · 101 , y = 0.993 · 100
=⇒ relativer Fehler 18% ≈ 35 eps
4. Mit ein, zwei zusätzlichen Stellen erreicht man exakte Rundung!
5. Diese Stellen heißen guard digits und werden auch in der Praxis (CPU) genutzt.
21
2 Fließkommazahlen
Tabellenmacher-Dilemma
Algebraische Funktionen:
√
z.B. Polynome, 1/x, x, rationale Funktionen,. . .
grob: endliche Kombination der Grundrechenarten und Wurzeln
Transzendente Funktionen: der Rest, z.B. exp(x), ln(x), sin(x), xy , . . .
Tabellenmacher-Dilemma:
Es lässt sich nicht im Voraus entscheiden, wie viele guard digits für eine bestimmte Kombination
von transzendenter Funktion f und Argument x benötigt werden, um f (x) exakt gerundet zu
berechnen.
√
IEEE754 garantiert exakte Rundung für ⊕, , , , ebenso x.
Weitere Probleme / Eigenschaften
Die folgenden Punkte sind zu beachten:
• Assoziativ- und Distributivgesetz gelten nicht, es kommt auf die Reihenfolge von Opera-
tionen an!
• Es gibt y ∈ F, y 6= 0, so dass x ⊕ y = x
• Beispiel: ( ⊕ 1) 1=1 1 = 0 6= = ⊕ 0 = ⊕ (1 1)
• Das Kommutativgesetz gilt jedoch:
x ~ y = y ~ x für ~ ∈ {⊕, }
• Weitere einfache gültige Regeln:
– (−x) y = −(x y)
– 1 x=x⊕0=x
– x y = 0 =⇒ x = 0 ∨ y = 0
– x z≤y z falls x ≤ y ∧ z > 0
22
2.4 Fehleranalyse
2.4 Fehleranalyse
Rundungsfehler pflanzen sich in Rechnungen fort.
F1 (x1 , . . . , xm )
• Sei F : Rm → Rn , in Komponenten F (x) = ..
.
Fn (x1 , . . . , xm )
• Zur Berechnung von F im Rechner nutze numerische Realisierung F 0 : Fm → Fn .
F 0 wird durch einen Algorithmus realisiert, d.h. aus
– endlich vielen (= Terminierung)
– elementaren (= bekannten, d.h. ⊕, , , )
Rechenoperationen zusammengesetzt:
F 0 (x) = ϕl (. . . ϕ2 (ϕ1 (x)) . . . )
Wichtig:
1. Zu einem F gibt es in der Regel viele Realisierungen, im Sinne unterschiedlicher Reihen-
folgen
a + b + c ≈ (a ⊕ b) ⊕ c 6= a ⊕ (b ⊕ c)!
2. Jeder Schritt ϕi steuert einen (unbekannten) Fehler bei.
3. Im Prinzip kann die Rechengenauigkeit beliebig gesteigert werden, d.h. eigentlich Folge
m n
(F 0 )(k) : (F(k) ) → (F(k) ) . Das machen wir aber nicht so formal.
F (x) − F 0 (rd(x)) = F (x) − F (rd(x)) + F (rd(x)) − F 0 (rd(x)) (2.5)
| {z } | {z }
Konditionsanalyse Rundungsfehleranalyse
• F (x): exaktes Ergebnis
• F 0 (rd(x)): numerische Auswertung
• F (rd(x)): exaktes Ergebnis für rd(x) ≈ x
Von hier aus:
• Analyse “in erster Näherung”
• absolute / relative Fehler
• Normen bilden, das lassen wir aber i.d.R. weg
23
2 Fließkommazahlen
Differentielle Konditionsanalyse
Wir nehmen an, dass F : Rm → Rn zweimal stetig differenzierbar.
Nach dem Satz von Taylor gilt für die Komponenten Fi :
m
X ∂Fi
Fi (x + ∆x) = Fi (x) + (x)∆xj + RiF (x; ∆x) i = 1, . . . , n.
∂xj
j=1
Für das Restglied gilt unter diesen Voraussetzungen
RiF (x; ∆x) = O k∆xk2 ,
d.h. der Fehler der Näherung ist quadratisch in ∆x.
Definition 2.10 (Landausche Symbole). Man schreibt:
g(t) = O(h(t)), (t → 0)
falls es t0 > 0 und c ≥ 0 gibt, so dass für alle t ∈ (0, t0 ] die Abschätzung |g(t)| ≤ c · |h(t)| gilt.
Sprechweise: “g(t) geht wie h(t) gegen 0”
(quantifiziert also, “wie schnell” eine Funktion gegen Null geht)
Weiter bedeutet
g(t) = o(h(t)), (t → 0)
dass es t0 > 0 und eine Funktion c(t), limt→0 c(t) = 0 gibt, so dass für alle t ∈ (0, t0 ] die
Abschätzung |g(t)| ≤ c(t) · |h(t)| gilt.
Bedeutung: “g(t) geht schneller als h(t) gegen 0” (falls h(t) → 0)
Somit können wir die Taylorformel umformen:
m
X ∂Fi
Fi (x + ∆x) − Fi (x) = (x)∆xj + RiF (x; ∆x)
∂xj | {z }
j=1 höhere Ordnungen
| {z }
führende (erste) Ordnung
.
Oft lässt man die Terme höherer Ordnung weg und schreibt dann “ =” statt “=”.
Sprechweise: “ist in erster Näherung gleich”
24
2.4 Fehleranalyse
Damit gilt:
m
Fi (x + ∆x) − Fi (x) . X ∂Fi ∆xj
= (x) (2.6)
Fi (x) ∂xj Fi (x)
j=1
m
. X ∂Fi xj ∆xj
= (x) · ,
∂xj Fi (x) xj
j=1 | {z } | {z }
Verstärkungsfaktor kij (x) ≤eps
∆xj
d.h. die Verstärkungsfaktoren kij (x) geben an, wie die (relativen) Eingabefehler xj in den
(relativen) Fehler der i-ten Komponente von F eingehen!
Kondition
Definition 2.11 (Kondition). Wir nennen die Auswertung y = F (x) “schlecht konditioniert”
im Punkt x, falls |kij (x)| 1, andernfalls “gut konditioniert”.
|kij (x)| < 1 heißt Fehlerdämpfung, |kij (x)| > 1 heißt Fehlerverstärkung.
Das Symbol “ ” steht für “viel größer als”. Meist bedeutet das, dass die eine Zahl mehrere
Größenordnungen größer als die andere ist (z.B. 1 Mio. 1).
Definition stellt ein Kontinuum dar: es gibt einen fließenden Übergang von “gut konditioniert”
zu “schlecht konditioniert”!
Warum relative Kondition?
Wegen Lemma 2.8 (Rundungsfehler) gilt
x − rd(x) 1
≤ eps = β 1−r ,
x 2
d.h. es gibt ∈ R, || ≤ eps, so dass
x − rd(x)
= ,
x
d.h. für die relativen Eingabefehler in Gl. (2.6) gilt gerade
∆xj
= j
xj
25
2 Fließkommazahlen
Beispiel
∂F ∂F
Beispiel 2.12. 1. Addition: F (x1 , x2 ) = x1 + x2 , ∂x1
= ∂x2 = 1.
Nach obiger Formel:
F (x1 + ∆x1 , x2 + ∆x2 ) − F (x1 , x2 )
F (x1 , x2 )
. x1 ∆x1 x2 ∆x2
=1· +1·
x1 + x2 x1 x1 + x2 x2
| {z } | {z }
=k1 =k2
Schlecht konditioniert für x1 → −x2 !
∂F ∂F
2. F (x1 , x2 ) = x21 − x22 , ∂x 1
= 2x1 , ∂x 2
= −2x2 .
F (x1 + ∆x1 , x2 + ∆x2 ) − F (x1 , x2 )
F (x1 , x2 )
. x1 ∆x1 x2 ∆x2
= 2x1 · 2 + (−2x2 ) · 2
x1 − x22 x1 x1 − x22 x2
| {z } | {z }
=k1 =k2
2x21 2x2
=⇒ k1 = , k2 = − 2 2 2 ,
x21− x22 x1 − x2
schlecht konditioniert für |x1 | ≈ |x2 |.
Rundungsfehleranalyse
Sogenannte “Vorwärtsrundungsfehleranalyse”, rückwärts: später.
Nach Fehler-Aufspaltung, Gl. (2.5): F (x) − F 0 (x) mit x ∈ Fm ,
F 0 “zusammengesetzt” aus Einzeloperationen ~ ∈ {⊕, , , }
Wegen Gl. (2.4) (exakt gerundete Arithmetik) und Lemma 2.8 (Rundungsfehler) gilt
(x ~ y) − (x ∗ y)
= mit || ≤ eps
(x ∗ y)
Vorsicht, hängt von x und y ab und ist daher für jede Operation verschieden!
=⇒ x ~ y = (x ∗ y) · (1 + ) für ein |(~, x, y)| ≤ eps
26
2.4 Fehleranalyse
Beispiel
Beispiel 2.13. F (x1 , x2 ) = x21 − x22 mit zwei Realisierungen:
1. Fa (x1 , x2 ) = (x1 x1 ) (x2 x2 )
2. Fb (x1 , x2 ) = (x1 x2 ) (x1 ⊕ x2 )
Erste Realisierung:
u = x1 x1 = (x1 · x1 ) · (1 + 1 )
v = x2 x2 = (x2 · x2 ) · (1 + 2 )
Fa (x1 , x2 ) = u v = (u − v) · (1 + 3 )
Fa (x1 , x2 ) − F (x1 , x2 ) . x2 x2
= 2 1 2 (1 + 3 ) + 2 2 2 (2 + 3 )
F (x1 , x2 ) x1 − x2 x2 − x1
Zweite Realisierung:
u = x1 x2 = (x1 − x2 ) · (1 + 1 )
v = x1 ⊕ x2 = (x1 + x2 ) · (1 + 2 )
Fb (x1 , x2 ) = u v = (u · v) · (1 + 3 )
Fb (x1 , x2 ) − F (x1 , x2 ) . x21 − x22
= (1 + 2 + 3 ) = 1 + 2 + 3
F (x1 , x2 ) x21 − x22
=⇒ Zweite Realisierung besser als erste.
Numerische Stabilität
Definition 2.14 (Numerische Stabilität). Wir nennen einen numerischen Algorithmus “nume-
risch stabil”, wenn die im Lauf der Rechnung akkumulierten Rundungsfehler den unvermeidbaren
Problemfehler aus der Konditionsanalyse nicht übersteigen.
Mit anderen Worten:
Verstärkungsfaktoren aus Rundungsfehleranalyse ≤ denen aus Konditionsanalyse =⇒ “nume-
risch stabil”
Beide Realisierungen a, b aus Bsp. 2.13 sind numerisch stabil.
27
2 Fließkommazahlen
2.5 Die quadratische Gleichung
Die Gleichung
y 2 − py + q = 0
hat für p2 /4 > q 6= 0 die beiden reellen und verschiedenen Lösungen
r
p p2
y1,2 = f± (p, q) = ± − q. (Definiert zwei f !)
2 4
q
p2
Die Konditionsanalyse liefert mit D := 4 − q:
f (p + ∆p, q + ∆q) − f (p, q)
f (p, q)
.
p p ∆p q ∆q
= 1± −
2D p ± 2D p D (p ± 2D) q
Das heißt:
p2
• Für 4 q und p < 0 ist
r
p p2
f− (p, q) = − −q
2 4
gut konditioniert
p2
• Für 4 q und p > 0 ist
r
p p2
f+ (p, q) = + −q
2 4
gut konditioniert
2
• Für p4 ≈ q sind f+ und f− beide schlecht konditioniert, dieses Problem kann nicht
vermieden werden
p2
Numerisch geschickte Auswertung für den Fall 4 q:
p < 0:
q
p p2 q
Berechne y2 = 2 − 4 − q, y1 = y2 nach Satz von Vieta (q = y1 · y2 ).
p > 0:
q
p p2 q
Berechne y1 = 2 + 4 − q, dann y2 = y1 .
=⇒ Jedes Problem ist ein Spezialfall!
28
2.6 Auslöschung
2.6 Auslöschung
Die obigen Beispiele enthalten das Phänomen der Auslöschung.
Dies tritt auf bei
• Addition x1 + x2 mit x1 ≈ −x2
• Subtraktion x1 − x2 mit x1 ≈ x2
Bemerkung 2.15. Bei der Auslöschung werden vor der entsprechenden Addition bzw. Sub-
traktion eingeführte Fehler extrem verstärkt.
Sind x1 , x2 ∈ F Maschinenzahlen, so gilt wie oben gezeigt:
(x1 x2 ) − (x1 − x2 )
≤ eps.
(x1 − x2 )
Das ist also kein Problem. Das Problem tritt erst ein, wenn x1 und x2 selbst schon mit Fehlern
behaftet sind.
Beispiel
Beispiel 2.16. Betrachte F = F(10, 4, 1).
x1 = 0.11258762 · 102 , x2 = 0.11244891 · 102
=⇒ rd(x1 ) = 0.1126 · 102 , rd(x2 ) = 0.1124 · 102
x1 − x2 = 0.13871 · 10−1 , aber rd(x1 ) − rd(x2 ) = 0.2 · 10−1
Ergebnis hat keine gültige Ziffer! Relativer Fehler:
0.2 · 10−1 − 0.13871 · 10−1 1
−1
≈ 0.44 ≈ 883 · · 10−3 !
0.13871 · 10 |2 {z }
=eps
Grundregel
Im Beispiel: Fehler durch Rundung der Eingangsgrößen.
Herkunft der Fehler ist egal, tritt ebenso auf, wenn x1 , x2 mit Fehlern vorhergehender Rechen-
schritte behaftet sind.
Regel 2.17. Setze potentiell gefährliche Operationen möglichst früh im Algorithmus ein. (Ver-
gleiche Bsp. 2.13).
29
2 Fließkommazahlen
2.7 Die Exponentialfunktion
Die Exponentialreihe
Die Funktion exp(x) = ex lässt sich für alle x ∈ R als Potenzreihe darstellen:
∞
X xk
exp(x) =
k!
k=0
Pn xk
Naheliegender Ansatz: breche zur Berechnung nach n Gliedern ab, exp(x) ≈ k=0 k! .
Verwende Rekursionsvorschrift:
y0 := 1, S0 := y0 = 1,
x
∀k > 0 : yk := · yk−1 , Sk := Sk−1 + yk
k
yn : Glieder der Reihe, Sn : Partialsummen
Fehler für verschiedene Werte von x
Ergebnisse für Rekursionsformel mit float, n = 100:
1*1040
relativer Fehler von exp(x)
1*1035
1*1030
1*1025
1*1020
1*1015
1*1010
1*105
1*100
1*10-5
1*10-10
-50 -40 -30 -20 -10 0 10 20 30 40 50
• Für negative Werte für x wird der relative Fehler beliebig groß
• Dieser Effekt hat nichts mit dem Abbruch der Reihe zu tun!
30
2.7 Die Exponentialfunktion
Abweichung für imaginäre Argumente
1 2
Einheitskreis Einheitskreis
Auswertung exp(z) Auswertung exp(z)
1.5
0.5 1
0.5
0 0
-0.5
-0.5 -1
-1.5
-1 -2
-1 -0.5 0 0.5 1 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2
10
Einheitskreis Einheitskreis
4 Auswertung exp(z) Auswertung exp(z)
5
2
0 0
-2
-5
-4
-10
-4 -2 0 2 4 -10 -5 0 5 10
Ergebnisse für das imaginäre Intervall [−50, 50] · i
• Für |z| ≤ π ist das Ergebnis noch halbwegs brauchbar
• Für |z| → 2π wächst der Fehler weiter an
• Danach entfernen sich die Werte vom Kreis (die Trajektorie nähert sich einer Geraden an
und kommt auch nicht zurück)
31
2 Fließkommazahlen
Visualisierung des Konvergenzverhaltens
exp(2πi · 0.1)
4 exp(2πi · 0.2)
exp(2πi · 0.3)
exp(2πi · 0.4)
exp(2πi · 0.5)
Einheitskreis
2
0
-2
-4
-4 -2 0 2 4
• gerade Potenzen tragen zum Realteil von exp(2πi · x) bei
• ungerade Potenzen tragen zum Imaginärteil bei
=⇒ Addieren eines Gliedes ändert abwechselnd Real- und Imaginärteil
exp(2πi · 0.6)
40 exp(2πi · 0.7)
exp(2πi · 0.8)
exp(2πi · 1.0)
exp(2πi · 1.0)
Einheitskreis
20
0
-20
-40
-40 -20 0 20 40
• Betrag der Zwischenergebnisse wächst exponentiell in x
32
2.8 Rekursionsformel für Integrale
• Form der Trajektorie nähert sich immer mehr Quadrat an
=⇒ Auslöschung
Konditionsanalyse für exp(x)
Für die Funktion exp gilt wegen exp0 = exp:
exp(x + ∆x) − exp(x) . x ∆x
= exp0 (x) · = ∆x
exp(x) exp(x) x
=⇒ absoluter Fehler von x wird zu relativem Fehler von exp(x)
(Vergleiche: exp ist Isomorphismus zwischen (R, +) und (R+ , ·).)
Wegen k = x ist exp gut konditioniert, falls x nicht zu groß
=⇒ betrachteter Algorithmus ist instabil für x < 0
Gibt es stabileren Algorithmus? Übungsaufgabe
2.8 Rekursionsformel für Integrale
Integrale der Form
Z 1
Ik = xk exp(x) dx
0
kann man mithilfe einer Rekursionsformel berechnen:
I0 = e − 1, ∀k > 0 : Ik = e − k · Ik−1
Für das erste Glied ist wegen exp0 (x) = exp(x) eine Stammfunktion bekannt, die weiteren
Glieder kann man dann über die obige Formel bestimmen.
Wie gut funktioniert das in der Praxis?
Die ersten 26 Werte für Ik , mit endlicher Genauigkeit berechnet:
33
2 Fließkommazahlen
k berechnetes Ik Fehler |∆Ik | k berechnetes Ik Fehler |∆Ik |
0 1.718281828459050 2.6 ·10−15 13 0.181983054536145 3.3 · 10−7
1 1 (Null) 14 0.170519064953013 4.6 · 10−6
2 0.718281828459045 1.5 · 10−15 15 0.160495854163853 7.0 · 10−5
3 0.563436343081910 5.5 · 10−16 16 0.150348161837404 1.1 · 10−3
4 0.464536456131406 1.0 · 10−15 17 0.162363077223183 1.9 · 10−2
5 0.395599547802016 6.0 · 10−15 18 -0.204253561558257 3.4 · 10−1
6 0.344684541646949 3.8 · 10−14 19 6.59909949806592 6.7 · 100
7 0.305490036930402 2.7 · 10−13 20 -129.263708132859 1.3 · 101
8 0.274361533015832 2.1 · 10−12 21 2717.25615261851 2.7 · 103
9 0.249028031316559 1.9 · 10−11 22 -59776.9170757787 6.0 · 104
10 0.228001515293454 1.9 · 10−10 23 1374871.81102474 1.4 · 106
11 0.210265160231056 2.1 · 10−9 24 -32996920.7463119 3.3 · 107
12 0.195099905686377 2.5 · 10−8 25 824923021.376079 8.2 · 108
Rekursionsformel Ik = e − k · Ik−1 führt zu Fehlerverstärkung um Faktor k im k-ten Schritt!
Bessere Optionen
1. Alle Ik sind von der Form a · e + b, mit a, b ∈ Z. Berechne diese Zahlen über die Rekursi-
onsformel, und werte erst ganz am Ende mit Fließkommazahlen aus.
2. Drehe die Rekursionsformel um: wenn Ik → Ik+1 den Fehler um k verstärkt, wird er
durch Ik+1 → Ik um den Faktor k reduziert!
Wegen 0 ≤ xk ≤ 1 und 0 ≤ exp(x) ≤ 3 auf [0, 1] muss 0 ≤ Ik ≤ 3 gelten. Wenn wir also z.B.
I50 := 1.5 setzen, kann der Fehler maximal 1.5 sein.
Benutze dann umgedrehte Rekursionsformel
Ik = (k + 1)−1 · (e − Ik+1 ).
Rekursionsformel für Integrale
Die Werte für Ik zwischen k = 25 und 50, rückwärts berechnet:
34
2.8 Rekursionsformel für Integrale
k berechnetes Ik Fehler |∆Ik | k berechnetes Ik Fehler |∆Ik |
50 1.5 1.4 ·100 37 0.0697442966294832 2.8 · 10−16
49 0.0243656365691809 2.9 · 10−2 36 0.0715820954548530 1.9 · 10−16
48 0.0549778814671401 5.9 · 10−4 35 0.0735194370278942 2.8 · 10−17
47 0.0554854988956647 1.2 · 10−5 34 0.0755646397551757 2.2 · 10−16
46 0.0566552410545400 2.6 · 10−7 33 0.0777269761383491 2.7 · 10−18
45 0.0578614475522718 5.7 · 10−9 32 0.0800168137066878 1.8 · 10−16
44 0.0591204529090394 1.3 · 10−10 31 0.0824457817110112 2.6 · 10−16
43 0.0604354858079547 2.9 · 10−12 30 0.0850269692499366 2.9 · 10−16
42 0.0618103800616533 6.7 · 10−14 29 0.0877751619736370 4.5 · 10−17
41 0.0632493201999379 1.8 · 10−15 28 0.0907071264305313 4.3 · 10−16
40 0.0647568904453441 1.4 · 10−16 27 0.0938419536438755 4.7 · 10−16
39 0.0663381234503425 2.1 · 10−16 26 0.0972014768450063 1.3 · 10−17
38 0.0679985565386847 1.5 · 10−16 25 0.1008107827543860 1.1 · 10−16
Trotz völlig unbrauchbarer Schätzung für den Startwert I50 führt die neue Rekursionsformel
Ik = (k + 1)−1 · (e − Ik+1 ) sehr schnell zu brauchbaren Ergebnissen!
Idee der Fehlerabschätzung
Nach Fehleranalyse gilt für Startwert Ik+m , m ≥ 1:
k! k!
|∆Ik | ≈ |∆Ik+m | ≤ · 1.5 ≤ (k + 1)−m · 1.5
(k + m)! (k + m)!
Idee: berechne benötigte Anzahl Schritte m aus gewünschter Genauigkeit |∆Ik | < tol.
tol
(k + 1)−m · 1.5 < tol =⇒ exp(−m · ln(k + 1)) <
1.5
tol
ln(tol) − ln(1.5)
=⇒ −m · ln(k + 1) < ln =⇒ m >
1.5 ln(k + 1)
Beispiel: k = 25, tol = 10−8 =⇒ m > 5.7
Ergebnis für m = 6:
k berechnetes Ik Fehler |∆Ik |
31 1.5 1.4 · 100
30 0.0392994138212595 4.6 · 10−2
29 0.0892994138212595 1.5 · 10−3
28 0.0906545660219926 5.3 · 10−5
27 0.0938438308013233 1.9 · 10−6
26 0.0972014073206564 7.0 · 10−8
25 0.1008107854284000 2.9 · 10−9
35
2 Fließkommazahlen
• Umgedrehte Rekursionsformel ist im Gegensatz zum naiven Ansatz stabil
• Fehlerabschätzung minimiert Aufwand für vorgegebene Genauigkeit
=⇒ stabil und effizient
36
3 Motivation linearer Gleichungssysteme
3.1 Strömung in Rohrleitungsnetzen
Netzwerk aus Röhren
Pumpe
Pumpe
Betrachte Netzwerk aus Röhren, in denen Flüssigkeit zirkuliert
• Um es einfach zu halten: nur durch Pumpen getrieben, keine Quellen und Senken im
Kreislauf vorhanden
• Mögliche Anwendungen: Trinkwasserversorgung, elektrische Schaltkreise, Blutkreislauf
(als Näherung)
Pumpe
Pumpe
Beschrieben durch gerichteten, (schwach) zusammenhängenden Graphen G = (V, E):
• Knotenmenge V = {v1 , . . . , vn }, |V | = n
• Kantenmenge E = {e1 , . . . , eM } ⊂ V × V, |E| = M
• gerichtet “von v nach w”: (v, w) ∈ E =⇒ (w, v) ∈
/E
• Röhren ER = {e1 , . . . , em }, Pumpen EP = {em+1 , . . . , eM }
37
3 Motivation linearer Gleichungssysteme
Gesetz von Hagen-Poiseuille
Fluss durch lange zylindrische Röhre e = (v, w):
πre4
qe = ∆pe
8ρde
qe : gerichteter Volumenstrom [m3 /s]
∆pe : gerichtete Druckdifferenz [Pa = N/m2 ]
re : Radius der Röhre e [m]
de : Länge der Röhre e [m]
ρ: dynamische Viskosität [Pa s]
qe > 0, ∆pe > 0 falls Fluss von v nach w, qe < 0, ∆pe < 0 sonst
3.2 Die Kirchhoffschen Gesetze
Quelle: wikipedia.org, gemeinfrei
Gustav Robert Kirchhoff (1824–1887):
• deutscher Physiker
• Beiträge zur Erforschung von
– elektrischen Schaltkreisen
38
3.2 Die Kirchhoffschen Gesetze
– Spektralanalyse
– Schwarzkörperstrahlung
• Kirchhoffsche Gesetze (Erhaltungsregeln für Schaltkreise)
Knotenregel (1. Kirchhoffsches Gesetz)
1
2
5 3
4
q1 + q4 + q5 = q2 + q3
Knotenregel (Erhaltung der Ladung):
Die Summe der zufließenden Ströme ist für jeden Knoten eines elektrischen Netzwerks gleich
der Summe der abfließenden Ströme.
Hier: Massenerhaltung
Flüssigkeit wird in den einzelnen Knoten weder vernichtet noch zwischengespeichert.
Definiere für Knoten v ∈ V eingehende und ausgehende Kanten:
Ev+ := {(u, w) ∈ E | w = v}, Ev− := {(u, w) ∈ E | u = v}
Die Trennung in ein- und ausgehend basiert auf der Orientierung der Kanten und hat nichts
mit dem physikalischen Fluss zu tun!
Knotenregel: X X
qe − qe = 0 ∀v ∈ V (3.1)
e∈Ev+ e∈Ev−
39
3 Motivation linearer Gleichungssysteme
Das definiert n − 1 = |V | − 1 linear unabhängige Bedingungen, da sich aufgrund globaler
Massenerhaltung die Erhaltung im letzten Knoten jeweils automatisch ergibt.
Beweis: Tafel
Maschenregel (2. Kirchhoffsches Gesetz)
1
2
5
3
4
∆p3 + ∆p2 = ∆p4 − ∆p5 + ∆p1
Maschenregel (Energieerhaltung):
Die Summe der Spannungsdifferenzen entlang eines geschlossenen Pfades ist Null. Anders ge-
sagt: Die Spannungsdifferenz zwischen zwei Knoten ist wegunabhängig.
(Sonst könnte man entlang des Pfades beliebig Energie erzeugen =⇒ Perpetuum Mobile)
Hier: Energieerhaltung
Es gibt keine spontanen (d.h. nicht durch Pumpen verursachte) Kreisbewegungen von Flüssig-
keit.
Betrachte C ⊂ E, c beliebiger geschlossener Pfad. Definiere:
C + := {e ∈ C | e genauso wie C orientiert}
C − := {e ∈ C | e entgegen C orientiert}
(Anmerkung: wir vernachlässigen in der Notation, dass Kanten in C mehrfach vorkommen
dürfen, evtl. mit geänderter Orientierung (!))
40
3.3 Das Knotenpotentialverfahren
Maschenregel: X X
∆pe − ∆pe = 0 ∀ geschl. Pfade C ⊂ E (3.2)
e∈C + e∈C −
3.3 Das Knotenpotentialverfahren
Grundidee des Knotenpotentialverfahrens:
1. Wegen der Maschenregel ist Potentialdifferenz zwischen je zwei Knoten wegunabhängig
und somit wohldefiniert
Druckdifferenz zu (beliebigem) Referenzknoten als Variablen
2. Da das Potential nur bis auf Konstante eindeutig ist, können wir die Referenz als Null
definieren
n − 1 Knotendrücke pv als Unbekannte
3. Wegen der Knotenregel gibt es n − 1 unabh. Nebenbedingungen
n − 1 lineare Gleichungen für die Unbekannten
Stelle System auf und löse nach Unbekannten (Druckdifferenzen) auf, daraus kann man dann
die Flüsse berechnen.
Wähle Referenzknoten r, ergibt Druckdifferenzen ∆pe der Röhren:
pw − pv
v 6= r ∧ w 6= r
∀e = (v, w) ∈ ER : ∆pe = −pv w = r (setze pw = 0)
pw v = r (setze pv = 0)
Wähle feste Anordnung ek , k ∈ {1, . . . , m}, vi , i ∈ {1, . . . , n} der Kanten und Knoten mit Wahl
r = vn
=⇒ p = (pv1 , . . . , pvn−1 )T , ∆p = (∆pe1 , . . . , ∆pem )
Als lineares Gleichungssystem:
[∆p]m = [B]m×(n−1) · [p]n−1
Berechne resultierende Flüsse über Leitfähigkeiten (Hagen-Poiseuille):
∀e ∈ ER : qe = Le · ∆pe
41
3 Motivation linearer Gleichungssysteme
Als lineares Gleichungssystem mit Diagonalmatrix L:
[q]m = [L]m×m · [∆p]m
Setze Umrechnung von p zu ∆p ein:
[q]m = [L]m×m · [B]m×(n−1) · [p]n−1 = [L · B]m×(n−1) · [p]n−1
Verwende Knotenregel Gl. (3.1):
X X X X
qe − qe = qe − qe
e∈Ev+ ∩ER e∈Ev− ∩ER e∈Ev− ∩EP e∈Ev+ ∩EP
Als lineares Gleichungssystem formuliert:
[B T ](n−1)×m · [q]m = [b]n−1 ,
wobei b = (bv1 , . . . , bvn−1 )T die Pumpströme enthält (die Einträge sind Null für innere Kno-
ten)
Es handelt sich um B T , weil eine Kante e genau dann zur Bilanz von v beisteuert, wenn v (mit
gleichem Vorzeichen) zur Druckdifferenz entlang e beiträgt
Kleines Beispiel
1 1 2
2 3 4
3 5 4
42
3.3 Das Knotenpotentialverfahren
Resultierende Matrix B für Beispielnetzwerk:
−1 1 0 0 −1 1 0
−1 0 1 0 −1 0 1
0 −1 1 0 B = 0 −1 1
0 −1 0 1 0 −1 0
0 0 −1 1 0 0 −1
• Zeilen: Knoten, die von i-ter Kante verbunden werden (mit Vorzeichen)
• Spalten: Kanten, die für j-ten Knoten ein- / ausgehend sind (mit Vorzeichen)
• Streiche letzte Spalte wegen p4 := 0 (für B) bzw. wegen redundanter Knotenregel (für
BT )
Knotenpotentialverfahren
Insgesamt:
[b]n−1 = [B T ](n−1)×m · [L · B]m×(n−1) · [p]n−1
= [B T · L · B](n−1)×(n−1) · [p]n−1
Eigenschaften des Systems: A := B T LB ist
• symmetrisch und positiv definit
• damit insbesondere invertierbar
• dünn besetzt (Anzahl Kanten meist O(n), nicht O(n2 ))
Beweis: Tafel
Anmerkungen
Das Ohmsche Gesetz besagt:
Ie = Re−1 (Uv − Uw )
mit Ie dem elektrischen Strom entlang Leiter e, Re dem Widerstand entlang e, und Uv dem
elektrischen Potential im Knoten v.
=⇒ die hergeleitete Methode lässt sich völlig analog auf elektrische Netzwerke anwenden
Falls das Netzwerk Spulen oder Kondensatoren enthält (RLC-Netzwerk), werden die Zustands-
größen komplex:
Ax = b mit x, b ∈ Cn , A ∈ Cn×n
43
3 Motivation linearer Gleichungssysteme
3.4 Partielle Differentialgleichungen
Grenzübergang zum Kontinuum
Wenn man das Röhren- / Leitungsnetz immer feiner und feiner macht, beschreibt man im
Limes ein unendlich feines Netzwerk, also ein Kontinuum. Die zugehörige Gleichung ist die
Poisson-Gleichung:
n n
X X ∂2u
−∆u := −∇ · (∇u) := − ∂x2i u := − =f
i=1 i=1
∂x2i
Das ist ebenfalls ein lineares Gleichungssystem, allerdings über einem unendlich-dimensionalen
Funktionenraum.
Anwendungen der Gleichung sind zum Beispiel:
• Elektrostatik im Vakuum (f : Ladungsdichte, u: resultierendes Potential)
• Gravitation (f : Massendichte, u: Gravitationspotential)
Poröse Medien
Falls die Leitfähigkeit κ (Inverse des Widerstandes) ortsabhängig ist, ergibt sich die Darcy-
Gleichung:
Xn
−∇ · (κ∇u) = − ∂xi (κ∂xi u) = f
i=1
Im Allgemeinen lässt sich diese Gleichung nicht geschlossen lösen, daher verwendet man numeri-
sche Methoden. Die Diskretisierung führt dabei auf beliebig große lineare Gleichungssysteme.
Mögliche Anwendungen:
• Stationäre Strömungen in porösen Medien (z.B. Grundwasser)
• Elektrostatik in Bauteilen und Materialien
• Wärmeleitung durch heterogene Festkörper
44
3.5 Radiosity-Methode in der Computergraphik
3.5 Radiosity-Methode in der Computergraphik
Raytracing
Aufgabe: stelle beleuchtete Szene graphisch dar (Computergraphik).
Grundidee des Raytracing-Algorithmus:
Image
Camera Light Source
View Ray
Shadow Ray
Scene Object
Quelle: Henrik, wikipedia.org, cc-by-sa
Beispiele von rekursivem Raytracing:
Quelle: Tim Babb, wikipedia.org, cc-by-sa Quelle: Greg L, wikipedia.org, cc-by-sa
45
3 Motivation linearer Gleichungssysteme
Raytracing vs. Radiosity-Methode
Idee von Raytracing:
• Verfolge “Sehstrahlen” vom Auge des Betrachters zu ihrem Ziel
• Berücksichtige dabei gerichtete Reflexion, diffuse Reflexion, Brechung, Wurf von Schatten
Vorteile / Nachteile von (rekursivem) Raytracing:
• Viele verschiedene optische Effekte darstellbar
• Aufwand wächst exponentiell mit der Anzahl der Rekursionen
• Ergebnis ist auf Blickwinkel festgelegt (!)
Idee der Radiosity-Methode:
• Berechne das Strahlungsgleichgewicht für die Szene
• Beleuchtete Oberflächen erzeugen passive Beleuchtung, indem sie Licht in den Raum
zurück werfen
Vorteile / Nachteile der Radiosity-Methode:
• Weniger Effekte berücksichtigt (z.B. nur diffuse Streuung)
• Aufwand ist auch bei dieser Methode sehr hoch
• Ergebnis ist unabhängig vom Betrachter ( Animationen)
Radiosity-Methode
Anwendung der Radiosity-Methode:
Quelle: Hugo Elias, wikipedia.org, cc-by-sa
46
3.5 Radiosity-Methode in der Computergraphik
Quelle: Hugo Elias, wikipedia.org, cc-by-sa
Definiere Oberfläche der Szene:
S = {x ∈ R3 | x auf Oberfläche eines Objekts}
(Lichtquellen sind dabei ebenfalls Objekte)
Bestimme “Energiedichte” B : S → R:
Z
B dA = von ω ⊂ S abgestrahlte Energie
ω
E(x): Eigenstrahlung (bei Lichtquellen), Senke (bei Absorption)
Benötigte Materialeigenschaften:
ρ(x): Reflexionsfaktor
Bestimmungsgleichung
Bestimmungsgleichung für B(x), x ∈ S:
Z
cos ϕx,x0 cos ϕx0 ,x
B(x) = E(x) + ρ(x) B(x0 ) 0 k2
V (x, x0 ) dx0
S kx − x
| {z }
“Kern”, “kernel” λ(x, x0 )
mit
ϕa,b = ](νa , b − a)
(
1 Licht aus Richtung νa
cos ϕa,b =
0 tangential zur Oberfläche
(
1 x0 von x aus sichtbar
V (x, x0 ) =
0 sonst
(Eigenschaften der beleuchteten Szene)
47
3 Motivation linearer Gleichungssysteme
νa
a b
Integralgleichung
Resultierende Integralgleichung:
Z
B(x) − ρ(x) B(x0 )λ(x, x0 ) dx0 = E(x)
S
Alternative Formulierung mit Dirac-Delta δ(x, x0 ):
Z
B(x0 ) δ(x, x0 ) − ρ(x)λ(x, x0 ) dx0 = E(x)
S
• Aufgabe: finde für gegebene Quellen / Senken E(x) resultierende Energiedichte B(x)
• System ist erneut linear in Unbekannten B(x), aber diesmal unendlichdimensional!
Kollokationsmethode
Zerlege die Oberfläche in eine endliche Menge von Teilgebieten Th = {t1 , . . . , tn }:
n
[
ti ⊂ S offen, ti ∩ tj = ∅ ∀i 6= j, ti = S
i=1
Dieser Vorgang heißt Diskretisierung und ist üblich bei Integral- und Differentialgleichungen.
Definiere außerdem xi jeweils als den Mittelpunkt von ti .
(Existenz eines solchen Punktes schränkt implizit Wahl der ti ein, aber wir sind sowieso nur an
“simplen” Teilgebieten ti interessiert!)
48
3.5 Radiosity-Methode in der Computergraphik
Approximiere B : S → R durch diskrete Funktion Bh : S → R:
m
X
Bh (x) = zj ϕj (x)
j=1
mit zj ∈ R Koeffizienten und ϕj : S → R Basisfunktionen. Insbesondere: ϕ1 , . . . , ϕm linear
unabhängig.
Wir verwenden stückweise konstante Funktionen (m = n):
(
1 falls x ∈ tj
ϕj (x) =
0 sonst
Erfülle Integralgleichung nur für x ∈ Xh = {x1 , . . . , xn }:
Z
Bh (xi ) − ρ(xi ) Bh (x0 )λ(xi , x0 ) dx = E(xi ) i = 1, . . . , n
S
n
X n
Z X
⇐⇒ zj ϕj (xi ) − ρ(xi ) zj ϕj (x0 )λ(xi , x0 ) dx = E(xi )
j=1 S j=1
Xn Z
0 0
⇐⇒ zj ϕj (xi ) − ρ(xi ) ϕj (x )λ(xi , x ) dx = E(xi )
j=1 S
| {z }
=:aij
⇐⇒ Az = b
mit A = [aij ]n,n n n
i,j=1 , z = [zj ]j=1 , b = [E(xi )]i=1 .
Anmerkungen
• Das Integral in aij wird i.d.R. ebenfalls numerisch berechnet (Stoff der Vorlesung).
• Man begeht dabei Diskretisierungsfehler
kBh − Bk = O(hα )
mit k·k Norm auf Funktionenräumen, h = maxt∈Th diam(t), und α “Konvergenzordnung”.
=⇒ bessere Approx. Bh von B bei feinerer Unterteilung Th
=⇒ Gleichungssysteme können beliebig groß werden
• Die Matrix A ist in diesem Fall nicht dünn besetzt, sondern voll besetzt! Das ist wichtig
bei der Wahl des Lösers.
49
4 Kondition linearer Gleichungssysteme
4 Kondition linearer Gleichungssysteme
4.1 Lösbarkeit linearer Gleichungssysteme
Gegeben seien Matrix A ∈ Km×n und Vektor b ∈ Km :
a11 · · · a1n b1
.. .. ,
A = [aij ]m,n
i,j=1 = .
..
. . b = [bi ]m = ...
i=1
am1 · · · amn bm
Gesucht ist [xj ]nj=1 = [x1 , . . . , xn ]T ∈ Kn , so dass
n
X
∀i : aij xj = bi ,
j=1
beziehungsweise in Kurzschreibweise:
Ax = b (4.1)
Im folgenden ist stets K = R oder K = C.
Das Gleichungssystem (4.1) heißt
• unterbestimmt, falls m < n
• quadratisch, falls m = n
• überbestimmt, falls m > n
Es gibt mindestens eine Lösung, falls
a11 ··· a1n b1
.. .. .. ..
rang(A) = rang([A|b]) = rang . . . .
am1 ··· amn bm
(Rang = Dimension Spaltenraum = Dimension Zeilenraum)
Für quadratische Matrizen sind folgende Aussagen äquivalent:
1. Ax = b ist für jedes b eindeutig lösbar
2. rang(A) = n (Vollrang)
3. det(A) 6= 0 (regulär)
4. A hat keinen Eigenwert Null
50
4.2 Vektornormen
4.2 Vektornormen
Die Quantifizierung von Fehlern erfordert Normen.
Definition 4.1 (Vektornorm). Eine Abbildung k · k : Kn → R+ heißt Norm, falls gilt:
(N1) Definitheit: kxk > 0 für x 6= 0
(N2) Positive Homogenität: kα · xk = |α| · kxk, x ∈ Kn , α ∈ K
(N3) Subadditivität: kx + yk ≤ kxk + kyk, x, y ∈ Kn (Dreiecksungleichung)
Häufig verwendete Normen
Beispiel 4.2 (Häufig verwendete Normen). Euklidische Norm, `2 -Norm:
n
!1/2
X
kxk2 = |xi |2
i=1
Maximumsnorm, `∞ -Norm:
kxk∞ = max |xi |
1,...,n
`1 -Norm:
n
X
kxk1 = |xi |
i=1
Folgerungen
Konvergenz in der Maximumsnorm entspricht komponentenweiser Konvergenz:
(t)
∀i = 1, . . . , n : lim xi = xi ⇐⇒ lim kx(t) − xk∞ = 0
t→∞ t→∞
Gilt auch im Unendlichdimensionalen, heißt dort aber aus naheliegenden Gründen Supremums-
norm.
Direkte und wichtige Folge von (N3) ist die inverse Dreiecksungleichung für Normen k · k:
kx − yk ≥ |kxk − kyk|
=⇒ Stetigkeit der Norm als Abbildung von Kn nach R+
51
4 Kondition linearer Gleichungssysteme
Äquivalenz aller Normen
Satz 4.3 (Äquivalenz aller Normen). Seien k · k, k · k0 Normen auf dem endlich-dimensionalen
(!) Kn . Dann gibt es Zahlen m, M > 0 aus R, so dass gilt:
m · kxk0 ≤ kxk ≤ M · kxk0 ∀x ∈ Kn
Anmerkung: die Zahlen m und M sind nicht nur von den betrachteten Normen, sondern auch
von der Dimension n abhängig. Insbesondere gilt die Aussage im Unendlichen nicht!
Beweis: Tafel
Schlussfolgerung
Im Endlichdimensionalen entspricht Konvergenz bezgl. einer beliebigen Norm k · k also Kon-
vergenz in der Maximumsnorm, so dass wir sagen können:
Konvergenz bzgl. einer Norm k · k entspricht komponentenweiser Konvergenz:
(t)
∀i = 1, . . . , n : lim xi = xi ⇐⇒ lim kx(t) − xk = 0
t→∞ t→∞
Es gibt also für endlichdimensionale Vektorräume nur einen Konvergenzbegriff. Das ist für nicht
endlich erzeugte Vektorräume i.A. falsch!
4.3 Matrixnormen
Normen auf Kn×m
Der Kn×m ist ein Vektorraum, der mit Kn·m identifiziert werden kann:
Kn×m ∼
= Kn·m
Daher definiert jede Norm auf Kn·m auch eine Norm auf Kn×m .
(t)
A(t) → A (t → ∞) meint dann aij → aij (t → ∞) ∀ 1 ≤ i, j ≤ n
Nicht jede dieser Normen ist mit der zusätzlichen Verknüpfung, der Matrixmultiplikation, ver-
träglich. Wir betrachten daher Normen mit zusätzlichen Eigenschaften.
52
4.3 Matrixnormen
Matrixnormen
Definition 4.4 (Matrixnormen). Eine Norm k · k auf Kn×n heißt verträglich mit einer Vek-
tornorm k · k auf Kn , falls gilt:
kAxk ≤ kAk · kxk A ∈ Kn×n , x ∈ Kn
Sie heißt Matrixnorm, wenn sie submultiplikativ ist:
kABk ≤ kAk · kBk A, B ∈ Kn×n
Die Submultiplikativität ergänzt die Subadditivität (Dreiecksungleichung), die für alle Normen
gilt.
Beispiel
Beispiel 4.5. Die Frobeniusnorm
1/2
n
X
kAkFr := |aij |2
i,j=1
ist eine Matrixnorm, die mit der euklidischen Norm verträglich ist.
Die Eigenschaften einer Matrixnorm (Definitheit, positive Homogenität, Subadditivität, Sub-
multiplikativität) und die Verträglichkeit kann man direkt nachrechnen ( Übung).
Zugeordnete Matrixnormen
Definition 4.6 (Zugeordnete Matrixnorm). Es sei k · k eine beliebige Vektornorm auf Kn .
Dann heißt
kAxk
kAk := sup = sup kAxk
x∈K \{0} kxk
n x∈Kn ,kxk=1
die k · k zugeordnete (oder natürliche) Matrixnorm. Sie ist verträglich und submultiplikativ.
Beweis: Tafel
53
4 Kondition linearer Gleichungssysteme
Hilfssatz 4.7. Die zugeordneten Matrixnormen zu k · k∞ und k · k1 sind:
n
X
kAk∞ := max |aij | (Zeilensummennorm)
i=1,...,n
j=1
Xn
kAk1 := max |aij | (Spaltensummennorm)
j=1,...,n
i=1
Beweis: [Rannacher, S. 104], Fall k · k∞ : Tafel
4.4 Eigenwerte und Eigenvektoren
Sei A ∈ Kn×n , λ ∈ K, und e ∈ Kn , so dass Ae = λe, dann heißt λ Eigenwert von A und e
zugehöriger Eigenvektor. (λ, e) heißt auch Eigenpaar.
Die Gleichung (A − λI)e = 0 ist für e 6= 0 nur für
P (λ) = det(A − λI) = 0
erfüllbar. Das charakteristische Polynom P (λ) hat genau n Nullstellen in C (Vielfachheit mit-
gezählt). Zu jedem Eigenwert λ gibt es mindestens einen Eigenvektor e.
Für alle Matrixnormen k · k und Eigenwerte λ von A gilt
|λ| ≤ kAk. (4.2)
Also: alle Eigenwerte λ ∈ C liegen innerhalb eines Kreises um Null mit Radius kAk. Mit kAk∞
erhält man z.B. eine konkrete Abschätzung für den betragsmäßig größten Eigenwert.
Folgerung aus (4.2):
kAk ≥ ρ(A) := max {|λ| | λ ist Eigenwert von A}
Die Zahl ρ(A) heißt Spektralradius von A (leider ist das keine Norm).
Beweis: Tafel
54
4.4 Eigenwerte und Eigenvektoren
Beispiele
Beispiel 4.8. • Vandermonde-Matrix für die Stützstellen 2, 3, 4:
1 2 4 (
det(A) = (3 − 2) · (4 − 2) · (4 − 3) = 2
A = 1 3 9
1 4 16 λ1 ≈ 18.67, λ2 ≈ 1.25, λ3 ≈ 0.0859
ρ(A) ≈ 18.67, kAkFr ≈ 19.62, kAk∞ = 21, kAk1 = 29
• Drehstreckung um / mit (θ = 45, r = 21/2 ):
(
cos( π4 ) − sin( π4 )
1/2 1 −1 λ1,2 = 1 ± i
A=2 · =
sin( π4 ) cos( π4 ) 1 1 e1,2 = (±i, 1)T
ρ(A) = 21/2 , kAkFr = kAk∞ = kAk1 = 2
• Nilpotente Matrix:
0 1
A= λ = 0 (geometrische Vielfachheit 1)
0 0
ρ(A) = 0, kAkFr = kAk∞ = kAk1 = 1
– A 6= 0, ρ(A) =⇒ ρ ist nicht definit
– ρ(A + AT ) = 1 > 0 = ρ(A) + ρ(AT ) =⇒ ρ ist nicht subadditiv
Idee: nutze Spektralradius als eine Art “Grenz-Maß” für Matrizen
Da ρ wie gezeigt keine Norm ist, funktioniert das leider nicht
Gerschgorin-Kreise
Man kann aus den Einträgen einer Matrix noch genauere Abschätzungen der Eigenwerte erhal-
ten:
Definition 4.9 (Gerschgorin-Kreise). Für A ∈ Cn×n heißen die Mengen
X
Si := z ∈ C | |z − aii | ≤ |aij |
j6=i
Gerschgorin-Kreise von A. Jeder komplexe Eigenwert λ von A liegt in mindestens einem Gerschgorin-
Kreis. (Analog kann man Kreise für die Spalten definieren.)
Beweis: Tafel
55
4 Kondition linearer Gleichungssysteme
Beispiel 4.10.
5 1 0
A = 1 10 1 λ1 ≈ 10.31, λ2 ≈ 4.82, λ3 ≈ 1.87
0 1 2
ρ(A) ≈ 10.31, kAkFr ≈ 11.53, kAk∞ = kAk1 = 12
Vorgehensweise: alle Eigenwerte liegen. . .
• Normen: . . . im Kreis um Null mit Radius 11.53
• Kreise: . . . in den Gerschgorin-Kreisen um 2, 5 und 10 (mit den Radien 1, 1 und 2)
• Symmetrie: . . . auf der reellen Achse
=⇒ λ1,2,3 ∈ [1, 3] ∪ [4, 6] ∪ [8, 11.53]
Spezielle Matrizen
Definition 4.11 (Spezielle Matrizen). Zu A ∈ Kn×m heißt AT mit (AT )ij = (A)ji die trans-
ponierte Matrix und A mit (A)ij = (A)ij die komplex konjugierte Matrix.
T
A ∈ Kn×n heißt hermitesch, falls A = A , d.h. aij = aji .
Reelle hermitesche Matrizen heißen symmetrisch (dann ist A = AT ).
T T
Eine komplexw. Matrix A ∈ Cn×n heißt normal, falls AA = A A.
T T T
Sie heißt unitär, falls zusätzlich AA = A A = I, also A−1 = A .
Reelle unitäre Matrizen heißen orthogonal (dann ist A−1 = AT ).
Skalarprodukte
Definition 4.12 (Skalarprodukte). Eine Abbildung (·, ·) : Kn × Kn → K heißt Skalarprodukt,
falls gilt:
(S1) Symmetrie: (x, y) = (y, x)
(S2) Linearität: (αx + βy, z) = α(x, z) + β(y, z), x, y, z ∈ Kn , α, β ∈ K
(S3) Definitheit: (x, x) > 0 für x ∈ Kn \ {0}
• (S3) nutzt, dass aus (S1) (x, x) ∈ R folgt
• Aus (S1) und (S2) folgt (z, αx + βy) = α(z, x) + β(z, y), also sind Skalarprodukte Biline-
arformen für K = R und Sesquilinearformen für K = C.
56
4.5 Die Spektralnorm
Ein Skalarprodukt erzeugt immer eine zugehörige Norm
p
kxk := (x, x), x ∈ Kn .
Es gilt die Cauchy-Schwarzsche Ungleichung
|(x, y)| ≤ kxk · kyk, x, y ∈ Kn .
Das Euklidische Skalarprodukt lautet
n
X
(x, y)2 = xi y i = xT y
i=1
(wird im folgenden oft verwendet). Damit gilt:
T
A=A ⇐⇒ (Ax, y)2 = (x, Ay)2 ∀x, y ∈ Kn
(andere Schreibweise für hermitesch, heißt auch selbstadjungiert).
4.5 Die Spektralnorm
Die der Euklidischen Norm zugeordnete Matrixnorm kAk2 heißt Spektralnorm.
Hilfssatz 4.13. Für die Spektralnorm hermitescher Matrizen gilt
kAk2 = max {|λ| | λ ist Eigenwert von A} = ρ(A)
Für jede Matrix A ∈ Kn×n gilt
n o
T
kAk2 = max |λ|1/2 | λ ist Eigenwert von A A
Beweis: Tafel
57
4 Kondition linearer Gleichungssysteme
4.6 Positiv definite Matrizen
Definition 4.14 (Positiv definite Matrizen). Eine Matrix A ∈ Kn×n heißt positiv definit (in
K), falls die beiden folgenden Eigenschaften erfüllt sind:
• (Ax, x)2 ∈ R ∀x ∈ Kn (trivial falls K = R)
• (Ax, x)2 > 0 ∀x ∈ Kn \ {0} (eigentliche Bedingung)
Das definiert eine wichtige Klasse von Matrizen mit vorteilhaften Eigenschaften. Ziel: Charak-
terisierung positiv definiter Matrizen.
Eigenschaft 4.15. Für A ∈ Cn×n gilt (Ax, x)2 ∈ R ∀x ∈ Cn genau dann, wenn A hermitesch
ist.
Beweis: Übung
Positiv definite Matrizen in Cn×n sind also immer hermitesch, aber: positiv definite Matrizen
in Rn×n sind nicht notwendigerweise symmetrisch.
Lemma 4.16. Eine hermitesche Matrix A ist genau dann positiv definit, wenn alle ihre (reellen)
Eigenwerte positiv sind. Alle Hauptdiagonalelemente sind (reell und) positiv.
Beweis: Tafel
Lemma 4.17. Sei A ∈ Rn×n symmetrisch und positiv definit. Dann liegt der betragsmäßig
größte Eintrag auf der Hauptdiagonalen.
Normalerweise beschränkt man sich auf symmetrisch positiv definite Matrizen (kurz auch
s.p.d.). Manche Autoren verlangen das daher direkt in der Definition. Wir fordern das aber
extra.
Beweis: Tafel
58
4.7 Störungstheorie
4.7 Störungstheorie
Seien A ∈ Kn×n regulär und x, b ∈ Kn . Dann ist F (A, b) = A−1 b “Lösungsoperator” der
Gleichung G(x) = Ax − b = 0.
Betrachte die relative Kondition von F , zunächst nur bzgl. Änderungen in b:
kF (A, b + δb) − F (A, b)k kδbk
≤ kAk · kA−1 k ·
kF (A, b)k kbk
Definition 4.18 (Konditionszahl). Die Zahl cond(A) = kAk · kA−1 k heißt Konditionszahl von
A (bzgl. einer verträglichen Matrixnorm k · k).
Als nächstes wollen wir Störungen in A selbst zulassen. Dabei stellt sich zunächst die Frage, für
welche Störungen überhaupt mit einer (eindeutigen) Lösung zu rechnen ist. Wenn A regulär
ist, wann ist dann A + δA regulär?
Hilfssatz 4.19. Sei B ∈ Kn×n mit Matrixnorm kBk < 1. Dann ist I + B regulär, und es gilt
k(I + B)−1 k ≤ c · (1 − kBk)−1
mit c > 0. Falls die Matrixnorm zugeordnet ist, kann c = 1 gewählt werden.
Beweis: Tafel
Störungssatz
−1
Satz 4.20 (Störungssatz). Sei A ∈ Kn×n regulär und kδAk < kA−1 k (für zugeordnete Ma-
trixnorm). Dann ist A
e = A + δA ebenfalls regulär, und für den relativen Fehler des gestörten
Systems
(A + δA) · (x + δx) = b + δb (oder kurz: Ae
ex = eb)
gilt:
kδxk cond(A) kδbk kδAk
≤ kδAk
· +
kxk 1 − cond(A) · kbk kAk
kAk
Beweis: Tafel
59
4 Kondition linearer Gleichungssysteme
Beispiele
kδAk kδbk
Beispiel 4.21. Es sei kAk ≈ 10−k und kbk ≈ 10−k , sowie cond(A) ≈ 10s (mit k, s > 0).
Weiter nehmen wir an, dass 10s · 10−k 1, also z.B. s − k ≤ −3.
Dann gilt
kδxk 10s
≤ · 2 · 10−k ≈ 10−k+s
kxk 1 − 10s · 10−k
Man verliert also bis zu s Stellen an Genauigkeit!
Beispiel 4.22 (Kondition und Determinante). • Betrachte Matrix B = 10−10 · I:
−10 10
10 0 −1 10 0
B= ,B =
0 10−10 0 1010
Es gilt cond∞ (B) = kBk∞ · kB −1 k∞ = 1, aber det(B) = 10−20 !
• Generell gilt cond(A) = cond(A) für 6= 0, aber det(A) = n det(A)
=⇒ Die Determinante einer Matrix sagt nichts über die Fehlerverstärkung beim Lösen des
zugehörigen Gleichungssystems aus
d1 0
Beispiel 4.23. Für A = ergibt sich cond∞ (A) = max(d1 , d2 )/ min(d1 , d2 ), also z.B.
0 d2
1010 für die Wahl d1 := 1010 , d2 := 1. Andererseits sind alle Gleichungen in A unabhängig und
xi = bi /di jeweils gut konditioniert. Was passiert hier?
• Bei Fehlern nur in b: das Problem liegt in der Verwendung von Normen, es handelt sich
nur um eine Abschätzung, und diese ist für die einzelnen Komponenten hier sehr schlecht.
• Bei Fehlern auch in A: kleine Schwankungen in den Nebendiagonaleinträgen führen dazu,
dass das gestörte System voll gekoppelt ist. Der Störungssatz behandelt die Auswirkungen.
60
5 Eliminationsverfahren zur Lösung linearer Gleichungssysteme
5.1 Dreieckssysteme (gestaffelte Gleichungssysteme)
Sei A ∈ Kn×n von oberer Dreiecksgestalt:
a11 · x1 + a12 · x2 + · · · + a1n · xn = b1
a22 · x2 + · · · + a2n · xn = b2
.. .. ..
. . .
ann · xn = bn
Dieses System ist genau dann eindeutig lösbar, wenn aii 6= 0, i = 1, . . . , n, gilt.
Aufgrund der einfachen Struktur lässt sich das System dann durch “Rückwärtseinsetzen” lö-
sen.
Lösung durch rückwärts einsetzen:
xn = bn /ann
xn−1 = bn−1 − a(n−1)n · xn /a(n−1)(n−1)
..
.
n
!
X
xi = bi − aik · xk /aii
k=i+1
Die Anzahl der benötigten Rechenoperationen ist:
n−1
X
N∆ (n) = (2i + 1) = n2
i=0
(Gilt für untere Dreieckssysteme natürlich analog.)
5.2 Gauß-Elimination
Sei nun A ∈ Kn×n regulär, aber mit beliebiger Struktur.
Ziel: Bringe A durch äquivalente Umformungen auf (obere) Dreiecksgestalt, um die Lösung
leicht berechnen zu können.
Dazu benutzt man:
61
5 Eliminationsverfahren zur Lösung linearer Gleichungssysteme
• Vertauschen zweier Gleichungen
• Addition eines Vielfachen einer Gleichung zu einer anderen
(Wir beschränken uns der Einfachheit halber im folgenden auf den Fall K = R, für die komple-
xen Zahlen funktioniert es aber analog.)
Wir schreiben
a11 · · · a1n b1
.. .. ..
[A, b] = . . .
an1 · · · ann bn
(erweiterte Koeffizientenmatrix)
Das Gaußsche Eliminationsverfahren erzeugt dann eine endliche Folge von Matrizen
h i h i h i
[A, b] = A(0) , b(0) , A(1) , b(1) , . . . , A(n−1) , b(n−1) ,
so dass nach n − 1 Schritten A(n−1) obere Dreiecksgestalt hat.
Erster Schritt: A(0) , b(0)
(1) (1)
A ,b
h i
Erster Schritt, Teil (a): A(0) , b(0) e(0) , eb(0)
A
(0)
• Bestimme ein r ∈ {1, . . . , n} mit ar1 6= 0
(0)
• Dieses Pivotelement ar1 existiert, weil A regulär ist
(0)
• Vertausche Zeilen 1 und r, dadurch e
a11 6= 0
Beispiel:
0 4 −2 −3 1 1 3 0 −2 2
−2 −6 3 1 2 −2 −6 3 1 2
3 7 4 −1 3 3 7 4 −1 3
1 3 0 −2 2 0 4 −2 −3 1
h i
Erster Schritt, Teil (b): Ae(0) , eb(0)
(1) (1)
A ,b
(0) (0)
• Berechne für i ∈ {2, . . . , n} den Faktor qi1 = e
ai1 /e
a11
• Ziehe das qi1 -fache der Zeile 1 von Zeile i ab
• Dadurch werden alle Einträge der ersten Spalte Null, bis auf den ersten (das Pivotelement)
62
5.2 Gauß-Elimination
(1) (0) (0)
∀ j ∈ {1, . . . , n} : aij = e
aij − qi1 e
a1j
(1) (0) (0)
bi = ebi − qi1eb1
1 3 0 −2 2 1 3 0 −2 2
−2 −6 3 1 2 0 0 3 −3 6
3 7 4 −1 3 0 −2 4 5 −3
0 4 −2 −3 1 0 4 −2 −3 1
k-ter Schritt: A(k−1) , b(k−1) A(k) , b(k)
h i
k-ter Schritt, Teil (a): A(k−1) , b(k−1) e(k−1) , eb(k−1)
A
• Die ersten k − 1 Zeilen und Spalten sind die einer oberen Dreiecksmatrix
(k−1)
• Finde ein neues r ∈ {k, . . . , n} mit Pivotelement ark 6= 0
• Tauschen der k-ten und r-ten Zeilen ändert die ersten k − 1 Zeilen und Spalten nicht
Beispiel für k = 2:
1 3 0 2 −2 1 3 0 −2 2
0 0
3 6 −3 0 −2
4 5 −3
0 −2 4 5 −3 0 0 3 −3 6
0 4 −2 −3 2 0 4 −2 −3 2
h i
k-ter Schritt, Teil (b): Ae(k−1) , eb(k−1) A(k) , b(k)
(k−1) (k−1)
• Berechne für i ∈ {k + 1, . . . , n} den Faktor qik = e
aik /e
akk
• Ziehe das qik -fache der Zeile k von Zeile i ab
• Diese Subtraktion ändert die k −1 ersten Zeilen und Spalten nicht (für die Spalten: wegen
(k−1)
akj
e = 0 für j < k)
(1) (0) (0)
∀ j ∈ {k, . . . , n} : aij = e
aij − qik e
akj
(1) (0) (0)
bi = ebi − qikeb1
1 3 0 −2 2 1 3 0 −2 2
0 −2 4 5 −3 0 −2 4 5 −3
0 0 3 −3 6 0 0 3 −3 6
0 4 −2 −3 2 0 0 6 7 −4
63
5 Eliminationsverfahren zur Lösung linearer Gleichungssysteme
Aufwand der Gauß-Elimination
• In jedem Schritt entsteht eine Zeile (und Spalte), die die korrekte Form hat und danach
nicht mehr verändert wird
=⇒ das Verfahren bricht nach einer endlichen Zahl Schritte ab
• Da es immer mindestens ein geeignetes Pivotelement gibt, lassen sich die einzelnen Teil-
schritte stets durchführen
=⇒ das Verfahren besteht aus einzelnen, berechenbaren Operationen
Formaler Algorithmus
Wie aufwändig ist die Gauß-Elimination?
Lemma 5.1. Der Aufwand zur Transformation von A auf obere Dreiecksgestalt durch Gauß-
Elimination beträgt
2
NGauß (n) = n3 + O(n2 )
3
Da der Aufwand des Rückwärtseinsetzens vernachlässigbar ist (N∆ (n) = n2 ), ist dies auch der
Aufwand für das Lösen eines linearen Gleichungssystems mit dem Gauß-Verfahren.
Beweis: Tafel
Vollständiger Algorithmus
Algorithmus 5.2 (Gauß-Elimination). In Pseudocode:
Eingabe : A ∈ Rn×n (wird überschrieben)
b ∈ Rn (wird überschrieben)
Ausgabe : x ∈ Rn
for ( k = 1; k < n; k = k+1 )
{
// Schritt (a)
// Finde ein geeignetes Pivotelement
Finde r in {k,...,n}
mit a_rk != 0 (sonst Fehler);
if (r != k)
{ // Tausche Zeilen r und k
for ( j = k; j <= n; j = j+1 )
64
5.2 Gauß-Elimination
{t = a_kj; a_kj = a_rj; a_rj = t;}
t = b_k; b_k = b_r; b_r = t;
}
// (... -->)
// (...)
// Schritt (b)
// Subtrahiere Vielfache der Zeile k
for ( i = k+1; i <= n; i = i+1 )
{
q_ik = a_ik / a_kk;
for (j = k; j <= n; j = j+1 )
a_ij = a_ij - q_ik * a_kj;
b_i = b_i - q_ik * b_k;
}
}
// A ist jetzt obere Dreiecksmatrix,
// löse daher durch Rückwärtseinsetzen
Vorsicht: In der Numerik zählen wir üblicherweise ab Eins, aber in C++ ist Null der kleinste
Index!
Beispiel
Beispiel 5.3. Ein weiteres Beispiel, diesmal ohne Vertauschungen:
2 4 6 8 40 2 4 6 8 40
16 33 50 67 330 0 1 2 3 10
4 15 31 44 167 0 7 19 28 87
10 29 63 97 350 0 9 33 57 150
2 4 6 8 40 2 4 6 8 40
0 1 2 3 10 0 1 2 3 10
0 0 5 7 17 0 0 5 7 17
0 0 15 30 60 0 0 0 9 9
Lösung durch Rückwärtseinsetzen: (4, 3, 2, 1)T
65
5 Eliminationsverfahren zur Lösung linearer Gleichungssysteme
5.3 LR-Zerlegung
Transpositionsmatrizen
Ziel: Finde möglichst kompakte Darstellung der Gauß-Elimination
abstrakte Formulierung mit Matrixmultiplikationen
Definition 5.4 (Transpositionsmatrizen). Eine Matrix Prs ∈ Rn×n , 1 ≤ r, s ≤ n, r 6= s, mit
1
i = j ∧ i 6= r ∧ i 6= s
[Prs ]ij = 1 (i = r ∧ j = s) ∨ (i = s ∧ j = r)
0 sonst
heißt Transpositionsmatrix.
Beispiel:
1 0 0 0
0 0 0 1
P24 =
0
0 1 0
0 1 0 0
Hilfssatz 5.5. Eine Transpositionsmatrix Prs hat die folgenden Eigenschaften:
1. A
e = Prs A ist A mit Zeilen r und s vertauscht
2. A
e = APrs ist A mit Spalten r und s vertauscht
T (symmetrisch)
3. Prs = Prs
−1 = P
4. Prs −1 T
rs =⇒ Prs = Prs (orthogonal)
Endliche Produkte von Transpositionsmatrizen heißen Permutationsmatrizen. Diese sind or-
thogonal, aber nicht symmetrisch.
Beweis: Tafel
66
5.3 LR-Zerlegung
Frobeniusmatrizen
Die Transpositionsmatrizen werden für Teil (a) der Eliminationsschritte verwendet, für Teil (b)
benötigen wir eine spezielle Form von Matrix.
Definition 5.6 (Frobeniusmatrizen). Eine Matrix Gk ∈ Rn×n heißt Frobeniusmatrix, wenn
Gk = I + G0k , wobei G0k ∈ Rn×n mit [G0k ]ij = 0 für j 6= k ∨ i ≤ k.
Beispiel:
1 0 0 0
0 1 0 0
G2 =
0
3 1 0
0 4 0 1
Hilfssatz 5.7. Eine Frobeniusmatrix Gk hat die folgenden Eigenschaften:
(
aij i≤k
1. Für Ae = Gk A gilt e
aij =
aij + gik akj i > k
(Das ist gerade der Eliminationsschritt, wenn gik := −qik )
2. G0k G0k = 0
3. G−1 0
k = (I − Gk )
Pk
und G−1 −1 −1
0
Pk 0
4. G1 G2 · · · Gk = I + j=1 Gj 1 G2 · · · Gk = I − j=1 Gj
5. ∀ k < r ≤ s : Prs Gk = G e k := I + Prs G0
e k Prs mit G
k
(Ändern der Reihenfolge führt zu neuer Frobeniusmatrix G
ek )
Beweis: Tafel
LR-Zerlegung
Für das Produkt von Matrizen Bi , a ≤ i ≤ b, definieren wir:
b
Y
Bi := Bb Bb−1 · · · Ba+1 Ba
i=a
(Achtung: von rechts nach links, da Multiplikation nicht kommutativ)
Für Transpositionsmatrizen ist dieses Produkt eine Permutationsmatrix, für Frobeniusmatrizen
nach dem eben gezeigten eine untere Dreiecksmatrix.
67
5 Eliminationsverfahren zur Lösung linearer Gleichungssysteme
Satz 5.8 (LR-Zerlegung). Sei A ∈ Rn×n regulär, dann gibt es eine Zerlegung P A = LR, wobei
1 0 ··· 0 r11 r12 · · · r1n
.. .. ..
l11 1 . . , R = 0 r22
. r2n
L= .. .. .. .. . . ..
,
. . . 0 . . . r(n−1)n
ln1 · · · ln(n−1) 1 0 ··· 0 rnn
n−1
Q
und P = Pα,sα , sα ≥ α, eine Permutationsmatrix ist. Im Fall P = I ist die Zerlegung
α=1
eindeutig.
Beweis: Tafel
Aus dem Beweis ergibt sich, dass
1 n−1
! !
Y Y
L= I− Pα,sα G0k
k=n−1 α=k+1
n−1 n−1
!
X Y
=I− Pα,sα G0k mit [G0k ]ik = −qik
k=1 α=k+1
• Die Nichtnullelemente −qik in G0k nehmen den selben Platz wie die eliminierten Einträge
(k−1)
aik
e ein
• Kein zusätzlicher Speicher notwendig, verwende einfach die linke untere Hälfte von A
(wird sowieso überschrieben)
• Die Einträge von G0k nehmen ebenfalls an späteren Zeilenvertauschungen teil (tausche
daher einfach ganze Zeile)
Wozu das Ganze?
Lösen eines Gleichungssystems Ax = b mittels LR-Zerlegung:
Ax = b
⇐⇒ P Ax = P b
⇐⇒ LRx = P b, y := Rx
(
Ly = P b (vorwärts einsetzen)
⇐⇒
Rx = y (rückwärts einsetzen)
Bei bekannter LR-Zerlegung lässt sich das Gleichungssystem also mit zwei sehr einfachen Schrit-
ten lösen.
Einzelne Schritte:
68
5.3 LR-Zerlegung
1. Berechne LR-Zerlegung und Matrix P für gegebenes A
2. berechne b0 = P b für gegebenes b
3. Löse Dreieckssystem Ly = b0
4. Löse Dreieckssystem Rx = y
Die LR-Zerlegung ist äquivalent zur Gauß-Elimination und hat daher den selben Aufwand
NLR (n) = 32 n3 + O(n2 ).
Wichtiger Unterschied:
Für eine neue rechte Seite Ae
x = eb kann die LR-Zerlegung wiederverwendet werden, die Gauß-
Elimination müsste erneut durchgeführt werden!
Algorithmus
Algorithmus 5.9 (LR-Zerlegung). In Pseudocode:
Eingabe : A ∈ Rn×n (wird überschrieben)
Ausgabe : L, R ∈ Rn×n (ohne Diagonale von L)
p : {1, . . . , n} → {1, . . . , n}
for ( k = 1; k < n; k = k+1 )
{
// Schritt (a)
// Finde ein geeignetes Pivotelement
Finde r in {k,...,n}
mit a_rk != 0 (sonst Fehler);
if (r != k)
{ // Tausche Zeilen r und k
// Tausche jetzt auch alte Spalten!
for ( j = 1; j <= n; j = j+1 )
{t = a_kj; a_kj = a_rj; a_rj = t;}
}
p_k = r; // merke Permutation
// (... -->)
// (...)
// Schritt (b)
// Subtrahiere Vielfache der Zeile k
for ( i = k+1; i <= n; i = i+1 )
{ // Verwende jetzt a_ik statt q_ik
69
5 Eliminationsverfahren zur Lösung linearer Gleichungssysteme
a_ik = a_ik / a_kk;
// Subtrahiere weiterhin ab Spalte k
// (sonst würden q_ij verändert!)
for (j = k; j <= n; j = j+1 )
a_ij = a_ij - a_ik * a_kj;
}
}
(Wir speichern nicht die kompletten Transpositionsmatrizen, sondern nur den Index, mit dem
jeweils getauscht werden muss. Die Diagonale von L ist konstant Eins und wird daher nicht
gespeichert.)
Beispiel
Beispiel 5.10 (LR-Zerlegung). Betrachte wieder das Beispiel bei der Gauß-Elimination.
Einträge von R sind rot markiert, Einträge von L blau, Permutationsindizes grün.
Tausche erste und vierte Zeile, erste Zeile ist damit fertig.
Notiere, dass die erste Zeile mit der vierten getauscht wurde:
0 4 −2 −3 1 3 0 −2
−2 −6 3 1 −2 −6 3 1
3 7 4 −1 3 7 4 −1
1 3 0 −2 0 4 −2 −3
(0, 0, 0)T (4, 0, 0)T
Eliminiere die verbleibenden Einträge der ersten Spalte, trage stattdessen die Faktoren qi1 ein:
1 3 0 −2 1 3 0 −2
−2 −6 3 1 −2 0 3 −3
3 7 4 −1 3 −2 4 5
0 4 −2 −3 0 4 −2 −3
(da das Pivotelement Eins ist, bleiben hier die Einträge einfach erhalten!)
Tausche die zweite Zeile mit der dritten:
1 3 0 −2 1 3 0 −2
−2 0 3 −3 3 −2 4 5
3 −2 4 5 −2 0 3 −3
0 4 −2 −3 0 4 −2 −3
(4, 0, 0)T (4, 3, 0)T
70
5.4 Rundungsfehleranalyse der LR-Zerlegung
(wichtig: dabei werden die Einträge von L mitgetauscht!)
Eliminiere die verbleibenden Einträge der zweiten Spalte und merke die Faktoren qi2 :
1 3 0 −2 1 3 0 −2
3
−2 4 5
3
−2 4 5
−2 0 3 −3 −2 0 3 −3
0 4 −2 −3 0 −2 6 7
Beobachtung: jeder Tausch-Schritt ergibt eine Zeile von R, jeder Eliminationsschritt eine Spalte
von L
Im letzten Schritt ist kein Tauschen nötig. Notiere dies durch den Index der aktuellen Zeile.
Eliminiere den letzen Eintrag und speichere dort stattdessen q43 :
1 3 0 −2 1 3 0 −2
3 −2 4 5 3 −2 4 5
−2 0 3 −3 −2 0 3 −3
0 −2 6 7 0 −2 2 13
(4, 3, 0)T (4, 3, 3)T
Prüfe nach, dass L · R = P · A:
1 0 0 0 1 3 0 −2
· 0 −2 4 5
3 1 0 0
−2 0 1 0 0 0 3 −3
0 −2 2 1 0 0 0 13
0 0 0 1 0 4 −2 −3
· −2 −6 3
0 0 1 0 1
=
0 1 0 0 3 7 4 −1
1 0 0 0 1 3 0 −2
5.4 Rundungsfehleranalyse der LR-Zerlegung
Absolutwertnotation
Sei A ∈ Rm×n , dann definieren wir die Matrix der Beträge B = |A| durch
bij := |aij |, 1 ≤ i ≤ m, 1 ≤ j ≤ n,
und analog |x| für x ∈ Rn (Vektor der Beträge).
Damit können wir die Abbildung rd auf Rm×n bzw. Rn erweitern:
rd(A) = A + A0 mit |A0 | ≤ |A| · eps
71
5 Eliminationsverfahren zur Lösung linearer Gleichungssysteme
(das sind m · n Ungleichungen), wegen
rd(aij ) = aij (1 + ij ) = aij + ij aij , |ij aij | ≤ |aij | · eps
fl-Notation
Sei E eine Formel, dann bezeichne fl(E) eine Berechnung der Formel in Gleitkommaarithmetik.
Die Reihenfolge wird dabei meist offensichtlich sein (von links nach rechts, Klammerung), sonst
muss sie angegeben werden.
Somit ist zum Beispiel:
fl(A + B) = (A + B) + H mit |H| ≤ eps · (|A + B|)
wegen
fl(aij + bij ) = aij ⊕ bij = (aij + bij )(1 + ij ) = (aij + bij ) + ij (aij + bij )
√
Man kann z.B. auch fl( x) und fl(sin(x)) schreiben, letzteres setzt Angabe von Algorithmus
voraus.
Rückwärtsanalyse
Bisher haben wir die “Vorwärtsanalyse” der Rundungsfehler betrieben, z.B. gilt für das Skalar-
produkt xT y:
| fl(xT y) − xT y| ≤ n · eps · |x|T |y| + O(eps2 ) (absolut)
oder
| fl(xT y) − xT y| |x|T |y|
≤ n · eps · + O(eps2 ) (relativ)
|xT y| |xT y|
Eine Alternative ist die sog. “Rückwärtsanalyse”. Dort versucht man, das Fließkommaergebnis
als exaktes Ergebnis eines modifizierten Ausdrucks zu schreiben, z.B.:
fl(xT y) = (x + f )T y mit |f | ≤ n · eps · |x| + O(eps2 )
Erinnerung: Ein Algorithmus ist stabil, wenn der akkumulierte Fehler der Rechnung den auf-
grund der Problemstellung unvermeidbaren Fehler (→ Kondition) nicht übersteigt (oder zu-
mindest in der selben Größenordnung ist).
Der tatsächlich auftretende Diskretisierungsfehler hat dann die gleiche Größenordnung wie der
kleinste für diese Aufgabe theoretisch mögliche Fehler.
72
5.4 Rundungsfehleranalyse der LR-Zerlegung
F(x)
k·∆x
x ∆x
x' F(x')
F'(x')-F(x')
F'(x')
Definition 5.11 (Rückwärtsstabilität). Ein Algorithmus F 0 für eine Funktion F heißt rück-
e mit F 0 (x) = F (e
wärtsstabil, wenn zu jeder Eingabe x eine modifizierte Eingabe x x) existiert,
für die |x − x
e| = O(|x − rd(x)|) = O(|x| · eps) gilt.
Interpretation: Der Algorithmus löst das gestellte Problem, aber für eine leicht gestörte Einga-
be.
Man kann zeigen, dass rückwärtsstabile Algorithmen stets stabil sind.
F(x)
k·∆x
x ∆x
x' F(x')
x~
F'(x')-F(x')
F'(x')
~
=F(x)
Beispiel
73
5 Eliminationsverfahren zur Lösung linearer Gleichungssysteme
Beispiel 5.12. Die Berechnungen F 0 (x, y) = x ⊕ y und G0 (x) = x ⊕ 1 sind jeweils stabile
Algorithmen für die Addition bzw. das Inkrement. F 0 ist rückwärtsstabil:
x ⊕ y = (x + y) · (1 + ) = x · (1 + ) + y · (1 + ) = x
e + ye
mit xe := x · (1 + ), |x − xe| = |x| · ≤ |x| · eps, ye analog
Aber G0 ist nicht rückwärtsstabil:
x ⊕ 1 = (x + 1) · (1 + ) = x · (1 + (1 + x−1 ) · ) + 1 = x
e+1
mit e := x · (1 + (1 + x−1 ) · ),
x e| = |x| · |1 + x−1 | ·
|x − x
Für x 1 wird der relative Fehler beliebig groß (vgl. ⊕ 1 = 1).
Skalarprodukt
Hilfssatz 5.13. Das Skalarprodukt xT y ∈ R, x, y ∈ Fn , ist rückwärtsstabil:
ŝ = fl(xT y) = (x + f )T y mit |f | ≤ n · eps · |x| + O(eps2 )
Bemerkung 5.14. Der Faktor n · eps ist der schlechteste Fall und sehr pessimistisch. Run-
dungsfehler von Operationen hängen von den Argumenten ab und können sich auch wegheben
(Vorzeichenwechsel!). Eine statistische Betrachtung würde auf einen deutlich geringeren durch-
schnittlichen Fehler führen.
Beweis: Tafel
Dyadisches Produkt
Beispiel 5.15. Das dyadische Produkt (äußere Produkt) x · y T ∈ Rn×n hingegen ist zwar stabil,
aber nicht rückwärtsstabil. Warum?
Ein dyadisches Produkt hat immer Rang Eins, fl(x · y T ) wird jedoch durch die Störungen i.A.
einen anderen Rang haben (und kann somit nicht das Ergebnis einer exakten Rechnung sein).
Lösen von Dreieckssystemen
Satz 5.16 (Lösen von Dreieckssystemen). Es seien x̂ und ŷ die numerischen Lösungen des
unteren bzw. oberen Dreieckssystems Lx = b und Ry = c. Dann gilt
(L + F )x̂ = b |F | ≤ n · eps · |L| + O(eps2 )
(R + G)ŷ = c |G| ≤ n · eps · |R| + O(eps2 )
Das Lösen von Dreiecksystemen durch Vorwärts- / Rückwärtseinsetzen ist also rückwärtsstabil.
Beweis: Übung
74
5.4 Rundungsfehleranalyse der LR-Zerlegung
LR-Zerlegung
Satz 5.17 (Rückwärtsanalyse der LR-Zerlegung). Sei A ∈ Fn×n . Es werde die LR-Zerlegung
von A ohne Zeilentausch berechnet (sofern möglich). Dann gilt für die numerisch berechneten
Matrizen L
b und R:
b
b·R
L b =A+H mit |H| ≤ 3(n − 1) · eps · |A| + |L| b + O(eps2 )
b · |R|
Beweis: Tafel
Satz 5.18 (Lösen mittels LR-Zerlegung). Seien L b und R
b die numerisch berechnete LR-Zerlegung
von A ∈ F n×n n
aus Satz 5.17. Sei weiter ŷ ∈ F die numerische Lösung von Ly
b = b, und schließ-
n
lich x̂ ∈ F die numerische Lösung von Rx = ŷ. Dann gilt für x̂ die Beziehung
b
(A + E)x̂ = b
mit
|E| ≤ n · eps · (3|A| + 5|L| b + O(eps2 )
b · |R|)
Beweis: Tafel
Folgerung
Nach Satz 5.18 ist x̂ exakte Lösung des modifizierten Systems (A + E)x̂ = b. Mit dem Stö-
rungssatz gilt dann in k · k∞ -Norm
( )
kx̂ − xk∞ kLk
b
∞ · kRk∞
b
≤ cond∞ (A) · 3n · eps + 5n · eps + O(eps2 )
kxk∞ kAk∞
• Der erste Term ist mit dem aus der Konditionsanalyse vergleichbar und daher gutartig.
• Der zweite Term kann problematisch werden, da L
b Einträge der Form e aii enthält (und
aij /e
aii beliebig klein werden kann).
e
Man kann zeigen, dass R
b genauso problematisch sein kann.
=⇒ Die Gauß-Elimination / LR-Zerlegung ist in dieser Form nicht numerisch stabil!
75
5 Eliminationsverfahren zur Lösung linearer Gleichungssysteme
Beispiel
Beispiel 5.19. Betrachte die Matrix A gegeben durch
1 −1 0 1
A= , A =
1 0 1 −
mit 0 < 1. cond∞ (A) = (1 + )2 , also ist A gut konditioniert, und Rückwärtseinsetzen
bietet sich als stabiler Algorithmus an.
Die (exakte) LR-Zerlegung ist aber in diesem Fall A = LR mit
1 0 1
L= , R=
−1 1 0 −−1
−1
=⇒ kLk
b
∞ ≈ kLk∞ = + 1 1, kRk
b
∞≈ kRk ∞ = −1 1
5.5 Pivotisierung
Die Rundungsfehleranalyse in Satz 5.18 führt auf den unvorteilhaften Term |L| b und wie
b · |R|,
wir gesehen haben, können diese beiden Faktoren prinzipiell beliebig groß werden.
Wählt man im Gauß-Algorithmus (bzw. bei der LR-Zerlegung) den Index r so, dass
(k) (k)
|ark | ≥ |aik | ∀ k ≤ i ≤ n,
dann gilt
|ˆlij | ≤ 1 und damit kLk
b
∞ ≤n
Diese Wahl des Pivotelements nennt man Spaltenpivotisierung.
Beispiel
Beispiel 5.20. Das Gauß-Verfahren für das lineare Gleichungssystem
−10−5 1
x 1
· 1 =
2 1 x2 0
führt in exakter Arithmetik auf
−10−5
1 x1 1
· =
0 1 + 2 · 105 x2 2 · 105
mit der Lösung (x1 , x2 )T = (−0.4999975, 0.999995)T ≈ (−0.5, 1)T .
76
5.5 Pivotisierung
In F(10, 4, 1) ergibt sich zwar ebenfalls ein Faktor q21 = 2 · 105 = 0.2 · 106 , das dadurch entste-
hende System ist aber
−0.1 · 10−4
1 x1 1
· = ,
0 0.2 · 106 x2 0.2 · 106
(1)
der Eintrag a22 ist also (als einziger!) leicht fehlerbehaftet.
Als Ergebnis erhält man das völlig inakzeptable (x1 , x2 )T = (0, 1)T .
Die berechnete Lösung hat also einen relativen Fehler von 50% in der Maximumsnorm, obwohl
die Matrix sehr gut konditioniert ist.
Rückwärtsanalyse:
Tatsächlich wurde das modifizierte System
−10−5 1
x 1
· 1 = ,
2 0 x2 0
exakt gelöst, welches aber eine völlig andere Lösung hat.
Spaltenpivotisierung (in F(10, 4, 1)) führt in diesem Fall auf das System
2 1 x 0
· 1 = ,
0 1 x2 1
mit der vollkommen akzeptablen Antwort (x1 , x2 )T = (−0.5, 1)T .
Beachte aber:
Multiplikation der ersten Zeile des ursprünglichen Systems mit −106 führt auf ein mathematisch
äquivalentes System, das auch bei Spaltenpivotisierung keine Vertauschungen braucht und daher
trotzdem das ursprüngliche, unbrauchbare, Ergebnis liefert.
Äquilibrierung
Solche Probleme kann man oft durch eine Skalierung (Äquilibrierung) des Gleichungssystems
vermindern:
n
X
−1 −1
Ax = b ⇐⇒ D Ax = D b mit dii = |aij |
j=1
e = eb mit kAk
⇐⇒ Ax ∞ =1
e
77
5 Eliminationsverfahren zur Lösung linearer Gleichungssysteme
(Vorsicht: keine Aussage über kA
e−1 k , es kann durchaus auch cond∞ (A)
∞
e > cond∞ (A) sein!)
Spezialfall der allgemeineren Idee der Vorkonditionierung:
Transformiere die ursprüngliche Aufgabenstellung in eine neue Formulierung, die besser kon-
ditioniert ist und sich stabil lösen lässt, und transformiere danach die Ergebnisse zurück (falls
nötig)
Beispiel
Beispiel 5.21. Betrachte die folgende vom Parameter c 6= 0 abhängige Matrix Ac :
1 2c−1 1
c 0 2 −1 2c −c −1
Ac = · = , Ac = · −1
0 1 −1 2 −1 2 3 c 2
Es gilt cond∞ (Ac ) = 1 + 2 · max{|c|, |c|−1 }.
Optimal ist also die Wahl c = 1 mit cond∞ (A1 ) = 3. Für c 6= 1 wird eine der Zeilen mehr betont
als die andere, was im Extremfall zu einem sehr unausgeglichenen Gleichungssystem führt.
Rundungsfehleranalyse bei Pivotisierung
Analog zu Satz 5.18 zeigt man, dass für die Lösung x̂ bei Spaltenpivotisierung gilt:
(A + E)x̂ = b
mit
|E| ≤ n · eps · (3|A| + 5P T · |L| b + O(eps2 )
b · |R|)
Nach Konstruktion ist kLk
b
∞ ≤ n, und mit der Definition
(k)
ρ := kAk−1
∞ · max |âij |
i,j,k
(sog. “Wachstumsfaktor”) zeigt man
kEk∞ ≤ 8n3 kAk∞ ρ · eps + O(eps2 )
Beweis: [Golub, Van Loan: Matrix Computations]
Man kann Beispiele mit sehr schlechtem Verhalten für ρ konstruieren:
Betrachte A ∈ Rn×n mit
1
j =i∨j =n
aij := −1 j < i
sonst
0
78
5.5 Pivotisierung
Beispiel:
1 0 0 0 1
−1 1 0 0 1
−1 −1 1
A= 0 1
−1 −1 −1 1 1
−1 −1 −1 −1 1
Für diesen Extremfall ist ρ = 2n−1 .
In der Praxis wachsen die Einträge aber normalerweise nur um bis zu eine Größenordnung
(ρ ≈ 10, [Golub, Van Loan]), daher können Gauss-Verfahren und LR-Zerlegung mit Spaltenpi-
votisierung als stabil angesehen werden.
Totale Pivotisierung
Man kann auch r, s ∈ {k, . . . , n} wählen, so dass
(k)
|a(k)
rs | ≥ |aij | ∀ k ≤ i, j ≤ n
und erreicht dann durch Zeilen- und Spaltenvertauschung, dass
(k)
akk = a(k)
e rs
Diese Wahl des Pivotelements nennt man totale Pivotisierung.
Man kann zeigen, dass damit
(k) 1/2
|aij | ≤ k 1/2 (2 · 31/2 · · · k 1/(k−1) ) · max |aij |,
also deutlich kleineres Wachstum der Einträge.
In Matrixform ergibt sich (mit Prk := Pk,rk und Qsk := Pk,sk )
Schritt 1: G1 Pr1 AQs1 · Qs1 x = G1 Pr1 b
Schritt 2: G2 Pr2 G1 Pr1 AQs1 Qs2 · Qs2 Qs1 x = G2 Pr2 G1 Pr1 b
.. ..
. .
"n−1 1
# n−1 n−1
Y Y Y Y
Schritt n − 1: (Gk Prk )A Qsk · Qsk x = (Gk Prk )b
k=1 k=n−1 k=1 k=1
| {z }
=R (obere Dreiecksmatrix)
Analog zur Spaltenpivotisierung zeigt man, dass dann
P AQT · Qx = P b
79
5 Eliminationsverfahren zur Lösung linearer Gleichungssysteme
mit der LR-Zerlegung P AQT = LR (beachte QT · Q = I).
Setze dazu
n−1
Y n−1
Y
P := Prk , Q := Qsk ,
k=1 k=1
n−1
!−1 n−1
Y Y
L := P · (Gk Prk ) , R := (Gk Prk ) · A · QT
k=1 k=1
R ist obere Dreiecksmatrix nach Konstruktion, und wie zuvor ist L untere Dreiecksmatrix (selbe
Argumentation, da unabhängig von Q).
Außerdem kann L wieder unten links in A gespeichert werden (Spaltenvertauschungen ändern
die bisher berechneten Einträge von L nicht, da sk ≥ k).
Lösen bei totaler Pivotisierung
Mit diesen Definitionen haben wir
P AQT · Qx = LR · Qx = P b
Löse mit den folgenden Schritten nach x auf:
b0 = P b (Permutation)
0
Ly = b (Dreieckssystem)
0
Rx = y (Dreieckssystem)
−1 0
x=Q x (Permutation)
Im letzten Schritt muss man dabei Q nicht explizit invertieren. Stattdessen werden die Vertau-
schungen einfach in umgekehrter Reihenfolge durchgeführt.
Anmerkungen
Aufwand der Pivotisierung:
• n2 /2 Vergleiche bei Spaltenpivotisierung
• n3 /3 Vergleiche bei totaler Pivotisierung
=⇒ totale Pivotisierung hat gleiche Komplexität wie eigentliche Zerlegung, während Spalten-
pivotisierung vernachlässigbar ist
Praktische Erfahrung zeigt trotz höheren Aufwands keine Vorteile für Rundungsfehler bei to-
taler Pivotisierung
Spaltenpivotisierung mit Zeilenskalierung ist in der Praxis effizient und numerisch stabil.
80
5.6 Spezielle Systeme
5.6 Spezielle Systeme
Symmetrisch positiv definite Matrizen
Satz 5.22. Eine symmetrisch positiv definite Matrix A ∈ Rn×n ist stets ohne Pivotisierung sta-
bil LR-zerlegbar. Für die Diagonalelemente der im Eliminationsprozess auftauchenden Matrizen
gilt
(k)
aii ≥ λmin (A), k ≤ i ≤ n.
Die Symmetrie der Matrix kann man ausnutzen, um die Matrix A mit geringerem Aufwand zu
zerlegen.
Beweis: Tafel, Stabilität: [Golub, Van Loan]
Cholesky-Zerlegung
Mit D = diag(R) gilt
A = LD(D−1 R) = LDU mit U := D−1 R,
und wegen der Symmetrie gilt U = LT , also A = LDLT . Da alle Diagonalelemente von D
positiv sind, ist die Matrix D1/2 mit
1/2
(D1/2 )ii = dii , (D1/2 )ij = 0 für i 6= j
wohldefiniert, und es gilt
A = LD1/2 · D1/2 LT = L
eLeT mit L
e := LD1/2
Diese spezielle Form der Zerlegung heißt Cholesky-Zerlegung.
Ein einzelner Schritt der Cholesky-Zerlegung lässt sich über das äußere Produkt darstellen:
α vT β β −1 · v T
β 0 1 0
A= = −1 · ·
v B β · v In−1 0 B − α−1 · vv T 0 In−1
mit β := α1/2 . Dieser Prozess kann dann mit B−α−1 ·vv T rekursiv fortgesetzt werden. Die dabei
auftretenden Matrizen sind verallgemeinerte Frobenius-Matrizen (β als Diagonaleintrag).
Aufgrund der Symmetrie ist dabei der Aufwand halb so groß wie bei der allgemeinen LR-
Zerlegung:
1
NChol (n) = n3 + O(n2 )
3
81
5 Eliminationsverfahren zur Lösung linearer Gleichungssysteme
Diagonaldominante Matrizen
Definition 5.23 (Diagonaldominanz). Eine Matrix A ∈ Rn×n heißt diagonaldominant, falls
X
|aij | ≤ |aii | i = 1, . . . , n
j6=i
und strikt diagonaldominant, falls die Summe echt kleiner ist.
(vgl. mit Kriterium der Gerschgorin-Kreise, was folgt dann für A?)
Satz 5.24. Diagonaldominante reguläre Matrizen erlauben eine LR-Zerlegung ohne Pivotisie-
rung.
Beweis: Tafel
Rangbestimmung und Determinante
Für eine allgemeine, nicht unbedingt reguläre Matrix A ∈ Rn×n lässt sich die LR-Zerlegung
ebenfalls durchführen, sie bricht aber evtl. mangels Pivotelement frühzeitig ab.
• Wenn dies bei totaler Pivotisierung nach k Schritten passiert, dann hat A Rang k.
• Da die Einträge aufgrund von Fehlern nie exakt Null werden, muss mit einer Toleranz
gearbeitet werden.
Falls A Vollrang haben sollte, kann man wegen
det(A) = det(P T LR) = det(P T ) · det(L) · det(R) = ± det(R)
über die LR-Zerlegung die Determinante von A berechnen
(das Vorzeichen hängt davon ab, ob P eine gerade oder ungerade Permutation repräsentiert).
82
5.6 Spezielle Systeme
Inversenberechnung
Die Inverse einer Matrix A ∈ Rn×n kann wie folgt berechnet werden:
• Berechne LR-Zerlegung (Aufwand 23 n3 + O(n2 ))
• Löse Axi = ei , i = 1, . . . , n, für die kanonische Basis ei
(Aufwand 2n · n2 : n Gleichungen, jeweils 2 Dreieckssysteme)
• Forme Matrix der Urbilder: A−1 = [x1 , . . . , xn ] (spaltenweise)
Der Gesamtaufwand einer solchen Matrixinversion ist daher
8
NInv (n) = n3 + O(n2 )
3
Die Matrixinversion spielt in der Praxis keine Rolle — stattdessen beschränkt man sich auf
Matrix-Vektor-Produkte mit A−1 !
Tridiagonalsysteme
Definition 5.25 (Tridiagonalmatrix). Eine Matrix A ∈ Rn×n heißt Tridiagonalmatrix, falls
aij = 0 für |i − j| > 1, 1 ≤ i ≤ n
Eine Tridiagonalmatrix ist Spezialfall einer Bandmatrix mit
aij = 0 für j < i − ml oder j > i + mr , 1 ≤ i ≤ n
Falls sich die LR-Zerlegung einer Tridiagonal- bzw. Bandmatrix ohne Pivotisierung durchführen
lässt, sind L und R ebenfalls tridiagonal bzw. Bandmatrizen der selben Bandbreite (ml , mr ).
Rechnung für Tridiagonalmatrizen: Tafel
83
5 Eliminationsverfahren zur Lösung linearer Gleichungssysteme
5.7 Nichtreguläre Systeme
Sei nun A ∈ Rm×n , b ∈ Rm , mit rang(A) beliebig. Dann hat Ax = b genau eine Lösung, gar
keine Lösung, oder unendlich viele.
Einige Grundbegriffe aus der linearen Algebra:
Bild: im(A) = {y ∈ Rm | ∃ x ∈ Rn : y = Ax}
Kern: ker(A) = {x ∈ Rn | Ax = 0}
Rang: rang(A) = dim(im(A)) = rang(AT ) = dim(im(AT ))
Orthogonales Komplement:
im(A)⊥ = {y ∈ Rm | ∀ y 0 ∈ im(A) : (y, y 0 )2 = 0}
Es gilt
im(A)⊥ = ker(AT )
und daher
dim(im(A)) + dim(ker(AT )) = m
dim(im(AT )) + dim(ker(A)) = n
Damit lässt sich der Lösungsbegriff für lineare Gleichungssysteme auf nichtreguläre Systeme
erweitern.
Least Squares
Satz 5.26 (Least Squares). Für allgemeines A ∈ Rm×n , b ∈ Rm , existiert ein x̄ ∈ Rn , so dass
kAx̄ − bk2 = minn kAx − bk2 .
x∈R
Diese Bedingung ist äquivalent dazu, dass x̄ Lösung der folgenden “Normalengleichung” ist:
AT Ax̄ = AT b, AT A ∈ Rn×n
Falls rang(A) = n (und damit zwingend m ≥ n), ist x̄ eindeutig, sonst hat jede weitere Lösung
die Form x̄ + y mit y ∈ ker(A).
84
5.7 Nichtreguläre Systeme
Beweis: Tafel
AT A ist symmetrisch positiv semidefinit, also könnte man Least Squares prinzipiell über die
Cholesky-Zerlegung lösen (mit symmetrischem Pivoting im nicht-definiten Fall).
Aber dieses Problem ist i.A. sehr schlecht konditioniert:
cond(AT A) ≈ cond(A)2 für m = n
Daher verwendet man üblicherweise spezielle Zerlegungen, die ähnlich zur Cholesky-Zerlegung
AT A = LLT = RT R
zur Verfügung stellen, aber ohne AT A explizit aufzustellen.
Householdermatrizen
Definition 5.27 (Householderreflektion). Sei v 6= 0 ∈ Rn gegeben, dann heißt die Matrix
vv T
Qv = I − 2
vT v
Householdermatrix oder Householderreflektion.
Dabei ist vv T das äußere Produkt, also eine Matrix mit Rang 1, und v T v das innere Produkt,
also eine Zahl. Beachte: jedes Vielfache αv, α 6= 0, von v führt auf die gleiche Matrix Qv
(Skaleninvarianz).
Hilfssatz 5.28. Eine Householdermatrix Qv hat die folgenden Eigenschaften:
• Qv ist orthogonal und symmetrisch, d.h. QTv = Q−1
v = Qv
• Es gilt für x ∈ Rn :
vT x
Qv x = x − 2βv, β= ∈R
vT v
• Speziell gilt x = αv =⇒ Qv x = −x
• (x, v)2 = 0 =⇒ Qv x = x
• Als Operator aufgefasst bewirkt Qv eine Reflektion an der Ebene mit Normale v
85
5 Eliminationsverfahren zur Lösung linearer Gleichungssysteme
Beweis: Tafel
Sei x 6= 0 ∈ Rn , und sei v := x ± kxk2 e1 , dann gilt
Qv x = ∓2kxk2 e1
Wenn x die erste Spalte einer Matrix A ist, dann eliminiert also eine Multiplikation von A mit
Qv alle Einträge ai1 , i > 1.
Das heißt, dass Householdermatrizen Qv genauso wie Frobeniusmatrizen Gk benutzt werden
können, um eine Matrix auf Dreiecksgestalt zu bekommen. Für die Rekursion beschränkt man
das obige Vorgehen auf die entsprechenden Unterräume.
QR-Zerlegung
Satz 5.29 (QR-Zerlegung). Zu jeder Matrix A ∈ Rm×n mit m ≥ n und rang(A) = n existiert
eine orthogonale Matrix Q ∈ Rm×m und eine obere Dreiecksmatrix R ∈ Rm×n , so dass
A = QR.
Diese Zerlegung kann durch n − 1 Multiplikationen mit Householdermatrizen berechnet werden.
Die ersten n Spalten von Q bilden eine orthonormale Basis von im(A).
Beweis: Tafel
Anmerkungen
Bemerkung 5.30 (Anmerkungen zur QR-Zerlegung). • Um Auslöschung zu vermeiden,
wird das Vorzeichen in v meist so gewählt, dass sich das Vorzeichen des Diagonalein-
trags nicht ändert. Falls dieser Eintrag Null ist, kann beliebig eines der beiden Vorzeichen
gewählt werden.
• Analog zu P bei der LR-Zerlegung wird Q nicht explizit aufgestellt, man skaliert die Nor-
malenvektoren vi jeweils so, dass der erste Eintrag Eins ist, und speichert sie dann in der
linken unteren Hälfte von A (so wie zuvor die Einträge von L).
• Für die Matrix Q−1 = QT existiert, wie wir gesehen haben, eine einfache und effizi-
ente Darstellung: sie ist das Produkt der gleichen Qv wie Q selbst, nur in umgekehrter
Reihenfolge.
• Die QR-Zerlegung benötigt den doppelten Aufwand der LR-Zerlegung:
4
NQR = n3 + O(n2 )
3
86
5.7 Nichtreguläre Systeme
• Man kann zeigen dass die QR-Zerlegung rückwärtsstabil ist, man erhält für die numerisch
berechneten Matrizen Q
b und R:
b
b·R
Q b =A+H mit kHk ≈ eps · kAk
Insbesondere treten hier keine problematischen Terme auf. Das gleiche gilt für das Lösen
mittels QR-Zerlegung.
Lösen nichtregulärer Systeme
1. Fall m > n (überbestimmt) und rang(A) = n:
Nach Satz 5.29 existiert QR-Zerlegung A = QR.
A = QR =⇒ AT A = (QR)T QR = RT QT QR = RT R
also
AT Ax̄ = AT b ⇐⇒ RT Rx̄ = RT QT b
Seien R
e und Qe die ersten n Zeilen von R, R = R , bzw. Spalten von Q, Q = [QQ
e ∗ ],
e
0
dann gilt (Rückwärtseinsetzen)
Rx̄
e =Qe T b ⇐⇒ AT Ax̄ = AT b
2. Fall m < n (unterbestimmt) und rang(A) = m:
Nach Satz 5.29 existiert QR-Zerlegung AT = QR.
AT = QR =⇒ AT A = QR(QR)T = QRRT QT
also
AT Ax = AT b ⇐⇒ QRRT QT x = QRb ⇐⇒ RRT QT x = Rb
Seien Re und Q
e die ersten m Zeilen von R bzw. Spalten von Q, dann gilt (Rückwärts-
einsetzen)
eT Q
R e T x̄ = b =⇒ AT Ax̄ = AT b
Die so bestimmte Lösung x̄ ist jene mit minimaler Norm kx̄k2 .
87
5 Eliminationsverfahren zur Lösung linearer Gleichungssysteme
3. Fall m = n (quadratisch) und rang(A) = n:
Nach Satz 5.29 existiert QR-Zerlegung A = QR. Da A Vollrang hat, gilt
AT Ax̄ = AT b ⇐⇒ Ax̄ = b ⇐⇒ QRx = b
Die Least-Squares-Lösung ist also die normale Lösung des Systems, und kann über die
QR-Zerlegung berechnet werden.
Überbestimmte, quadratische und unterbestimmte Systeme mit maximalem Rang haben
also immer eine eindeutig bestimmte Least-Squares-Lösung x̄, für die kx̄k2 minimal ist
(in den ersten beiden Fällen, weil es nur eine solche Lösung gibt). Diese lässt sich jeweils
mittels der QR-Zerlegung berechnen.
4. Fall rang(A) < min{m, n} (rangdefizitär):
Eine solche Lösung mit minimaler Norm gibt es immer, auch wenn der Rang von A nicht
maximal ist, allerdings hängt sie nicht mehr stetig von der Problemstellung ab.
Kleine Störungen in A oder b führen dann zu beliebig großen Änderungen in x̄ (da die ge-
störte Matrix praktisch garantiert regulär ist). So etwas nennt man ein schlecht gestelltes
Problem.
Moore-Penrose-Inverse
Das QR-Verfahren berechnet jeweils die Lösung der Gleichung x̄ = A+ b mit der Moore-Penrose-
Inversen A+ von A, gegeben durch
−1 T
falls m > n und rang(A) = n
T
(A A) A
−1
A+ = AT (AAT ) falls m < n und rang(A) = m
falls m = n und rang(A) = n
−1
A
Für die verbleibenden Fälle (rangdefizitär) ist die Pseudoinverse A+ zwar ebenfalls eindeutig,
es gibt aber keine einfache geschlossene Formel, und sie hängt nicht stetig von Änderungen in
A ab.
88
5.7 Nichtreguläre Systeme
Anwendung: Ausgleichsrechnung
Problemstellung der Gaußschen Ausgleichsrechnung:
Gegeben:
1. n Funktionen u1 , . . . , un : R → R
2. m ≥ n Datenpunkte (xi , yi ) ∈ R2 , 1 ≤ i ≤ m
Gesucht: n Koeffizienten c1 , . . . , cn , so dass
n
X m
X
y = u(x) = cj · uj (x) und (u(xi ) − yi )2 minimal
j=1 i=1
“Modell” für die Abhängigkeit der Größe y vom Parameter x, mit {u1 , . . . , un } als Ansatz-
raum.
Formuliere im Rahmen der Linearen Algebra:
c = (c1 , . . . , cn )T , y = (y1 , . . . , ym )T , A mit aij = uj (xi )
Dann gilt
2
m
X m
X Xn
(u(xi ) − yi )2 = cj · uj (xi ) − yi = kAc − yk22
i=1 i=1 j=1
Optimale Lösung ist also durch Normalengleichung gegeben:
n
X
u(x) := c̄j · uj (x) mit AT Ac̄ = AT y
j=1
• Für m > n, rang(A) = n gibt es eine optimale Kombination der uj , z.B. eine Ausgleichs-
gerade
y = u(x) = c1 + c2 · x,
die die gegebenen Datenpunkte so gut wie möglich repräsentiert (aber i.A. nicht durch
die Punkte geht).
• Für m = n, rang(A) = n kann das Modell die Datenpunkte exakt darstellen. Im Rahmen
der Polynominterpolation werden wir diesen Fall genauer betrachten.
• Die Rangdefizienz rang(A) < n entspricht hier einer Redundanz im Modell, d.h. die
gleichen Datenpunkte lassen sich auf völlig verschiedene Weisen darstellen.
89
5 Eliminationsverfahren zur Lösung linearer Gleichungssysteme
• Ein recht künstliches, aber leicht interpretierbares Beispiel ist
u1 (x) = x + 1, u2 (x) = x − 1, u3 (x) = 1
mit einer größeren Zahl Datenpunkten (xi , xi ), 1 ≤ i ≤ n ∈ N. Die Vektoren (1, 0, −1)T
und (0, 1, 1)T sind beides gültige Lösungen, und kleine Fehler in den Daten sorgen dafür,
dass jeweil einer dieser beiden völlig verschiedenen Vektoren gewählt wird.
Singulärwertzerlegung (SVD)
Satz 5.31 (Singulärwertzerlegung). Sei A ∈ Rm×n , dann gibt es orthogonale Matrizen U ∈
Rm×m und V ∈ Rn×n und eine Diagonalmatrix
D = diag(σ1 , . . . , σp ) ∈ Rm×n
mit σ1 ≥ σ2 ≥ · · · ≥ σp ≥ 0 und
U T AV = D (und daher A = U DV T )
Die Einträge von D heißen Singulärwerte, sie sind bis auf die Reihenfolge eindeutig. Die Spalten
von U und V heißen (linke und rechte) Singulärvektoren.
Beweis: Tafel
Anmerkungen
Anmerkungen zur Singulärwertzerlegung (SVD):
• Die SVD kann als Verallgemeinerung der Diagonalisierung von symmetrisch positiv se-
midefiniten Matrizen betrachtet werden.
• Die SVD zerlegt eine beliebige Matrix A in eine sehr einfache und elegante Kette von drei
Operationen:
Drehung / Spiegelung — Streckung — Drehung / Spiegelung
• Trotz ihrer praktischen Relevanz kann hier auf die Berechnung der SVD aus Zeitgründen
nicht weiter eingegangen werden.
90
5.7 Nichtreguläre Systeme
Vergleich der Verfahren
Alle vorgestellten Verfahren sind zur Lösung regulärer Systeme Ax = b geeignet. Welches sollte
man verwenden?
Hier eine Gegenüberstellung:
Methode Einschränkung Stabilität Aufwand T (n)
LR-Zerlegung A quadratisch mit Pivot. 2/3 · n3
Cholesky A s.p.d. stabil 1/3 · n3
QR-Zerlegung — sehr stabil 4/3 · n3
SVD — sehr stabil 12 · n3
LR- und Cholesky-Zerlegung sind am effizientesten, aber die QR-Zerlegung kann sich aus Grün-
den der Stabilität anbieten. Die SVD wird bei speziellen Anwendungen verwendet, wo der hohe
Aufwand gerechtfertigt ist.
91
6 Interpolation und Approximation
6 Interpolation und Approximation
6.1 Einführung
Ziel: Darstellung und Auswertung von Funktionen im Rechner.
Anwendungen:
• Rekonstruktion eines funktionalen Zusammenhangs aus “gemessenen Funktionswerten”,
Auswertung an Zwischenstellen
• Teuer auszuwertende Funktionen effizienter auswerten
• Darstellung von Fonts (2D), Körpern (3D) im Rechner
• Lösung von Differential- und Integralgleichungen
• Datenkompression
Wir beschränken uns auf Funktionen in einer Variablen, z.B.
f ∈ C r [a, b]
Dies ist ein unendlichdimensionaler Funktionenraum. Im Rechner betrachten wir Funktionen-
klassen, die durch endlich viele Parameter bestimmt sind (nicht unbedingt lineare Unterräume),
z.B.:
p(x) = a0 + a1 x + · · · + an xn (Polynome)
a0 + a1 x + · · · + an xn
r(x) = (rationale Funktionen)
b0 + b1 x + · · · + bm xm
n
1 X
t(x) = a0 + (ak cos(kx) + bk sin(kx)) (trigonom. Polynome)
2
k=1
n
X
e(x) = ak exp(bk x) (Exponentialsumme)
k=1
Approximation
Grundaufgabe der Approximation:
Gegeben eine Menge von Funktionen P (Polynome, rationale Funktionen, . . . ) und eine Funk-
tion f (z.B. f ∈ C[a, b]), finde g ∈ P , so dass der Fehler f − g in geeigneter Weise minimiert
wird.
92
6.2 Polynominterpolation
Beispiele:
Z b 1/2
2
(f − g) dx → min (2-Norm)
a
max |f (x) − g(x)| → min (∞-Norm)
a≤x≤b
max |f (xi ) − g(xi )| → min für a ≤ xi ≤ b, i = 0, . . . , n
i∈{0,...,n}
Interpolation
Man spricht von Interpolation, wenn g durch
g(xi ) = yi := f (xi ) i = 0, . . . , n
festgelegt wird.
Dies ist ein Spezialfall der Approximationsaufgabe:
• Der Fehler f − g wird nur in den Stützstellen xi , i = 0, . . . , n, betrachtet.
• In diesen endlich vielen Punkten muss die Abweichung nicht nur minimal, sondern Null
sein.
6.2 Polynominterpolation
Pn sei die Menge der Polynome über R vom Grad kleiner gleich n ∈ N0 :
n
X
Pn := {p(x) = ai xi | ai ∈ R}
i=0
Pn ist ein n + 1-dimensionaler Vektorraum.
Die Monome 1, x, x2 , . . . , xn bilden eine Basis von Pn .
Zu gegebenen (paarweise verschiedenen) n + 1 Stützstellen x0 , x1 , . . . , xn ist die Interpolations-
aufgabe dann
Finde p ∈ Pn : p(xi ) = yi := f (xi ), i = 0, . . . , n
93
6 Interpolation und Approximation
Das ist äquivalent zum linearen Gleichungssystem
x0 x20 · · · xn0
1 a0 y0
1 x1 x21 · · · xn1 a1 y1
.. .. .. . . .. · .. = ..
. . . . . . .
1 xn x2n · · · xnn an yn
| {z }
=:V [x0 ,...,xn ]
Die Matrix V [x0 , . . . , xn ]
• heißt Vandermonde-Matrix
• ist genau dann regulär, wenn alle xi paarweise verschieden sind
• ist schlecht konditioniert: cond∞ (V ) = O(2n ) für xi positiv!
• benötigt Rechenaufwand in O(n3 ) zum Lösen des LGS
Problem:
Aufstellen der Vandermonde-Matrix V [x0 , . . . , xn ] und anschließendes Lösen mit den vorgestell-
ten Verfahren ist wegen der sehr schlechten Kondition nicht zu empfehlen.
Geht das auch besser und einfacher?
Ursache ist die Wahl der Monombasis 1, x, x2 , . . . , xn , die auf eine sehr ungünstige Formu-
lierung der Interpolationsaufgabe führt. Im folgenden werden wir uns mögliche Alternativen
anschauen.
Langrange-Interpolation
Definition 6.1 (Lagrange-Polynome). Zu den paarweise verschiedenen Stützstellen xi , i =
0, . . . , n, definiere die sog. Lagrange-Polynome:
(n)
Y x − xj
Li (x) := , i = 0, . . . , n (n + 1 Stück)
xi − xj
j6=i
(n)
Die Li haben Grad n, es gilt
(
(n) 1 i=k
Li (xk ) = δik =
6 k
0 i=
(n)
und die Li bilden eine Basis von Pn .
Beweis: Tafel
94
6.2 Polynominterpolation
Existenz und Eindeutigkeit
In der Lagrange-Basis sind die gesuchten Koeffizienten ai gerade die Werte in den Stützstellen:
ai = yi . Die Interpolationsaufgabe ist so also besonders einfach zu lösen.
Satz 6.2 (Eindeutigkeit der Polynominterpolation). Zu gegebenen paarweise verschiedenen
Stützstellen x0 , . . . , xn gibt es genau ein Polynom p vom Grad n mit
p(xi ) = yi i = 0, . . . , n, yi ∈ R
Die Interpolationsaufgabe ist also eindeutig lösbar.
Beweis: Tafel
Newton-Darstellung
Nachteil der Lagrange-Polynome:
Hinzufügen einer Stützstelle ändert alle bisherigen Basispolynome. Eignet sich daher nicht für
die “inkrementelle” Erstellung des Interpolationspolynoms.
Besser ist hier die Newton-Darstellung mit Basispolynomen
i−1
Y
N0 (x) = 1; i = 1, . . . , n : Ni (x) = (x − xj )
j=0
Ni hat jeweils Grad i, es gilt
∀k < i : Ni (xk ) = 0
und die N0 , . . . , Nn bilden eine Basis von Pn .
Gestaffelte Berechnung
In x0 sind alle Newton-Basispolynome bis auf N0 Null, in x1 alle bis auf N0 und N1 , usw.
Die Interpolationsaufgabe
n
X k
X
p(xk ) = ai Ni (xk ) = ai Ni (xk ) = yk k = 0, . . . , n
i=0 i=0
führt daher auf die folgende gestaffelte Berechnung:
k−1
" #
X
a0 = y0 ; k = 1, . . . , n : ak = yk − ai Ni (xk ) /Nk (xk )
i=0
95
6 Interpolation und Approximation
Hintergrund
Die Polynominterpolation in der Sprache der Linearen Algebra:
• Die Monombasis xi , i = 0, . . . , n, führt auf eine Matrix V [x0 , . . . , xn ] (Vandermonde-
Matrix), die dicht besetzt und sehr schlecht konditioniert ist.
(n)
• Die Lagrange-Basis Li , i = 0, . . . , n, führt stattdessen auf eine Einheitsmatrix, so dass
die Lösung trivial ist. Eine Erweiterung der Stützstellen ändert aber alle Basisfunktionen.
• Die Newton-Basis Ni , i = 0, . . . , n, wiederum führt auf eine untere Dreiecksmatrix. Die
gestaffelte Berechnung entspricht daher gerade dem Vorwärtseinsetzen. Das Lösen ist
aufwendiger als bei der Lagrange-Basis, weitere Stützstellen führen aber nur zu einer
Erweiterung der unteren Dreiecksmatrix.
Dividierte Differenzen
Satz 6.3 (Dividierte Differenzen). Man definiert rekursiv die sog. “dividierten Differenzen”
∀i = 0, . . . , n : y[xi ] := yi (Werte an den Stützstellen)
∀k = 1, . . . , n − i :
y[xi+1 , . . . , xi+k ] − y[xi , . . . , xi+k−1 ]
y[xi , . . . , xi+k ] :=
xi+k − xi
Dann gilt
n
X
p(x) = y[x0 , . . . , xi ]Ni (x)
i=0
Beweis: [Rannacher]
Die dividierten Differenzen ordnet man in einem Tableau an:
y0 = y[x0 ] → y[x0 , x1 ] → y[x0 , x1 , x2 ] → y[x0 , x1 , x2 , x3 ]
% % %
y1 = y[x1 ] → y[x1 , x2 ] → y[x1 , x2 , x3 ]
% %
y2 = y[x2 ] → y[x2 , x3 ]
%
y3 = y[x3 ]
In der ersten Zeile finden sich dann die gesuchten Koeffizienten ai , i = 0, . . . , n, der Basispoly-
nome.
Diese Form der Berechnung benötigt die Werte der Ni (xk ) nicht und ist zudem stabiler.
96
6.2 Polynominterpolation
Beispiel
Beispiel 6.4. Gegeben seien die folgenden Paare von Stützstellen und Werten:
(x0 = 0, y0 = 1), (x1 = 1, y1 = 4), (x2 = 2, y2 = 3)
Die Monombasis führt auf das Gleichungssystem
1 0 0 a0 1
1 1 1 · a1 = 4
1 2 4 a2 3
mit der Lösung [1, 5, −2]T , also
p(x) = 1 + 5 · x − 2 · x2
Die Lagrange-Basis führt bei diesen Werten auf
(2) (2) (2)
p(x) = 1 · L0 (x) + 4 · L1 (x) + 3 · L2 (x)
x−1 x−2 x−0 x−2 x−0 x−1
=1· · +4· · +3· ·
0−1 0−2 1−0 1−2 2−0 2−1
1 3
= · (x − 1) · (x − 2) − 4 · x · (x − 2) + · x · (x − 1)
2 2
Die einzelnen Basispolynome können für eine gegebene Menge von Stützstellen x0 , . . . , xn einmal
ausmultipliziert werden, danach sind weitere Interpolationsaufgaben leicht zu lösen.
Für die Newton-Basis ergibt sich das Tableau
4−1 (−1)−3
y0 = a0 = 1 → a1 = 1−0 = 3 → a2 = 2−0 = −2
% %
3−4
y1 = 4 → 2−1 = −1
%
y2 = 3
und daher
p(x) = 1 · N0 (x) + 3 · N1 (x) − 2 · N2 (x)
= 1 + 3 · x − 2 · x · (x − 1)
Eine weitere Stützstelle kann hier leicht hinzugefügt werden.
97
6 Interpolation und Approximation
Auswertung von Polynomen
Um ein Polynom der Form
n
X
p(x) = ai xi = a0 + a1 x + · · · + an xn
i=0
an einer Stelle x auszuwerten, verwendet man aus Effizienz- und Stabilitätsgründen das sog.
Horner-Schema:
bn := an ; k = n − 1, . . . , 0 : bk := ak + x · bk+1 ; p(x) = b0
Für ein Polynom in Newton-Darstellung
n
X
p(x) = ai Ni (x) = a0 + a1 N1 (x) + · · · + an Nn (x)
i=0
ergibt sich analog die Rekursion
bn := an ; k = n − 1, . . . , 0 : bk := ak + (x − xk ) · bk+1 ; p(x) = b0
Interpolationsfehler
Sei yi = f (xi ), i = 0, . . . , n, die Auswertung einer Funktion f an n + 1 paarweise verschiedenen
Stützstellen, und p(x) das Polynom vom Grad n, das die zugehörige Interpolationsaufgabe
löst.
Die Differenz
e(x) := f (x) − p(x)
erfüllt nach Konstruktion die Bedingung
e(xi ) = 0 für i = 0, . . . , n
Frage: Wie groß kann die Differenz an anderen Stellen werden?
Satz 6.5 (Interpolationsfehler). Sei f (x) n + 1-mal stetig differenzierbar auf [a, b] und
a ≤ x0 < x1 < · · · < xn ≤ b.
Dann gibt es zu jedem x ∈ [a, b] ein ξx ∈ (x0 , . . . , xn , x) (kleinstes Intervall, das alle diese
Punkte enthält), so dass
n
f (n+1) (ξx ) Y
f (x) − p(x) = (x − xj )
(n + 1)!
j=0
Beweis: Tafel
98
6.2 Polynominterpolation
Anmerkungen
Für den Spezialfall äquidistanter Stützstellen, d.h.
xk+1 − xk = h für k = 0, . . . , n − 1,
gilt damit
|f (x) − p(x)| ≤ |f (n+1) (ξx )| · hn+1
Für |f (n+1) | beschränkt und n → ∞ gilt also |f (x) − p(x)| → 0.
Allerdings sind die höheren Ableitungen auch einfacher Funktionen für n → ∞ oft nicht be-
schränkt, sondern wachsen sehr schnell.
Runges Gegenbeispiel
1
0.5
0
-0.5
-1
-1 -0.5 0 0.5 1
−1
Polynominterpolation von Runges Funktion f (x) = (1 + 25x2 ) mit äquidistanten Stützstellen
(3, 5, 9, 17 bzw. 33 Wertepaare).
Die Minima / Maxima der letzten beiden Polynome sind −14.35/1.40 bzw. −5059/2.05 (!).
Anmerkungen
Bemerkung 6.6. Laut Approximationssatz von Weierstraß lässt sich jede Funktion aus C 0 ([a, b])
gleichmäßig durch Polynome approximieren.
Die beobachteten Phänomene sind kein Widerspruch, denn:
• Die Approximation muss nicht durch Interpolation erfolgen (der Beweis nutzt Bernstein-
Polynome).
99
6 Interpolation und Approximation
• Mit nicht-äquidistanten Stützstellen erhält man bereits deutlich bessere Ergebnisse (wenn
man weiß, wie sie zu wählen sind. . . ).
Bemerkung 6.7. Allgemein gilt für “Methoden hoher (Polynom-) Ordnung”, dass entsprechen-
de Differnzierbarkeit vorliegen muss.
Konditionierung
Wenn p(x; y) das Interpolationspolynom zu den Ordinatenwerten (y0 , . . . , yn )T und festen Ab-
szissenwerten (x0 , . . . , xn )T ist, gilt
n n
(n) (n)
X X
p(x; y + ∆y) − p(x; y) = (yi + ∆yi )Li (x) − yi Li (x)
i=0 i=0
n
(n)
X
= ∆yi Li (x)
i=0
Somit ist
n (n)
p(x; y + ∆y) − p(x; y) X Li (x)yi ∆yi
= ·
p(x; y) p(x; y) yi
i=0
(n)
Für große n kann Li sehr groß werden, die Interpolationsaufgabe ist dann schlecht konditio-
niert!
6.3 Anwendungen der Polynominterpolation
Numerische Differentiation
Problem:
Berechne die Ableitung (der Ordnung k) einer tabellarisch gegebenen Funktion oder einer im
Rechner gegebenen Funktion (im Sinne der Informatik).
Idee:
Erstelle Interpolationspolynom zu bestimmten Stützstellen, leite dieses ab und werte aus.
Dabei sei zunächst Ableitungsordnung = Polynomgrad.
Die Lagrange-Polynome sind
(n)
Y x − xj
(xi − xj )−1 ·xn + αn−1 xn−1 + · · · + α0 ,
Y
Li (x) = =
xi − xj
j6=i j6=i
| {z }
=:λi ∈R
100
6.3 Anwendungen der Polynominterpolation
also liefert n-faches Differenzieren:
dn (n)
L (x) = n! · λi
dxn i
Somit gilt für die n-te Ableitung eines Interpolationspolynoms vom Grad n:
n n
!
dn X (n)
X
yi Li (x) = n! · yi λi (unabh. von x)
dxn
i=0 i=0
Eine Aussage über den entstehenden Fehler liefert:
Satz 6.8. Sei f ∈ C n ([a, b]) und a = x0 < x1 < · · · < xn = b. Dann gibt es ξ ∈ (a, b), so dass
n
X
f (n) (ξ) = n! · yi λ i
i=0
Die über das Interpolationspolynom berechnete Ableitung stimmt also in mindestens einem
Punkt mit der von f überein.
Beweis: Tafel
Für äquidistante Stützstellen, xi = x0 + ih, 0 ≤ i ≤ n, erhält man speziell
n
−n
X n−i n y1 − y0
f (n)
(x) ≈ h (−1) yi , z.B. f (1) (x) ≈ ,
i h
i=0
y2 − 2y1 + y0 y3 − 3y2 + 3y1 − y0
f (2) (x) ≈ , f (3) (x) ≈
h2 h3
Mittels Taylorentwicklung zeigt man
f (x + h) − f (x − h)
f 0 (x) = + O(h2 ) für f ∈ C 3 ,
2h
f (x + h) − 2f (x) + f (x − h)
f 00 (x) = + O(h2 ) für f ∈ C 4
h2
101
6 Interpolation und Approximation
Extrapolation zum Limes
Eine Größe a(h) sei im Rechner für h > 0 berechenbar, nicht jedoch für h = 0. Man möchte
a(0) = lim a(h)
h→0
mit guter Genauigkeit berechnen.
Beispiel 6.9. Mögliche Anwendungen:
1. Regel von de l’Hospital:
cos(x) − 1
a(0) = lim (= 0)
h→0 sin(x)
2. Numerische Differentiation:
f (x + h) − f (x)
f 0 (x) = lim
h→0 h
(für kleine h tritt Auslöschung auf )
3. Numerische Integration:
Z b n
X 1
f (x) = lim h·f a+ i− (b − a) · h
a h→0 2
i=1
4. Numerische Lösung des Anfangswertproblems
y 0 (t) = f (t, y(t)) auf [0, T ]; y(0) = y0
Setze
h = N −1 ; yn = yn−1 + h · f (t, yn−1 ); y(T ) ≈ yN
Hier ist h → 0 gleichbedeutend mit N → ∞ und entsprechend ansteigendem Aufwand.
Grundlegende Idee
Idee der Extrapolation:
Zu h0 > h1 > · · · > hn > 0 bestimme Interpolationspolynom
p(hi ) = a(hi ) i = 0, . . . , n
und berechne
a(0) ≈ p(0)
(Extrapolation statt Interpolation, da 0 6∈ [hn , . . . , h0 ])
102
6.3 Anwendungen der Polynominterpolation
Beispiel
Beispiel 6.10. Für a(h) = (cos(h) − 1) · (sin(h))−1 erhält man für
h0 = 1/8 : a(h0 ) = −6.258151 · 10−2
h1 = 1/16 : a(h1 ) = −3.126018 · 10−2
h2 = 1/32 : a(h2 ) = −1.562627 · 10−2
(halbieren von h halbiert also a(h)), und bei Extrapolation mit einem Polynom p2 vom Grad 2:
a(0) ≈ p2 (0) = −1.02 · 10−5
also sehr viel besser als die ursprünglich berechneten Näherungen oder eine mögliche direkte
Auswertung für h 1 (Auslöschung)!
Warum ist das so gut?
Sei hi = h · ri mit r < 1 (geometrische Verteilung), z.B. r = 1/2, und sei p das Interpolations-
polynom von a zu den Stützstellen hi . Dann gilt für den Fehler der Extrapolation
|a(n+1) (ξ)|
|p(0) − a(0)| ≤ kV −T k∞ hn+1
(n + 1)! · (1 − rn+1 )
mit der Vandermonde-Matrix V und ξ ∈ (0, h), sofern a hinreichend differenzierbar ist.
Beweis: Tafel
Extrapolation bei Ableitungen
Entscheidend ist die Taylor-Entwicklung von a in Null (Mclaurin-Entwicklung). Für den be-
kannten Differenzenquotienten der zweiten Ableitung erhält man
f (x + h) − 2f (x) + f (x − h)
a(h) =
h2
h 2 h2n
= f 00 (x) + f (4) (x) + · · · + f (2n+2) (x)
2 · 4! 2 · (2n + 2)!
h2n+2
+ [f (2n+4) (ξ+ ) + f (2n+4) (ξ− )]
2 · (2n + 4)!
= px (h2 ) + O(h2(n+1) )
Man kann also (falls f hinreichend differenzierbar ist) mit jeder Auswertung sogar zwei h-
Potenzen gewinnen!
103
6 Interpolation und Approximation
6.4 Bernstein-Polynome und Kurvendarstellung
Bernstein-Polynome
Wir gehen nun von der Interpolation zur Approximation über.
Speziell für Kurven, d.h. Funktionen u(t) : [a, b] → Rd , haben sich die Bernstein-Polynome
bewährt.
Definition 6.11 (Bernstein-Polynome). Die Polynome
(n) n
βi (t) := (1 − t)n−i ti i = 0, . . . , n
i
heißen Bernstein-Polynome auf dem Interval [0, 1].
Mittels der affin-linearen Transformation
t−a
φ : [a, b] → [0, 1], p(t) =
b−a
definiert man zusätzlich Bernstein-Polynome auf einem allgemeinen Intervall [a, b]:
(n) (N ) n
βi,[a,b] (t) := βi (φ(t)) = (b − a)−n (b − t)n−i (t − a)i
i
Eigenschaften
Satz 6.12 (Eigenschaften der Bernstein-Polynome). Die Bernstein-Polynome haben die fol-
genden Eigenschaften:
1. Zerlegung der Eins:
n
(n)
X
βi (t) = 1
i=0
(n)
2. t = 0 ist i-fache Nullstelle von βi
(n)
3. t = 1 ist (n − i)-fache Nullstelle von βi
4. Symmetrie:
(n) (n)
βi (t) = βn−i (1 − t)
5. Positivität:
(n) (n)
0 ≤ βi (t) ≤ 1 für t ∈ [0, 1], 0 < βi (t) für t ∈ (0, 1)
(n)
6. βi , n 6= 0, hat in [0, 1] genau ein Maximum in i/n.
104
6.4 Bernstein-Polynome und Kurvendarstellung
(n)
7. Die βi sind lin. unabh. und bilden eine Basis von Pn .
8. Die Bernstein-Polynome erlauben eine rekursive Darstellung:
(n−1)
(1 − t)β0
(t) i=0
(n) (n−1) (n−1)
βi (t) = tβi−1 (t) + (1 − t)βi (t) 0<i<n
(n−1)
tβn−1 (t) i=n
9. Für die erste Ableitung gilt die Rekursionsformel:
(n−1)
d (n) −nβ0
(t) i=0
β = n[βi−1 (t) − βi(n−1) (t)]
(n−1)
0<i<n
dt i (n−1)
nβn−1 (t) i=n
Achtung: mit dem Superskript sind hier jeweils die Grade der Polynome gemeint, nicht
höhere Ableitungen!
Zusammenfassend bilden die Bernstein-Polynome eine “Konvexkombination” (positive Zerle-
gung der Eins) mit zusätzlicher Symmetrieeigenschaft und rekursiver Darstellung.
Beweis: Tafel
Graphische Darstellung
Bernstein-Polynome vom Grad 5:
1
0.8
0.6
0.4
0.2
0
0 0.2 0.4 0.6 0.8 1
105
6 Interpolation und Approximation
Approximation mit Bernstein-Polynomen
Approximation des Datenvektors y = (0, 8, 6, 2, 0, 8)T mit Bernstein-Polynomen:
8
7
6
5
4
3
2
1
0
0 0.2 0.4 0.6 0.8 1
Die Daten werden nirgends im Intervall interpoliert, stattdessen hat jeder Datenpunkt einen
“Bereich des Einflusses”.
Bézier-Kurven
Kurvendarstellung mittels Bernstein-Polynomen:
Definition 6.13 (Bézier-Kurven). Für gegebene Punkte b0 , . . . , bn ∈ Rd heißt das vektorwertige
Polynom
n
(n)
X
B(t) := bi βi (t)
i=0
Bézier-Kurve. Bézier war Ingenieur bei Renault, und die nach ihm benannten Kurven sind die
Grundlage von Vektorgraphiken.
Beispiel
Beispiel 6.14. Betrachte die Punkte b0 = (0, 0)T , b1 = (1/3, −1)T , b2 = (2/3, −1)T und
b3 = (1, 1)T , die zugehörige Bézier-Kurve ist:
3
X (3) t
B(t) = bi βi (t) =
t + 3t2 − 3t
3
i=0
106
6.4 Bernstein-Polynome und Kurvendarstellung
(1,1)
(0,0)
(1/3,-1) (2/3,-1)
Eigenschaften
Es gelten folgende Eigenschaften:
1. Die Bézier-Kurve B(t) liegt in der konvexen Hülle der Bézier-Punkte b0 , . . . , bn .
2. Es ist B(0) = b0 und B(1) = bn .
3. Die Ableitung (Tangente der Kurve) hat an den Endpunkten die Richtung (b1 − b0 ) bzw.
(bn − bn−1 ).
Dadurch sind also glatte, lokalisierte Kurven mit klar definierten Enden gegeben.
Bernstein-Polynome erlauben eine gute Kontrolle über die Werte eines Polynoms durch dessen
Koeffizienten.
Vergleich mit anderen Polynombasen:
• Monom- und Newton-Basis: keine einfache Kontrolle der Werte
• Langrange-Basis: exakte Kontrolle an den Stützstellen, dazwischen aber (sehr) große Ab-
weichungen
Pn (n)
• Bernstein-Basis: min{bi } ≤ i=0 bi βi ≤ max{bi }
Das ist eine direkte Folge der Konvexkombination.
107
6 Interpolation und Approximation
Algorithmus von de Casteljau
(0) (0)
Für gegebene Bézier-Punkte b0 , . . . , bn ergibt sich mittels der Rekursionsformel aus Satz 6.12:
n n−1
(0) (n) (0) (0) (n−i)
X X
B(t) = bi βi (t) = [bi (1 − t) + bi+1 (t)] βi (t)
i=0 i=0
| {z }
(1)
=:bi
Rekursiv setzt man daher die (von t abhängigen) “Koeffizienten”
(k)
bi := b(k−1)
i
(k−1)
(1 − t) + bi+1 t 0 < k ≤ n, 0 ≤ i ≤ n − k,
(n)
um B(t) = b0 zu berechnen.
Dieser Algorithmus ist nach de Casteljau, Physiker und Mathematiker bei Citroën, benannt.
Die Struktur der Auswertung ist mit den dividierten Differenzen vergleichbar:
(0) (1) (2) (3)
b0 = b0 → b0 → b0 → b0 = B(t)
% % %
(0) (1) (2)
b1 = b1 → b1 → b1
% %
(0) (1)
b2 = b2 → b2
%
(0)
b3 = b3
Bei fester Wahl von t ist damit auch ein Verfahren gegeben, um eine Bézier-Kurve in zwei
Teilkurven zu zerlegen. Die neuen Bézier-Punkte sind
(0) (1) (n) (n) (n−1)
b0 , b0 , . . . , b0 bzw. b0 , b1 , . . . , b(0)
n
(die ersten bzw. letzten Punkte der Rekursionsebenen).
(0) (n) (0)
Offensichtlich liegen die Start- / Endpunkte b0 , b0 und bn auf der Kurve. Dass die Kurven
(j)
sich auch tatsächlich überlagern, ist eine Folge der Berechnung der bi und der rekursiven
Darstellung der Bernstein-Polynome.
Beispiel
Am besten lässt sich das an der graphischen Darstellung erkennen:
Beispiel 6.15. Zerlegung der Kurve aus Beispiel 6.14 an der Stelle t = 1/3, bzw. Auswertung
an dieser Stelle: rekursives Dritteln der Verbindungslinien, resultierende neue Punkte in rot
und blau.
108
6.5 Splines
6.5 Splines
Wir kehren jetzt wieder zur Interpolation zurück.
Bis jetzt galt:
• Anzahl Stützstellen = Polynomgrad +1
• Viele Stützstellen =⇒ hoher Polynomgrad =⇒ evtl. starke Abweichung zwischen den
Stützstellen
Idee: Benutze stückweise Polynome niedrigen Grades.
Definition 6.16 (Spline-Raum). Sei X = (x0 , x1 , . . . , xn ) mit a = x0 < x1 < · · · < xn = b
eine Zerlegung des Intervalls [a, b] und sei m ∈ N. Die Menge
S m (X) := {s ∈ C m−1 ([a, b]) : s|[xi ,xi+1 ] ∈ Pm , 0 ≤ i < n}
heißt Spline-Raum vom Grad m über der Zerlegung X.
Beispiel
Beispiel 6.17. s ∈ S 1 (X) bedeutet:
• s ist auf jedem Teilintervall [xi , xi+1 ] ein Polynom vom Grad 1
• S 1 (X) ⊂ C 0 ([a, b]), also ist s stetig
109
6 Interpolation und Approximation
Dieser Vektorraum ist relevant beim Lösen partieller Differentialgleichungen, und wird dort als
Finite-Elemente-Raum P1 bezeichnet.
Kubische Splines
In der Praxis ist S 3 (X) sehr beliebt. S 3 (X) heißt Raum der kubischen Splines.
Geschichte: “Straklatte” zur Konstruktion glatter Kurven im Schiffs- und Flugzeugbau.
Ein dünnes Balsaholz biegt sich unter Energieminimierung
b
|y 00 (t)|2 b
Z Z
2
= dt ≈ |y 00 (t)| dt → min
a 1 + |y 0 (t)| a
Quelle: Pearson Scott Foresman, gemeinfrei
110
6.5 Splines
Eine Funktion s ∈ S 3 (X) besteht stückweise aus n Polynomen:
(
pi (x) x ∈ [xi−1 , xi ), i ∈ {1, . . . , n}
s(x) =
pn (xn ) x = xn (Intervalabschluss)
Für die Polynome pi gelten folgende Bedingungen:
• Interpolationsbedingung (Stetigkeit) in allen Punkten:
(
pi (xi−1 ) = yi−1
2n Bedingungen: i = 1, . . . , n :
pi (xi ) = yi
• Stetigkeit von erster und zweiter Ableitung an inneren Punkten:
(
p0i (xi ) = p0i+1 (xi )
2(n − 1) Bedingungen: i = 1, . . . , n − 1 :
p00i (xi ) = p00i+1 (xi )
=⇒ zusammen 4n − 2 Bedingungen
Man hat 4n Freiheitsgrade, pro Polynom pi vier (da Grad 3).
Die fehlenden 2 Bedingungen erhält man als Randbedingungen an den Stellen x0 und xn . Dabei
gibt es verschiedene Varianten:
1. “Natürliche” Randbedingungen, “Krümmung Null”:
p001 (x0 ) = p00n (xn ) = 0
2. Hermite-Randbedingungen:
p01 (x0 ) = f 0 (x0 ), p0n (xn ) = f 0 (xn )
(f ist die zu interpolierende Funktion)
3. Periodische Randbedingungen:
p01 (x0 ) = p0n (xn ), p001 (x0 ) = p00n (xn )
Wir behandeln im folgenden nur die natürlichen Randbedingungen.
111
6 Interpolation und Approximation
Kubische Splines
Kubische Splines mit natürlichen Randbedingungen lassen sich über eine Minimierungseigen-
schaft charakterisieren:
Zu gegebenen Funktionswerten f (xi ) = yi ist der natürliche kubische Spline unter den die
Werte interpolierenden Funktionen aus C 2 ([a, b]) jene, die das Funktional
Z b
2
|f 00 (x)| dx
a
minimiert, und somit durch eine besonders geringe Oszillation ausgezeichnet ist.
Beweis: [Rannacher]
Berechnung kubischer Splines
Wir schreiben die Teilpolynome eines kubischen Splines in der Form
(i) (i) (i) (i)
pi (x) = a0 + a1 (x − xi ) + a2 (x − xi )2 + a3 (x − xi )3 , i = 1, . . . , n
und nehmen natürliche Randbedingungen an, also
p001 (x0 ) = p00n (xn ) = 0
Dadurch ist der Spline wie eben gezeigt eindeutig festgelegt, aber:
(i)
Wie bestimmt man die lokalen Koeffizienten aj ?
Satz 6.18 (Berechnung kubischer Splines). Definiere die lokale Maschenweite hi := xi − xi−1
(0) (n) (i)
und a2 = a2 = 0. Die Koeffizienten a2 der gewählten Darstellung des kubischen Splines
sind dann die Lösung des linearen Gleichungssystems
(i−1) (i) (i+1) yi+1 − yi yi − yi−1
hi a2 + 2(hi + hi+1 )a2 + hi+1 a2 =3 −
hi+1 hi
der Dimension n − 1. Die restlichen Koeffizienten ergeben sich zu:
(i) (i−1)
(i) (i) yi − yi−1 hi (i) (i−1) (i) a2 − a2
a0 = yi , a1 = + (2a2 + a2 ), a3 =
hi 3 3hi
Beweis: Tafel
112
6.5 Splines
Lösung des Tridiagonalsystems
Zur Lösung des entstehenden Tridiagonalsystems:
• Die Gaußelimination / LR-Zerlegung hat in diesem Fall Aufwand O(n) (Bandstruktur)
• Das Gleichungssystem ist symmetrisch und strikt diagonaldominant:
X
|aij | < |aii |
j6=i
Die Koeffizienten können also effizient und stabil ohne Pivotisierung berechnet werden.
Beispiel
Beispiel 6.19. Approximationseigenschaften der verschiedenen Ansätze:
Fehlerabschätzung
Satz 6.20 (Fehlerabschätzung für kubische Splines). Sei f ∈ C 4 ([a, b]). Erfüllt der kubische
Spline s zu f
s00 (a) = f 00 (a) und s00 (b) = f 00 (b)
(Hermitesche Randbedingungen), so gilt
1
max |f (x) − s(x)| ≤ h4 max |f (4) (x)|
a≤x≤b 2 a≤x≤b
Selbst unter noch wesentlich schwächeren Voraussetzungen konvergiert die Spline-Interpolation
gleichmäßig gegen f .
Beweis: siehe [Rannacher]
113
6 Interpolation und Approximation
Beispiel
Beispiel 6.21. Interpolation der Funktion exp(−x2 ) auf dem Intervall [−10, 10] mit n äquidi-
stanten Stützstellen, maximaler Punktfehler:
n S1 S3
4 6.0 · 10−1 7.4 · 10−1
8 3.0 · 10−1 3.9 · 10−1
16 1.1 · 10−1 2.8 · 10−2
32 6.9 · 10−2 7.1 · 10−3
64 2.2 · 10−2 3.3 · 10−4
128 6.0 · 10−3 1.9 · 10−5
256 1.5 · 10−3 1.2 · 10−6
512 3.8 · 10−4 7.3 · 10−8
n Pn
4 8.0 · 10−1
6 1.0 · 100
8 2.3 · 100
10 5.9 · 100
Interpolation mit Splines konvergiert mit Ordnung h2 für S 1 bzw. h4 für S 3 .
6.6 Trigonometrische Interpolation
Aufgabe:
Interpolation “periodischer” Funktionen, d.h. es gebe ein ω ∈ R, ω > 0, so dass
f (x + ω) = f (x) ∀x ∈ R
Es bietet sich die Interpolation mit “trigonometrischen Summen”
m
a0 X
ak cos(2πkω −1 x) + bk sin(2πkω −1 x)
tn (x) = +
2
k=1
an, denn jeder der Summanden hat Periode ω.
Es gibt 2m + 1 Parameter, also setze n = 2m (konsistent mit Definition von n bei Polynomin-
terpolation).
114
6.6 Trigonometrische Interpolation
Der Einfachheit halber betrachten wir nur ω = 2π, also:
m
a0 X
tn (x) = + [ak cos(kx) + bk sin(kx)]
2
k=1
Die Interpolation führen wir dann mit den äquidistanten Stützstellen
2πj
xj = j = 0, . . . , n
n+1
durch. (Für den allgemeinen Fall transformiere man das Argument und die Stützstellen).
Es stellt sich heraus, dass die Interpolationsaufgabe
tn (xj ) = f (xj ) = yj j = 0, . . . , n
zunächst im Körper C einfacher zu lösen ist, d.h. betrachte das komplexe trigonometrische
Polynom
Xn
∗
tn (x) = ck eikx
k=0
mit i = (−1)1/2 der imaginären Einheit und ck ∈ C komplexen Koeffizienten. Für φ ∈ R gilt
die Eulersche Identität:
eiφ = cos(φ) + i sin(φ)
Komplexe Einheitswurzeln
Hilfssatz 6.22 (Komplexe Einheitswurzeln). Die (n + 1)-ten komplexen Einheitswurzeln
2πk
wk := eixk = ei n+1
haben für alle j, k ∈ Z die folgenden Eigenschaften:
1. wkn+1 = 1, wkj = wjk und wk−j = wj−k
j mod(n+1) j mod(n+1)
2. wkj = wk = wkj mod(n+1) = wk mod(n+1)
3. Es gilt
n
(
n+1 k mod(n + 1) = 0
wkj
X
=
j=0
0 sonst
Beweis: Tafel
Die komplexen Einheitswurzeln für n = 7 (also n + 1 = 8):
115
6 Interpolation und Approximation
Im
i w1
1
0 w0 Re
Komplexe Multiplikation ist eine Drehstreckung, und für Werte auf dem Einheitskreis entfällt
die Streckung (der Faktor ist Eins).
Es handelt sich also um eine Drehung, und die Einheitswurzeln drehen jeweils um ein rationales
Vielfaches von 360 Grad.
Komplexe Trigonometrische Interpolation
Satz 6.23 (Komplexe Trigonometrische Interpolation). Zu gegebenen Zahlen y0 , . . . , yn ∈ C
gibt es genau eine Funktion der Gestalt
n
X
∗
tn (x) = ck eikx ,
k=0
2πj
die den n + 1 Interpolationsbedingungen t∗n (xj ) = yj für xj = n+1 entspricht. Die komplexen
Koeffizienten sind bestimmt durch
n n
ck = (n + 1)−1 yj e−ijxk = (n + 1)−1 yj wk−j k = 0, . . . , n
X X
j=0 j=0
Beweis: Tafel
Aufwand und reeller Fall
Mit der Formel aus Satz 6.23 ergibt sich ein Aufwand von O(n2 ) für die Berechnung der
Koeffizienten ck .
Damit wurde auch die reelle Interpolationsaufgabe gelöst, da man in Satz 6.23 yj ∈ R annehmen
kann: aufgrund der Eulerschen Identität
eiφ = cos(φ) + i sin(φ)
ist dann auch eine passende reelle Funktion gegeben.
Wie berechnet man nun die reellen Koeffizienten ak und bk ?
116
6.6 Trigonometrische Interpolation
Diskrete Fourier-Analyse
Satz 6.24 (Diskrete Fourier-Analyse). Für n ∈ N0 gibt es zu gegebenen reellen Zahlen y0 , . . . , yn
genau ein trigonometrisches Polynom der Form
m
a0 X θ
tn (x) = + [ak cos(kx) + bk sin(kx)] + am+1 cos((m + 1)x)
2 2
k=1
das die Interpolationsaufgabe löst, wobei (m = n/2, θ = 0) für n gerade, und (m = (n−1)/2, θ =
1) für n ungerade. Es gilt
n n
ak = 2(n + 1)−1 bk = 2(n + 1)−1
X X
yj cos(jxk ), yj sin(jxk )
j=0 j=0
Beweisskizze: Tafel, genauer: [Rannacher]
Diskrete Fouriertransformation (DFT)
Die schnelle Fouriertransformation (FFT):
Einer der berühmtesten Algorithmen der angewandten Mathematik und Informatik! Entwickelt
1965 von James Cooley und John Tukey.
Ergebnis der bisherigen Überlegungen sind zwei Transformationen:
N −1
X 2πjk
−1
ck = N yj e−i N (Hintransformation)
j=0
N −1
X 2πjk
yj = ck ei N (Rücktransformation)
k=0
mit N := n + 1. Diese beiden Gleichungen bezeichnet man als diskrete Fouriertransformation
(DFT).
Setze
c := (c0 , . . . , cN −1 )T , y := (y0 , . . . , yN −1 )T ,
dann ist die diskrete Fouriertransformation äquivalent zu den Gleichungssystemen
c = N −1 W y, y = U c mit
−i 2πjk
(W )k,j = e N = wk−j und (U )j,k = wjk
Es gilt hier also W −1 = N −1 U , so dass man eine explizite Formel für die Inverse hat. Da
beide Matrizen dicht besetzt sind, ist der Aufwand für Hin- und Rücktransformation jeweils
O(N 2 ).
Grundidee:
117
6 Interpolation und Approximation
Quelle: Phonical, wikipedia.org, cc-by-sa
Die DFT zerlegt ein (periodisch fortgesetztes) diskretes Signal in sein Frequenzspektrum.
Schnelle Fouriertransformation (FFT)
Falls N gerade ist, kann man die DFT der Länge N durch zwei DFT der Länge N/2 berechnen
(Teile-und-herrsche-Verfahren, divide and conquer).
Setze c̃k := N ck (d.h. ignoriere Vorfaktor), dann gilt
( g 2πk
c̃k + e−i N c̃uk 0 ≤ k < N/2
c̃k = 2πk
c̃gk−N/2 + e−i N c̃uk−N/2 N/2 ≤ k < N
mit den geraden und ungeraden Anteilen
N/2−1 N/2−1
−i 2πjk −i 2πjk
c̃gk :=
X X
y2j e N/2 , c̃uk := y2j+1 e N/2
j=0 j=0
Wenn N eine Zweierpotenz ist, N = 2d , lässt sich die DFT somit rekursiv auf den Fall der
Länge 2 zurückführen. Dieser Fall lässt sich sehr einfach direkt lösen.
Im Gegensatz zum direkten Ansatz mittels Matrixmultiplikation ergibt sich damit der Auf-
wand
TFFT (N ) = 2 · TFFT (N/2) + Tmerge (N ) = O(N · log(N )),
also deutlich schneller für große N . Dies ist ein Spezialfall des Master-Theorems für rekursiv
definierte Funktionen.
118
6.6 Trigonometrische Interpolation
Beispiel
Beispiel der rekursiven Zerlegung für N = 8:
Quelle: Virens, wikipedia.org, cc-by-sa
Berechnung der DCT für N/2 = 4 für die geraden bzw. ungeraden Indizes, danach Rekonstruk-
tion der Ausgabe durch Zusammenführen von je zwei Zwischenergebnissen.
Anwendungen
Wichtige Anwendungen der FFT sind:
• Spektralanalyse (Zerlegung eines Signals in Frequenzen)
• Datenkompression (Vernachlässigung hochfrequenter Anteile)
• Lösen partieller Differentialgleichungen (Spektralmethoden)
Benutzt wird dies unter anderem für:
• Telekommunikation (z.B. WLAN, LTE-Mobilfunk)
• Multimedia (Datenformate wie JPEG, MPEG, MP3)
• Medizin (Kernspin- und Computertomographie)
• Forschung (z.B. digitale Oszilloskope)
119
6 Interpolation und Approximation
6.7 Gauß-Approximation von Funktionen
Prähilberträume
Wir betrachten nun die Approximation von Funktionen in Prähilberträumen:
Definition 6.25 (Prähilbertraum). Ein Vektorraum von Funktionen über R oder C mit Ska-
larprodukt (·, ·) heißt Prähilbertraum. Das Skalarprodukt induziert eine zugehörige Norm kf k =
(f, f )1/2 .
Im folgenden sei H ein Prähilbertraum und S ⊂ H ein endlichdimensionaler Teilraum von
H.
Beispiele
Beispiel 6.26. Beispiele für Prähilberträume:
• Raum der stetigen Funktionen C 0 ([a, b]) mit Skalarprodukt
Z b
(f, g) = f (x)g(x) dx
a
• Gegeben Ψ = {ψ1 , . . . , ψN }, ψi ∈ C 0 ([a, b]), ist
N
X
S = span(Ψ) = {f | f = ci ψi }
i=1
ein endlichdimensionaler Prähilbertraum. Analog für andere Prähilberträume statt C 0 .
• Raum der quadratintegrierbaren Funktionen
Z b
2
L ((a, b)) = {f | |f (t)|2 dt < ∞}
a
Rb
Dabei ist mit a das “Lebesgue-Integral” gemeint. Im Gegensatz zu C 0 ([a, b]) ist L2 ((a, b))
vollständig bzgl. der Norm
kf k = (f, f )1/2
und damit sogar ein Hilbertraum. L2 ((a, b)) ist das Analogon des Rn für Funktionenräume.
120
6.7 Gauß-Approximation von Funktionen
Allgemeine Gauß-Approximation
Wir betrachten die folgende Aufgabe:
Zu f ∈ H finde g ∈ S, so dass
kf − gk → min
wobei kf k = (f, f )1/2 die durch das Skalarprodukt induzierte Norm bezeichne.
Satz 6.27 (Allgemeine Gauß-Approximation). Die obige Aufgabe hat genau eine Lösung g ∈ S.
Diese ist charakterisiert durch
(g, φ) = (f, φ) ∀φ ∈ S
Beweis: Tafel
Berechnung der Approximation
Wähle Basis Ψ = {ψ1 , . . . , ψN } des endlichdimensionalen Raums S, N = dim(S). Stelle damit
die sogenannte Massenmatrix (Gramsche Matrix) A und rechte Seite b mit
[A]ij = (ψi , ψj ), [b]i = (f, ψi )
auf. A ist symmetrisch positiv definit bzw. hermitesch. Bestimme die gesuchten Koeffizienten
αi durch Lösen des linearen Gleichungssystems
Aα = b.
Approximation mit Orthonormalbasen
Besonders einfach wird die Lösung der Approximationsaufgabe, wenn Ψ eine Orthonormalbasis
ist, d.h. (ψi , ψj ) = δij . Dann gilt
N
X
αj (ψj , ψi ) = αi = (f, ψi ) ∀i = 1, . . . , N
j=1
und somit
N
X N
X
g= αj ψj = (f, ψj )ψj
j=1 j=1
121
6 Interpolation und Approximation
Beispiel
Beispiel 6.28 (Fourierreihe). Für N = 2m + d, m ∈ N ist
ΨF = π −1/2 · {2−1/2 , cos(x), . . . , cos(mx), sin(x), . . . , sin(mx)}
eine Orthonormalbasis von Funktionen auf dem Interval [−π, π]. Damit gilt dann
m
a0 X
g(x) = + [ak cos(kx) + bk sin(kx)]
2
k=1
Z π Z π
mit ak = f (x) cos(kx) dx, bk = f (x) sin(kx) dx
−π −π
Vorsicht: die ak und bk sind nicht die αi von oben, sondern unterscheiden sich davon um einen
Faktor π 1/2 !
Für unendlich viele Glieder nennt man die entstehende Reihe, also
∞
a0 X
g(x) = + [ak cos(kx) + bk sin(kx)],
2
k=1
Fourier-Reihe von f . Diese konvergiert gegen ein Element aus L2 ((−π, π)).
Die DFT lässt sich als näherungsweie Auswertung der Integrale mittels “Trapezregel” verstehen.
Fehlerkontrolle
Bisher war S ⊂ H fest gewählt, damit wurde das “optimale” g bestimmt. Der Fehler kf − gk
wurde einfach akzeptiert.
Daher nun Verfeinerung der Approximationsaufgabe:
Finde S ⊂ H mit dim(S) möglichst klein, so dass
kf − gk < TOL,
wobei TOL eine vorgegebene Zahl (Toleranz) sei.
Ein Spezialfall ist die “Kompression”: hier ist H selbst bereits endlichdimensional. Man versucht,
das “Signal” f möglichst platzsparend zu speichern.
Mittels Orthonormalbasen lässt sich der Fehler recht einfach messen:
122
6.7 Gauß-Approximation von Funktionen
Sei SN ⊂ H, N ∈ N, eine Folge von Approximationsräumen mit dim(SN ) = N , und gN die
Bestapproximation von f in SN . Sei zusätzlich ΨN eine Orthonormalbasis von SN .
Dann gilt
N
X
0 ≤ kf − gN k2 = (f, f ) − (f, ψi )2
i=1
Oft gilt sogar SN ⊂ SN +1 mit ΨN ⊂ ΨN +1 , das setzen wir aber nicht voraus.
Wir möchten
N
X
2
kf − gN k = (f, f ) − (f, ψi )2 ≤ TOL ·(f, f ),
i=1
und das ist äquivalent zu
N
X
(f, ψi )2 ≥ (1 − TOL) · (f, f ) (6.1)
i=1
(Frage: “Welche Anteile von f werden von SN bereits repräsentiert?”)
Achtung: hier wurde eine relative Toleranz TOL gewählt, diese lässt sich jedoch in eine absolute
umrechnen, und umgekehrt.
Anmerkungen
Bemerkung 6.29. Zur Fehlerkontrolle:
• Der Wert (f, f ) wird als (zumindest hinreichend genau) berechenbar angenommen. Dafür
kann man z.B. numerische Integration verwenden.
• Falls ΨN ⊂ ΨN +1 , sind einfach weitere Basisfunktionen hinzuzufügen bis (6.1) erreicht
ist. Das gilt zum Beispiel für die Fourierreihe, Bsp. 6.28.
• Bei ΨN ⊂ ΨN +1 folgt aus der Fehlerdarstellung unmittelbar
kf − gN +1 k2 ≤ kf − gN k2 ,
der Fehler ist also monoton fallend!
123
6 Interpolation und Approximation
Approximation mit Haar-Wavelets
Es bezeichne tr(f ) = {x | f (x) 6= 0} den Träger einer Funktion. Bei der ONB aus Sinus- und
Cosinus-Funktionen haben alle ψi im wesentlichen globalen Träger, d.h. tr(f ) = [a, b].
Oft variiert die zu approximierende Funktion aber nur lokal sehr stark. In diesem Fall möchte
man ein Funktionensystem mit folgenden Eigenschaften:
• (ψi , ψj ) = δij (Orthogonalität)
• ΨN ⊂ ΨN +1 (Geschachteltheit)
• diam(tr(ψi )) → 0 für i → ∞
Diese Eigenschaften besitzten sogenannte Waveletfunktionen, von denen wir das sogenannte
Haar-Wavelet untersuchen.
Definition 6.30. Wir definieren das “mother wavelet” ψ(x) und die Abschneidefunktion χ(x)
auf (0, 1]:
1 0<x≤1 (
1 0<x≤1
ψ(x) := −1 1<x≤2, χ(x) :=
0 sonst
0 sonst
Für l ∈ N0 (Stufe) und 0 ≤ i < 2l−1 (Index) ist das Haar-Wavelet
l−1
ψil := max(2 2 , 1) · ψ(2l x − 2i) · χ(x).
S2j−1 −1
{ψij }.
Sl
Die Waveletbasis der Stufe l ist Ψl := {ψ00 } ∪ j=1 i=0
Eigenschaften der Wavelets
Wir fassen einige Eigenschaften der Haar-Wavelets zusammen:
• Für l > 0 gilt
i i+1
0
x≤ 2l−1
∨x> 2l−1
l−1 i+ 12
ψil = 2 2 i
2l−1
<x≤ 2l−1
i+ 12
−2 l−1
2 <x≤ i+1
2l−1 2l−1
• Speziell für l = 1 gibt es ein einzelnes Wavelet
ψ01 = max(1, 1) · ψ(2x) · χ(x) = ψ(2x)
• Und speziell für l = 0 gibt es ein einzelnes Wavelet
ψ00 = max(2−1/2 , 1) · ψ(x) · χ(x) = χ(x)
124
6.7 Gauß-Approximation von Funktionen
Graphische Darstellung
0
ψ0
1
-1 0 1
1
ψ0
1
0 1
-1
2 2
1/2 ψ0 ψ1
2
1/2
0 1
-2
3 3 3 3
ψ0 ψ1 ψ2 ψ3
2
0 1
-2
Darstellung der Haar-Wavelets für l = 0, . . . , 3.
Für l = 0 und l = 1 gibt es jeweils nur ein Wavelet (konstante Einsfunktion bzw. gestauchtes
“mother wavelet”). Die Wavelets aller weiteren Level sind Kopien des Wavelets ψ01 , die in x-
Richtung gestaucht und zur Kompensation in ihrem Wert gestreckt wurden.
Weitere Eigenschaften
Weiterhin gilt:
• Die Wavelets bilden eine Baumstruktur bzgl. ihrer Träger:
l+1
tr(ψil ) = tr(ψ2i l+1
) ∪ tr(ψ2i+1 ), tr(ψil ) ∩ tr(ψjl ) = ∅ für i 6= j
• Die Wavelets sind paarweise orthogonal:
(
1 i=j∧l =k
(ψil , ψjk ) =
0 sonst
• Offensichtlich erfüllt die Konstruktion Ψl ⊂ Ψl+1 , da immer nur Basisfunktionen hinzu-
genommen werden.
125
6 Interpolation und Approximation
Damit sind alle geforderten Eigenschaften erfüllt.
Beweis: Tafel
Aufgespannter Raum
Um die durch Haar-Wavelets darstellbaren Funktionen etwas anschaulicher zu machen, definie-
ren wir die Funktionen (
l 1 2il < x ≤ i+1
2l
φi (x) :=
0 sonst
S2l −1
für l ∈ N0 und 0 ≤ i < 2l , sowie Φl := i=0 {φli }.
S l := span(Φl ) ist der Raum der stückweise konstanten Funktionen auf dem Intervall (0, 1]
bzgl. der äquidistanten Unterteilung in 2l Teilintervalle. Es gilt
span(Ψl ) = span(Φl ) = S l .
Beweisskizze: Tafel
Darstellung von Funktionen
Da die Haar-Wavelets eine Orthonormalbasis bilden, ist die Approximationsaufgabe aus Satz 6.27
analog zur Fourier-Reihe durch Skalarprodukte der Form (f, ψij ) berechenbar:
−1
l 2X j
(f, ψij ) · ψij
X
g(x) =
j=0 i=0
Im Unterschied zur Fourier-Reihe lassen sich die gesuchten Koeffizienten aus den Mittelwerten
von f auf den Trägern der ψij berechnen, was eine stark lokalisierte und zugleich durch Rekursion
besonders effiziente Berechnung ermöglicht.
126
6.7 Gauß-Approximation von Funktionen
Datenkompression mit Wavelets
Anwendung von (6.1), der Formel für die Fehlerkontrolle:
Gegeben: f ∈ S l in Basis Ψl
Gesucht: Teilraum Se ⊂ S l , so dass kf − f˜k ≤ TOL ·(f, f ) und Se möglichst klein.
Algorithmus:
• M ← {ψ00 } (Basis von S)
e
2
• ρ ← (f, ψ00 ) (repräsentierter Anteil)
• while (ρ < (1 − TOL) · (f, f )) do:
– wähle Kind ψkl+1 eines ψil ∈ M mit (f, ψkl+1 ) maximal
– M ← M ∪ {ψkl+1 }
2
– ρ ← ρ + (f, ψkl+1 )
Ausblick
Eine leichte Abänderung der Gaußschen Approximationsaufgabe:
Finde u ∈ S (endlich-dimensionaler Funktionenraum), so dass
Z Z
∇u · ∇φ = f φ ∀φ ∈ S
Ω Ω
Dies führt auf die (numerische) Lösung der “Poisson-Gleichung”
n
X ∂2u
−∆u = − =f in Ω ⊂ Rd
i=1
∂x2i
(vgl. das Einführungskapitel über lineare Gleichungssysteme)
127
6 Interpolation und Approximation
6.8 Tschebyscheff-Approximation
Tschebyscheff-Approximation
Zentrales Hilfmittel bei der Gauß-Approximation war die Existenz des Skalarproduktes, das
zur Projektion auf endlichdimensionale Teilräume verwendet wurde. Man kann jedoch auch für
allgemeine Normen ohne zugehöriges Skalarprodukt eine solche Bestapproximation einführen.
Satz 6.31 (Allgemeine Tschebyscheff-Approximation). Sei E ein normierter Vektorraum mit
Norm k · k und S ⊂ E ein endlichdimensionaler Teilraum. Dann gibt es zu jedem f ∈ E eine
(nicht unbedingt eindeutige) beste Approximation g ∈ S:
kf − gk = min kf − φk
φ∈S
Beweis: Tafel
Wir haben für die Polynominterpolation auf [a, b] die Fehlerformel
n
f (n+1) (ξx ) Y
f (x) − p(x) = l(x), l(x) := (x − xj )
(n + 1)!
j=0
hergeleitet. Eine “optimale” Wahl der Stützstellen wäre jene, die
max |l(x)| = klk∞
a≤x≤b
minimiert, den Teil der Formel, der durch die Stützstellen beeinflussbar ist und nicht von f
abhängt.
(“Optimal” mit Anführungszeichen, da das implizit die Position ξx beeinflusst!)
Es gilt l(x) = xn+1 − p(x), p ∈ Pn , daher lässt sich die Minimierung von klk∞ als Tschebys-
cheffsche Approximationsaufgabe auffassen:
Finde p ∈ Pn : max kxn+1 − p(x)k → min
a≤x≤b
Dies entspricht dann einer entsprechenden Wahl der Stützstellen {x0 , . . . , xn }, die zu diesem
Polynom p ∈ Pn führt. Man kann zeigen, dass die Tschebyscheffsche Approximationsaufgabe
in diesem Fall sogar eine eindeutige Lösung hat.
128
6.8 Tschebyscheff-Approximation
Optimale Stützstellenwahl
Satz 6.32 (Optimale Stützstellen). Auf dem Intervall [a, b] = [−1, 1] ist die Tschebyscheff-
Approximation g ∈ Pn zu f (x) = xn+1 bzgl. k · k∞ gegeben durch
g(x) = xn+1 − 2−n Tn+1 (x)
mit dem (n + 1)-ten Tschebyscheff-Polynom
Tn+1 (x) = cos((n + 1) arccos(x)) ∈ Pn+1 .
Dass dies tatsächlich ein Polynom (n+1)-ten Grades in x ist, kann man über Winkelidentitäten
zeigen.
Die Nullstellen
π 2k + 1
xk = cos , k = 0, . . . , n
2 n+1
von Tn+1 sind die “optimalen” Stützstellen auf [−1, 1]. Für andere Intervalle [a, b] transformiert
man diese affin-linear mittels
1
x̃k = [(b − a) · xk + (a + b)].
2
Für das Interpolationspolynom p ∈ Pn einer Funktion f ∈ C n+1 [a, b] erhält man dann die
Fehlerabschätzung
1 (b − a)n+1 (n+1)
kf − pn k∞ ≤ kf k∞ .
(n + 1)! 22n+1
Beweis: [Rannacher]
Runges Funktion
Wir kehren zu Runges Beispielfunktion zurück. Die Stützstellen sind die Projektion äquidistan-
ter Punkte auf dem Einheitskreis:
Quelle: Stephen G. Johnson, wikipedia.org, cc-by-sa
129
6 Interpolation und Approximation
Anschaulich werden dadurch an den Intervallgrenzen, wo sonst große Über- und Unterschwinger
auftreten würden, zusätzliche Stützstellen zur Kompensation platziert.
1
0.8
0.6
0.4
0.2
0
-0.2
-0.4
-1 -0.5 0 0.5 1
−1
Polynominterpolation von Runges Funktion f (x) = (1 + 25x2 ) , diesmal mit “optimalen”
Stützstellen (wieder 3, 5, 9, 17 bzw. 33 Wertepaare). Das letzte Polynom verdeckt dabei die
eigentliche Funktion. Offensichtlich deutlich bessere Resultate!
130
7 Numerische Integration
7.1 Einführung
Numerische Integration (Quadratur)
Wir behandeln die numerische Berechnung bestimmter Integrale in einer Raumdimension:
Z b
I[a,b] (f ) = f (x) dx
a
Alle von uns behandelten Verfahren führen auf die Form
n
(n)
X
I[a,b] (f ) = wi f (xi ) + Fehler,
i=0
wobei wi ∈ R die Gewichte und xi ∈ R die Stützstellen der Integrationsformel sind.
Gründe für Numerische Integration
Es gibt mehrere Gründe, warum die Anwendung von numerischer Integration (numerischer
Quadratur) sinnvoll bzw. notwendig sein kann:
• Oft ist eine Funktion f nur an endlich vielen Punkten bekannt, z.B. da sie aus Messreihen
stammt. Häufig ist man neben der Darstellung der Funktion, z.B. mittels Least Squares
oder Splines, auch an ihrer Ableitung oder ihrem Integral interessiert.
• Es gibt Funktionen, z.B. f (x) = exp(−x2 ), deren Integral wichtig ist, aber es steht keine
geschlossene Formel für dessen Berechnung zur Verfügung.
7.2 Newton-Cotes-Formeln
Die Newton-Cotes-Formeln sind sogenannte interpolatorische Quadraturformeln.
Idee: Stelle das Interpolationspolynom p zu gewissen Stützstellen auf und berechne das Integral
von p exakt.
Formal: Stützstellen und Werte (xi , f (xi )), i = 0, . . . , n, Lagrange-Darstellung:
n Y x − xj
(n) (n)
X
pn (x) = f (xi )Li (x), Li (x) =
xi − xj
i=0 j6=i
und somit Z b n Z b
(n) (n)
X
I[a,b] (f ) ≈ I[a,b] (f ) = pn (x) dx = f (xi ) Li (x) dx
a i=0 a
131
7 Numerische Integration
Ordnung einer Quadraturformel
Definition 7.1 (Ordnung einer Quadraturformel). Eine Quadraturformel I (n) (f ) hat mindes-
tens die Ordnung m, wenn sie Polynome vom Grad m − 1 exakt integriert.
Eine Formel zweiter Ordnung integriert beispielsweise lineare Funktionen exakt.
Die Newton-Cotes-Formeln benutzen Polynominterpolation und haben daher für n + 1 Stütz-
stellen mindestens Ordnung n + 1. Wir werden aber später sehen, dass das noch besser geht.
Newton-Cotes-Formeln
Die Newton-Cotes-Formeln nutzen äquidistante Stützstellen. Es gibt davon zwei Varianten:
Abgeschlossene Formeln:
Die Intervallgrenzen a und b sind Stützstellen, es gilt
b−a
xi = a + iH, i = 0, . . . , n, mit H =
n
Offene Formeln:
Die Grenzen a und b sind keine Stützstellen, und es gilt
b−a
xi = a + (i + 1)H, i = 0, . . . , n, mit H =
n+2
Generell muss für die Summe der Gewichte einer Quadraturformel
X n Xn Z b
wi = wi · 1 = 1 dx = b − a
i=0 i=0 a
gelten. Um Gewichte zu erhalten, die unabhängig vom gewählten Intervall sind, wird dieser
Term ausgeklammert, und man erhält für die geschlossenen Newton-Cotes-Formeln
n
(n)
X
I[a,b] (f ) = (b − a) · w̃i f (xi )
i=0
mit
b nY
s−j
Z Z
(n)
w̃i = (b − a)−1 Li dx = n−1 ds
a 0 i−j
j6=i
Beweis: Tafel
132
7.2 Newton-Cotes-Formeln
Beispiele
Abgeschlossene Formeln für n = 1, 2, 3 und H = (b − a)/n sind:
Die Trapezregel bzw. Sehnen-Trapezregel:
b−a
I (1) (f ) = · [f (a) + f (b)]
2
Die Simpsonregel bzw. Keplersche Fassregel:
(2) b−a a+b
I (f ) = · [f (a) + 4f + f (b)]
6 2
Die 3/8-Regel bzw. Pulcherrima:
b−a
I (3) (f ) = · [f (a) + 3f (a + H) + 3f (b − H) + f (b)]
8
Offene Formeln für n = 0, 1, 2 und H = (b − a)/(n + 2) sind:
Die Mittelpunktregel, Tangenten-Trapezregel bzw. Rechteckregel:
(0) a+b
I (f ) = (b − a) · f
2
Die zweite offene Regel (kein besonderer Name):
b−a
I (1) (f ) = · [f (a + H) + f (b − H)]
2
Die dritte offene Regel (ebenfalls kein Name):
(2) b−a a+b
I (f ) = · [2f (a + H) − f + 2f (b − H)]
3 2
133
7 Numerische Integration
Anmerkungen
Bemerkung 7.2. Ab n = 7 für die abgeschlossenen bzw. n = 2 für die offenen Formeln treten
negative Gewichte auf. Das ist ungünstig, denn:
• Für strikt nicht-negative Funktionen f kann I (n) (f ) < 0 sein (Stoffkonzentration, Masse-
nerhaltung,. . . ).
• Es herrscht erhöhte Gefahr der Auslöschung.
• Die Kondition kann schlechter werden, während sie für rein positive Gewichte beschränkt
ist.
Restgliedabschätzung
Satz 7.3 (Restglieder). Den begangenen Fehler kann man folgendermaßen abschätzen:
1. Trapezregel: n = 1, Ordnung 2, Es gilt
b−a (b − a)3 00
I(f ) − · [f (a) + f (b)] = − f (ξ), ξ ∈ [a, b]
2 12
für f ∈ C 2 ([a, b]). Es werden also Polynome bis zum Grad 1 exakt integriert, denn für
diese ist f 00 (x) = 0 auf [a, b].
Allgemein gilt: die Ordnung der ungeraden Formeln entspricht der Anzahl an Stützstellen,
während die Ordnung der geraden Formeln um Eins höher ist.
2. Simpsonregel: n = 2, Ordnung 4, für f ∈ C 4 ([a, b]) gilt
(b − a)5 (4)
b−a a+b
I(f ) − · [f (a) + 4f + f (b)] = − f (ξ)
6 2 2880
3. Mittelpunktregel: n = 0, Ordnung 2, für f ∈ C 2 ([a, b]) gilt
(b − a)3 00
a+b
I(f ) − (b − a) · f = f (ξ)
2 24
also “halber Fehler der Trapezregel” bei nur einer Auswertung!
Beweis: Tafel
134
7.3 Summierte Quadraturformeln
7.3 Summierte Quadraturformeln
Das Erhöhen des Polynomgrads ist wenig sinnvoll, da
• früh negative Gewichte auftreten
• die Lagrange-Interpolation mit äquidistanten Stützstellen nicht punktweise konvergiert
• man entsprechende Differenzierbarkeit von f für die Abschätzung benötigt
Idee der summierten Quadraturformeln:
• Unterteile das Interval [a, b] in N Teilintervalle
b−a
[xi , xi+1 ], xi = a + ih, i = 0, . . . , N − 1, h =
N
• Wende eine der obigen Formeln auf jedem Teilintervall an und summiere die Ergebnisse.
Beispiele
Beispielsweise ergeben sich bei N Teilintervallen der Schrittweite h:
Die summierte Trapezregel:
N −1 N −1
" #
(1)
X h f (a) X f (b)
Ih (f ) = [f (xi ) + f (xi+1 )] = h · + f (xi ) +
2 2 2
i=0 i=1
Die summierte Simpsonregel:
N −1 N −1
" #
(2) h f (a) X X xi + xi+1 f (b)
Ih (f ) = · + f (xi ) + 2 f +
3 2 2 2
i=1 i=0
Die summierte Mittelpunktregel:
N −1
(0)
X xi + xi+1
Ih (f ) =h· f
2
i=0
135
7 Numerische Integration
Restgliedabschätzung
Satz 7.4 (Restglieder für summierte Quadraturen). Für die je Teilintervall verwendete Qua-
draturformel gelte
(n)
I[xi ,xi+1 ] (f ) − I[xi ,xi+1 ] (f ) = αn hm+1 f (m) (ξi ), ξi ∈ [xi , xi+1 ].
Dann gilt für die summierte Quadraturformel
(n)
I[a,b] (f ) − Ih (f ) = αn (b − a)hm f (m) (ξ), ξ ∈ [a, b].
Beweis: Tafel
Wie wir sehen, bleibt bei der Summation die Ordnung der Formel erhalten (da in der Abschät-
zung die gleiche Ableitung auftritt).
Konkret bedeutet das:
• Die Trapezregel und Mittelpunktregel haben Ordnung 2, also haben es die summierten
Regeln ebenfalls. Es gilt
(1 bzw. 0) b−a
|I[a,b] (f ) − Ih (f )| = | max f (2) | · h2 ∈ O(h2 ),
α [a,b]
mit α = 12 bzw. 24, falls die zweite Ableitung beschränkt ist.
• Die Simpsonregel hat Ordnung 4, also gilt das auch für die summierte Simpsonregel, und
(2) b−a
|I[a,b] (f ) − Ih (f )| = | max f (4) | · h4 ∈ O(h4 ),
2880 [a,b]
falls die vierte Ableitung beschränkt ist.
Fehlerkontrolle
Es gelten
(2) 1 (1) 2 (0)
Ih (f ) = Ih (f ) + Ih (f )
3 3
sowie für die halbierte Schrittweite h/2
(1) 1 (1) 1 (0)
Ih/2 (f ) = Ih (f ) + Ih (f ).
2 2
Diese Zusammenhänge kann man zur Fehlerkontrolle gemäß der Form
|I(f ) − Ih (f )| ≤ TOL
verwenden.
136
7.4 Quadraturen höherer Ordnung
7.4 Quadraturen höherer Ordnung
Wie wir gesehen haben, sind die Newton-Cotes-Formeln nur bis n = 7 bzw. n = 2 brauchbar.
Die summierten Quadraturformeln erlauben zwar beliebig feine Berechnungen, haben aber im-
mer die selbe Ordnung wie die zugrundeliegenden Formeln.
Wie erreicht man eine höhere Ordnung?
Zwei Ansätze:
• Extrapolation auf summierte Formeln anwenden
• Nicht-äquidistante Stützstellen wählen (Gauß-Quadratur)
Romberg-Integration
(n)
Bei den summierten Formeln ist man an limh→0 Ih (f ) interessiert. Dafür ist Extrapolation
zum Limes anwendbar.
Herzstück des Ansatzes ist die Euler-Maclaurinsche Summenformel:
m−1
(1)
X B2k (2k−1)
I(f ) − Ih (f ) = h2k (f (b) − f (2k−1) (a))
(2k)!
k=1
B2m
+ h2m (b − a)f (2m) (ξ), ξ ∈ [a, b]
(2m)!
Speziell die Trapezsumme hat also eine Fehlerentwicklung in geraden Potenzen, was Extra-
polation besonders effizient macht. Die in der Formel auftretenden Konstanten Bi sind die
sogenannten Bernoulli-Zahlen.
Gauß-Integration
Frage: Lässt sich durch nicht-äquidistante Stützstellen die Genauigkeit von Quadraturformeln
verbessern?
Idee: Wähle xi und wi so, dass Polynome von möglichst hohem Grad exakt integriert werden
(Ordnungsmaximierung).
Satz 7.5. Die maximal erreichbare Ordnung einer Quadraturformel mit n + 1 Stützstellen ist
2n + 2 (d.h. Polynome vom Grad 2n + 1 werden exakt integriert).
137
7 Numerische Integration
Beweis: Tafel
Satz 7.6. Es gibt genau eine interpolatorische Quadraturformel zu n + 1 paarweise versch.
Stützstellen in [−1, 1] mit der Ordnung 2n + 2.
Ihre Stützstellen sind die Nullstellen λ0 , . . . , λn ∈ (−1, 1) des (n + 1)-ten Legendre-Polynoms
Ln+1 , wobei
(2k + 1)x k
L0 (x) = 1, L1 (x) = x, Lk+1 (x) = Lk (x) − Lk−1 (x)
k+1 k+1
Die Gewichte erhält man mittels
Z 1 Y x − λj 2
wi = dx > 0, 0, . . . , n
−1 j6=i λi − λj
Beweis: [Rannacher]
Beispiele
Durch Auswerten der Nullstellen der Legendre-Polynome ergeben sich mit
b−a a+b
h= , c=
2 2
zum Beispiel die folgenden Gauß-Quadraturen:
Eine Stützstelle, Ordnung 2:
I (0) = (b − a) · f (c) (Mittelpunktsregel)
Zwei Stützstellen, Ordnung 4:
b−a
I (1) = · [f (c − (1/3)1/2 h) + f (c + (1/3)1/2 h)]
2
Drei Stützstellen, Ordnung 6:
b−a
I (2) = · [5f (c − (3/5)1/2 h) + 8f (c) + 5f (c + (3/5)1/2 h)]
18
Natürlich lassen sich mithilfe der Gauß-Quadraturen wieder summierte Formeln definieren,
indem man das Intervall unterteilt und die Ergebnisse summiert.
138
7.5 Ausblick
7.5 Ausblick
Adaptive Quadratur
Oftmals ist eine äquidistante Aufteilung des Intervals wenig sinnvoll, da die zu integrierende
Funktion lokal unterschiedliches Verhalten zeigt.
Idee der adaptiven Quadratur:
I2
I1 I3 (negativ!)
I0
• Berechne sukzessive Updates für die Näherung durch feinere Unterteilung (Prinzip des
Archimedes)
• Breche Verfeinerung ab, wenn lokales Update |Ij | klein genug
Mehrdimensionale Quadratur
Im Eindimensionalen gilt:
Intervall = zusammenhängende Menge
= wegzusammenhängende Menge
= einfach zusammenhängende Menge
Daher ist es einfach, auf beliebigen “Gebieten” Ω ⊂ R zu integrieren. Im Mehrdimensionalen
werden die Intervalle zu (Hyper-) Quadern, die oben genannten Begriffe sind aber nicht mehr
äquivalent:
Quader 6= zusammenhängende Menge
6= wegzusammenhängende Menge
6= einfach zusammenhängende Menge
139
7 Numerische Integration
Ω
• Die vorgestellten Formeln lassen sich leicht auf Quader verallgemeinern.
• Im allgemeinen muss man diese Quader und die Integrale darauf transformieren (Trans-
formationssatz, Verallgemeinerung der Substitutionsregel).
• Man kann auch Dreiecke für die Unterteilung verwenden.
• Bei komplizierten Geometrien entsteht zusätzlich ein Geometriefehler durch Approxima-
tion der Geometrie.
Fluch der Dimension
Wenn die Dimension d sehr groß ist, sind die hier behandelten Methoden nicht brauchbar.
Betrachte Ω = [0, 1]d . Zerlegt man [0, 1] in zwei Teilintervalle je Dimension, so hat man den
d-dimensionalen Würfel in 2d Teilwürfel zerlegt. Der Aufwand steigt somit exponentiell in d
an. Dies bezeichnet man als “Fluch der Dimension”.
Eine Möglichkeit ist dann die Monte-Carlo-Integration
N
C X
I(f ) ≈ f (ξi )
N
i=1
mit zufälligen Punkten ξi ∈ Ω.
140
8 Iterative Lösung von Gleichungssystemen
8.1 Einführung
In diesem Abschnitt betrachten wir die Lösung von algebraischen Gleichungen
f (x) = 0 mit f : Rn → Rn .
Dabei beschränken wir uns zunächst auf den Fall n = 1 (skalar).
Idee der Intervallschachtelung:
Angenommen man kennt ein Teilintervall [a0 , b0 ], so dass f (a0 ) · f (b0 ) < 0 (Vorzeichenwechsel),
und f sei stetig. Dann hat f nach dem Zwischenwertsatz mindestens eine Nullstelle in [a0 , b0 ].
Bilde Folge von Teilintervallen durch Bisektion und Prüfen der Bedingung.
Algorithmus: Intervallschachtelung
Gegeben: I0 = [a0 , b0 ] mit f (a0 ) · f (b0 ) < 0 und > 0
Algorithmus:
t = 0;
while (bt − at > ) :
xt = (at + bt )/2;
if (f (xt ) == 0) :
at+1 = xt − ; bt+1 = xt ;
else if (f (at ) · f (xt ) < 0) :
at+1 = at ; bt+1 = xt ;
else :
at+1 = xt ; bt+1 = bt ;
t = t + 1;
Idee:
• Halbiere das aktuelle Intervall
• Falls zufällig in der Mitte eine Nullstelle liegt, beende das Verfahren
• Falls die Bedingung für das linke Intervall erfüllt ist, wähle es als nächstes Intervall
• Andernfalls muss die Nullstelle im rechten Intervall sein, setze die Suche dort fort
141
8 Iterative Lösung von Gleichungssystemen
Analyse
Es gilt
at ≤ at+1 < bt+1 ≤ bt
und (sofern nicht f (xt ) = 0)
t+1
1 1
|bt+1 − at+1 | ≤ |bt − at | = |b0 − a0 |
2 2
Anmerkungen:
• Die Konvergenz ist linear mit Rate 1/2
• Gut für monotone Funktionen geeignet (globale Konvergenz)
• Nur für reelle Funktionen im R1 geeignet
8.2 Newton-Verfahren
Die Funktion f sei (mindestens) einmal stetig differenzierbar.
Für gegebenes xt gibt es die “Tangente”
Tt (x) = f 0 (xt )(x − xt ) + f (xt )
mit der Nullstelle
f (xt )
Tt (x) = 0 ⇐⇒ x = xt − .
f 0 (xt )
Das führt auf die Iterationsvorschrift
f (xt )
xt+1 = xt − .
f 0 (xt )
xt x t+1
Offensichtlich ist |f 0 (xt )| > 0 erforderlich, d.h. wir setzen voraus, dass die Nullstelle einfach
ist.
142
8.2 Newton-Verfahren
Newton-Verfahren im Mehrdimensionalen
Das Newton-Verfahren lässt sich auf Systeme f : Rn → Rn erweitern:
Es existiere die Taylorentwicklung von f :
n
X ∂fi
fi (x) = fi (xt + ∆x) = fi (xt ) + (xt )∆xi + Ri (xt , ∆x) i = 1, . . . , n
∂xj
j=1
oder in vektorieller Schreibweise
f (xt + ∆x) = f (xt ) + J(xt )∆x + R(xt , ∆x)
mit der “Jacobimatrix”
∂fi
[J(xt )]ij = (xt )
∂xj
Das Ignorieren des Restglieds entspricht einer “Linearisierung von f ”.
Finde näherungsweise eine Nullstelle von f :
f (x) ≈ f (xt ) + J(xt )∆x = 0 ⇐⇒ ∆x = −J −1 (xt )f (xt )
Dies führt auf die Iteration
xt+1 = xt − J −1 (xt )f (xt )
Dabei erfordert jeder Schritt die Lösung eines linearen Gleichungssystems mit der Jacobima-
trix!
Konvergenz des Newton-Verfahrens
Satz 8.1 (Newton-Verfahren). Die Funktion f ∈ C 2 ([a, b]) habe in (a, b) (Inneres!) eine Null-
stelle z, und es sei
m := min |f 0 (x)| > 0, M := max |f 00 (x)|.
a≤x≤b a≤x≤b
Es sei ρ > 0 so gewählt, dass
M
q := ρ < 1, Kρ (z) := {x ∈ R | |x − z| ≤ ρ} ⊂ [a, b]
2m
Dann sind für jeden Startwert x0 ∈ Kρ (z) die Newton-Iterierten xt ∈ Kρ (z) definiert und
konvergieren gegen die Nullstelle z.
Dabei gilt die a-priori Fehlerabschätzung
2m (2t )
|xt − z| ≤ q , t∈N
M
143
8 Iterative Lösung von Gleichungssystemen
und die a-posteriori Fehlerabschätzung
M
|xt − xt−1 |2 , t ∈ N.
|xt − z| ≤ m−1 |f (xt )| ≤
2m
(a-priori: nutzt nur die Voraussetzungen, a-posteriori: nutzt auch die bereits berechneten Ite-
rierten)
Beweis: Tafel
Beispiel: Wurzelberechnung
Beispiel 8.2 (Wurzelberechnung mit Newton-Verfahren). Sei a > 0 und n ≥ 1. Löse xn = a,
d.h.
f (x) = xn − a = 0, f 0 (x) = n · xn−1 .
Das führt auf die Iterationsvorschrift
xt+1 = n−1 · [(n − 1) · xt + a · x1−n
t ].
Nach Satz 8.1 konvergiert die Iteration, falls x0 nahe genug an a1/n ist. Hier konvergiert sie
jedoch sogar global, d.h. für alle x0 > 0 (aber nicht unbedingt von Anfang an quadratisch).
Anmerkungen
Bemerkung 8.3. Zur Konvergenz des Newton-Verfahrens:
• Das Newton-Verfahren konvergiert nur lokal, d.h. wenn |x0 − z| ≤ ρ (“Einzugsbereich”).
Dabei ist ρ i.d.R. unbekannt und möglicherweise sehr klein.
• Das Newton-Verfahren konvergiert quadratisch,
|xt − z| ≤ c · |xt−1 − z|2 ,
im Gegensatz zur Intervallschachtelung, die nur linear konvergiert.
• Gedämpftes Newton-Verfahren:
Verbesserung der Konvergenz außerhalb des Einzugsbereichs:
f (xt )
xt+1 = xt − λt
f 0 (xt )
mit der Wahl einer Folge λt ∈ (0, 1] als “Dämpfungsstrategie”.
• Mehrfache Nullstellen:
Falls z eine p-fache Nullstelle ist, mit p > 1, konvergiert das Newton-Verfahren dennoch,
aber nur noch linear. Man kann zeigen, dass die alternative Vorschrift
f (xt )
xt+1 = xt − p ·
f 0 (xt )
in diesem Fall quadratisch konvergiert.
144
8.3 Sukzessive Approximation (Fixpunktiteration)
Sekanten-Methode
Unter Umständen ist die Berechnung der Ableitungen sehr teuer.
Idee der Sekanten-Methode:
Ersetze die Tangente durch eine Sekante
f (xt ) − f (xt−1 )
s(x) = f (xt ) + (x − xt ) ,
xt − xt−1
dies führt auf die Iteration
xt − xt−1
xt+1 = xt − f (xt ) .
f (xt ) − f (xt−1 )
Die Sekanten-Methode konvergiert lokal mit
2m γt
|xt − z| ≤ q , t ∈ N,
M
wobei γt die Fibonacci-Zahlen bezeichnet:
γ0 = γ1 = 1, γt+1 = γt + γt−1
Es gilt γt ≈ 0.723 · (1.618)t (“goldener Schnitt”), also liegt die Konvergenzordnung zwischen 1
und 2. Problematisch ist die bei der Auswertung auftretende Auslöschung.
8.3 Sukzessive Approximation (Fixpunktiteration)
Fixpunktiteration
f (x)
Mit g(x) = x − f 0 (x) hat das Newton-Verfahren die Form xt+1 = g(xt ). Da die Nullstelle z
wegen
f (z) = 0 =⇒ g(z) = z
ein Fixpunkt der Iteration xt+1 = g(xt ) ist, nennt man das auch Fixpunktiteration. Hier un-
tersuchen wir allgemeine Iterationen dieser Art.
Falls die Berechnung von f 0 (x) sehr teuer ist, könnte man z.B. f 0 nur einmal “in der Nähe von
z” auswerten, und wie folgt iterieren:
f (x)
xt+1 = xt −
f 0 (c)
Frage: Wann konvergiert so eine Iteration? Insbesondere wollen wir dabei auch f : Rn → Rn , n >
1, zulassen.
145
8 Iterative Lösung von Gleichungssystemen
Kontraktionen
Antwort gibt der sogenannte “Banachsche Fixpunktsatz”.
Dafür zunächst eine Definition:
Definition 8.4 (Kontraktion). Sei G ⊂ Rn eine nichtleere, abgeschlossene Punktmenge und
g : G → G Lipschitz-stetig mit Konstante q < 1, d.h.
kg(x) − g(y)k ≤ q · kx − yk
mit einer Vektornorm k · k im Rn . Ein solches g nennt man eine Kontraktion auf G.
Sukzessive Approximation
Satz 8.5 (Sukzessive Approximation). Sei g eine Kontraktion auf G ⊂ Rn . Dann existiert
genau ein Fixpunkt z ∈ G von g und für jeden Startpunkt x(0) ∈ G konvergiert die Folge der
Iterierten x(t+1) = g(x(t) ) gegen z.
Es gelten die a-priori und a-posteriori Fehlerabschätzungen
q qt
kx(t) − zk ≤ kx(t) − x(t−1) k ≤ kx(1) − x(0) k.
1−q 1−q
(wir schreiben den Index oben in Klammern, um bei Vektoren unten Platz für den Komponen-
tenindex zu haben.)
Beweis: Tafel
8.4 Iterationsverfahren zur Lösung linearer Gleichungssysteme
Wir kehren zurück zur Lösung von linearen Gleichungssystemen
Ax = b, A ∈ Rn×n regulär, b ∈ Rn .
Definition 8.6 (dünn besetzte Matrizen). Eine Folge von Matrizen {A(n) | n ∈ N} heißt dünn
besetzt, falls
(n) (n)
|{aij | aij 6= 0}| =: nnz(A(n) ) = O(n)
(nnz = “number of non-zeros”).
146
8.4 Iterationsverfahren zur Lösung linearer Gleichungssysteme
Aufgrund von “fill in” ist die Gauß-Elimination für dünn besetzte Matrizen oft schlecht geeignet.
Außerdem: für große n nicht mehr tragbarer Aufwand in O(n3 ).
Es gilt
Lösen von Ax = b ⇐⇒ “Nullstellensuche” f (x) = b − Ax = 0.
Definiere mit (zunächst beliebiger) Matrix C die Iteration
x(t+1) = g(x(t) ) = x(t) + C −1 f (x(t) )
= x(t) + C −1 (b − Ax(t) )
= (I − C −1 A) x(t) + C −1 b
| {z }
=:B
mit “Iterationsmatrix” B.
Für x := A−1 b gilt
g(x) = (I − C −1 A)x + C −1 b = x − (C −1 A) · (A−1 b) + C −1 b = x,
also ist x Fixpunkt von g.
Für die Lipschitz-Konstante der Funktion g gilt
kg(x) − g(y)k = kB(x − y)k ≤ kBk · kx − yk,
d.h. falls kBk < 1 ist (für verträgliche Matrixnorm k · k), ist g Kontraktion auf Rn .
Beispiele für Iterationsverfahren
Zerlege A = L + D + R mit L strikte untere Dreiecksmatrix, D Diagonalmatrix und R strikte
obere Dreiecksmatrix.
Jacobi-Verfahren:
Setze C = D und daher B = −D−1 (L + R), also
x(t+1) = x(t) + D−1 (b − Ax(t) )
Das Jacobi-Verfahren berechnet für die einzelnen linearen Gleichungen jeweils den Defekt, also
die Diskrepanz zwischen linker und rechter Seite, und ändert dann die zugehörige Kompo-
nente von x so, dass die Gleichung erfüllt wäre, wenn die anderen Komponenten unverändert
blieben.
Gauß-Seidel-Verfahren:
147
8 Iterative Lösung von Gleichungssystemen
Setze C = L + D und daher B = −(L + D)−1 R, also
x(t+1) = x(t) + (L + D)−1 (b − Ax(t) ) (Vorwärtseinsetzen)
Das Gauß-Seidel-Verfahren operiert ähnlich, verwendet jedoch die bereits berechneten aktua-
lisierten Werte aus vorhergehenden Zeilen, statt wie beim Jacobi-Verfahren alle Korrekturen
simultan auszurechnen.
Solche Iterationsverfahren konvergieren normalerweise nur für bestimmte Klassen von Matrizen
(da C regulär und kBk < 1 sein muss).
Konvergenz des Jacobi-Verfahrens
Erinnerung:
Eine Matrix heißt strikt diagonaldominant, falls
X
|aij | < |aii | ∀ i = 1, . . . , n.
j6=i
Beispiele: LGS für kubische Splines, Radiosity-Methode.
Satz 8.7. Das Jacobi-Verfahren und das Gauß-Seidel-Verfahren konvergieren für strikt diago-
naldominante Matrizen.
Beweis: Tafel, Gauß-Seidel: [Rannacher]
Aufwand
Es gibt viele ähnliche Aussagen für symmetrisch positiv definite Matrizen, schwach diagonal-
dominante Matrizen, sogenannte M-Matrizen, . . .
Der Aufwand für eine Iteration sei α(n), typischerweise mit α(n) = O(n). Wegen
kx(t) − xk ≤ kBkt kx(0) − xk
log()
sind für eine Reduktion des Fehlers um den Faktor 1 insgesamt t ≥ log(kBk) Iterationen
nötig, also ist der Gesamtaufwand
log()
Tfix (n) = α(n).
log(kBk)
Problem: hohe Kosten wenn kBk nahe Eins, kBk ist problemabhängig und wächst häufig mit
n.
148
8.4 Iterationsverfahren zur Lösung linearer Gleichungssysteme
Relaxationsverfahren
Man nutzt daher Relaxationsverfahren: Für die Wahl Cω = (L + ω −1 D) mit Relaxationspara-
meter ω 6= 0 ergibt sich
Bω = −(ωL + D)−1 ((ω − 1)D + ωR).
Für ω = 1 erhält man wieder das Gauß-Seidel-Verfahren, es stellt also einen Spezialfall dar.
Wie sollte man ω wählen?
Satz 8.8. Es gilt kBω k ≥ |ω − 1|, daher muss zwingend ω ∈ (0, 2) sein.
Beweis: Tafel
Für s.p.d. Matrizen kann man zeigen, dass ω ∈ (0, 2) nicht nur notwendig, sondern auch hinrei-
chend ist. Für solche Matrizen konvergiert also das Relaxationsverfahren, und das Gauß-Seidel-
Verfahren als Spezialfall.
Die Grundidee des SOR-Verfahrens (Successive Over-Relaxation) ist es, ein ω > 1 zu wählen,
das möglichst nah am (in der Praxis unbekannten) optimalen Wert liegt, um ein minimales
kBω k und somit schnellere Konvergenz als bei Jacobi/Gauß-Seidel zu erhalten.
Generell kann man durch Relaxation oft Fixpunktiterationen modifizieren und die Konvergen-
zeigenschaften verbessern. Statt
x(t+1) = g(x(t) ) = x(t) + (g(x(t) ) − x(t) )
verwendet man dabei
x(t+1) = (1 − ω)x(t) + ωg(x(t) )
= x(t) + ω(g(x(t) ) − x(t) )
mit einem Relaxationsparameter ω > 0. Für 0 < ω < 1 erhält man ein robusteres, lang-
sameres Verfahren (vgl. Newton), für ω > 1 eine Beschleunigung der Konvergenz langsamer
Verfahren.
Abstiegsverfahren
Für symmetrisch positiv definite Matrizen gibt es eine weitere Klasse von Verfahren, die übli-
cherweise deutlich effizienter ist, die sogenannten Abstiegsverfahren.
149
8 Iterative Lösung von Gleichungssystemen
Satz 8.9 (Minimierungseigenschaft). Sei A ∈ Rn×n symmetrisch positiv definit. Mit dem qua-
dratischen Funktional
1
Q(y) := (Ay, y) − (b, y)
2
ist die Lösung x des Gleichungssystems Ax = b durch die Eigenschaft
∀ y ∈ Rn , y 6= x : Q(x) < Q(y)
charakterisiert.
Beweis: Tafel
Abstiegsverfahren bestehen aus zwei Komponenten:
• Einer Folge von Suchrichtungen d(t) , die garantieren, dass Q sukzessive kleiner wird, man
setzt
x(t+1) = x(t) + αt · d(t) .
• Zugehörige Schrittweiten αt , die jeweils Q entlang der Suchrichtung d(t) minimieren:
Q(x(t+1) ) = min(x(t) + αd(t) )
α∈R
Die Iterierten x(t) konvergieren dann gegen das Minimum von Q, und daher gegen die Lösung
von Ax = b.
Bestimmung der Schrittweite
Der Gradient von Q ist das Residuum des Gleichungssystems,
1
∇Q(y) = (A + AT )y − b = Ay − b,
2
und es gilt (Kettenregel)
d
Q(x(t) + αd(t) ) = ∇Q(x(t) + αt · d(t) ) · d(t)
dα
= (Ax(t) − b, d(t) ) + α · (Ad(t) , d(t) )
Daher sind die Schrittweiten eindeutig durch die Suchrichtungen festgelegt:
(r(t) , d(t) )
αt = − mit r(t) := ∇Q(x(t) ) = Ax(t) − b
(Ad(t) , d(t) )
150
8.4 Iterationsverfahren zur Lösung linearer Gleichungssysteme
Allgemeines Abstiegsverfahren
Definition 8.10 (Allgemeines Abstiegsverfahren). Es sei x(0) ein vorgegebener Startwert mit
zugehörigem Gradienten r(0) = Ax(0) − b. Ausgehend von einer vorgegebenen oder sich während
des Verfahrens ergebenden Folge von Suchrichtungen d(t) berechnet das Abstiegsverfahren in
jeder Iteration t die folgenden Größen:
(t)
(r ,d ) (t)
• Schrittweite αt = − (Ad (t) ,d(t) )
• Neue Iterierte x(t+1) = x(t) + αt d(t)
• Neuer Gradient r(t+1) = r(t) + αt Ad(t)
Das Verfahren kann abgebrochen werden, sobald r(t) klein genug ist.
Beispiele von Abstiegsverfahren
Koordinatenrelaxation:
Die wohl einfachste Wahl der Suchrichtungen ist
d(t) = e(t mod n)+1
(zyklisch die einzelnen Raumdimensionen), dies ist im wesentlichen äquivalent zum Gauß-Seidel-
Verfahren (warum?).
Bei dieser Wahl kann es passieren, dass in einzelnen Richtungen d(t) keine Verbesserung möglich
ist; nach einem vollen Zyklus ist dies aber garantiert der Fall.
Gradientenverfahren:
Das Gradientenverfahren (Verfahren des steilsten Abstiegs) wählt den negativen Gradienten als
Suchrichtung:
d(t) = −∇Q(x(t) )
Im Gegensatz zur Koordinatenrelaxation ist dies immer eine Abstiegsrichtung.
Satz 8.11 (Konvergenz des Gradientenverfahrens). Sei A ∈ Rn×n symmetrisch positiv defi-
nit. Dann konvergiert das Gradientenverfahren für jeden Startwert x(0) gegen die Lösung des
Gleichungssystems Ax = b.
Beweis: Tafel
151
8 Iterative Lösung von Gleichungssystemen
Konvergenzgeschwindigkeit
Eine genauere Aussage über die Konvergenzgeschwindigkeit des Gradientenverfahrens liefert
der folgende Satz:
Satz 8.12 (Fehlerabschätzung Gradientenverfahren). Für das Gradientenverfahren gilt die Feh-
lerabschätzung
1 − 1/κ t (0)
(t)
kx − xkA ≤ kx − xkA
1 + 1/κ
mit der Spektralkonditionszahl κ = cond2 (A) und der “A-Norm”
kykA := (y, Ay)1/2 .
Beweis: Tafel
Wie wir sehen können, ist die Spektralnorm der Matrix entscheidend für die Konvergenzrate.
Je größer der Unterschied zwischen kleinstem und größtem Eigenwert ist, desto langsamer
konvergiert das Verfahren.
Im Zweidimensionalen kann man dies besonders gut veranschaulichen: sukzessive Richtungen
sind zwar stets orthogonal, was lokal den steilsten Abstieg garantiert, global entsteht jedoch
ein “Zickzack-Kurs”, der sich extrem langsam dem Ziel nähert.
2
1.5 Gradientenverfahren
1 Isolinien
0.5
0
-0.5
-1
-1.5
-2
-10 -8 -6 -4 -2 0 2 4 6 8 10
Beobachtung: Für die Einheitsmatrix konvergiert das Verfahren in einem Schritt. Die Äquipo-
tentialflächen sind (Hyper-) sphären, die um die Lösung zentriert sind, und der Gradient zeigt
stets ins Zentrum.
Können wir Problemstellung und Verfahren so transformieren, dass man dem möglichst nahe
kommt?
Zwei Ansätze:
152
8.4 Iterationsverfahren zur Lösung linearer Gleichungssysteme
• Verwendung eines problemspezifischen Skalarproduktes statt der üblichen Orthogonalität
(führt auf das Verfahren konjugierter Gradienten), und/oder
• explizite Transformation des Vektorraums (Konzept der Vorkonditionierung)
Konjugierte Gradienten
Definition 8.13 (CG-Verfahren). Das CG-Verfahren ist ein Abstiegsverfahren, es setzt also in
jeder Iteration
(r(t) , d(t) )
αt = − , x(t+1) = x(t) + αt d(t) , r(t+1) = r(t) + αt Ad(t) .
(Ad(t) , d(t) )
Man setzt d(0) = r(0) und definiert für t > 0 die einzelnen Suchrichtungen über
(r(t+1) , r(t+1) )
βt = , d(t+1) = −r(t) + βt d(t) .
(r(t) , r(t) )
Das Verfahren der konjugierten Gradienten (CG-Verfahren) benutzt Informationen aus vorher-
gehenden Iterationen, und lässt sich daher nicht als Fixpunktiteration auffassen.
Relevant für die Analyse des Verfahrens sind die Krylow-Räume
Kt (y; A) := span{y, Ay, A2 y, . . . , At−1 y}
bezüglich eines Vektors y, sowie das “A-Skalarprodukt”, dessen zugehörige “A-Norm” wir bereits
kennengelernt haben:
(y1 , y2 )A := (Ay1 , y2 ), kykA := (Ay, y)1/2
Satz 8.14 (Konvergenz des CG-Verfahrens). Der t-te Schritt des CG-Verfahrens minimiert Q
auf dem affinen Unterraum x(0) + Kt (r(0) , A). Die bisherigen Suchrichtungen bilden dabei eine
A-orthogonale Basis des Krylow-Raums Kt (r(0) , A), es gilt
∀i ≤ t, j ≤ t : (d(i) , d(j) )A = 0.
Falls Kt+1 (r(0) , A) = Kt (r(0) , A) gelten sollte, folgt r(t) = 0, daher konvergiert das Verfahren
in exakter Arithmetik nach spätestens n Schritten.
Beweis: [Rannacher]
Während das Gradientenverfahren Orthogonalität immer nur für den nächsten Schritt garan-
tieren kann, baut das CG-Verfahren also systematisch eine A-orthogonale Basis auf.
Wie wir gesehen haben, konvergiert das CG-Verfahren in einer endlichen Anzahl von Schritten,
und zählt somit eigentlich zu den direkten Verfahren. In der Praxis wird es jedoch als iteratives
Verfahren eingesetzt, und zwar aus zwei Gründen:
153
8 Iterative Lösung von Gleichungssystemen
• Bei Fließkommaarithmetik ist die A-Orthogonalität nicht exakt erfüllt, so dass das Ver-
fahren nicht abbricht.
• Bei sehr großen Matrizen erreicht man oft bereits nach vergleichsweise wenigen Schritten
eine brauchbare Näherungslösung.
Satz 8.15 (Fehlerabschätzung CG-Verfahren). Für das CG-Verfahren gilt die Fehlerabschät-
zung
!t
1 − 1/κ1/2
kx(t) − xkA ≤ 2 kx(0) − xkA
1 + 1/κ1/2
mit der Spektralkonditionszahl κ = cond2 (A). Zur Reduktion des Anfangsfehlers um den Faktor
sind höchstens
1
t() ≤ κ1/2 ln(2/) + 1
2
Schritte erforderlich.
Beweis: [Rannacher]
Die Konditionszahl geht im Gegensatz zum Gradientenverfahren also nur über ihre Wurzel ein,
so dass bei hohen Konditionszahlen (breites Eigenspektrum, “schwierige” Probleme) mit einer
deutlich schnelleren Konvergenz zu rechnen ist.
Besonders drastisch sieht dies im Zweidimensionalen aus, da hier das CG-Verfahren nach Kon-
struktion stets in zwei Schritten konvergiert.
2
1.5 Gradientenverfahren
1 CG-Verfahren
0.5 Isolinien
0
-0.5
-1
-1.5
-2
-10 -8 -6 -4 -2 0 2 4 6 8 10
154
8.4 Iterationsverfahren zur Lösung linearer Gleichungssysteme
Ausblick
Die Abstiegsverfahren lassen sich auch zur Minimierung allgemeiner Funktionen verwenden,
denen nicht die Lösung eines linearen Gleichungssystems zugrunde liegt (Nichtlineare Optimie-
rung). Zu beachten ist dabei:
• Es können mehrere lokale Minima vorliegen, die Verfahren konvergieren dann quasi zu-
fällig gegen eines davon, oder auch gar nicht (Sattelpunkte!).
• Die Herleitung von Fehlerabschätzungen ist deutlich schwieriger.
• Bei allgemeinen Funktionen gibt es keine geschlossene Formel für die optimale Schritt-
weite αt . Stattdessen muss man mittels Liniensuche (eindimensionale Optimierung) eine
hinreichend gute Schrittweite bestimmen.
155
9 Schlussbemerkungen
9 Schlussbemerkungen
Top 10 der wichtigsten Algorithmen
Die 10 Algorithmen mit “dem größten Einfluss auf Wissenschaft und Technik” im 20. Jahrhun-
dert, laut Computing in Science and Engineering (2000):
• Metropolis algorithm for Monte Carlo
• Simplex method for linear programming
• Krylov subspace iteration methods
• The decompositional approach to matrix computations
• The Fortran optimizing compiler
• QR algorithm for computing eigenvalues
• Quicksort algorithm for sorting
• Fast Fourier transform
• Integer relation detection
• Fast multipole method
In rot: wurde im Rahmen der Vorlesung behandelt / angesprochen.
Die Top 10 laut Princeton Companion to Applied Mathematics (2015, größerer Fokus auf
angewandte Mathematik):
• Newton and quasi-Newton methods
• Matrix factorizations (LU, Cholesky, QR)
• Singular value decomposition, QR and QZ algorithms
• Monte-Carlo methods
• Fast Fourier transform
• Krylov subspace methods (conjugate gradients, Lanczos, GMRES, minres)
• JPEG
• PageRank
• Simplex algorithm
• Kalman filter
In rot: wurde im Rahmen der Vorlesung behandelt / angesprochen.
156
Weiterführende Veranstaltungen
Auswahl an weiterführenden Veranstaltungen in den Bereichen Numerik und Wissenschaftliches
Rechnen:
• Numerik 1 (MD1)
• Numerische Lineare Algebra (MH5)
• Numerik gewöhnlicher Differentialgleichungen (MH6)
• Numerik partieller Differentialgleichungen (MH7)
• Object-Oriented Prog. for Scientific Computing (IOPWR)
• Paralleles Höchstleistungsrechnen (IPHR)
• Parallele Lösung großer Gleichungssysteme (IPLGG)
157