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
Authors Dr. Ole Klein
License CC-BY-SA-4.0