Einführung in die Numerik Axel Dreves Institut für Mathematik und Rechneranwendung Fakultät für Luft- und Raumfahrttechnik Universität der Bundeswehr München Werner-Heisenberg-Weg 39 85577 Neubiberg/München E-Mail: [email protected] Version: 19. März 2015 Vorwort Dies ist ein Skript zur Vorlesung Einführung in die Numerik“ im Bachelorstudiengang ” Mathematical Engineering an der Universität der Bundeswehr München. Wesentliche Teile des Skripts sind aus dem gleichnamigen Skript von Prof. Dr. Matthias Gerdts von der Universität der Bundeswehr München, sowie dem Skript zu der Vorlesung Numerische ” Mathematik“ von Prof. Dr. Christian Kanzow an der Universität Würzburg entstanden. Empfohlene Literatur: • Deuflhard, P. und Hohmann, A. Numerische Mathematik . de Gruyter, Berlin, 3. Auflage 2002. • Hämmerlin, G. und Hoffmann, K.-H. Numerische Mathematik . Springer, BerlinHeidelberg-New York, 4. Auflage, 1994. • Hanke-Bourgois, M. Grundlagen der Numerischen Mathematik und des Wissenschaftlichen Rechnens. Teubner, 2002. • Kincaid, D. und Cheney, W. Numerical Analysis: Mathematics of Scientific Computing. Brooks/Cole Series in Advanced Mathematics, Thomson Learning, 3rd Edition, 2002. • Opfer, G. Numerische Mathematik für Anfänger . Vieweg, 5. Auflage, 2008. • Plato, R. Numerische Mathematik kompakt. Vieweg, 2000. • Schaback, R. und Werner, H. Numerische Mathematik . Springer, Berlin-HeidelbergNew York, 4. Auflage, 1992. • Stoer, J. Numerische Mathematik 1 . Springer, Berlin-Heidelberg-New York, 9. Auflage, 2004. i Inhaltsverzeichnis 1 Einleitung 1 1.1 Was ist Numerische Mathematik? . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Zahldarstellung und Gleitpunktarithmetik . . . . . . . . . . . . . . . . . . 2 1.3 Rundung und Auslöschung . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.4 Taylorentwicklung und Landau-Symbole . . . . . . . . . . . . . . . . . . . 10 2 Lineare Gleichungssysteme 12 2.1 Grundlagen und Lösbarkeit . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.2 Gauß’sches Eliminationsverfahren und LR-Zerlegung . . . . . . . . . . . . 17 2.2.1 Einfach lösbare Gleichungssysteme . . . . . . . . . . . . . . . . . . 18 2.2.2 Das Gauß’sche Eliminationsverfahren . . . . . . . . . . . . . . . . . 19 2.2.3 Pivotisierung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.2.4 Die LR-Zerlegung einer Matrix . . . . . . . . . . . . . . . . . . . . 25 2.2.5 Pivotisierung in Matrixnotation . . . . . . . . . . . . . . . . . . . . 29 2.3 Cholesky-Zerlegung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 2.4 Fehlereinfluss und Kondition . . . . . . . . . . . . . . . . . . . . . . . . . . 33 2.5 Iterative Lösung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 2.5.1 2.6 Spezielle Verfahren . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 Das konjugierte Gradientenverfahren . . . . . . . . . . . . . . . . . . . . . 48 3 Lineare Ausgleichsprobleme 54 3.1 Gauß’sche Normalengleichungen . . . . . . . . . . . . . . . . . . . . . . . . 55 3.2 QR-Zerlegung einer Matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 3.2.1 QR-Zerlegung zur Lösung von linearen Gleichungssystemen . . . . . 58 3.2.2 Lösung des linearen Ausgleichsproblems mittels QR-Zerlegung . . . 59 3.2.3 QR-Zerlegung nach Householder . . . . . . . . . . . . . . . . . . . . 60 4 Nichtlineare Gleichungen 63 4.1 Bisektion (Intervallschachtelungsverfahren) . . . . . . . . . . . . . . . . . . 65 4.2 Fixpunktiteration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 4.3 Newton-Verfahren . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 4.4 Sekantenverfahren . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 ii 5 Interpolation 5.1 Polynominterpolation . . . . . . . . . . . . . . . . . . . . . 5.1.1 Newton’sche dividierte Differenzen . . . . . . . . . 5.1.2 Fehlerdarstellung interpolierender Polynome . . . . 5.2 Splineinterpolation . . . . . . . . . . . . . . . . . . . . . . 5.2.1 Fehlerdarstellung interpolierender kubischer Splines . . . . . . . . . . . . . . . 102 . 105 . 109 . 112 . 118 . 120 7 Eigenwert-Probleme 7.1 Matrixzerlegungen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2 Potenz-Methoden . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.3 QR-Algorithmus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 . 124 . 128 . 134 6 Numerische Integration 6.1 Interpolatorische Quadraturformeln 6.2 Extrapolationsverfahren . . . . . . 6.3 Gauß’sche Quadraturformeln . . . . 6.4 Adaptive Verfahren am Beispiel der 6.5 Mehrdimensionale Integrale . . . . Literaturverzeichnis . . . . . . . . . . . . . . . . . . . . . . . . . . . Simpson-Regel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 84 86 91 94 101 . . . . . 140 iii Kapitel 1 Einleitung 1.1 Was ist Numerische Mathematik? Die Numerische Mathematik (kurz: Numerik) beschäftigt sich mit der Entwicklung, Analyse und Implementierung von Algorithmen zur Lösung von mathematischen Problemen. Die Problemstellungen sind dabei häufig durch konkrete Anwendungen motiviert. Ein (numerischer) Algorithmus ist dabei das Werkzeug, mit dem ein mathematisches Problem, z.B. das Lösen eines linearen Gleichungssystems, die Approximation einer Funktion, die Berechnung von Eigenwerten, etc., gelöst werden kann. Natürlich kann es verschiedene Algorithmen zur Lösung des gleichen Problems geben, die in bestimmten Situationen jeweils ihre Vor- und Nachteile haben. Die meisten Probleme, die in der Numerik behandelt werden, können entweder nicht analytisch gelöst werden (d.h. es gibt keine geschlossene Lösungsdarstellung) und erfordern daher approximative Lösungsverfahren, oder sie sind zu rechenintensiv, um Lösungen per Hand ausrechnen zu können. Beispiel 1.1.1 Für die folgenden Probleme kann die Lösung nicht explizit angegeben werden: • Die nichtlinearen Gleichungen x − cos(x) = 0, x5 − 3x4 − πx3 + x2 − 43x + exp(x) + x − 2 = 0, können nicht geschlossen gelöst werden. • Die Werte der Integrale Z 5 exp(x) dx, x 1 Z 2 1 sin(x) dx, x Z 1 0 können nicht analytisch berechnet werden. • Die nichtlineare Diffentialgleichung (Pendelgleichung) x00 (t) = − kann nicht geschlossen gelöst werden. 1 3g sin(x(t)) 2` 1 √ dx 3 x +1 √ 2=0 Es gibt sehr viele solcher Beispiele. Die meisten der in der Praxis auftretenden Beispiele sind nicht geschlossen lösbar. Das Gebiet der Numerik bekam insbesondere durch die fortschreitende Entwicklung von leistungsfähigen Computern eine immer stärkere Bedeutung, da diese die Lösung immer komplexerer Aufgabenstellungen erlauben, z.B. • Wettervorhersage mittels numerischer Berechnungsverfahren • Simulation und Optimierung von technischen Systemen (Crashtests, Strömungssimulation, Baustatik, . . . ) • Signal- und Bildverarbeitung (jpeg) Im Rahmen der Vorlesung werden verschiedene mathematische Problemstellungen und Lösungsalgorithmen behandelt. Dabei geht es nicht nur um die Herleitung und Formulierung von Algorithmen, sondern insbesondere auch um deren Eigenschaften. Zentrale Begriffe und Fragestellungen in der Numerik sind die folgenden: • Stabilität: Wie reagiert ein Algorithmus auf Rundungsfehler, Störungen oder Parameteränderungen? • Konvergenz und Konvergenzgeschwindigkeit: Konvergiert ein Algorithmus gegen die Lösung eines Problems und wenn ja, wie schnell? • Approximationsfehler und Approximationsordnung: Wie gut approximiert ein Algorithmus die Lösung in Abhängigkeit von Diskretisierungsparametern? Darüber hinaus besteht ein wesentlicher Teil der Numerischen Mathematik in der konkreten Implementierung von Algorithmen auf dem Computer. 1.2 Zahldarstellung und Gleitpunktarithmetik Ein Computer hat nur eine begrenzte Anzahl von Speicherplätzen zur Verfügung. Die Anzahl der auf dem Computer darstellbaren Zahlen ist damit begrenzt auf solche Zahlen, die sich mit endlich vielen Ziffern vor oder nach dem Komma beschreiben lassen. Mathematisch entspricht dies einer Teilmenge der rationalen Zahlen Q. Man nennt diese Zahlen auch Maschinenzahlen. Die reelle Zahl 12345.678 kann geschrieben werden als 12345.678 = 1 · 104 + 2 · 103 + 3 · 102 + 4 · 101 + 5 · 100 + 6 · 10−1 + 7 · 10−2 + 8 · 10−3 . 2 Allgemeiner kann eine Zahl (im Dezimalsystem) geschrieben werden als an an−1 · · · a1 a0 . b1 b2 b3 · · · = | {z } | {z } ganzzahliger Anteil Bruchanteil n X k ak · 10 + |k=0 {z } ganzzahliger Anteil ∞ X bk · 10−k . } |k=1 {z Bruchanteil Tatsächlich benötigen wir mitunter einen unendlichen String von Ziffern rechts des Dezimalpunkts, um eine beliebige reelle Zahl darstellen zu können, z.B. 1 = 0.3333333 . . . , 3 π = 3.1415926 . . . , e = 2.7182818 . . . , √ 2 = 1.4142135 . . . . Die obigen Zahlen sind im Dezimalsystem zur Basis B = 10 dargestellt. Jedes andere B ∈ N, B > 1, ist ebenfalls zulässig und man erhält die Darstellung (an an−1 · · · a1 a0 .b1 b2 b3 · · · )B = n X k=0 k ak · B + ∞ X bk · B −k , k=1 wobei ai , bj ∈ {0, 1, . . . , B − 1}. Im alltäglichen Leben sind wir es gewohnt im Dezimalsystem mit B = 10 zu rechnen. Bei der internen Darstellung im Computer wird hingegen zumeist eines der folgenden Systeme bevorzugt: • Binärsystem (B = 2). Die Binärzahl 10010.1001 entspricht der Dezimalzahl 18.5625, denn (10010.1001)2 = 1·24 +0·23 +0·22 +1·21 +0·20 +1·2−1 +0·2−2 +0·2−3 +1·2−4 = 18.5625. • Oktalsystem (B = 8). • Hexadezimalsystem (B = 16). Das Hexadezimalsystem verwendet die Zahlen“ ai , bj ∈ {0, 1, . . . , 9, A, B, C, D, E, F }. ” Es kann ebenfalls passieren, dass die Darstellung einer Zahl in einem Zahlensystem endlich ist, also nur endlich viele Ziffern benötigt, während sie unendliche viele Ziffern zur Darstellung in einem anderen Zahlensystem benötigt, z.B. (0.1)10 = (0.00011001100110011 . . .)2 . 3 Soll die reelle Zahl x im Computer realisiert werden, muss auf Grund der begrenzten Speicherkapazität die Anzahl der Ziffern zur Darstellung der Mantisse und des Exponenten auf eine gewisse Länge beschränkt werden, etwa in der Form x = m · Be, (1.1) m = ±(m1 · B −1 + m2 · B −2 + . . . + m` · B −` ) ∈ R (1.2) wobei die Mantisse, B ∈ N, B > 1 die Basis und e = ±(e1 · B n−1 + e2 · B n−2 + . . . + en−1 · B + en ) ∈ Z (1.3) den Exponenten bezeichnen. Hierbei sind ` und n die jeweilige Länge der Mantisse bzw. des Exponenten. Für die Ziffern mi , i = 1, . . . , `, und ei , i = 1, . . . , n, gilt die Zusatzforderung mi , ei ∈ {0, 1, 2, . . . , B − 1}. Damit ist schon eine gewisse Normierung vorgenommen worden. Da 0 ≤ m1 ≤ B − 1 gefordert wurde, ist |m| wegen ` ` ` X X B−1 X 1 1 1 −i |m| = mi · B ≤ = − <1 = 1 − Bi B i−1 B i B` i=1 i=1 i=1 stets kleiner als 1. Definition 1.2.1 (Gleitpunktdarstellung) Die Darstellung einer Zahl x in der Form x = ±0.m1 m2 . . . m` · B ±e1 e2 ...en (1.4) heißt Gleitpunktdarstellung von x. x heißt dann auch Gleitpunktzahl. Die Gleitpunktdarstellung der Zahl x ist jedoch noch nicht eindeutig, wie das folgende Beispiel zeigt: 12345.678 = 0.12345678 · 105 = 0.012345678 · 106 = 0.0012345678 · 107 .. . Der Dezimalpunkt verschiebt sich durch Multiplikation mit entsprechenden Zehnerpotenzen. Um Zahlen auf dem Computer effizient verwalten zu können, benötigt man eindeutige Zahldarstellungen. 4 Eine normierte Gleitpunktdarstellung von x entsteht durch die Zusatzforderung, dass für x 6= 0 die erste Ziffer m1 der Mantisse m ungleich 0 sei. Die normierte Gleitpunktdarstellung von x ist eindeutig. Definition 1.2.2 (Normierte Gleitpunktdarstellung) Die normierte Gleitpunktdarstellung der reellen Zahl x 6= 0 ist definiert als x = ±0.m1 m2 m3 · · · m` . . . · B e = ±m · B e wobei m1 6= 0 gelte. Die Zahl e = ±e1 e2 . . . en ∈ Z heißt Exponent, die Zahl m heißt normierte Mantisse. Die Menge aller normierten Gleitpunktdarstellungen liefert die auf dem Computer darstellbaren Zahlen, die sogenannten Maschinenzahlen. Definition 1.2.3 (Maschinenzahlen) Die Menge ±e1 e2 ...en MB , m1 6= 0} ∪ {0}, `,n := {x ∈ R | x = ±0.m1 m2 . . . m` · B wobei die Länge des Bruchanteils der Mantisse auf ` ∈ N Ziffern und die Länge des Exponenten auf n ∈ N Ziffern beschränkt sind, heißt Menge der Maschinenzahlen zur Basis B (kurz: Maschinenzahlen). Auf Grund der beschränkten Länge des Exponenten gibt es einen größtmöglichen Exponenten emax und einen kleinstmöglichen Exponenten emin und e ∈ [emin , emax ]. Die Zahl µ = B emin −1 ist die kleinste positive Maschinenzahl und ν = B emax · (1 − B −` ) ist die größte Maschinenzahl. Die Maschinenzahlen sind nicht gleichmäßig verteilt, wie das folgende Beispiel zeigt. Beispiel 1.2.4 (Verteilung der Maschinenzahlen) Betrachte die Maschinenzahlen MB `,n mit ` = 3, n = 2 und Basis B = 2. Die folgende Tabelle listet die verfügbaren Mantissen und Exponenten auf: Mantissen (0.100)2 (0.101)2 (0.110)2 (0.111)2 = = = = Exponenten 1 2 5 8 3 4 7 8 (00)2 = 0 ±(01)2 = ±1 ±(10)2 = ±2 ±(11)2 = ±3 Der größtmögliche Exponent ist emax = (11)2 = 21 +20 = 3, der kleinstmögliche Exponent 1 ist emin = −(11)2 = −3. Die kleinste positive Maschinenzahl ist µ = 2−4 = 16 , die größte 3 −3 Maschinenzahl ist ν = 2 (1 − 2 ) = 7. Die Verteilung der positiven Maschinenzahlen sieht wie folgt aus: 5 0 −2 2 1.3 −1 2 0 2 1 2 2 2 7 Rundung und Auslöschung Rundungsfehler treten zwangsläufig auf, da z.B. 1 = 0.3333333 . . . , 3 π = 3.1415926 . . . , √ e = 2.7182818 . . . , 2 = 1.4142135 . . . , keine Maschinenzahlen sind und somit nicht korrekt auf einem Computer dargestellt werden können. Die Umwandlung einer gegebenen reellen Zahl x in eine Maschinenzahl wird Rundung genannt. Ist x bereits eine Maschinenzahl, so soll sie invariant bzgl. Rundung sein. Mathematisch ist die Rundung eine Abbildung rd : R → MB `,n , x 7→ rd(x), (1.5) die eine reelle Zahl x = ±0.m1 m2 . . . mk mk+1 . . . m` m`+1 . . . · B e gemäß der üblichen Vorschrift ±0.m1 m2 . . . m` · B e , ±0.m1 m2 . . . (mk + 1)0 . . . 0 · B e , rd(x) = ±0.10 . . . 0 · B e+1 , falls 0 ≤ m`+1 < B/2, falls B/2 ≤ m`+1 ≤ B − 1, mk 6= B − 1, mk+1 = . . . = m` = B − 1, falls B/2 ≤ m`+1 ≤ B − 1, m1 = . . . = m` = B − 1 auf die nächstgelegene Maschinenzahl rundet. Die bei der Rundung auftretenden Fehler können abgeschätzt werden: Es gilt 1 · B e−` 2 für den absoluten Fehler und, falls x 6= 0 ist, gilt |rd(x) − x| ≤ |rd(x) − x| B ≤ · B −` |x| 2 6 für den relativen Fehler. Die Zahl rd(x) lässt sich zudem als rd(x) = x(1 + ε) mit |ε| ≤ B 2 · B −` darstellen. Die Schranke eps = B · B −` 2 für den relativen Rundungsfehler heißt Maschinengenauigkeit. Sie hängt von der Computerarchitektur ab. Es können zwei Ausnahmesituationen auftreten: • Ein Overflow tritt auf, falls rd(x) = ±0.m1 . . . m` . . . · B e mit e > emax dargestellt werden soll. In diesem Fall generiert der Computer üblicherweise eine Fehlermeldung und bricht das aktuelle Programm ab. • Ein Underflow tritt auf, falls rd(x) = ±0.m1 . . . m` . . . · B e mit e < emin dargestellt werden soll. Der Computer behandelt diesen Fall üblicherweise automatisch, indem x auf Null gesetzt wird. War x 6= 0, so ist der relative Fehler immer 1, während der absolute Fehler kleiner als µ (kleinste positive Maschinenzahl) und damit sehr klein ist. Beispiel 1.3.1 (Overflow) Wir machen einen Test mit MATLAB. Die Befehlszeile X = 1000.0; for i=1:7 X=X*X, end liefert die Ausgabe 1000000 1.0000e+12 1.0000e+24 1.0000e+48 1.0000e+96 1.0000e+192 Inf Bemerkung 1.3.2 (Runden durch Abschneiden) Ein alternativer Weg des Rundens ist das Abschneiden (chopping). Für eine Zahl x = ±0.m1 m2 . . . m` m`+1 . . .·B e wird rd(x) definiert durch die Maschinenzahl rd(x) = ±0.m1 m2 . . . m` · B e , d.h. die Ziffern m`+1 . . . werden einfach vernachlässigt. Absoluter Fehler: |rd(x) − x| ≤ B e−` . Relativer Fehler für x 6= 0: |rd(x) − x| ≤ B −`+1 . |x| 7 Beachte, dass diese Fehler doppelt so groß sind wie bei der exakten Rundung. Für Maschinenzahlen x und y werden die Grundrechenarten +, −, ·, / durch Operationen ⊕, , , gemäß x ∗ y = rd(x ∗ y), ∗ ∈ {+, −, ·, /} realisiert. Für diese Realisierung gilt x ∗ y = (x ∗ y)(1 + ε∗ ), |ε∗ | ≤ eps, ∗ ∈ {+, −, ·, /}. Insbesondere die Subtraktion zweier Zahlen x und y kann kritisch sein, wenn x ≈ y. Dies kann wie folgt eingesehen werden: Der Ansatz z = rd(rd(x) + rd(y)) = (x(1 + ε1 ) + y(1 + ε2 ))(1 + ε3 ) liefert in erster Näherung (d.h. durch Vernachlässigung der Fehlerterme höherer Ordnung in der Taylorentwicklung) xε1 + yε2 |(x + y) − z| · = ε3 + . |(x + y)| (x + y) Speziell für x ≈ −y ist der relative Fehler im schlimmsten Fall unbeschränkt. Diese Beobachtung ist der Grund für schwerwiegende numerische Probleme, die in Berechnungen auftreten können, da ein Fehler in einem Zwischenergebnis durch weitere Rechenschritte vergrößert werden kann und der Gesamtfehler sich somit aufschaukelt. Dies kann zu völlig unbrauchbaren Ergebnissen führen. Dieses Phänomen heißt Auslöschung. Beispiel 1.3.3 Seien x = 0.3721448693 und y = 0.3720214371 gegeben. Angenommen, unser Computer besitzt 5 Stellen Genauigkeit. Dann gelten rd(x) = 0.37214 und rd(y) = 0.37202. Subtraktion liefert z = rd(rd(x) − rd(y)) = rd(0.37214 − 0.37202) = 0.00012. Das korrekte Resultat ist x − y = 0.0001234322. Der relative Fehler beträgt somit |(x − y) − z| 0.0000034322 = ≈ 3 · 10−2 . |(x − y)| 0.0001234322 Dieser Fehler ist sehr groß im Vergleich zu den relativen Fehlern von rd(x) und rd(y), die jeweils durch 0.5 · 10−4 beschränkt sind. Durch Subtraktion zweier in etwa gleich großer Zahlen haben wir also 3 Stellen an Genauigkeit verloren! Gelegentlich ist es möglich, den Verlust an Genauigkeit zu vermeiden, indem alternative, analytisch äquivalente Rechenvorschriften verwendet werden. Beispiel 1.3.4 √ Betrachte die Funktion f (x) = x2 + 1 − 1. Für x ≈ 0 entsteht bei der Subtraktion ein 8 √ x2 + 1 ≈ 1 gilt. Jedoch kann f dargestellt werden als √ x2 + 1 + 1 √ x2 x2 + 1 − 1 · √ =√ . f (x) = x2 + 1 + 1 x2 + 1 + 1 Verlust an Genauigkeit, da dann Wenn wir wiederum mit 5 Stellen Genauigkeit rechnen und x = 10−3 wählen, dann wird f (x) gemäß der ersten Formel fälschlicherweise zu 0 berechnet, während die zweite Formel 0.5 · 10−6 liefert. Der exakte Wert ist nahe bei 0.4999998750000625 · 10−6 . Dies führt auf einen relativen Fehler von 1 für die erste Formel bzw. von ungefähr 0.25 · 10−6 für die zweite Formel. Der IEEE Standard Die Zahldarstellung auf Computern orientiert sich in der Regel am IEEE Standard. Dieser Standard besitzt folgende Merkmale: • Binärsystem, korrekte Rundung • größte ganze Zahl: 231 − 1 = 2147483647 • Gleitpunktdarstellung in einfacher Genauigkeit (single precision) x = ±(1.f )2 ·2c−127 ± c im Exponenten f in normierter Mantisse (1.f )2 |{z} | 1 bit {z }| 8 bits {z } 23 bits Es gilt 0 < c < (11111111)2 = 255 und −127 < c − 127 < 128. Die Werte c = 0 und c = 255 sind reserviert für die Darstellung von u.a. ±0 und ±∞. • Gleitpunktdarstellung in 2c−1023 ± c im Exponenten |{z} | 1 bit {z 11 bits doppelter Genauigkeit (double precision) x = ±(1.f )2 · f in normierter Mantisse (1.f )2 }| {z 52 bits } Es gilt 0 < c < (11111111111)2 = 2047 und −1023 < c − 1023 < 1024. Die Werte c = 0 und c = 2047 sind reserviert für die Darstellung von u.a. ±0 und ±∞. Länge Exponent Länge Mantisse Maschinengenauigkeit kleinste positive Zahl größte Zahl Single Precision Double Precision 8 Bits 23 Bits −23 2 ≈ 1.192 · 10−7 2−126 ≈ 1.175 · 10−38 (2 − 2−23 )2127 ≈ 3.403 · 1038 11 Bits 52 Bits −52 2 ≈ 2.220 · 10−16 2−1022 ≈ 2.225 · 10−308 (2 − 2−52 )21023 ≈ 1.798 · 10308 9 1.4 Taylorentwicklung und Landau-Symbole In der Numerik werden wir an vielen Stellen den aus der Analysis bekannten Satz von Taylor benötigen. Daher wollen wir ihn hier wiederholen. Satz 1.4.1 (Satz von Taylor) f : [a, b] → R sei (n+1)-mal stetig differenzierbar in [a, b]. Dann gilt für jedes x, x0 ∈ [a, b] f (x) = n X f (k) (x0 ) k=0 k! (x − x0 )k + Rn+1 (x; x0 ). Für das Restglied gilt die Restglieddarstellung nach Lagrange: Rn+1 (x; x0 ) = f (n+1) (ξ) (x − x0 )n+1 , (n + 1)! wobei ξ = ξ(x, x0 ) zwischen x und x0 liegt. Falls f also (n + 1)-mal stetig differenzierbar in [a, b] ist, kann f in x ∈ [a, b] durch das n-te Taylorpolynom in x0 ∈ (a, b) Tn (x; x0 ) := n X f (k) (x0 ) k! k=0 (x − x0 )k approximiert werden. Der Approximationsfehler in x kann abgeschätzt werden durch f (n+1) (ξ) |f (x) − Tn (x; x0 )| = (x − x0 )n+1 (n + 1)! f (n+1) (ξ) ≤ max · |x − x0 |n+1 ξ∈[a,b] (n + 1)! = C · |x − x0 |n+1 . Beispiel 1.4.2 Approximation von f (x) = exp(x) in x0 = 0 durch Tn (x; 0): T0 (x; 0) = 1, T1 (x; 0) = 1 + x, T2 (x; 0) = 1 + x + .. . Tn (x; 0) = n X xi i=0 10 i! . x2 , 2! (1.6) 20 exp(x) n = 4 : T4 (x; 0) 15 10 n = 2 : T2 (x; 0) 5 n = 1 : T1 (x; 0) n = 0 : T0 (x; ; 0) 0 –2 –1 0 x 1 2 3 Der Satz von Taylor gilt nicht nur im eindimensionalen, sondern auch im Rn . Hier verwendet man allerdings nur Approximationen bis zur Ordnung 2, mit dem Gradienten als erste und der Hesse-Matrix als zweite Ableitung. Um Konvergenzverhalten zu untersuchen, ist es notwendig zwei Funktionen f (x) und g(x) miteinander zu vergleichen, wenn x gegen einen Wert x0 strebt. Hierzu verwendet man die Landau-Symbole O und o. Definition 1.4.3 (Landau-Symbole) Seien f, g : D ⊆ Rn → R gegeben. (a) Wir schreiben f (x) = O(g(x)) für x → x0 , falls lim sup x→x0 |f (x)| < ∞. |g(x)| (b) Wir schreiben f (x) = o(g(x)) für x → x0 , falls lim x→x0 f (x) = 0. g(x) Mit Hilfe der O-Notation können wir (1.6) schreiben als f (x) − Tn (x; x0 ) = O(|x − x0 |n+1 ) für x → x0 . Ist f stetig differenzierbar, so gilt f (x) = f (x0 ) + f 0 (x0 )(x − x0 ) + o(kx − x0 k) mit limx→x0 o(kx−x0 k) kx−x0 k = 0. 11 Kapitel 2 Lineare Gleichungssysteme Das Folgende gilt sowohl für reellwertige als auch für komplexwertige Vektoren und Matrizen. Daher verwenden wir das Symbol K, welches für R oder C steht. Gesucht ist die Lösung x = (x1 , . . . , xn )> ∈ Kn des linearen Gleichungssystems a11 x1 + a12 x2 + . . . + a1n xn = b1 , a21 x1 + a22 x2 + . . . + a2n xn = b2 , .. . an1 x1 + an2 x2 + . . . + ann xn = bn , wobei aij , i, j = 1, . . . , n, und bi , a11 a12 a21 a22 A= .. .. . . an1 an2 i = 1, . . . , n, gegebene Zahlen sind. Mit x1 b1 · · · a1n · · · a2n x2 b2 .. .. , b = .. , x = .. . . . . xn bn · · · ann (2.1) (2.2) lautet das lineare Gleichungssystem in Matrixschreibweise (A ∈ Kn×n , x, b ∈ Kn ). Ax = b (2.3) Lineare Gleichungssysteme treten in nahezu allen Anwendungen und Verfahren auf. Beispiel 2.0.1 (Elektrische Schaltkreise, Kirchhoffsches Gesetz und lineare Gleichungssysteme) Gegeben sei das folgende stationäre elektronische Netzwerk mit gegebenen ohmschen Widerständen R1 , . . . , R5 und gegebenem Eingangsstrom I0 . I4 R3 I1 I0 R1 R2 4 1 R4 I2 R2 3 R5 I8 I10 I5 R4 2 I6 I3 R1 I9 R3 I7 12 Das erste Kirchhoffsche Gesetz (Knotensatz) An jedem Knotenpunkt ist die Summe aller zu- (positiven) und abfließenden (negativen) Ströme unter Beachtung der durch die Pfeile angegebenen Richtungen in jedem Zeitpunkt gleich Null liefert die Beziehungen −I1 + I3 = −I0 , I1 − I2 − I4 + I6 = 0, I2 − I3 − I7 = 0, I4 − I5 + I8 = 0, I5 − I6 + I7 − I9 = 0, −I8 + I9 − I10 = 0. Aus dem Ohmschen Gesetz U =R·I und dem zweiten Kirchhoffschen Gesetz (Maschensatz) In einer Masche ist die Summe aller Teilspannungen in jedem Zeitpunkt gleich Null folgen die Gleichungen R1 · I1 + R4 · I2 + R2 · I3 = 0, R2 · I8 + R4 · I5 + R1 · I9 = 0, R4 · I2 + R5 · I6 + R3 · I7 = 0, R3 · I4 + R4 · I5 + R5 · I6 = 0. Diese 10 Gleichungen für die Ströme I1 , . . . , I10 können als lineares Gleichungssystem geschrieben werden: −1 0 1 0 0 0 0 0 0 0 I1 −I0 1 −1 0 −1 0 1 0 0 0 0 I2 0 0 1 −1 0 0 0 −1 0 0 0 I3 0 0 0 1 −1 0 0 1 0 0 I4 0 0 0 0 0 0 1 −1 1 0 −1 0 · I5 = 0 . 0 0 0 0 0 0 0 −1 1 −1 I6 0 R1 R4 R2 0 0 0 0 0 0 0 I7 0 0 I8 0 0 0 0 R 0 0 R R 0 4 2 1 0 0 0 R5 R3 0 0 0 I9 0 0 R4 0 0 0 R3 R4 R5 0 0 0 0 I10 0 13 Unser Ziel ist es, Verfahren zur Bestimmung einer Lösung x von (2.3) zu entwickeln und zu analysieren. 2.1 Grundlagen und Lösbarkeit Im folgenden werden wir verschiedene Typen von Matrizen und einige Eigenschaften verwenden. Eine Matrix A ∈ Kn×n heißt • Diagonalmatrix, falls aij = 0 für alle i 6= j gilt, d.h. ∗ .. A= . ∗ • Tridiagonalmatrix, falls aij = 0 für alle i, j ∈ {1, . . . , n} mit |i − j| > 1 gilt, d.h. ∗ ∗ ∗ ∗ ∗ .. . A= ∗ ∗ .. .. . . ∗ ∗ ∗ • untere Dreiecksmatrix, falls aij = 0 für alle i < j gilt, d.h. A= ∗ ∗ .. . ∗ ∗ .. . . . . ∗ ··· ∗ • obere Dreiecksmatrix, falls aij = 0 für alle i > j gilt, d.h. ∗ ∗ ··· ∗ ∗ ··· ∗ A= . . .. . . ∗ • normierte untere (obere) Dreiecksmatrix, falls A eine untere (obere) Dreiecksmatrix ist, deren Diagonaleinträge alle gleich Eins sind. 14 Desweiteren heißt A ∈ K n×n • hermitesch, wenn A∗ = A mit A∗ = Ā> gilt, wobei A> die transponierte Matrix von A bezeichnet und Ā die konjugiert komplexe Matrix von A ist. • symmetrisch, wenn K = R und A> = A gilt. • unitär, falls A∗ A = I gilt. • orthogonal, falls K = R und A> A = I gilt. • positiv semidefinit, falls x∗ Ax ≥ 0 für alle x ∈ Kn gilt. • positiv definit, falls x∗ Ax > 0 für alle x ∈ Kn , x 6= 0, gilt. Definition 2.1.1 Sei A ∈ Km×n eine Matrix. Dann heißen Kern(A) := {x ∈ Kn | Ax = 0} Bild(A) := {y ∈ Km | ∃x ∈ Kn : y = Ax} Rang(A) := dim(Bild(A)) Kern von A, Bild von A, Rang von A. Kern(A) ist ein Unterraum des Kn und Bild(A) ist ein Unterraum von Km . Rang(A) ist die maximale Anzahl von linear unabhängigen Spalten von A, was gleich der Anzahl der linear unabhängigen Zeilen von A ist. Es stellt sich die Frage, welche Voraussetzungen an A und b gestellt werden müssen, damit das lineare Gleichungssystem (2.3) überhaupt eine Lösung x besitzt. Zunächst ist klar, dass das lineare Gleichungssystem genau eine Lösung besitzt, falls die Matrix A invertierbar (regulär) ist. Die eindeutige Lösung des Gleichungssystems ist dann durch x = A−1 b (2.4) gegeben. Es kann aber auch der Fall eintreten, dass (2.3) mehrere Lösungen oder gar keine Lösung besitzt. Offensichtlich kann die Matrix A dann nicht regulär sein. Beispiel 2.1.2 Das folgende lineare Gleichungssystem besitzt keine Lösung: ! ! ! 1 0 x1 1 = 0 0 x2 1 Das folgende lineare Gleichungssystem besitzt genau eine Lösung: ! ! ! 1 0 x1 1 = 0 1 x2 1 15 Das folgende lineare Gleichungssystem besitzt unendlich viele Lösungen: ! ! ! 1 0 x1 1 = 0 0 x2 0 Es gilt das folgende Kriterium, welches die Frage der Lösbarkeit von (2.3) abschließend beantwortet. Satz 2.1.3 Das lineare Gleichungssystem (2.3) besitzt genau dann mindestens eine Lösung, falls Rang(A) = Rang(A|b), mit der um die rechte Seite b erweiterte a11 a21 (A|b) = .. . an1 Matrix a12 · · · a1n b1 a22 · · · a2n b2 .. .. .. .. . . . . . an2 · · · ann bn Beweis: Offenbar besitzt das lineare Gleichungssystem genau dann eine Lösung, wenn b ∈ Bild(A) gilt. Sei nun b ∈ Bild(A). Dann gilt auch Bild(A) = Bild(A|b) und somit Rang(A) = Rang(A|b). Gilt nun Rang(A) = Rang(A|b), so lässt sich b als Linearkombination der Spalten von A schreiben. Also gilt b ∈ Bild(A) und somit besitzt das Gleichungssystem eine Lösung. 2 Die Struktur der Lösungsmenge eines linearen Gleichungssystems wird im folgenden Satz behandelt. Satz 2.1.4 Seien A ∈ Kn×n und b ∈ Kn gegeben. Sei x̂ eine Lösung von Ax = b. Die gesamte Lösungsmenge von Ax = b ist gegeben durch x̂ + Kern(A). Beweis: Sei y ∈ Kern(A), d.h. Ay = 0. Dann gilt A(x̂ + y) = Ax̂ + Ay = b und x̂ + y ist Lösung. Sei x Lösung. Dann gilt A(x − x̂) = b − b = 0 und somit x − x̂ ∈ Kern(A) bzw. x ∈ x̂ + Kern(A). 2 Wir fassen für A ∈ Km×n noch einige aus der linearen Algebra bekannte Resultate zusammen: 16 • Für A ∈ Km×n gelten n = dim(Kern(A)) + dim(Bild(A)) m = dim(Kern(A> )) + Rang(A) • Sind A ∈ Kn×n und B ∈ Kn×n invertierbar, so gelten (A · B)−1 = B −1 · A−1 und (A∗ )−1 = (A−1 )∗ . • Sei A ∈ Kn×n . Die folgenden Aussagen sind äquivalent: – – – – – – – – A ist regulär (invertierbar) Kern(A) = {0} Bild(A) = Kn Rang(A) = n Ax = 0 ⇔ x = 0 Ax = b hat für jedes b genau eine Lösung det(A) 6= 0 alle Eigenwerte von A sind von Null verschieden • Für A ∈ Km×n gelten die folgenden Beziehungen – Kern(A) = Bild(A∗ )⊥ und Bild(A∗ ) = Kern(A)⊥ – Kern(A) = Kern(A∗ A) – Bild(A∗ ) = Bild(A∗ A) Hierin bezeichnet V ⊥ = {w ∈ Kn | w∗ v = 0 ∀v ∈ V } das orthogonale Komplement des Unterraums V ∈ Kn . Um das lineare Gleichungssystem numerisch lösen zu können, gehen wir im Folgenden davon aus, dass die Matrix A regulär ist und somit genau eine Lösung existiert. Die numerische Berechnung einer nicht eindeutigen Lösung im Falle einer singulären Matrix A ist komplizierter und wird hier nicht betrachtet. 2.2 Gauß’sches Eliminationsverfahren und LR-Zerlegung Es werden Algorithmen zur numerischen Lösung des linearen Gleichungssystems (2.3) entwickelt. Es wird vorausgesetzt, dass die Matrix A regulär ist. Formal kennen wir die Lösung dann schon, sie ist durch (2.4) gegeben. Dort wird aber die Inverse von A benötigt, die im Allgemeinen nicht bekannt ist und für große Matrizen auch nicht in einfacher Weise berechnet werden kann. 17 2.2.1 Einfach lösbare Gleichungssysteme Zunächst widmen wir uns linearen Gleichungssystemen, die einfach“ lösbar sind. ” Wir betrachten insbesondere lineare Gleichungssysteme, in denen die Matrix A Dreiecksgestalt hat. Diese Gleichungssysteme werden später im Gauß’schen Eliminationsverfahren bzw. der LR-Zerlegung eine wichtige Rolle spielen. Zunächst betrachten wir das lineare Gleichungssystem Ly = b (2.5) für y ∈ Kn mit einem Vektor b ∈ Kn und einer linken unteren Dreiecksmatrix `11 `21 `22 ∈ Kn×n . L= . .. . . . `n1 `n2 · · · `n,n In Komponenten lautet es `11 y1 = b1 , `21 y1 + `22 y2 = b2 , .. . `n1 y1 + `n2 y2 + . . . + `nn yn = bn . Dieses Gleichungssystem können wir von oben beginnend sukzessive nach den Variablen yi , i = 1, . . . , n auflösen: Algorithmus 2.2.1 (Vorwärtssubstitution) Input: untere Dreiecksmatrix L, rechte Seite b Output: Lösung y von Ly = b. y1 = b1 /`11 ; for i = 2:n, y i = bi ; for j = 1:(i-1), yi = yi − `ij · yj ; end yi = yi /`ii ; end Die Vorwärtssubstitution ist genau dann durchführbar, wenn `ii 6= 0 für alle i = 1, . . . , n gilt, also wenn L regulär ist. Der Algorithmus heißt dann wohldefiniert. Der Aufwand der Vorwärtssubstitution beträgt O(n2 ). 18 Entsprechend kann man auch bei der Lösung des linearen Gleichungssystems Rx = c für x ∈ Kn mit rechter oberer Dreiecksmatrix r11 r12 · · · r22 · · · R= .. . (2.6) r1n r2n .. . rnn vorgehen. In diesem Fall kann man die Gleichungen r11 x1 + r12 x2 + . . . + r1n xn = c1 , .. . rn−1,n−1 xn−1 + rn−1,n xn = cn−1 , rnn xn = cn beginnend mit der letzten Zeile nach xi , i = n, n − 1, . . . , 1, auflösen: Algorithmus 2.2.2 (Rückwärtssubstitution) Input: obere Dreiecksmatrix R, rechte Seite c Output: Lösung x von Rx = c. xn = cn /rnn ; for i = (n-1):-1:1, x i = ci ; for j = (i+1):n, xi = xi − rij · xj ; end xi = xi /rii ; end Der Algorithmus ist wiederum genau dann wohldefiniert, wenn rii 6= 0 für alle i = 1, . . . , n gilt, also genau dann, wenn R regulär ist. Der Aufwand der Rückwärtssubstitution beträgt wiederum O(n2 ). 2.2.2 Das Gauß’sche Eliminationsverfahren Wir haben im vorhergehenden Abschnitt gesehen, dass lineare Gleichungssysteme mit Dreiecksstruktur direkt gelöst werden können. Die Idee des Gauß’schen Eliminationsverfahrens besteht darin, die Matrix A durch elementare Zeilenumformungen schrittweise in eine rechte obere Dreiecksmatrix R und die rechte Seite b in einen Vektor c zu überführen, so dass ein Gleichungssystem der Form (2.6) entsteht. Dieses kann dann 19 mittels Rückwärtssubstitution gelöst werden. Wichtig ist hierbei, dass die Lösung des Ausgangsproblems mit der des transformierten Problems übereinstimmt. Dies ist bei der Verwendung elementarer Zeilenumformungen gewährleistet. Elementare Zeilenumformungen sind • die Multiplikation einer Zeile mit einem Wert ungleich Null, • die Addition bzw. Substraktion zweier Zeilen, • sowie das Vertauschen zweier Zeilen. Schematisch läuft der Gauß’sche Eliminationsalgorithmus wie folgt ab. Gauß’scher Eliminationsalgorithmus: (i) Setze A(1) := A und b(1) := b. (ii) Beginnend mit A(1) x = b(1) führe n − 1 Transformationsschritte durch bis ein äquivalentes lineares Gleichungssystem mit oberer Dreiecksstruktur erreicht ist: A(1) x = b(1) A(2) x = b(2) ⇔ ⇔ ... ⇔ A(n) x = b(n) (iii) Löse das lineare Gleichungssystem A(n) x = b(n) durch Rückwärtssubstitution. Beispiel 2.2.3 Betrachte | 6 −2 2 4 12 −8 6 10 3 −13 9 3 −6 4 1 −18 {z } 12 34 27 −38 {z 6 −2 2 4 12 12 −8 6 10 34 3 −13 9 3 27 −6 4 1 −18 −38 x1 x2 x3 x4 =A(1) = | =b(1) } Gauß’sche Elimination liefert folgendes Ergebnis: Start: h A(1) b(1) i = 20 Subtraktion des 2-fachen der 1. Zeile von der 2. Zeile und des 0.5-fachen der 1. Zeile von der 3. Zeile und des −1-fachen der 1. Zeile von der 4. Zeile liefert: 6 −2 2 4 12 i 0 −4 2 h 2 10 A(2) b(2) = 0 −12 8 1 21 0 2 3 −14 −26 Subtraktion des 3-fachen der 2. Zeile von der 3. von der 4. Zeile liefert: 6 −2 i 0 −4 h A(3) b(3) = 0 0 0 0 Zeile und des −0.5-fachen der 2. Zeile 12 2 4 2 2 10 2 −5 −9 4 −13 −21 Subtraktion des 2-fachen der 3. Zeile von der 4. Zeile liefert: 6 −2 2 4 12 i 0 −4 2 h 2 10 A(4) b(4) = 0 0 2 −5 −9 0 0 0 −3 −3 Rückwärtssubstitution: x = (1, −3, −2, 1)> . Dieses Vorgehen lässt sich zum Gauß’sche Eliminationsverfahren verallgemeinern: Algorithmus 2.2.4 (Gauß’scher Eliminationsalgorithmus) Input: Matrix A ∈ Kn×n , rechte Seite b ∈ Kn Output: Lösung x mit Ax = b, A und b werden überschrieben. (Elimination) for k = 1:n-1 for i = (k+1):n piv = aik /akk ; for j = (k+1):n aij = aij − piv ∗ akj ; end bi = bi − piv ∗ bk ; end end 21 (Rückwärtssubstitution) Führe Rückwärtssubstitution für A und b durch. Man beachte, dass in der Schleife für j der Index k weggelassen wurde, da hierdurch die Nullen unterhalb der Hauptdiagonalen erzeugt werden. 2.2.3 Pivotisierung Der Gauß’sche Eliminationsalgorithmus ist durchführbar, solange die in jedem Schritt überschriebenen Pivotelemente akk 6= 0 für alle k = 1, . . . , n − 1 erfüllen. Es kann gezeigt werden, dass dies für strikt diagonaldominante Matrizen und symmetrische, positiv definite Matrizen ! der Fall ist. Allerdings ist das Verfahren für sehr einfache Matrizen, wie 0 1 A= nicht durchführbar. Darüber hinaus entstehen numerische Probleme be1 0 dingt durch Rundungsfehler und Fehlerfortpflanzung bereits für akk ≈ 0 wie das folgende Beispiel zeigt. Beispiel 2.2.5 Betrachte ! ! ! ε 1 x1 1 = 1 1 x2 2 wobei ε ≈ 0. Die exakte Lösung ist gegeben durch x1 = 1 ≈ 1, 1−ε x2 = Der Gauß’sche Eliminationsalgorithmus liefert ! ! ε 1 x1 = 0 1 − ε−1 x2 1 − 2ε ≈ 1. 1−ε 1 2 − ε−1 ! . Die Lösung dieses Systems lautet x2 = 2 − ε−1 , 1 − ε−1 x1 = (1 − x2 )ε−1 . Für betragsmäßig kleine Werte von ε erhält man auf Grund von Rundungsfehlern x2 ≈ 1 und anschließend x1 ≈ 0, also ein völlig falsches Ergebnis. Falls wir die Gleichungen vertauschen, d.h. ! ! 1 1 x1 = ε 1 x2 2 1 ! , und anschließend den Gauß’schen Algorithmus anwenden, erhalten wir ! ! ! 1 1 x1 2 = 0 1−ε x2 1 − 2ε 22 und die (beinahe) exakte Lösung x2 = 1 − 2ε ≈ 1, 1−ε x1 = 2 − x2 ≈ 1. Zeilenvertauschungen haben im obigen Beispiel zum Erfolg geführt. Dies bedeutet, dass die natürliche Reihenfolge bei der Wahl der Pivotelemente geändert wurde. Ziel dabei ist es, in jedem Schritt i die Zeile i mit einer der verbleibenden Zeilen k ≥ i zu vertauschen, so dass wir ein von Null verschiedenes Diagonalelement akk haben. Ein solches Element existiert stets, da A regulär ist. Um die numerisch problematische Division durch ein betragsmäßig sehr kleines Pivotelement zu vermeiden, wird unter den in Frage kommenden Elementen der Spalte i das betragsmäßig größte Element als neues Pivotelement ausgewählt. Der Gauß’sche Algorithmus zusammen mit der Pivotisierung ist für reguläre Matrizen A bei exakter Rechnung stets durchführbar. Natürlich müssen Zeilenvertauschungen auch in b(i) berücksichtigt werden. In einer effizienten Implementierung werden Zeilenvertauschungen nicht explizit durchgeführt, um zeitaufwändiges (und unnötiges) Umspeichern zu vermeiden. Stattdessen werden die Zeilenvertauschungen in einem Indexvektor l = [l1 , l2 , . . . , ln ] protokolliert. Die folgende Variante des Gauß’schen Eliminationsverfahrens verwendet eine skalierte Spaltenpivotsuche zur Bestimmung des Pivotelements. Algorithmus 2.2.6 (Gauß’sche Elimination mit skalierter Pivotsuche) (i) Setze A(1) := A, b(1) := b, l := [l1 , l2 , . . . , ln ] = [1, 2, . . . , n] und k := 1. (ii) Berechne die Skalierungsfaktoren (k) sli := max |ali ,j | 1≤j≤n (k ≤ i ≤ n). (k) (iii) Bestimme das Pivotelement alj ,k mit k ≤ j ≤ n gemäß (k) alj ,k sl j (k) a l ,k := max i k≤i≤n sli (iv) Vertausche Werte lk und lj in l. (k) (v) Verwende Zeile lk von A(k) als Pivotzeile und das Element alk ,k als Pivotelement und berechne A(k+1) und b(k+1) . (vi) Falls k < n − 1, setze k := k + 1 und gehe zu (iii). (vii) Löse A(n) x = b(n) durch Rückwärtssubstitution in der Reihenfolge ln , ln−1 , . . . , l1 . 23 Bemerkung 2.2.7 • Es existiert eine Variante des Algorithmus, bei dem die Skalierungsfaktoren si in jedem Durchlauf neu berechnet werden. • Es ist wichtig zu bemerken, dass lediglich die Einträge in dem Indexvektor l vertauscht werden. Dies vermeidet das kostenintensive Umspeichern von Zeilen. • Die Berechnung von A(k+1) und b(k+1) in Schritt (v) erfolgt wie in Algorithmus (k) (k) (2.2.4), wobei der Zeilenindex i durch li ersetzt wird, d.h. ali ,k /alk ,k multipliziert mit der Zeile lk wird von den Gleichungen li , i = k + 1, . . . , n, subtrahiert, wohingegen die Zeilen l1 , . . . , lk unverändert bleiben. • Der Aufwand berechnet sich zu 13 n3 + O(n2 ). Beispiel 2.2.8 Betrachte | 3 −13 9 3 x1 −6 4 1 −18 x2 6 −2 2 4 x3 12 −8 6 10 x4 {z } =A(1) = | −19 −34 16 26 {z =b(1) . } Anwendung der Gauß’schen Elimination mit skalierter Spaltenpivotsuche liefert: Init: l = [1, 2, 3, 4], s = [13, 18, 6, 12]. Schritt 1: k = 1. Pivotelement (nicht eindeutig): (1) (1) a 6 a3,1 3 6 6 12 li ,1 max , , , = = = max . 1≤i≤4 sli 13 18 6 12 6 s3 Zeile 3 ist Pivotzeile. Vertausche Zeilenindizes: l = [3, 2, 1, 4]. Elimination: 0 −12 8 1 −27 0 2 3 −14 (2) −18 A(2) = , b = . 6 −2 2 16 4 0 −4 2 2 −6 Schritt 2: k = 2. Pivotelement: (2) (2) a a 2 12 4 12 1,2 l ,2 , , = = max i = max . 2≤i≤4 sli 18 13 12 13 s1 24 Zeile 1 ist Pivotzeile. Vertausche Zeilenindizes: l = [3, 1, 2, 4]. Elimination: A (3) = 0 −12 8 1 −27 0 0 13/3 −83/6 (3) −45/2 , b = . 6 −2 2 4 16 0 0 −2/3 5/3 3 Schritt 3: k = 3. Pivotelement: (3) (3) a a 2/3 13/3 13/3 l ,3 2,3 max i = max , = = . 3≤i≤4 sli s2 18 12 18 Zeile 2 ist Pivotzeile und l = [3, 1, 2, 4]. Elimination: A(4) = 0 −12 8 1 0 0 13/3 −83/6 6 −2 2 4 0 0 0 −6/13 (4) , b = −27 −45/2 16 −6/13 . Rückwärtssubstitution: x4 = 1, x2 = 1, x1 = 3, x3 = −2. 2.2.4 Die LR-Zerlegung einer Matrix In diesem Abschnitt werden wir sehen, dass das Gauß’sche Eliminationsverfahren nicht nur ein oberes Dreieckssystem produziert, sondern tatsächlich eine Faktorisierung A=L·R (2.7) der Matrix A erzeugt, wobei L eine untere Dreiecksmatrix und R eine obere Dreiecksmatrix ist. Das lineare Gleichungssystem Ax = b kann dann gelöst werden, indem zunächst y := Rx definiert wird und mit Hilfe der Vorwärtssubstitution Ly = b gelöst wird. Damit erhält man y. Im Anschluss wird mit der Rückwärtssubstitution Rx = y gelöst (y ist darin bereits bekannt). Dies liefert schließlich die Lösung x von Ax = b. Dieser Vorgang wird Vorwärts-Rückwärtssubstitution genannt. Diese Vorgehensweise ist besonders dann von Vorteil, falls das lineare Gleichungssystem Ax = b für viele verschiedene rechte Seiten b gelöst werden muss. Die LR−Zerlegung muss nur einmal berechnet werden. 25 Beispiel 2.2.9 (Berechnung der Inversen A−1 ) Die Inverse X = A−1 einer invertierbaren Matrix A ist gegeben durch A · X = I. Falls x(i) die i-te Spalte von X und e(i) den i-ten Einheitsvektor bezeichnen, müssen wir zur Bestimmung von X die folgenden n linearen Gleichungssysteme lösen: Ax(i) = e(i) , i = 1, 2, . . . , n. Falls die LR−Zerlegung von A bekannt ist, müssen wir lediglich n Vorwärts-RückwärtsSubstitutionen durchführen. Der Gesamtaufwand zur Invertierung der Matrix beträgt damit O(n3 ) (LR: O(n3 ); VRS: O(n · n2 )). Würde das Gauß’sche Eliminationsverfahren für jedes der n Gleichungssysteme von Neuem angewendet werden, würde der Aufwand O(n4 ) betragen. Bemerkung 2.2.10 Es ist niemals notwendig, die Inverse A−1 einer Matrix explizit zu berechnen, wenn wir die Lösung x = A−1 b eines linearen Gleichungssystems berechnen wollen. Es genügt, die LR-Zerlegung von A zu berechnen und Vorwärts-Rückwärtssubstitution durchzuführen. Beispiel 2.2.11 A= 6 −2 2 4 12 −8 6 10 3 −13 9 3 −6 4 1 −18 . Jeder Schritt des Gauß’schen Eliminationsverfahrens entspricht der Multiplikation von links mit einer Matrix Lk , so dass L3 · L2 · L1 · A = R | {z } =A(2) | {z } =A(3) | {z } ⇒ −1 A = L−1 · L−1 2 · L3 ·R |1 {z } =:L =A(4) wobei L1 = 1 −2 − 12 1 0 1 0 0 0 0 1 0 0 0 0 1 , L2 = 1 0 0 0 1 0 0 −3 1 1 0 0 2 26 0 0 0 1 , L3 = 1 0 0 0 0 0 1 0 0 1 0 −2 0 0 0 1 . Zusammenfassend erzeugt das Gauß’sche Eliminationsverfahren die LR-Zerlegung von A: L= 0 0 0 1 0 1 0 0 2 1 = 1 1 2 0 3 1 0 2 1 −1 − 2 2 1 −1 0 {z | 1 2 0 0 1 0 0 0 0 1 · } | =L−1 1 1 0 0 1 0 3 0 − 21 {z 0 0 1 0 =L−1 2 0 0 0 1 · } | 1 0 0 0 0 0 1 0 0 1 0 2 {z 0 0 0 1 } =L−1 3 und R= 6 −2 2 4 0 −4 2 2 0 0 2 −5 0 0 0 −3 . Allgemein lässt sich zeigen, dass der Transformationsschritt i 7→ i + 1 im Gauß’schen Eliminationsverfahren mit Hilfe von linksseitigen Matrixmultiplikationen beschrieben werden kann. Es gilt dabei A(i+1) = Li · A(i) mit Li = 1 .. . 1 li+1,i 1 .. .. . . lni 1 , (i) lji = − aji (i) , j = i + 1, . . . , n. (2.8) aii Die Inverse L−1 zu Li lässt sich leicht angeben (und durch Nachrechnen überprüfen): i L−1 i = 1 .. . 1 −li+1,i 1 .. .. . . −lni 1 27 . Das Gauß’sche Eliminationsverfahren liefert also nach n − 1 Schritten die Beziehung R = A(n) = Ln−1 · A(n−1) = Ln−1 · Ln−2 · A(n−2) .. . = Ln−1 · Ln−2 · · · L1 · A(1) . (2.9) −1 −1 (1) Mit L := L−1 = A folgt 1 · · · Ln−2 · Ln−1 und A A = L · R. Da das Produkt von zwei normierten unteren Dreiecksmatrizen wieder eine normierte untere Dreiecksmatrix ergibt, ist L normierte untere Dreiecksmatrix mit L= 1 −l21 1 .. .. .. . . . −ln1 · · · −ln,n−1 1 , (2.10) wobei die lij ’s durch (2.8) gegeben sind, was einfach zu verifizieren ist. Um Speicherplatz zu sparen, speichern effiziente Computerimplementierungen L und R in A, d.h. die Elemente lij werden unterhalb der Diagonalen gespeichert und die Elemente von R werden auf und oberhalb der Diagonalen gespeichert. Algorithmus 2.2.12 (LR-Zerlegung ohne Pivotisierung) Input: Matrix A ∈ Kn×n , rechte Seite b ∈ Kn Output: Lösung x ∈ Kn des LGS Ax = b, A wird überschrieben und enthält unterhalb der Diagonalen die Matrix L und auf und oberhalb der Diagonalen die Matrix R (LR Zerlegung) for k = 1:n-1 for i = (k+1):n piv = aik /akk ; for j = (k+1):n aij = aij − piv · akj ; end aik = piv; end end ( Vorwärtssubstitution Ly = b (Ergebnis y steht in x)) 28 x 1 = b1 ; for i=2:n x i = bi ; for j=1:(i-1) xi = xi − aij · xj ; end end ( Rückwärtssubstitution Rx=y ) xn = xn /ann ; for i=(n-1):1 for j=(i+1):n xi = xi − aij · xj ; end xi = xi /aii ; end 2.2.5 Pivotisierung in Matrixnotation Die Pivotsuche lässt sich ebenfalls durch Matrixmultiplikationen darstellen. Die Multiplikation einer Matrix von links mit der Permutationsmatrix P = 1 .. Zeile i Zeile j . 1 0 1 1 .. . 1 1 0 1 .. . 1 Spalte i Spalte j vertauscht die Zeilen i und j, während die Multiplikation von rechts die Spalten i und j vertauscht. Zusätzlich gilt P > P = I und P = P > . Die Pivotsuche im Schritt i entspricht daher der Multiplikation von A(i) von links mit einer Permutationsmatrix Pi bevor die 29 Subdiagonalelemente eliminiert werden: A(i+1) = Li · Pi · A(i) . Analog zu (2.9) erhält man R = Ln−1 · Pn−1 · Ln−2 · Pn−2 · · · L1 · P1 · A. (2.11) Nun definieren wir für i = 1, . . . , n − 1 L̃i := Pn−1 · . . . · Pi+1 · Li · Pi+1 · . . . · Pn−1 . Man kann zeigen dass L̃i die gleiche Struktur wie Li hat, und das lediglich von Null verschiedenen Elemente unterhalb der Diagonalen vertauscht werden. Wegen P · P = I folgt R = Ln−1 · Pn−1 · Ln−2 · Pn−2 · . . . · L2 · P2 · L1 · P1 · A = L̃n−1 · . . . · L̃2 · L̃1 · Pn−1 · . . . · P2 · P1 · A. −1 −1 mit P = Pn−1 · . . . · P2 · P1 und L̃ = L̃−1 1 · L̃2 · . . . · Ln−1 erhält man P · A = L̃ · R. Satz 2.2.13 (LR-Zerlegung mit Pivotisierung) Sei A ∈ Kn×n invertierbar. Dann definiert das Gauß’sche Eliminationsverfahren mit Pivotisierung eine Zerlegung P A = LR mit einer rechten oberen Dreiecksmatrix R, einer normierten linken unteren Dreiecksmatrix L und einer Permutationsmatrix P . Beweis: Bricht das Gauß’sche Eliminationsverfahren mit Pivotisierung nicht zusammen, (i) d.h. sind alle Pivotelemente aki 6= 0 für alle i = 1, . . . , n, so liefern die vorangegangenen Betrachtungen die Behauptung. Zu klären ist, dass das Verfahren für invertierbare Matrizen tatsächlich nicht zusammenbricht, d.h. dass alle Pivotelemente von Null verschieden sind. Angenommen, das Pivotelement im i-ten Teilschritt wäre Null, dann gilt gemäß Definition des Pivotelements (i) aji = 0 für alle i ≤ j ≤ n und A(i) hätte folgende Gestalt: (i) (i) (i) a11 · · · a1i . . . a1n . . . .. .. .. (i) 0 . . . ain . A(i) = (i) 0 . . . ai+1,n .. .. . . (i) 0 . . . ann Die Determinante des rechten unteren Teilblocks ist folglich 0 und somit ist det(A(i) ) = 0. Wegen ! i−1 Y det(A(i) ) = det(Li−1 · Pi−1 · · · L1 · P1 · A) = det(Lj ) det(Pj ) det(A) j=1 30 und det(Lj ) = 1, det(Pj ) = ±1, j = 1, . . . , i − 1, folgt dann det(A) = 0 im Widerspruch zur Invertierbarkeit von A. 2 Beispiel 2.2.14 Wir lösen das lineare Gleichungssystem 2 18 9 21 Ax = b mit A = 4 4 4 , b = 4 , 1 33 9 9 indem wir die LR-Zerlegung mit Pivotisierung berechnen. 1.Schritt: Pivotelement a21 = 4. 4 4 4 0 1 0 1 0 0 P1 = 1 0 0 , L1 = − 12 1 0 , A(2) = 0 16 7 0 32 8 0 0 1 − 14 0 1 (2) 2.Schritt: Pivotelement a32 = 32. 1 0 0 1 0 0 P2 = 0 0 1 , L2 = 0 1 0 , 0 − 21 1 0 1 0 4 4 4 = 0 32 8 0 0 3 A(3) Insgesamt: 0 1 0 P = 0 0 1 , 1 0 0 1 0 0 L̃ = 41 1 0 , 1 1 1 2 2 4 4 4 R = 0 32 8 0 0 3 Lösung: Es ist b̃ := P · b = P · Ax = L̃ · Rx. Es wird also L̃ · Rx = b̃ = (4, 9, 21)> mit Vorwärts-Rückwärtssubstitution gelöst. Vorwärtssubstitution: L̃y = b̃ ⇒ y = (4, 8, 15)> Rückwärtssubstitution: Rx = y ⇒ x = (−3, −1, 5)> 2.3 Cholesky-Zerlegung Treten bei einer Anwendung ausschließlich lineare Gleichungssysteme mit hermitescher, positiv definiter Matrix A ∈ Kn×n auf, sollte die Cholesky-Zerlegung an Stelle der LRZerlegung verwendet werden, da diese nur etwa den halben Aufwand erfordert. Definition 2.3.1 (Cholesky-Zerlegung) Sei A ∈ Kn×n . Eine Zerlegung der Form A = L · L∗ 31 mit linker unterer (nicht normierter) Dreiecksmatrix L und positiven Diagonalelementen heißt Cholesky-Zerlegung von A. Existiert zu einer hermiteschen Matrix A ∈ Kn×n eine Cholesky-Zerlegung A = L · L∗ , so ist die Matrix L nach Definition regulär, und es gilt x∗ Ax = x∗ L · L∗ x = kL∗ xk22 > 0 ∀x ∈ Kn \ {0}. Es gilt sogar folgende Äquivalenz, die man mit Hilfe vollständiger Induktion beweisen kann: Satz 2.3.2 Sei A ∈ Kn×n hermitesch. A besitzt genau dann eine Cholesky-Zerlegung, wenn A positiv definit ist. Diese Äquivalenz von Cholesky-Zerlegung und positiver Definitheit einer hermiteschen Matrix liefert einen praktisch durchführbaren Algorithmus, um eine gegebene hermitesche Matrix auf positive Definitheit zu testen. Man muss lediglich die Cholesky-Zerlegung berechnen. Falls sie existiert (d.h. es gilt `kk > 0 für alle k), ist A positiv definit, andernfalls nicht. Die konkrete Berechnung von L erfolgt durch zeilenweisen Koeffizientenvergleich von `¯11 `¯21 · · · `¯n1 a11 a12 · · · a1n `11 `¯22 · · · `¯n2 a21 a22 · · · a2n `21 `22 . .. .. = .. .. . . .. .. .. . . · . . . . . . . . . `¯nn an1 an2 · · · ann `n1 `n2 · · · `nn Es folgt (beachte: für reelle Matrizen, soll eine reelle Zerlegung entstehen): √ a11 = `11 `¯11 = |`11 |2 ⇒ `11 = a11 , a21 = `21 `¯11 ⇒ `21 = a21 /`¯11 , p a22 = `21 `¯21 + `22 `¯22 ⇒ `22 = a22 − |`21 |2 , a31 = `31 `¯11 ⇒ `31 = a31 /`¯11 , a33 a32 = `31 `¯21 + `32 `¯22 ⇒ `32 = (a32 − `31 `¯21 )/`¯22 , p = `31 `¯31 + `32 `¯32 + `33 `¯33 ⇒ `33 = a33 − |`31 |2 − |`32 |2 , .. .. .. . . . ! j−1 j−1 X X aij = `ik `¯jk + `ij `¯jj ⇒ `ij = aij − `ik `¯jk /`¯jj , k=1 aii = i−1 X k=1 `ik `¯ik + `ii `¯ii v u i−1 X u t ⇒ `ii = aii − `ik `¯ik k=1 k=1 32 j = 1, . . . , i − 1, Insgesamt entstehen die Formeln `ii v u i−1 X u = taii − |`ik |2 , i = 1, . . . , n, k=1 j−1 `ij = aij − X ! `ik `¯jk /`¯jj , j = 1, . . . , i − 1. k=1 Diese Formeln sind für positiv definite Matrizen wohldefiniert, da die Existenz einer Cholesky-Zerlegung bereits durch Satz 2.3.2 garantiert ist. Die obigen Formeln sind auch eindeutig unter der Forderung, dass für relle Matrizen eine reelle Zerlegung entstehen soll. Beispiel 2.3.3 Cholesky-Zerlegung: 3.1623 0.6325 −0.3162 3.1623 0 0 10 2 −1 A= 2 5 0 2.1448 0.0933 0 0 ≈ 0.6325 2.1448 0 0 1.9726 −0.3162 0.0933 1.9726 −1 0 4 Da A symmetrisch ist und die Cholesky-Zerlegung existiert (lkk > 0 für alle k), folgt, dass A positiv definit ist. Der Aufwand der Cholesky-Zerlegung beträgt 16 n3 + O(n2 ) und ist damit in etwa halb so groß wie bei der LR-Zerlegung. 2.4 Fehlereinfluss und Kondition Gegenstand dieses Abschnitts ist die Untersuchung der Kondition eines linearen Gleichungssystems. Beispiel 2.4.1 (Vandermonde-Matrix) Wir versuchen, das lineare Gleichungssystem n−1 n 1 2 4 8 ··· 2 x1 (2 − 1)/1 n−1 1 3 9 27 ··· 3 (3n − 1)/2 x2 1 x n−1 n 4 16 64 · · · 4 (4 − 1)/3 3 = . . .. .. .. .. . .. .. . .. . . . . . . 2 3 n−1 n 1 n + 1 (n + 1) (n + 1) · · · (n + 1) xn ((n + 1) − 1)/n für verschiedene Werte von n mittels Gauß’scher Elimination zu lösen. Die auftretende Matrix heißt Vandermonde-Matrix. Die exakte Lösung ist durch xi = 1 für i = 1, 2, . . . , n gegeben. Auf einem Computer mit ca. 16 Stellen Genauigkeit erhalten wir das exakte Resultat für n ≤ 14. Für n = 15 jedoch erhalten wir plötzlich folgendes Ergebnis (gerundet 33 auf 2 Stellen): x1 ≈ −894.97, x2 ≈ 2074.34, x3 ≈ −2138.14, x4 ≈ 1309.97, x5 ≈ −531.35, x6 ≈ 153.62, x7 ≈ −30.88, x8 ≈ 5.94, x9 ≈ 0.43, x10 ≈ 1.05, x11 ≈ . . . ≈ x15 ≈ 1.00 Dieses Ergebnis ist völlig wertlos und hat nichts mit der Lösung zu tun. Für n > 15 wird es noch schlimmer. Was ist der Grund für dieses seltsame Verhalten in Beispiel 2.4.1? Offenbar spielen Rundungsfehler eine Rolle, da ansonsten bei exakter Rechnung die exakte Lösung gefunden werden müsste. Wir interessieren uns nun für die Frage, wie Fehler in A und b (etwa Rundungsfehler, die beim Gauß’schen Eliminationsverfahren auftreten) die Lösung x von Ax = b beeinflussen. Wir werden sehen, dass die sogenannte Konditionszahl einer Matrix eine entscheidende Rolle spielt. Definition 2.4.2 (Vektornorm) Sei V ein Vektorraum über K. Die Abbildung k · k : V → R heißt Norm, wenn sie die folgenden Eigenschaften besitzt: 1. kxk = 0 gilt genau dann, falls x = 0 ist, 2. kcxk = |c| kxk für alle x ∈ V und c ∈ K, 3. kx + yk ≤ kxk + kyk für alle x, y ∈ V . Beispiel 2.4.3 (Vektornormen) Häufig verwendete Vektornormen auf Kn mit x = (x1 , . . . , xn )> sind (i) die Maximumnorm: kxk∞ := max |xi |; 1≤i≤n (ii) die Betragssummennorm: kxk1 := n X |xi |; i=1 v u n uX (iii) die Euklidische Norm: kxk2 := t |xi |2 . i=1 (iv) p-Norm: kxkp := n X !1/p |xi |p für p ∈ N. i=1 Bemerkung 2.4.4 In endlichdimensionalen Vektorräumen sind alle Normen äquivalent, d.h. für beliebige 34 Normen k · k(1) und k · k(2) über demselben Vektorraum V gibt es Konstanten c1 , c2 , so dass für alle x ∈ V gilt c1 kxk(2) ≤ kxk(1) ≤ c2 kxk(2) . Für unendlichdimensionale Vektorräume gilt dies nicht! Formal kann man einer Vektornorm k · kV eine passende“ Matrixnorm k · kM zuordnen. ” Definition 2.4.5 (Matrixnorm) Sei k · kV eine Vektornorm. Die durch k · kV induzierte Matrixnorm k · kM auf dem Raum aller linearen Abbildungen A : Kn → Km ist definiert als kAxkV = sup kAxkV = sup kAxkV , kxkV 6=0 kxkV kxkV =1 kxkV ≤1 kAkM := sup (2.12) wobei x ∈ Kn und A ∈ Km×n gilt. Beachte: In Definition 2.4.5 wird in Kn und Km dieselbe Vektornorm verwendet! Allgemeiner könnte man auch unterschiedliche Vektornormen verwenden, was wir aus Gründen der Übersichtlichkeit hier nicht tun. Die Wohldefiniertheit der Matrixnorm ergibt sich aus der Stetigkeit der Abbildung x 7→ kAxkV und der Kompaktheit der Menge {x ∈ Kn | kxkV = 1}. Dass k · kM tatsächlich eine Norm ist, ergibt sich direkt aus den Normeigenschaften von k · kV und kann leicht nachgerechnet werden. Die so definierte Matrixnorm hat zwei schöne Eigenschaften: Satz 2.4.6 Sei k · kV eine Vektornorm und k · kM die induzierte Matrixnorm. Dann gelten: (a) Verträglichkeit, d.h. es gilt kAxkV ≤ kAkM · kxkV für alle x ∈ Kn , A ∈ Km×n . (b) Submultiplikativität, d.h. es gilt kA · BkM ≤ kAkM · kBkM für alle A ∈ Km×n , B ∈ Kn×q . Beweis: (a) Für x = 0 ist die Behauptung klar. Für kxkV 6= 0 gilt kAxkV = A kxkx V · kxkV V und somit ! ! x kAxk V · kxkV = sup kAxkV ≤ sup · kxkV = kAkM · kxkV . A kxkV kxkV 6=0 kxkV kxkV 6=0 V (b) Nach (a) gilt kA · BxkV ≤ kAkM · kBxkV und somit kA · BxkV kBxkV ≤ kAkM · sup = kAkM · kBkM . kxkV kxkV 6=0 kxkV 6=0 kxkV kA · BkM = sup 35 2 Die konkrete Berechnung einer induzierten Matrixnorm für eine gegebene Vektornorm ist mitunter kompliziert. Wir fassen gänge Matrixnormen zusammen. Satz 2.4.7 (Zeilensummennorm, Spaltensummennorm, Spektralnorm) (a) Die durch die Vektornorm k · k∞ induzierte Matrixnorm ist die sogenannte Zeilenn X summennorm kAk∞ := max |aij |. 1≤i≤m j=1 (b) Die durch die Betragssummennorm k · k1 induzierte Matrixnorm ist die sogenannte m X Spaltensummennorm kAk1 := max |aij |. 1≤j≤n i=1 (c) Die durch die euklidische Vekotrnorm k · k2 induzierte Matrixnorm ist die Spektralnorm p kAk2 = ρ(A∗ A), wobei ρ(A∗ A) den Spektralradius von A∗ A bezeichnet, d.h. ρ(A∗ A) = max{|λ| | λ ist Eigenwert von A∗ A}. Bemerkung 2.4.8 Es gibt weitere Matrixnormen, die aber im Allgemeinen. nicht durch Vektornormen induziert sind: • Die Matrixnorm kAk = max 1≤i≤m,1≤j≤n |aij | ist durch keine Vektornorm induziert, da sie nicht submultiplikativ ist wie folgendes Beispiel zeigt: ! ! 1 1 2 2 A=B= , A·B = , 2 = kA · Bk > kAk · kBk = 1. 1 1 2 2 • Die Frobenius-Norm s X kAkF := |aij |2 , i=1,...,m,j=1,...n ist verträglich und submultiplikativ bezüglich der euklidischen Vektornorm, aber sie wird nicht durch diese induziert. Mit diesen Vorbereitungen untersuchen wir den Fehlereinfluss in linearen Gleichungssystemen. Hierbei sei x die Lösung von Ax = b. 36 Störungen in der rechten Seite b: Wie wirken sich Störungen ∆b in der rechten Seite b in (2.3) auf die Lösung x + ∆x des gestörten Gleichungssystems A(x + ∆x) = b + ∆b aus? Für die Abweichung der gestörten Lösung von der exakten Lösung gilt bezüglich einer Vektornorm k · kV und der zugeordneten Matrixnorm k · kM die Beziehung k∆xkV = kx + ∆x − xkV = kA−1 (b + ∆b) − A−1 bkV = kA−1 (b + ∆b − b)kV ≤ kA−1 kM · k∆bkV . Mit Hilfe der absoluten Konditionszahl κa = kA−1 kM ergibt sich für den absoluten Fehler die Abschätzung k∆xkV ≤ κa · k∆bkV . (2.13) Mit (2.13) und der Verträglichkeit folgt für kxkV , kbkV 6= 0 die Abschätzung kA−1 kM · k∆bkV k∆bkV k∆xkV k∆xkV = kA−1 kM · ≥ ≥ . kbkV kAxkV kAxkV kAkM · kxkV Umformen ergibt eine Abschätzung des relativen Fehlers k∆xkV k∆bkV ≤ kAkM · kA−1 kM · . kxkV kbkV (2.14) Definition 2.4.9 (Kondition einer Matrix) Die Zahl κM (A) := kAkM · kA−1 kM heißt Konditionszahl der Matrix A bzgl. der Matrixnorm k · kM . A heißt schlecht konditioniert, falls κM (A) sehr groß ist. Die Kondition erfüllt wegen 1 = kIkM = kA−1 AkM ≤ kA−1 kM · kAkM = κM (A) stets κM (A) ≥ 1 . Mit dieser Definition lässt sich (2.14) schreiben als k∆xkV k∆bkV ≤ κM (A) · . kxkV kbkV 37 Die Kondition κM (A) einer Matrix A macht also eine Aussage darüber, wie stark sich Störungen in der rechten Seite auf das Ergebnis auswirken können. Insbesondere bei schlecht konditionierten Matrizen mit κM (A) >> 1 können sich bereits kleine Störungen in b sehr stark auf das Ergebnis auswirken. Störungen in der Matrix A und der rechten Seite b: Es gilt folgender Satz, der hier nicht bewiesen werden soll. Satz 2.4.10 Es gelte Ax = b und (A + ∆A)(x + ∆x) = b + ∆b. Dann gilt k∆xkV κM (A) k∆AkM k∆bkV ≤ · + , kxkV kAkM kbkV 1 − κM (A) · k∆AkM kAkM falls κM (A) · k∆AkM kAkM < 1 und kxkV 6= 0 ist. Beispiel 2.4.11 (Vandermonde-Matrix) Nun können wir den Grund für das seltsame Verhalten in Beispiel 2.4.1 angeben. Die Vandermonde-Matrix ist mit κM (A) > 1021 eine extrem schlecht konditionierte Matrix. Faustregel: Falls κM (A) = 10k , dann muss man mit dem Verlust von mindestens k Stellen Genauigkeit bei der Lösung von Ax = b rechnen. Wir untersuchen noch einige Eigenschaften der Kondition bzgl. der Spektralnorm k · k2 . Diese Norm bietet einige Vorteile, da für unitäre Matrizen Q ∈ Kn×n für beliebige x ∈ Kn die Beziehung kQxk22 = x∗ Q∗ Qx = x∗ x = kxk22 bzw. kQxk2 = kxk2 gilt. Da mit Q auch Q−1 unitär ist, folgt ebenso kQ−1 xk2 = kxk2 . Aus der Definition einer Matrixnorm folgt damit kQk2 = kQ−1 k2 = 1 und somit κ2 (Q) = 1. Unitäre Matrizen haben also eine bzgl. der Spektralnorm optimale Kondition. Dies ist der Grund, warum wir später mit orthogonalen Householder-Matrizen arbeiten werden, vgl. Abschnitt 3.2.3. 2.5 Iterative Lösung Der Aufwand zur Berechnung der LR-Zerlegung und der Cholesky-Zerlegung (und später der QR-Zerlegung) für quadratische Matrizen ist jeweils von der Ordnung O(n3 ). Der Aufwand wächst also kubisch mit der Dimension n des Gleichungssystems. Für sehr große 38 Gleichungssysteme mit teilweise mehr als einer Milliarde Gleichungen, wie sie in praktischen Anwendungen durchaus entstehen können, ist der Aufwand damit zu groß. Es werden alternative Methoden benötigt, die weniger aufwendig sind. Beispiel 2.5.1 (Diskretisierung einer partiellen Differentialgleichung) Große, dünn besetzte Gleichungssysteme entstehen bei Diskretisierungsverfahren zur approximativen Lösung partieller Differentialgleichungen. Als Beispiel sei die 2-dimensionale partielle Differentialgleichung −uxx (x, y) − uyy (x, y) = f (x, y), u(x, y) = 0, (x, y) ∈ Ω := (0, 1) × (0, 1), (x, y) ∈ ∂Ω genannt. Hierbei ist eine Funktion u : Ω → R gesucht, die die obige partielle Differentialgleichung erfüllt. Eine gängige Vorgehensweise besteht in der Diskretisierung der Differentialgleichung auf dem äquidistanten Gitter G := {(xi , yj ) | xi = ih, yj = jh, 0 ≤ i, j ≤ N }, h = 1/N. Approximation der zweiten partiellen Ableitungen durch ui+1,j − 2ui,j + ui−1,j , h2 ui,j+1 − 2ui,j + ui,j−1 uyy (xi , yj ) ≈ , h2 uxx (xi , yj ) ≈ 1 ≤ i, j ≤ N − 1, 1 ≤ i, j ≤ N − 1, mit ui,j ≈ u(xi , yj ) liefert das lineare Gleichungssystem Au = b für den Vektor 2 u = (u1,1 , u1,2 , . . . , u1,N −1 , u2,1 , u2,2 , . . . , u2,N −1 , . . . , uN −1,1 , uN −1,2 , . . . , uN −1,N −1 )> ∈ R(N −1) , mit 2 b = h2 (f1,1 , f1,2 , . . . , f1,N −1 , f2,1 , f2,2 , . . . , f2,N −1 , . . . , fN −1,1 , fN −1,2 , . . . , fN −1,N −1 )> ∈ R(N −1) , fi,j := f (xi , yj ), M1 D2 D2 M2 D3 .. .. .. A= . . . DN −2 MN −2 DN −1 DN −1 MN −1 , 4 −1 −1 4 −1 2 .. .. .. Mi = ∈ R(N −1) . . . −1 4 −1 −1 4 und Di = −I ∈ R(N −1)×(N −1) . Die Matrix A ist sehr groß (Dimension wächst quadratisch mit N ), aber sehr dünn besetzt (pro Zeile sind maximal 5 Einträge 6= 0). Zudem kann man zeigen, dass sie positiv definit und symmetrisch ist. 39 In der Praxis werden nur die von Null verschiedenen Einträge der Matrix A und die zugehörigen Spalten- und Zeilenpositionen gespeichert. Damit kann man Multiplikationen der Form A · d mit einem Vektor d sehr effizient berechnen (der Aufwand hängt linear von der Anzahl der nichtverschwindenden Einträge in A ab). Wir werden im Folgenden sehen, dass die sogenannten iterativen Verfahren lediglich Matrix-Vektor Multiplikationen benötigen. Die Idee der iterativen Verfahren zur Lösung von linearen Gleichungssystemen (2.3) basiert auf der geschickten Zerlegung der Matrix A in A=M −N (2.15) mit Matrizen M und N , die noch genauer spezifiziert werden müssen. Es wird angenommen, dass M regulär ist und leicht“ invertiert werden kann. Dann lässt sich das lineare ” Gleichungssystem b = Ax = (M − N )x = M x − N x umformen zu x = M −1 N x + M −1 b. (2.16) Mit T := M −1 N und c := M −1 b ist dies eine Fixpunktgleichung für x: x = T x + c. (2.17) Damit ist jede Lösung der Fixpunktgleichung (2.17) auch Lösung des linearen Gleichungssystems (2.3) und umgekehrt. Definition 2.5.2 (Fixpunkt) x̂ heißt Fixpunkt der Abbildung g : Kn → Kn , x 7→ g(x), wenn x̂ = g(x̂) gilt. Eine einfache Idee zur Berechnung eines Fixpunkts x̂ von g besteht darin, ausgehend von einem Startwert x(0) eine Fixpunktiteration x(i+1) = g(x(i) ), i = 0, 1, 2, . . . durchzuführen, indem das Ergebnis immer wieder in die Funktion eingesetzt wird. Auf diese Weise erhält man eine Folge x(i) . Es stellt sich natürlich die Frage, ob die Folge auch tatsächlich gegen den gesuchten Fixpunkt x̂ konvergiert. Beispiel 2.5.3 40 (i) Wir wollen einen Fixpunkt x̂ mit x̂ = cos(x̂) bestimmen. Wir führen die Fixpunktiteration ausgehend von x(0) = 1 für die Funktion g(x) = cos(x) durch. Die Fixpunktiteration lässt sich ganz einfach mit dem Taschenrechner durchführen, indem 1 eingetippt wird und anschließend immer wieder auf cos gedrückt wird. Nach 35-maligem Drücken der cos-Taste erhält man mit sechs Nachkommastellen x̂ = 0.739085. Weiteres Drücken der cos-Taste ändert das Ergebnis nicht mehr. Die Folge konvergiert offensichtlich. (ii) Jetzt wollen wir einen Fixpunkt von g(x) = exp(x) mit Startwert x(0) = 1 berechnen. Hier das Ergebnis: x(1) = 2.71828183 x(2) = 15.15426223 x(3) = 0.38142791 · 107 x(4) = 0.22316774 · 101656521 Der nächste Iterationsschritt liefert einen Exponentenüberlauf. Die Folge konvergiert offensichtlich nicht. Beispiel 2.5.3 zeigt, dass die Fixpunktiteration nicht immer sinnvolle Resultate liefert. Offensichtlich hängt die Konvergenz von Eigenschaften der Funktion g ab. Ein hinreichendes Kriterium für die Konvergenz der Fixpunktfolge liefert Satz 2.5.4 (Banachscher Fixpunktsatz) Sei eine Selbstabbildung g : D → D mit abgeschlossener Menge D ⊆ Kn gegeben. Für alle x, y ∈ D sei die Kontraktionsbedingung kg(x) − g(y)k ≤ qkx − yk mit q < 1 (2.18) für eine geeignete Norm k · k erfüllt. Dann gilt: (a) g besitzt genau einen Fixpunkt x̂ in D. (b) Die Fixpunktiteration x(i+1) = g(x(i) ), i = 0, 1, 2, . . . konvergiert für jeden Startwert x(0) aus D gegen den Fixpunkt x̂. (c) Es gilt die a-priori Fehlerabschätzung kx(i) − x̂k ≤ qi kx(1) − x(0) k. 1−q 41 Beweis: Auf Grund der Selbstabbildungseigenschaft ist die Fixpunktiteration wohldefiniert. (i) Existenz eines Fixpunkts: Sei x(0) ∈ D beliebig. Wir zeigen, dass die Folge eine Cauchyfolge ist. Die Kontraktionseigenschaft von g liefert für jedes i ∈ N die Beziehung kx(i+1) − x(i) k ≤ qkx(i) − x(i−1) k ≤ . . . ≤ q i kx(1) − x(0) k. Mit Hilfe der Dreiecksungleichung folgt dann für jedes k ∈ N die Abschätzung kx(i+k) − x(i) k ≤ kx(i+k) − x(i+k−1) k + kx(i+k−1) − x(i+k−2) k + . . . + kx(i+1) − x(i) k ≤ q i+k−1 + q i+k−2 + . . . + q i kx(1) − x(0) k qi kx(1) − x(0) k. ≤ 1−q Wegen q < 1 ist {x(i) } eine Cauchyfolge in D. Da D abgeschlossen ist, konvergiert die Folge dann gegen ein x̂ ∈ D. Gemäß der Kontraktionsbedingung ist g Lipschitz-stetig und insbesondere stetig. Damit gilt x̂ = lim x(i+1) = lim g(x(i) ) = g( lim x(i) ) = g(x̂), i→∞ i→∞ i→∞ d.h. x̂ ist Fixpunkt. (ii) Eindeutigkeit: Angenommen, es gäbe einen weiteren Fixpunkt x̃ 6= x̂ in D. Dann gilt kx̃ − x̂k = kg(x̃) − g(x̂)k ≤ qkx̃ − x̂k. Wegen q < 1 kann dies nur für kx̃ − x̂k = 0 gelten im Widerspruch zu x̃ 6= x̂. (iii) Die Fehlerabschätzung folgt mit (i) aus kx(i) − x̂k ≤ kx(i+1) − x̂k + kx(i+1) − x(i) k = kg(x(i) ) − g(x̂)k + kx(i+1) − x(i) k (i) ≤ qkx(i) − x̂k + q i kx(1) − x(0) k. 2 Anwendung der Fixpunktiteration auf die Fixpunktgleichung (2.17) liefert das folgende iterative Verfahren zur Lösung des linearen Gleichungssystems Ax = b: Algorithmus 2.5.5 (Iteratives Lösungsverfahren für Ax = b ) 42 (i) Zerlege die Matrix A gemäß (2.15) mit einer invertierbaren Matrix M , wähle einen beliebigen Startvektor x(0) und setze k := 0. (ii) Falls rk = kb − Ax(k) k2 hinreichend klein: STOP (iii) Berechne M x(k+1) = N x(k) + b (2.19) (iv) Setze k := k + 1 und gehe zu (ii). Bemerkung 2.5.6 • (2.19) ist unter den bisherigen Annahmen äquivalent mit (2.17). • Konvergiert die Folge {x(k) } gegen x̂, dann gilt M x̂ = N x̂ + b = (M − A)x̂ + b = M x̂ − Ax̂ + b und somit löst x̂ das lineare Gleichungssystem Ax̂ = b. Unter Anwendung des Banachschen Fixpunktsatzes lautet ein hinreichendes Konvergenzkriterium wie folgt. Satz 2.5.7 (Hinreichendes Konvergenzkriterium) Es gelte kT kM < 1 (2.20) bezüglich einer mit der Vektornorm k · kV verträglichen Matrixnorm k · kM . Dann existiert genau ein Fixpunkt x̂ der Funktion g(x) = T x + c und die Fixpunktiteration x(i+1) = T x(i) + c, i = 0, 1, 2, . . . (2.21) konvergiert für jeden Startwert x(0) ∈ Kn gegen x̂. Beweis: Aus der Ungleichung kT x + c − (T y + c)kV = kT (x − y)kV ≤ kT kM · kx − ykV folgt mit kT kM < 1 die Kontraktionseigenschaft (2.18). Da die Abbildung x 7→ T x+c affin linear ist, liegt eine Selbstabbildung mit D = Kn vor. Aus dem Banachschen Fixpunktsatz folgt die Behauptung. 2 Die bisherige Konvergenzaussage hängt noch von den gewählten Normen k · kV und k · kM ab und setzt voraus, dass solche Normen mit (2.20) gefunden werden können. Ein notwendiges und hinreichendes Kriterium für die Konvergenz von (2.21) liefert 43 Satz 2.5.8 Sei A = M − N invertierbar und T = M −1 · N . Dann konvergiert (2.21) genau dann für jeden Startwert x(0) ∈ Rn gegen die Lösung x̂ von Ax = b, wenn für den Spektralradius ρ(T ) von T die Ungleichung ρ(T ) < 1 gilt. Es stellt sich noch die Frage, wie schnell die so konstruierten Näherungen x(i) gegen den Fixpunkt x̂ konvergieren. Wegen kx(i+1) − x̂kV = kT x(i) + c − (T x̂ + c)kV = kT (x(i) − x̂)kV ≤ kT kM · kx(i) − x̂kV mit kT kM < 1 erhalten wir im Allgemeinen nur eine sogenannte lineare Konvergenzgeschwindigkeit für die Fixpunktiteration, da sich der Fehler im schlechtesten Fall in jeder Iteration nur um den Faktor kT kM verringert. Insbesondere kann die Konvergenz für kT kM ≈ 1 sehr langsam sein. 2.5.1 Spezielle Verfahren Es werden spezielle Verfahren durch geeignete Wahl von M und N konstruiert. Die Fixpunktiteration (2.21) ist besonders dann günstig, falls M leicht invertiert werden kann und die auftretenden Multiplikationen M −1 N und M −1 b billig sind, wie es häufig bei dünn besetzten Matrizen der Fall ist. Wir zerlegen A in A=D−L−R wobei D die Diagonalelemente von A, L die Elemente unterhalb der Diagonalen von −A und R die Elemente oberhalb der Diagonalen von −A enthalten. 2.5.1.1 Gesamtschrittverfahren (Jacobi-Verfahren) Wir wählen in (2.15) M = D, N = L + R und erhalten den Gesamtschrittoperator T = D−1 (L + R). Die Fixpunktiteration (2.21) für das Gesamtschrittverfahren lautet somit x(i+1) = D−1 (L + R)x(i) + b . (2.22) Die Inverse der Diagonalmatrix D = diag(a11 , a22 , . . . , ann ) ist sehr leicht zu berechnen. Es ist D−1 = diag(1/a11 , 1/a22 , . . . , 1/ann ). Um das Verfahren durchführen zu können, müssen die Diagonalelemente ungleich Null sein. Algorithmus 2.5.9 (Gesamtschrittverfahren (Jacobi-Verfahren)) 44 (i) Wähle M = diag(a11 , a22 , . . . , ann ). (ii) Löse (2.19) bzgl. x(i+1) : (i+1) xj = 1 ajj ! bj − X (i) ajk xk , j = 1, . . . , n. k6=j Ein einfach überprüfbares Kriterium zur Konvergenz des Verfahrens liefert der folgende Satz. Satz 2.5.10 (Konvergenz des Gesamtschrittverfahrens) Die Matrix A ∈ Kn×n sei strikt diagonaldominant, d.h. es gelte |aii | > n X |aij |, i = 1, . . . , n. j=1 j6=i Dann konvergiert das Gesamtschrittverfahren für jeden Startvektor x(0) ∈ Rn gegen die Lösung x̂ von Ax = b und es gilt die Fehlerabschätzung kx(i) − x̂k∞ ≤ (kT k∞ )i kx(1) − x(0) k∞ , 1 − kT k∞ i = 1, 2, . . . , wobei T = D−1 (L + R) den Gesamtschrittoperator bezeichnet. Beweis: Die strikte Diagonaldominanz impliziert kT k∞ < 1, wobei k·k∞ die Zeilensummennorm bezeichnet. Die Konvergenz folgt dann aus Satz 2.5.7. Die Fehlerabschätzung folgt aus dem Banach’schen Fixpunktsatz. 2 Ein analoger Satz gilt für die Spaltensummennorm k · k1 , wenn gefordert wird, dass |ajj | > n X |aij |, j = 1, . . . , n. i=1 i6=j 2.5.1.2 Einzelschrittverfahren (Gauß-Seidel-Verfahren) Wir wählen in (2.15) M = D − L, N = R und erhalten den Einzelschrittoperator T = (D − L)−1 R. Die Fixpunktiteration (2.21) für das Einzelschrittverfahren lautet somit x(i+1) = (D − L)−1 Rx(i) + b . Algorithmus 2.5.11 (Einzelschrittverfahren (Gauß-Seidel-Verfahren)) 45 (2.23) (i) Wähle M als den unteren Dreiecksanteil von A, einschließlich Diagonale. (ii) Löse (2.19) bzgl. x(i+1) mittels Vorwärtssubstitution: (i+1) xj 1 = ajj bj − j−1 X (i+1) ajk xk − k=1 ! n X (i) ajk xk , j = 1, . . . , n. k=j+1 (i+1) Man beachte, dass die auf der rechten Seite auftretenden Komponenten xk , k < j der (i+1) neuen Iterierten bereits zur Berechnung von xj verwendet werden. Ein einfach überprüfbares Kriterium zur Konvergenz des Verfahrens liefert der folgende Satz. Satz 2.5.12 (Konvergenz des Einzelschrittverfahrens) Die Matrix A ∈ Kn×n erfülle das Sassenfeld-Kriterium: (i) Sämtliche Diagonalelemente von A sind von Null verschieden. (ii) Es gilt β = max βi < 1 mit i=1,...,n n β1 = βi 1 X |a1j |, |a11 | j=2 1 = |aii | i−1 X |aij |βj + j=1 n X ! |aij | , i = 2, . . . , n. j=i+1 Dann konvergiert das Einzelschrittverfahren für jeden Startvektor x(0) ∈ Kn gegen die Lösung x̂ von Ax = b und es gilt die Fehlerabschätzung kx(i) − x̂k∞ ≤ βi kx(1) − x(0) k∞ , 1−β i = 1, 2, . . . . Beweis: Wir zeigen, dass kT k∞ ≤ β für den Einzelschrittoperator T = (D − L)−1 R gilt. Sei x ∈ Kn mit kxk∞ ≤ 1 und y = T x, d.h. (D − L)y = Rx. Dann folgt a11 y1 = − n X j=2 n a1j xj ⇒ 1 X |a1j | = β1 . |y1 | ≤ |a11 | j=2 Die Behauptung folgt nun per Induktion. Sei also |yi | ≤ βi für i = 1, 2, . . . , j − 1 bewiesen. Dann folgt ! j−1 j−1 n n X X X X 1 |ajk |βk + |ajk | = βj . ajk yk + ajj yj = − ajk xk ⇒ |yj | ≤ |ajj | k=1 k=1 k=j+1 k=j+1 46 Somit gilt also |yi | ≤ βi für i = 1, 2, . . . , n. Damit gilt kT k∞ = sup kT xk∞ ≤ β. Die kxk∞ ≤1 Konvergenz folgt dann aus Satz 2.5.7. Die Fehlerabschätzung folgt aus dem Banach’schen Fixpunktsatz. 2 Beispiel 2.5.13 (Diskretisierung einer partiellen Differentialgleichung, Teil II) Wir betrachten wiederum Beispiel 2.5.1 mit N = 50 und f (x, y) = 1000 · sin((x − 0.5) · (y − 0.5)). Für N = 50 hat A im Gleichungssystem Au = b die Dimension 492 = 2401 und von den 24012 = 5764801 Einträgen sind lediglich 11809 Einträge von Null verschieden. Wir lösen das Gleichungssystem mit dem Gauß-Seidel-Verfahren, wählen u(0) = 0 als Startwert und brechen die Iteration ab, wenn das relative Residuum folgende Bedingung erfüllt: kb − Au(k) k2 rk = ≤ 10−5 . r0 kb − Au(0) k2 Nach 874 Iterationen ist dies der Fall und wir erhalten die numerische Lösung u der partiellen Differentialgleichung: 1.5 1.5 1 0.5 1 0 0.5 -0.5 0 -1 -0.5 -1.5 -1 -1.5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 Beispiel 2.5.14 Wir lösen das lineare Gleichungssystem 2 −1 0 x1 2 1 6 −2 x2 = −4 4 −3 8 x3 5 47 0.1 0 mit dem Gauß-Seidel-Verfahren und erhalten nach 22 Iterationen die Näherungslösung x = (0.62, −0.76, 0.03)> . Man erkennt sehr gut, dass der Fortschritt in den ersten Iterationen sehr groß ist, die Konvergenz des Verfahrens am Ende allerdings sehr langsam ist. k 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 rk = kb − Ax(k) k22 6.7e-00 9.1e-01 4.2e-01 8.5e-02 3.3e-02 7.9e-03 2.6e-03 7.2e-04 2.0e-04 6.5e-05 1.5e-05 5.8e-06 1.1e-06 5.0e-07 8.2e-08 4.3e-08 6.3e-09 3.7e-09 5.2e-10 3.0e-10 4.5e-11 2.5e-11 4.1e-12 1 "error.dat" u 1:2 0.9 0.8 0.7 rk 0.6 0.5 0.4 0.3 0.2 0.1 0 0 5 10 15 20 Iteration Häufig beobachtet man in der Praxis, dass das Einzelschrittverfahren schneller konvergiert als das Gesamtschrittverfahren. 2.6 Das konjugierte Gradientenverfahren Ein gänzlich anderer Ansatz zur Lösung des linearen Gleichungssystems (2.3) wird beim konjugierten Gradientenverfahren verfolgt. Wir beschränken uns hier auf den reellen Fall K = R. 48 25 Lemma 2.6.1 Seien A ∈ Rn×n symmetrisch und positiv definit und b ∈ Rn . Dann gilt: x̂ löst das lineare Gleichungssystem Ax = b genau dann, wenn x̂ die Funktion 1 f (x) = x> Ax − b> x 2 minimiert. Beweis: Setzt man 1 1 g(x) = f (x) + b> A−1 b = (Ax − b)> A−1 (Ax − b), 2 2 so folgt: x̂ löst Ax = b genau dann, wenn x̂ die Funktion g minimiert, da A−1 positiv definit ist. Da sich f und g nur um eine Konstante unterscheiden, sind die Minima von g und f gleich, und die Behauptung ist bewiesen. 2 Das konjugierte Gradientenverfahren berechnet nun ausgehend von einer Startschätzung x(0) die nächste Iterierte x(k+1) = x(k) + αk d(k) , k = 0, 1, 2, . . . mit Hilfe einer geeignet gewählten Suchrichtung d(i) und einer optimalen Schrittweite αi . Satz 2.6.2 (Konjugierte Gradienten Verfahren) Seien f (x) = 12 x> Ax − b> x mit einer symmetrisch, positiv definiten Matrix A ∈ Rn×n , b ∈ Rn und x(0) ∈ Rn . Seien d(0) , d(1) , . . . , d(n−1) ∈ Rn \ {0} A-konjugiert, d.h. (d(i) )> Ad(j) = 0 ∀i, j ∈ {0, 1, . . . , n − 1}, i 6= j. Sei {x(k) } berechnet aus dem Verfahren der sukzessiven eindimensionalen Minimierung, d.h. x(k+1) = x(k) + αk d(k) , k = 0, 1, 2, . . . wobei die Schrittweite aus f (x(k) + αk d(k) ) = min f (x(k) + αd(k) ) α∈R bestimmt wird. Dann bricht das Verfahren nach spätestens n Schritten mit x(n) als dem Minimum von f ab. 2 Beweis: Setzt man φ(α) := f (x(k) + αd(k) ), so minimiert αk die quadratische Funktion φ(α) genau dann, wenn φ0 (αk ) = 0 ist. Also wenn 0 = ∇f (x(k) + αk d(k) )> d(k) = (Ax(k) − b + αk Ad(k) )> d(k) 49 (2.24) gilt. Definieren wir (den Gradienten) g (k) := Ax(k) − b. (2.25) Durch Umstellen erhalten wir, da (d(k) )> Ad(k) wegen der positiven Definitheit von A und d(k) 6= 0 positiv ist, −(g (k) )> d(k) αk = (k) > (k) . (2.26) (d ) Ad Mit Hilfe der Definition von x(k+1) erhalten wir (g (k+1) )> d(k) = (Ax(k+1) − b)> d(k) = (Ax(k) − b + αk Ad(k) )> d(k) = 0. Für alle i 6= j erhalten wir, unter Ausnutzung der A-konjugierten Richtungen, (g (i+1) − g (i) )> d(j) = (Ax(i+1) − Ax(i) )> d(j) = (αi d(i) )> Ad(j) = 0. Daraus folgt für alle j = 0, 1, . . . , k (g (k+1) )> d(j) = (g (j+1) )> d(j) + k X (g (i+1) − g (i) )> d(j) = 0 (2.27) i=j+1 Da die Richtungen d(0) , d(1) , . . . , d(n−1) ∈ Rn \ {0} A-konjugiert sind, und die Matrix A symmetrisch positiv definit ist, sind sie orthogonal bezüglich des Skalarprodukts hx, yiA := x> Ay. Orthogonale Vektoren sind linear unabhängig. Wegen (2.27) ist g (n) orthogonal (bezüglich des Euklidischen Skalarprodukts) zu d(0) , d(1) , . . . , d(n−1) . Damit muss aber 0 = g (n) = Ax(n) − b 2 gelten und x(n) minimiert f. Der Satz liefert mit anderen Worten die Darstellung x̂ = x (0) + n−1 X αk d(k) k=0 für eine Lösung x̂ von Ax = b. Nun bleibt die Frage, wie man A-konjugierte Richtungen erzeugt. Dabei sollten die Richtungen möglichst zu einem Abstieg im Funktionswert von f führen, d.h. sie sollten ∇f (x(k) )> d(k) < 0 50 für alle k erfüllen. Wir setzen daher d(0) = −∇f (x(0) ). Sind nun bereits ` + 1 Vektoren d(0) , d(1) , . . . , d(`) erzeugt, so dass (d(i) )> Ad(j) = 0 ∀i, j = 0, 1, . . . , `, i 6= j, ∇f (x(i) )> d(i) < 0 i = 0, 1, . . . , `, so ist g (`+1) nach (2.27) linear unabhängig zu d(0) , d(1) , . . . , d(`) . Wir bestimmen nun d(`+1) mit dem Gram-Schmidt-Orthogonalisierungsansatz: d (`+1) = −g (`+1) + ` X βi` d(i) , (2.28) i=0 wobei die Koeffizienten βi` aus der Forderung (d(`+1) )> Ad(j) = 0 für alle j = 0, 1, . . . , ` zu βi` = (g (`+1) )> Ad(i) (d(i) )> Ad(i) bestimmt werden. Damit ist die Konstruktion der A-konjugierten Richtungen abgeschlossen. Es lassen sich jedoch noch effizientere update-Formeln herleiten. Aus (2.28) folgt unter Verwendung von (2.27) (g (`+1) )> d(`+1) = −kg (`+1) k22 < 0 und somit auch α`+1 = (2.29) kg (`+1) )k22 −(g (`+1) )> d(`+1) = > 0. (d(`+1) )> Ad(`+1) (d(`+1) )> Ad(`+1) Ferner folgt aus (2.28) (mit j statt ` + 1) (g (`+1) )> g (j) = (g (`+1) )> j−1 X ! βij−1 d(i) − d(j) ∀j = 0, 1, . . . , `, i=0 und somit unter Verwendung von (2.27) (g (`+1) )> g (j) = 0 ∀j = 0, 1, . . . , `. Die Orthogonalität impliziert für alle j = 0, 1, . . . , ` − 1 1 (`+1) > (j+1) (g ) (g − g (j) ) αj 1 (`+1) > (g ) (Ax(j+1) − Ax(j) ) = αj 0 = = (g (`+1) )> Ad(j) , und nach Definition der Koeffizienten βj` erhalten wir βj` = 0 ∀j = 0, 1, . . . , ` − 1. 51 (2.30) Für den verbleibenden Koeffizienten nutzen wir g (`+1) − g (`) = Ax(`+1) − Ax(`) = α` Ad(`) und erhalten β`` = (2.26) = = (2.29),(2.30) = (g (`+1) )> Ad(`) (d(`) )> Ad(`) (g (`+1) )> Ad(`) α` −(g (`) )> d(`) (g (`+1) )> (g (`+1) − g (`) ) −(g (`) )> d(`) kg (`+1) k22 =: β` . kg (`) k22 Damit folgt für die neue A-konjugierte Richtung d(`+1) = −g (`+1) + β` d(`) . Nun können wir das konjugierte Gradienten-Verfahren (CG-Verfahren) angeben: Algorithmus 2.6.3 (CG-Verfahren) (i) Wähle einen Startvektor x(0) ∈ Rn , ε ≥ 0, berechne g (0) := Ax(0) − b und setze d(0) := −g (0) , k := 0. (ii) Falls kg (k) k2 ≤ ε, STOP. (iii) Berechne kg (k) k22 αk := x(k+1) , > (d(k) ) Ad(k) := x(k) + αk d(k) , g (k+1) := g (k) + αk Ad(k) , kg (k+1) k22 , βk := kg (k) k22 d(k+1) := −g (k+1) + βk d(k) . (iv) Setze k := k + 1, und gehe zu (ii). In der Praxis führt man häufig aufgrund von auftretenden Rundungsfehlern nach einer gewissen Anzahl von Iterationen einen Restart des Verfahrens durch. Obwohl das CG-Verfahren (zumindest theoretisch) nach maximal n Schritten terminiert, ist es in der Praxis aus Aufwandsgründen oft nicht möglich, n Iterationen tatsächlich durchzuführen (n kann sehr groß sein, etwa 1 Milliarde). 52 Daher ist es sinnvoll, das CG-Verfahren als iteratives Verfahren zu betrachten. Es stellt sich dann die Frage nach der Konvergenzgeschwindigkeit. Es lässt sich folgende Abschätzung beweisen (vgl. Kelley [Kel95, S. 17]): !i p κ(A) − 1 (i) kx − x̂kA ≤ 2 p kx(0) − x̂kA , (2.31) κ(A) + 1 wobei κ(A) = λmax /λmin die Kondition der Matrix A bzgl. der Norm k · k2 bezeichnet. Die Fehlerabschätzung zeigt, dass das Verfahren umso schneller konvergiert, je näher die Kondition κ(A) bei Eins liegt. Umgekehrt ist die Konvergenz umso langsamer, je größer die Kondition von A ist. Es ist daher erstrebenswert, eine möglichst gut konditionierte Matrix A zu haben, was durch die sogenannte Vorkonditionierung (oder Präkonditionierung) erreicht werden kann. Beispiel 2.6.4 (Diskretisierung einer partiellen Differentialgleichung, Teil III) Wir betrachten wiederum Beispiel 2.5.1 mit N = 100 und f (x, y) = 1000 · sin((x − 0.5) · (y − 0.5)). Die partielle Differentialgleichung wird wie zuvor diskretisiert und führt auf das lineare Gleichungssystem Au = b. Anschließend wird das äquivalente Problem 1 > u Au − b> u → min 2 mit dem CG-Verfahren gelöst. Dabei werden nur die von Null verschiedenen Einträge von A mit einer entsprechenden Indizierung gespeichert, so dass Multiplikationen Ad sehr effizient durchgeführt werden können. Für N = 100 hat A die Dimension 992 = 9801 und von den 98012 = 96059601 Einträgen sind lediglich 48609 Einträge von Null verschieden. Als Startwert wählen wir u(0) = 0. Für N = 100 und ε = 10−8 ist das Abbruchkriterium kg (k) k2 ≤ ε nach 153 Iterationen erfüllt. Die folgende Abbildung zeigt die numerische Lösung u der partiellen Differentialgleichung: 53 Kapitel 3 Lineare Ausgleichsprobleme Häufig gilt es, über- oder unterbestimmte Gleichungssysteme zu lösen. Beispiel 3.0.1 (Lineares Ausgleichsproblem) Die Messung eines physikalischen Vorgangs hat an den Ortspunkten pi = (p1,i , p2,i ), i = 1, . . . , 9, die folgenden Messwerte yi , i = 1, . . . , 9, ergeben: 2 3 4 5 6 7 8 9 Messpunkt i 1 p1,i 6 7 9 11 13 15 17 18 19 p2,i 0 1 2 3 4 5 6 7 8 yi 96 189 283 373 467 553 647 733 832 Es wird ein affin-linearer Zusammenhang zwischen den Messpunkten und den Messwerten vermutet und folgender Ansatz gewählt y(p1 , p2 ) = x1 + x2 p1 + x3 p2 , wobei x = (x1 , x2 , x3 )> ∈ R3 unbekannt ist. Ein naheliegender Versuch, x zu bestimmen, besteht darin, das lineare Gleichungssystem yi = y(p1,i , p2,i ) = x1 + x2 p1,i + x3 p2,i , zu lösen. In Matrixschreibweise lautet dies Ax = b mit 96 1 6 189 1 7 283 1 9 373 1 11 , A = 1 13 b= 467 553 1 15 647 1 17 733 1 18 832 1 19 i = 1, . . . , 9, 0 1 2 3 4 5 6 7 8 . Offenbar ist dieses lineare Gleichungssystem überbestimmt, da die Anzahl der Gleichungen die Anzahl der Unbekannten übersteigt. Leider stellt sich bei näherer Betrachtung heraus, dass das Gleichungssystem keine Lösung besitzt. Der Versuch, den Vektor x über die Lösung des Gleichungssystems Ax = b zu bestimmen, schlägt also fehl. 54 Ein sinnvoller Zugang ist es daher, den Vektor x so zu wählen, dass der Fehler im Gleichungssystem kAx−bk22 möglichst klein wird. Dies führt auf das Minimierungsproblem kAx − bk22 → min . Dieses Problem ist auch als lineares Ausgleichsproblem oder Least-Squares-Problem bekannt, welches bereits auf Gauß zurück geht. In diesem Kapitel beschäftigen wir uns mit dem folgenden Problem. Problem 3.0.2 (Lineares Ausgleichsproblem) Gegeben seien eine Matrix A ∈ Rm×n und ein Vektor b ∈ Rm . Bestimme die Lösung x = (x1 , . . . , xn )> ∈ Rn des Minimierungsproblems Minimiere 1 f (x) := kAx − bk22 2 (3.1) Bemerkung 3.0.3 • Im Fall m = n und A invertierbar, ist die Lösung eindeutig durch das lineare Gleichungssystem Ax = b bestimmt. • Interessant ist insbesondere der Fall m > n, der in praktischen Anwendungen besonders häufig auftritt, zumeist ist m dort wesentlich grösser als n. • Der Fall m < n kann ebenfalls auftreten. In diesem Fall ist das Gleichungssystem Ax = b unterbestimmt, es kann jedoch unlösbar sein, wenn Rang(A) 6= Rang(A|b). • Anstatt der Norm k · k2 können auch andere Normen verwendet werden, etwa k · k∞ oder k · k1 , jedoch führt dies im Allgemeinen auf andere Lösungen und die Berechnung von Lösungen gestaltet sich mitunter schwierig, da die Funktion f nicht mehr differenzierbar sein muss. 3.1 Gauß’sche Normalengleichungen Wir betrachten das lineare Ausgleichsproblem von einem geometrischen Standpunkt. Dazu definieren wir den linearen Vektorraum V := {y ∈ Rm | y = Ax, x ∈ Rn } = Bild(A). V ist ein Teilraum des Rm mit Dimension dim(V ) = Rang(A). Das lineare Ausgleichsproblem lautet somit: 55 Finde ŷ ∈ V , so dass kb − ŷk2 ≤ kb − yk2 ∀y ∈ V. (3.2) Definition 3.1.1 Sei (X, k · k) ein normierter Vektorraum, V ⊂ X ein Teilraum und b ∈ X. ŷ ∈ V heißt Bestapproximierende an b in V , wenn (3.2) gilt. b ŷ V Abbildung 3.1: Bestapproximierende ŷ an b in V . Aus der geometrischen Anschauung ist klar, dass ŷ genau dann Bestapproximierende an b in V ist, wenn b − ŷ senkrecht auf V steht, vgl. Abbildung 3.1. Dies besagt der folgende Satz, der in recht allgemeiner Form gilt. Satz 3.1.2 (Projektionssatz) Sei X ein Vektorraum über K = R oder K = C und h·, ·i ein Skalarprodukt auf X × X. Sei b ∈ X und V ⊂ X ein Teilraum. Dann gelten: (i) Es existiert höchstens eine Bestapproximierende ŷ an b in V . (ii) ŷ ∈ V ist Bestapproximierende an b in V genau dann, wenn hb − ŷ, vi = 0 ∀v ∈ V. (3.3) Beweis: Wir zeigen zunächst (ii). Teil (i) ergibt sich unterwegs. ⇒: Sei ŷ Bestapproximierende an b in V . Angenommen, (3.3) gilt nicht. Dann gibt es v ein v ∈ V mit hb − ŷ, vi = a 6= 0. Insbesondere ist v 6= 0. Setze w = ŷ + a kvk 2 ∈ V . Es folgt av av |a|2 2 2 , b − ŷ − = kb − ŷk − < kb − ŷk2 . kb − wk = b − ŷ − 2 2 2 kvk kvk kvk Dies ist ein Widerspruch dazu, dass ŷ Bestapproximierende ist. 56 ⇐: Es gelte (3.3). Dann gilt für alle v ∈ V mit v 6= ŷ, kb − vk2 = hb − v, b − vi = hb − ŷ + ŷ − v, b − ŷ + ŷ − vi = kb − ŷk2 + kŷ − vk2 + hb − ŷ, ŷ − vi + hŷ − v, b − ŷi = kb − ŷk2 + kŷ − vk2 > kb − ŷk2 . Damit ist ŷ Bestapproximierende. (i) ist hiermit ebenfalls gezeigt, da obige Abschätzung auf einen Widerspruch führt, wenn man annimmt, dass es zwei unterschiedliche Bestapproximierende gibt. 2 Anwendung des Projektionssatzes 3.1.2 auf das lineare Ausgleichsproblem liefert folgenden Satz: Satz 3.1.3 (Gauß’sche Normalengleichungen) x̂ ∈ Rn löst das lineare Ausgleichsproblem genau dann, wenn die Gauß’schen Normalengleichungen A> Ax̂ = A> b gelten. Beweis: Mit ŷ = Ax̂ und v = Ax lautet (3.3): 0 = hb − Ax̂, Axi = (b − Ax̂)> Ax = A> b − A> Ax̂ > x. Da diese Gleichung für alle x ∈ Rn gelten muss, muss zwangsläufig A> Ax̂ − A> b = 0 gelten. 2 Die Gauß’schen Normalengleichungen charakterisieren also die Lösung des linearen Ausgleichproblems. Es gilt zu klären, wann die Normalengleichungen (bzw. das lineare Ausgleichsproblem) eine eindeutige Lösung besitzen. Offenbar ist A> A symmetrisch und positiv semidefinit. Satz 3.1.4 Für m ≥ n besitzen die Gauß’schen Normalengleichungen genau dann eine eindeutige Lösung, wenn Rang(A) = n ist. Beweis: A> A ist positiv definit, genau dann, wenn 0 < x> A> Ax = kAxk22 für alle x 6= 0 gilt. kAxk2 > 0 für alle x 6= 0 ist aber äquivalent dazu, dass Rang(A) = n. 2 Beispiel 3.1.5 (vgl. Beispiel 3.0.1) Wir betrachten wiederum Beispiel 3.0.1 und lösen das lineare Ausgleichsproblem (3.1) mit 57 Hilfe der Normalengleichungen A> Ax = A> b. Es gilt 9 115 36 4173 A> A = 115 1655 565 , A> b = 62916 , 36 565 204 22176 x1 x = x2 . x3 Die Lösung der Normalengleichungen lautet x1 = 533 = 106.6, 5 x2 = − 96 ≈ −1.4769, 65 x3 = 6109 ≈ 93.9846. 65 Ein Nachteil der Gauß’schen Normalengleichungen liegt in der möglicherweise schlechten Kondition der Matrix A> A im Vergleich zur Matrix A, da für quadratische Matrizen κ2 (A> A) = κ2 (A)2 gilt. Aus diesem Grund erfolgt die numerische Lösung des Ausgleichsproblems in der Regel nicht über die Normalengleichungen, sondern wird mit Hilfe einer QR-Zerlegung der Matrix A durchgeführt. Eigenschaften und Verfahren zur Berechnung einer QR-Zerlegung werden im folgenden Abschnitt untersucht. 3.2 QR-Zerlegung einer Matrix Ähnlich zur LR-Zerlegung beim Gauß’schen Eliminationsverfahren ist es unser Ziel, eine nicht notwendig quadratische Matrix A ∈ Rm×n zu faktorisieren in A = Q · R. Definition 3.2.1 (QR-Zerlegung) Sei A ∈ Rm×n mit m ≥ n gegeben. Eine Zerlegung der Form A = Q · R mit Q ∈ Rm×m und R ∈ Rm×n , wobei Q orthogonal ist und R die Gestalt ! R̄ R= 0 mit einer invertierbaren rechten oberen Dreiecksmatrix R̄ ∈ Rn×n hat, heißt QR-Zerlegung von A. Wir werden sehen, dass eine solche QR-Zerlegung stets möglich ist, wenn Rang(A) = n gilt. Wir beschränken die Darstellung auf den Fall K = R. Im komplexen Fall gelten die Beziehungen analog. 3.2.1 QR-Zerlegung zur Lösung von linearen Gleichungssystemen Ist eine QR-Zerlegung von A bekannt, so kann das lineare Gleichungssystem (2.3) im Fall n = m gelöst werden, indem die Orthogonalität von Q ausgenutzt wird: b = Ax = Q · Rx ⇔ 58 c := Q> b = Rx. Das lineare Gleichungssystem Rx = c kann wiederum mit Rückwärtssubstitution gelöst werden. Der Hauptvorteil der QR-Zerlegung im Vergleich zur LR-Zerlegung liegt in der Tatsache begründet, dass orthogonale Transformationen invariant bzgl. der Spektralnorm sind. Für einen beliebigen Vektor x gilt kQxk2 = kxk2 , kQk2 = 1. Speziell gilt kQ> bk2 = kbk2 und somit werden durch Multiplikation mit Q> möglicherweise vorhandene Fehler in b nicht verstärkt. Für den relativen Fehler in (2.14) erhält man für die LR-Zerlegung wegen κ2 (A) = κ2 (L · R) ≤ κ2 (L) · κ2 (R) die im Vergleich zu (2.14) schlechtere“ Abschätzung ” k∆bk2 k∆xk2 ≤ κ2 (L) · κ2 (R) · . kxk2 kbk2 Bei Verwendung der QR-Zerlegung gilt κ2 (A) = κ2 (Q · R) = κ2 (R) und somit folgt k∆bk2 k∆xk2 ≤ κ2 (R) · . kxk2 kbk2 Fazit: Die QR-Zerlegung ist numerisch stabiler als die LR-Zerlegung. Der Nachteil der QR-Zerlegung liegt im etwa doppelt so hohen Aufwand im Vergleich zur LR-Zerlegung. Bemerkung 3.2.2 Die rechte obere Dreicksmatrix R in der LR-Zerlegung von A ist nicht zu verwechseln mit der rechten oberen Dreiecksmatrix R, die bei der QR-Zerlegung berechnet wird. Die jeweiligen Matrizen sind im Allgemeinen verschieden. 3.2.2 Lösung des linearen Ausgleichsproblems mittels QR-Zerlegung Wir wenden uns wieder dem linearen Ausgleichsproblem (3.1) zu und möchten die Funktion f (x) = 21 kAx − bk22 minimieren, wobei A ∈ Rm×n mit m ≥ n und b ∈ Rm gegeben sind. Wir setzen voraus, dass eine QR-Zerlegung von A gegeben ist gemäß ! R̄ A = Q · R, R= , Q ∈ Rm×m , R ∈ Rm×n . 0 Hierin sei R̄ ∈ Rn×n eine invertierbare rechte obere Dreiecksmatrix und Q eine orthogonale Matrix. Wegen der Invarianz der Euklidischen Norm gilt kAx − bk22 = kQ> (Ax − b)k22 = kRx − Q> bk22 . Definiere c= c1 c2 ! := Q> b, c1 ∈ Rn , c2 ∈ Rm−n . 59 Damit folgt kRx − Q> bk22 = kR̄x − c1 k22 + kc2 k22 . Dieser Ausdruck wird genau dann minimal, wenn R̄x = c1 gilt, welches ein lineares Gleichungssystem für x darstellt und mittels Rückwärtssubstitution gelöst werden kann. Die Lösung dieses Gleichungssystems x̂ = R̄−1 c1 ist auch die Lösung der Minimierungsaufgabe und somit des linearen Ausgleichsproblems. Der Wert des Zielfunktionals kann ebenfalls direkt abgelesen werden und lautet 1 f (x̂) = kc2 k22 . 2 3.2.3 QR-Zerlegung nach Householder In diesem Abschnitt wird die QR-Zerlegung nach Householder besprochen, die eine numerisch stabile Methode zur Berechnung einer QR-Zerlegung A = Q·R gemäß Definition 3.2.1 darstellt. Definition 3.2.3 (Householder Spiegelung) Eine Matrix H ∈ Rm×m heißt Householder-Spiegelung oder Householder-Matrix, wenn es ein u ∈ Rm mit kuk2 = 1 gibt, so dass H = I − 2uu> . Eine Householder-Matrix ist offensichtlich symmetrisch und wegen H > H = HH = (I − 2uu> )(I − 2uu> ) = I − 2uu> − 2uu> + 4u |{z} u> u u> = I =1 orthogonal. Geometrisch beschreibt die Abbildung v 7→ Hv, wie im Bild dargestellt, eine Spiegelung von v an einer Ebene durch 0 mit Normalenvektor u. 60 Man bestimmt die Househoder-Spiegelungen so, dass ein beliebiger Vektor v ∈ Rn auf ein Vielfaches des ersten Einheitsvektors e1 abgebildet wird: ! Hv = (I − 2uu> )v = v − 2uu> v = ρe1 . Wegen kvk2 = kHvk2 = kρe1 k2 = |ρ| folgt ρ = ±kvk2 . Setzt man γ := 2u> v so folgt γu = v − ρe1 und damit muss γ = kv − ρe1 k2 gelten. Es gilt also Lemma 3.2.4 Sei v ∈ Rn . Definiert man ρ := ±kvk2 , γ := kv − ρe1 k2 , u := γ1 (v − ρe1 ), so gilt Hv = ρe1 für die Householder-Spiegelung H := I − 2uu> . Numerisch bevorzugt man folgende Wahl von ρ zur Vermeidung von Auslöschung bei v1 : ( −kvk2 , falls v1 ≥ 0, ρ := kvk2 , falls v1 < 0. Zur Anwendung: Seien bereits orthogonale Matrizen Q1 , . . . , Qi ∈ Rm×m bekannt, mit ! Ri Ci Qi Qi−1 · · · Q1 A = , 0 Mi wobei Ri ∈ Ri×i eine obere Dreiecksmatrix ist. Sei v ∈ Rm−i der erste Spaltenvektor von Mi . Wir bestimmen die Householder-Spiegelung Hi+1 ∈ R(m−i)×(m−i) so, dass Hi+1 v ein Vielfaches der ersten Einheitsvektors in Rm−i ist. Dann setzen wir ! I 0 Qi+1 := . 0 Hi+1 Qi+1 ist eine orthogonale Matrix und es gilt Qi+1 Qi Qi−1 · · · Q1 A = =: I 0 0 Hi+1 ! Ri+1 Ci+1 0 Mi+1 Ri Ci 0 Mi ! ! = Ri Ci 0 Hi+1 Mi ! , wobei Ri+1 ∈ R(i+1)×(i+1) wieder eine obere Dreiecksmatrix ist, da in Hi+1 Mi die erste Spalte ein Vielfaches des Einheitsvektors ist. 61 Nach min{m−1, n} Schritten hat man eine obere Dreiecksmatrix erreicht. Da das Produkt orthogonaler Matrizen wieder eine orthogonale Matrix ist, hat man so eine QR-Zerlegung erhalten. Im Pseudocode erhalten wir den folgenden Algorithmus: Algorithmus 3.2.5 (QR-Zerlegung nach Householder) Gegeben: Matrix A ∈ Rm×n , m ≥ n. for i = 1 : min{m − 1, n} v := A(i : m, i); if v1 ≥ 0 ρ(i) := −kvk2 ; else ρ(i) := +kvk2 ; end ui := v − ρ(i)e1 ; ui := ui /kui k2 ; A(i : m, i : n) := A(i : m, i : n) − 2ui (ui )> A(i : m, i : n); end Der Algorithmus ist wohldefiniert, falls Rang(A) = n ist. Die Diagonalelemente von R werden im Vektor ρ abgespeichert, der Rest steht in der Matrix A oberhalb der Diagonalen. Die Householder-Matrizen können in Form der Vektoren ui im unteren Dreieck von A abgespeichert werden. Der Aufwand beträgt mn2 − 31 n3 + O(mn + n2 ). Im Fall von quadratischen Matrizen ist daher der Aufwand mit 32 n3 + O(n2 ) in etwa doppelt so hoch wie bei der LR-Zerlegung. Beispiel 3.2.6 Mit gerundete Werten erhält man die QR-Zerlegung 5 7 −0.700 0.099 0.707 −7.140 −9.660 0 2.380 −5 −7 = A = Q · R = 0.700 −0.099 0.707 −1 1 0.140 0.990 0 0 0 62 Kapitel 4 Nichtlineare Gleichungen In vielen Problemstellungen ist es erforderlich, eine nichtlineare Gleichung F (x) = 0, F : R → R, x ∈ R, (4.1) F : Rn → Rn , x = (x1 , . . . , xn )> ∈ Rn (4.2) oder ein nichtlineares Gleichungssystem F1 (x1 , . . . , xn ) .. = 0, F (x) = . Fn (x1 , . . . , xn ) zu lösen. Offenbar ist (4.1) als Spezialfall mit n = 1 in (4.2) enthalten. Es gibt jedoch spezielle Verfahren für den eindimensionalen Fall n = 1. Nichtlineare Gleichungssysteme können meist nicht analytisch gelöst werden, sondern erfordern numerische Verfahren, wie das Beispiel F (x) = x − cos(x) zeigt. In industriellen Anwendungen treten häufig Probleme mit mehreren tausend Variablen auf. Beispiel 4.0.1 Bei der Simulation eines Autos tritt das Problem auf, den Kontaktpunkt (x, y) eines Autoreifens mit Mittelpunkt (xR , yR ) und Radius r auf einer unebenen Fahrbahn zu bestimmen, vgl. Abbildung. Das Fahrbahnprofil sei durch die nichtlineare Funktion f : R → R beschrieben. 63 (xR,yR) α r f(x) (x,y) y x Der Kontaktpunkt ist implizit gegeben durch die nichtlineare Gleichung f (xR + r sin α) = yR − r cos α, die den Winkel α festlegt. Weitere Anwendungen für nichtlineare Gleichungssysteme treten bei der Diskretisierung von Anfangswertproblemen mittels impliziter Integrationsverfahren, oder der Diskretisierung von nichtlinearen partiellen Differentialgleichungen und Randwertproblemen auf. Ferner treten sie in Optimierungsverfahren der restringierten (Lagrange-Newton-Verfahren) und unrestringierten Optimierung (Newtonverfahren) auf, oder bei der Bestimmung von Eigenwerten mittels Polynomnullstellen. Im Folgenden werden Algorithmen zur approximativen Bestimmung von zumindest einer Nullstelle von F untersucht, wobei vorausgesetzt wird, dass eine Nullstelle existiert. Definition 4.0.2 (Nullstelle) x̂ ∈ Rn heißt Nullstelle von (4.2), wenn F (x̂) = 0 gilt. Die Verfahren in diesem Kapitel haben die Gemeinsamkeit, dass sie iterativ arbeiten, d.h. es wird eine Folge {x(k) } mit x(k+1) = g(x(k) ), k = 0, 1, 2, . . . berechnet, wobei g eine geeignete Funktion ist. In jedem Iterationsschritt muss dabei ein einfach zu lösendes Problem gelöst werden. Das Hauptaugenmerk wird auf folgenden Fragestellungen liegen: • Unter welchen Voraussetzungen konvergierte die Folge {x(k) } gegen x̂? 64 • Wie schnell konvergiert {x(k) } gegen x̂? Die Konvergenzgeschwindigkeit einer Folge ist wie folgt definiert, wobei eine Folge umso schneller konvergiert, je höher die Konvergenzordnung ist. Definition 4.0.3 (lineare Konvergenz, quadratische Konvergenz, superlineare Konvergenz, Konvergenzordnung) Die Folge {x(i) } konvergiere gegen x̂. (a) Die Folge {x(i) } heißt linear konvergent, falls es eine Konstante 0 ≤ C < 1 und eine Zahl i0 gibt mit kx(i+1) − x̂k ≤ C · kx(i) − x̂k ∀i ≥ i0 . (b) Die Folge {x(i) } heißt quadratisch konvergent, falls es eine Konstante C und eine Zahl i0 gibt mit kx(i+1) − x̂k ≤ C · kx(i) − x̂k2 ∀i ≥ i0 . (c) Die Folge {x(i) } heißt konvergent mit Ordnung p > 1, falls es eine Konstante C und eine Zahl i0 gibt mit kx(i+1) − x̂k ≤ C · kx(i) − x̂kp ∀i ≥ i0 . (d) Die Folge {x(i) } heißt superlinear konvergent, falls es eine Folge {Ci } mit limi→∞ Ci = 0 und eine Zahl i0 gibt mit kx(i+1) − x̂k ≤ Ci · kx(i) − x̂k ∀i ≥ i0 . Beachte, dass die Konvergenz der Folge {x(i) } in der Definition explizit vorausgesetzt wurde. Im Fall der linearen und superlinearen Konvergenz ist dies nicht nötig, da lineare und superlineare Konvergenz automatisch die Konvergenz der Folge nach sich ziehen. Für die übrigen Konvergenzbegriffe gilt dies nicht, so dass beim Beweis der Konvergenz mit Ordnung p > 1 stets zunächst die Konvergenz der Folge gezeigt werden muss. 4.1 Bisektion (Intervallschachtelungsverfahren) Im eindimensionalen Fall n = 1 wird unter der Voraussetzung, dass F stetig ist, häufig das Bisketionsverfahren zur Bestimmung einer Nullstelle verwendet. Ausgehend von einem Startintervall [a, b], a < b, mit F (a) · F (b) < 0 (4.3) 65 F b(3) b(2) b(1) b(0) a(2) a(3) a(1) a(0) x Abbildung 4.1: Bisketion: Iterierte Intervallhalbierung. berechnet das Bisketionsverfahren in jedem Schritt ein neues Intervall mit halber Länge, in dem eine Nullstelle von F enthalten ist. Man erhält so eine beliebig genaue Einschachtelung einer Nullstelle. Das Bisektionsverfahren kann u.a. angewendet werden, um Eigenwerte zu berechnen, indem eine Nullstelle des charakteristischen Polynoms gesucht wird. Bedingung (4.3) besagt, dass F (a) und F (b) unterschiedliche Vorzeichen besitzen. In diesem Fall garantiert der Zwischenwertsatz1 , dass im Intervall (a, b) (mindestens) eine Nullstelle von F existiert. Durch iterierte Intervallschachtelung, vgl. Abbildung 4.1, erlaubt das Bisektionsverfahren die beliebig genaue Approximation einer Nullstelle in [a, b] (es ist aber nicht gesagt, welche Nullstelle man bekommt). Algorithmus 4.1.1 (Bisektionsverfahren) Wähle a < b mit F (a) · F (b) < 0, eine Toleranz tol. while |b − a| > tol c := (a + b)/2; if F (a) · F (c) > 0 a := c; else b := c; end 1 Sei f : [a, b] → R stetig in dem abgeschlossenen Intervall [a, b] und sei y eine beliebige Zahl mit f (a) ≤ y ≤ f (b) oder f (b) ≤ y ≤ f (a). Dann existiert ein Punkt c ∈ [a, b] mit f (c) = y. 66 end Nullstelle := (a + b)/2; Aus der Intervallhalbierung in jedem Iterationsschritt erhält man sofort Satz 4.1.2 Sei F : R → R stetig in [a, b], a < b, mit F (a) · F (b) < 0. Nach k Iterationen des Bisektionsverfahrens ist die berechnete Nullstellenapproximation als Intervallmittelpunkt von jeder Nullstelle in [a(k) , b(k) ] höchstens 2b−a k+1 entfernt. Beispiel 4.1.3 √ Berechnung von 2 als (positive) Nullstelle von F (x) = x2 − 2 mit dem Bisektionsverfahren mit Startintervall [1, 1.5] liefert x̂ ≈ 1.4142. Iteration 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 a 1.0000000000 1.2500000000 1.3750000000 1.3750000000 1.4062500000 1.4062500000 1.4140625000 1.4140625000 1.4140625000 1.4140625000 1.4140625000 1.4140625000 1.4141845703 1.4141845703 1.4141845703 1.4141998291 1.4142074585 b maximaler Fehler 1.5000000000 0.2500000000 1.5000000000 0.1250000000 1.5000000000 0.0625000000 1.4375000000 0.0312500000 1.4375000000 0.0156250000 1.4218750000 0.0078125000 1.4218750000 0.0039062500 1.4179687500 0.0019531250 1.4160156250 0.0009765625 1.4150390625 0.0004882812 1.4145507812 0.0002441406 1.4143066406 0.0001220703 1.4143066406 0.0000610352 1.4142456055 0.0000305176 1.4142150879 0.0000152588 1.4142150879 0.0000076294 1.4142150879 0.0000038147 Das Bisketionsverfahren bietet folgende Vorteile: • Es werden nur Funktionsauswertungen benötigt und das Verfahren ist sehr einfach zu implementieren. • Die Funktion muss lediglich stetig sein. Differenzierbarkeit wird nicht benötigt. • Es ist sehr genau und konvergiert stets, wenn die Anfangswerte richtig gewählt wurden. 67 Leider hat das Bisektionsverfahren auch einige Nachteile: • Das Verfahren lässt sich nur für n = 1 anwenden. • Das Verfahren liefert nur eine langsame Konvergenz. • Es muss ein Startintervall [a, b] mit F (a) · F (b) < 0 bestimmt werden, was z.B. für F (x) = x2 nicht möglich ist. 4.2 Fixpunktiteration Wir haben im Zusammenhang mit iterativen Verfahren zur Lösung von linearen Gleichungssystemen bereits den Banachschen Fixpunktsatz 2.5.4 zur iterativen Berechnung eines Fixpunkts der Gleichung x = g(x) kennengelernt. Die Aufgabenstellung (4.2) lässt sich nun leicht auf Fixpunktgestalt“ ” transformieren, etwa in der Form F (x) = 0 ⇔ x = x + ρ(x) · F (x) =: g(x), wobei ρ(x) ∈ Rn×n eine reguläre, problemabhängige Matrix ist, die dazu dienen soll, die Voraussetzungen des Banachschen Fixpunktsatzes zu erfüllen. Der naive Ansatz ρ(x) = ±I für alle x führt nur selten zum Erfolg. Anstelle des Nullstellenproblems (4.2) wird also ein Fixpunkt x̂ der Funktion g gesucht. Unter der Annahme, dass ρ(x̂) regulär ist, ist der Fixpunkt x̂ wegen x̂ = x̂ + ρ(x̂) · F (x̂) ⇔ ρ(x̂) · F (x̂) = 0 ⇔ F (x̂) = 0 auch Nullstelle von F . Für einen geeigneten Startwert x(0) wird die Fixpunktiteration x(i+1) = g(x(i) ) = x(i) + ρ(x(i) ) · F (x(i) ), i = 0, 1, . . . zur Approximation von x̂ durchgeführt, vgl. Abbildung 4.2. Diese Methode ist auch dann anwendbar, wenn F nicht differenzierbar ist. Um das Konvergenzresultat im Banachschen Fixpunktsatz anwenden zu können, müssen folgende Voraussetzungen erfüllt sein, vgl. Satz 2.5.4: • Selbstabbildungseigenschaft, d.h. es muss eine abgeschlossene Menge D ⊆ Rn existieren, so dass g(D) ⊆ D gilt. • g muss die Kontraktionsbedingung kg(x) − g(y)k ≤ qkx − yk für ein q < 1 und alle x, y ∈ D erfüllen. 68 g g(x0) g(x2) g(x1) x0 x2 x3x1 x Abbildung 4.2: Schematische Darstellung der Fixpunktiteration. Kennt man die Kontraktionszahl q < 1, so folgt aus dem Banachschen Fixpunktsatz die Fehlerabschätzung qi kx(i) − x̂k ≤ kx(1) − x(0) k. 1−q Der Nachteil der Fixpunktiteration liegt in der im Allgemeinen nur linearen Konvergenzgeschwindigkeit. Ist g stetig differenzierbar, so kann die Kontraktionsbedingung (2.18) mit Hilfe der Ableitung von g überprüft werden. Satz 4.2.1 Es sei D ⊆ Rn abgeschlossen und konvex und g : D → D stetig differenzierbar. Gibt es eine Konstante q < 1 mit kg 0 (x)k ≤ q für alle x ∈ D, so erfüllt g in D die Kontraktionsbedingung (2.18). Beweis: Seien x, y ∈ D gegeben. Anwendung des Mittelwertsatzes in Integralform auf die stetig differenzierbare Funktion g liefert Z 1 g(x) − g(y) = g 0 (x + t(y − x))(y − x)dt. 0 Wir schätzen mit einer geeigneten Vektornorm und einer verträglichen Matrixnorm ab: Z 1 kg(x) − g(y)k ≤ kg 0 (x + t(y − x))(y − x)kdt Z0 1 ≤ kg 0 (x + t(y − x))k · ky − xkdt 0 Z 1 0 ≤ max kg (z)k · ky − xkdt z∈D 0 ≤ qkx − yk. 69 2 Bemerkung 4.2.2 Es ist möglich, die Kontraktionszahl q numerisch zu schätzen, sobald einige Iterierte zur Verfügung stehen. Man betrachte die Quotienten q≈ kx(i+1) − x(i) k , kx(i) − x(i−1) k i = 1, 2, . . . . Hintergrund: Aus der Kontraktivität von g folgt kx(i+1) − x(i) k = kg(x(i) ) − g(x(i−1) )k ≤ qkx(i) − x(i−1) k. Das folgende Beispiel zeigt einerseits, dass die Fixpunktiteration nicht immer konvergiert, andererseits wird der Einfluss von ρ verdeutlicht. Beispiel 4.2.3 √ Berechnung von 2 als (positive) Nullstelle von F (x) = x2 − 2. Ausgehend vom Startwert x(0) = 2 wird zunächst mit ρ(x) = 1 für alle x gerechnet. i x(i) Fehler 1 4.0000000000 2.5857864376 2 18.0000000000 16.5857864376 3 340.0000000000 338.5857864376 4 115938.0000000000 115936.5857864376 Die Fixpunktiteration divergiert. Hier ist g(x) = x + x2 − 2 und g 0 (x) = 1 + 2x. Damit ist |g 0 (x)| ≥ 1 für x 6∈ (−1, 0). In einer Umgebung des Startwerts x(0) = 2 ist die Kontraktionsbedingung nicht erfüllt. Wird bei gleichem Startwert hingegen ρ(x) = −1/4 gewählt, so erhält man: 70 i 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 x(i) 1.5000000000 1.4375000000 1.4208984375 1.4161603451 1.4147828143 1.4143802114 1.4142623658 1.4142278560 1.4142177488 1.4142147886 1.4142139215 1.4142136676 1.4142135932 1.4142135714 1.4142135650 Fehler Konvergenzrate 0.0857864376 0.1464466094 0.0232864376 0.2714466094 0.0066848751 0.2870716094 0.0019467827 0.2912220000 0.0005692520 0.2924065231 0.0001666490 0.2927509058 0.0000488034 0.2928515566 0.0000142936 0.2928810180 0.0000041864 0.2928896454 0.0000012262 0.2928921722 0.0000003591 0.2928929122 0.0000001052 0.2928931292 0.0000000308 0.2928931913 0.0000000090 0.2928932135 0.0000000026 0.2928932004 Die Fixpunktiteration konvergiert offensichtlich. Es ist g(x) = x − (x2 − 2)/4 und g 0 (x) = 1 − x/2. Damit ist |g 0 (x)| < 1 für x ∈ (0, 4) und g ist kontrahierend in einer Umgebung des Startwerts x(0) = 2. Zu beachten ist noch, dass die Iterationsfolge im Intervall (0, 4) verbleibt. Somit sind die Voraussetzungen des Banachschen Fixpunktsatzes erfüllt und es folgt die Konvergenz. Die Konvergenzkonstante beträgt hier ca. 0.29289. Die Wahl von ρ ist entscheidend und leider problemabhängig. Bemerkung 4.2.4 (Wahl von ρ) Wir werden im nächsten Abschnitt sehen, dass die Wahl ρ(x) = −F 0 (x)−1 sehr geeignet ist und auf das sogenannte Newtonverfahren x(i+1) = x(i) − F 0 (x(i) )−1 F (x(i) ), i = 0, 1, 2, . . . führt. Insbesondere kann das Newtonverfahren also als Fixpunktiteration mit g(x) = x − F 0 (x)−1 F (x) interpretiert werden. 4.3 Newton-Verfahren Wir suchen eine Nullstelle des nichtlinearen Gleichungssystems F (x) = 0 für eine stetig differenzierbare Funktion F : Rn → Rn . Die Idee beim Newton-Verfahren besteht darin, 71 die Funktion F um die Iterierte x(i) zu linearisieren, d.h., wir betrachten das Taylorpolynom erster Ordnung Ti (x) = F (x(i) ) + F 0 (x(i) )(x − x(i) ), und bestimmen die neue Iterierte x(i+1) als Nullstelle der linearen Funktion Ti . Dabei ist ∂F1 ∂F1 (x) · · · (x) ∂xn ∂x1. .. .. .. F 0 (x) = . . ∂Fn ∂Fn (x) · · · ∂xn (x) ∂x1 die Jacobimatrix von F and der Stelle x = (x1 , . . . , xn )> . Für den eindimensionalen Fall ist das Vorgehen in Abbildung 4.3 veranschaulicht. F x(i+1) x x(i) Abbildung 4.3: Newtonverfahren: Lokale Approximation von F durch Tangente im Punkt x(i) . Die Nullstelle der Tangente definiert die neue Iterierte x(i+1) . Es führt auf die Darstellung x(i+1) = x(i) − F 0 (x(i) )−1 · F (x(i) ). (4.4) Man kann das Newtonverfahren daher, wie in Bemerkung 4.2.4 bereits angedeutet als Fixpunktiteration interpretieren. In der Praxis wird die Inverse der Jacobimatrix F 0 (x)−1 niemals explizit berechnet, da hierfür n lineare Gleichungssysteme der Dimension n × n gelöst werden müssten. Stattdessen wird (4.4) in der Form F 0 (x(i) ) · d(i) = −F (x(i) ), x(i+1) = x(i) + d(i) , i = 0, 1, 2, . . . gelöst. Diese Variante benötigt nicht mehr die Inverse der Jacobimatrix und erfordert nur die Lösung eines linearen Gleichungssystems für d(i) . Zusammenfassend erhalten wir Algorithmus 4.3.1 (Newton-Verfahren) 72 (i) Wähle x(0) ∈ Rn , tol > 0 und setze i := 0. (ii) Falls kF (x(i) )k ≤ tol, STOP. (iii) Löse (z.B. mit der LR-Zerlegung) das lineare Gleichungssystem F 0 (x(i) ) · d(i) = −F (x(i) ). (iv) Setze x(i+1) := x(i) + d(i) , i := i + 1 und gehe zu (ii). Wir wollen die Konvergenz des Verfahren untersuchen und benötigen dazu drei Lemmata. Ziel ist es, die lokal superlineare Konvergenz bzw. sogar die lokal quadratische Konvergenz zu zeigen. Zunächst zeigen wir mit zwei Lemmata, dass die Jacobimatrix F 0 (x) für x hinreichend nahe an x̂ invertierbar bleibt, sofern F 0 (x̂) invertierbar ist. Lemma 4.3.2 (Banach-Lemma) Seien A, B ∈ Rn×n mit kI − B · Ak < 1 für eine induzierte Matrixnorm k · k. Dann sind A und B invertierbar, und es gilt die Abschätzung kA−1 k ≤ kBk . 1 − kI − B · Ak Beweis: Wir zeigen zunächst eine allgemeine Abschätzung. Sei M ∈ Rn×n mit kM k < 1 gegeben. Dann folgt für alle x ∈ Rn k(I − M )xk = kx − M xk ≥ kxk − kM xk ≥ kxk − kM k · kxk = (1 − kM k) · kxk. Daher ist (I − M )x = 0 nur für x = 0 möglich und somit ist I − M invertierbar. Setzen wir x := (I − M )−1 y für ein beliebiges y ∈ Rn ein, so erhalten wir kyk ≥ (1 − kM k) · k(I − M )−1 yk und daher k(I − M )−1 yk 1 ≤ . k(I − M ) k := max y6=0 kyk 1 − kM k −1 Wähle wir nun speziell M := I − B · A, so ist nach Voraussetzung kM k = kI − B · Ak < 1 und, wie eben gezeigt, ist I − M = B · A invertierbar und somit auch A und B. Es folgt kA−1 k = k(I − M )−1 Bk ≤ k(I − M )−1 k · kBk kBk ≤ 1 − kI − B · Ak 73 2 Lemma 4.3.3 Seien F : Rn → Rn stetig differenzierbar, x̂ ∈ Rn und F 0 (x̂) invertierbar. Dann existiert ein r > 0, so dass auch F 0 (x) für alle x ∈ Ur (x̂) = {x ∈ Rn | kx − x̂k < r} invertierbar ist. Weiter existiert eine Konstante C > 0 mit kF 0 (x)−1 k ≤ C ∀x ∈ Ur (x̂). Beweis: Da F 0 in x̂ stetig ist, existiert ein r > 0 mit kF 0 (x̂) − F 0 (x)k ≤ 1 2kF 0 (x̂)−1 k ∀x ∈ Ur (x̂). Also ist kI − F 0 (x̂)−1 F 0 (x)k ≤ kF 0 (x̂)−1 k · kF 0 (x̂) − F 0 (x)k ≤ 1 2 ∀x ∈ Ur (x̂). Das Banach-Lemma 4.3.2 mit A := F 0 (x) und B := F 0 (x̂)−1 liefert kF 0 (x)−1 k ≤ kF 0 (x̂)−1 k ≤ 2kF 0 (x̂)−1 k =: C. 0 −1 0 1 − kI − F (x̂) F (x)k | {z } ≤1/2 2 Nun wollen wir noch ein Lemma, welches wir für die Konvergenzgeschwindigkeit brauchen, beweisen. Lemma 4.3.4 Sei {x(i) } ⊂ Rn eine gegen x̂ ∈ Rn konvergente Folge. Dann gilt: (a) Ist F : Rn → Rn stetig differenzierbar, so ist kF (x(i) ) − F (x̂) − F 0 (x(i) )(x(i) − x̂)k = o(kx(i) − x̂k). (b) Ist F : Rn → Rn stetig differenzierbar und ferner F 0 lokal Lipschitz-stetig in x̂, d.h. ∃r > 0, L > 0 : kF 0 (x) − F 0 (y)k ≤ Lkx − yk ∀x, y ∈ Ur (x̂) so gilt kF (x(i) ) − F (x̂) − F 0 (x(i) )(x(i) − x̂)k = O(kx(i) − x̂k2 ). Beweis: 74 (a) Aus der stetigen Differenzierbarkeit von F folgt kF (x(i) ) − F (x̂) − F 0 (x(i) )(x(i) − x̂)k ≤ kF (x(i) ) − F (x̂) − F 0 (x̂)(x(i) − x̂)k + kF 0 (x(i) ) − F 0 (x̂)k k(x(i) − x̂)k | {z } | {z } →0, da F’ stetig o(kx(i) −x̂k), da F differenzierbar in x̂ = o(kx(i) − x̂k). (b) Mit dem Mittelwertsatz in Integralform und der lokalen Lipschitz-Stetigkeit von F 0 erhalten wir = = ≤ = kF (x(i) ) − F (x̂) − F 0 (x(i) )(x(i) − x̂)k Z 1 0 (i) (i) 0 (i) (i) F (x̂ + t(x − x̂))(x − x̂)dt − F (x )(x − x̂) Z0 1 0 (i) 0 (i) (i) F (x̂ + t(x − x̂)) − F (x ) · (x − x̂)dt 0 Z 1 L(1 − t)kx(i) − x̂k2 dt 0 Z 1 (i) 2 (1 − t)dt L · kx − x̂k · 0 L · kx(i) − x̂k2 = O(kx(i) − x̂k2 ). = 2 2 Nun können wir den Konvergenzsatz für das Newton-Verfahren beweisen. Satz 4.3.5 (Konvergenzsatz für das Newton-Verfahren) Seien F : Rn → Rn stetig differenzierbar, x̂ eine Nullstelle von F und F 0 (x̂) sei invertierbar. Dann gibt ein r > 0, so dass für alle Startwerte x(0) ∈ Ur (x̂) gilt: (a) Das Newton-Verfahren ist wohldefiniert und erzeugt eine Folge {x(i) }, die gegen x̂ konvergiert. (b) Die Konvergenzrate ist superlinear. (c) Ist F 0 zusätzlich noch lokal Lipschitz-stetig um x̂, so ist die Folge {x(i) } sogar quadratisch konvergent. Beweis: Nach Lemma 4.3.3 existiert ein r1 > 0 und ein C > 0 mit kF 0 (x)−1 k ≤ C ∀x ∈ Ur1 (x̂). 75 Nach Lemma 4.3.4(a) existiert ein r2 > 0 mit kF (x(i) ) − F (x̂) − F 0 (x(i) )(x(i) − x̂)k ≤ 1 kx(i) − x̂k 2C ∀x ∈ Ur2 (x̂). Sei r := min{r1 , r2 }. Wähle x(0) ∈ Ur (x̂). Dann ist x(1) wohldefiniert und es gilt kx(1) − x̂k = kx(0) − x̂ − F 0 (x(0) )−1 F (x(0) )k ≤ kF 0 (x(0) )−1 k · kF 0 (x(0) )(x(0) − x̂) − F (x(0) ) + F (x̂) k | {z } =0 1 kx(0) − x̂k ≤ C· 2C 1 (0) = kx − x̂k. 2 Induktiv folgt nun die Wohldefiniertheit der Folge {x(i) } mit 1 1 kx(i) − x̂k ≤ kx(i−1) − x̂k ≤ . . . ≤ i kx(0) − x̂k. 2 2 Dies beweist, dass die Folge {x(i) } gegen x̂ konvergiert und somit Teil (a). Analog erhält man die Abschätzung kx(i+1) − x̂k ≤ C · kF 0 (x(i) )(x(i) − x̂) − F (x(i) ) + F (x̂)k, woraus mit Lemma 4.3.4 die Aussagen (b) und (c) unmittelbar folgen. 2 Bemerkung 4.3.6 Beachte, dass das Newtonverfahren nur für solche Anfangswerte x(0) mit kx̂ − x(0) k ≤ r konvergiert. Für Anfangswerte außerhalb dieser Umgebung kann das Verfahren divergieren. Man spricht hier von lokaler Konvergenz und vom lokalen Newton-Verfahren. Beispiel 4.3.7 √ Wir betrachten erneut das Problem, 2 als (positive) Nullstelle von F (x) = x2 − 2 zu berechnen. Das Newtonverfahren liefert mit Startwert x(0) = 2 die folgende Ausgabe: i x(i) Fehler 1 1.5000000000 0.0857864376 2 1.4166666667 0.0024531043 3 1.4142156863 0.0000021239 4 1.4142135624 0.0000000000 Hier wird die Stärke des Newtonverfahrens im Vergleich zum Bisektionsverfahren deutlich. Während das Bisektionsverfahren 17 Iterationen benötigte, konvergiert das Newtonverfahren in nur 4 Iterationen, wobei die Anzahl der führenden Nullen im Fehler sich von 76 Iteration zu Iteration in etwa verdoppeln. Dies ist die quadratische Konvergenz des Newtonverfahrens. Beispiel 4.3.8 (Funktion von Himmelblau) Wir betrachten ein Beispiel aus der unrestringierten Optimierung. Gesucht ist ein lokales Minimum der Funktion f (x, y) = (x2 + y − 11)2 + (x + y 2 − 7)2 . Aus der Theorie der nichtlinearen Optimierung ist bekannt, dass ein lokales Minimum (x̂, ŷ) von f notwendig die Bedingung ∇f (x̂, ŷ) = 0 erfüllt, d.h. wir suchen eine Nullstelle der nichtlinearen Funktion F (x, y) := ∇f (x, y) = 4x(x2 + y − 11) + 2(x + y 2 − 7) 2(x2 + y − 11) + 4y(x + y 2 − 7) ! , welche die Jacobimatrix F 0 (x, y) = 12x2 + 4y − 42 4(x + y) 4(x + y) 4x + 12y 2 − 26 ! besitzt. Ein Blick auf den Graphen von f zeigt, dass es in diesem Beispiel mehrere Nullstellen gibt, nämlich • 4 lokale Minimalstellen (zugleich global) mit Funktionswert 0; darunter der Punkt (3, 2)> . • 4 Sattelpunkte • ein lokales Maximum in (−0.270845, −0.923039)> 77 800 600 4 400 2 200 y0 0 –4 –2 –2 x0 –4 2 –4 –2 0 x 2 4 4 –4 –2 0 y 2 4 Als Abbruchkriterium verwenden wir kF (x, y)k ≤ 10−13 und als Startpunkt (x(0) , y (0) ) = (4, 2.5). Das Newtonverfahren liefert quadratische Konvergenz gegen das Minimum (3, 2)> . (i) (i) i x1 x2 kF (x(i) )k 0 4.000000 2.500000 1.35e+02 1 3.281417 2.056664 2.62e+01 2 3.035131 1.988137 2.42e+00 3 3.000634 1.999744 4.20e-00 4 3.000000 2.000000 1.41e-05 5 3.000000 2.000000 1.50e-12 6 3.000000 2.000000 0.00e+00 Für andere Startwerte konvergiert das Newtonverfahren mitunter gegen andere Nullstellen von F . In der Praxis stellt sich die Frage, wann das Newtonverfahren angehalten werden kann. Das Kriterium kF (x(i) )k ≤ tol ist zwar naheliegend, aber nicht skalierungsinvariant, d.h. durch Multiplikation von F mit einer Konstanten, lässt sich die Bedingung stets erfüllen und gibt mitunter keine Auskunft über die Güte der Iterierten x(i) . Als geeigneteres Konvergenzkriterium hat sich beispielsweise der natürliche Monotonietest bewährt, vgl. Deuflhard und Hohmann [DH91]. In jedem Iterationsschritt wird dabei geprüft, ob kd¯(i+1) k ≤ Θkd(i) k gilt, wobei meist Θ = 0.5 gewählt wird. Darin bezeichnet d(i) die Newtonrichtung und 78 d¯(i+1) ist die Lösung des linearen Gleichungssystems F 0 (x(i) )d¯(i+1) = −F (x(i+1) ). Die Lösung dieses Gleichungssystems erfordert nur eine Vorwärts-Rückwärtssubstitution, da die LR-Zerlegung von F 0 (x(i) ) bereits bei der Berechnung von x(i+1) berechnet wurde. Bemerkung 4.3.9 (Varianten des Newtonverfahrens) In der Praxis werden häufig Varianten des Newtonverfahrens verwendet, die hauptsächlich versuchen, die Vorteile des Newtonverfahrens (superlineare bzw. quadratische Konvergenzgeschwindigkeit, Anwendbarkeit für n > 1) zu erhalten und die Nachteile des Newtonverfahrens zu umgehen. Die Hauptnachteile des Verfahrens bestehen in der im Allgemeinen nur lokalen Konvergenz des Verfahrens, d.h. sobald der Startwert nicht hinreichend nahe an der Nullstelle liegt, besteht die Gefahr der Divergenz des Verfahrens. Desweiteren wird die Jacobimatrix von F benötigt. Die Berechnung derselben kann sehr aufwändig sein und ist oft nur numerisch möglich, mitunter ist sie auch nicht einmal invertierbar. • Vereinfachtes Newtonverfahren: Um den Rechenaufwand zu reduzieren wird die Jacobimatrix F 0 nicht in jedem Iterationsschritt neu berechnet, sondern es wird vereinfachend die konstante Matrix A0 := F 0 (x(0) ) für alle Iterationen (bzw. für eine bestimmte Anzahl von Iterationen) verwendet. Dadurch entsteht ein linear konvergentes Verfahren. • Quasi-Newton-Verfahren: Quasi-Newton-Verfahren ersetzen die Jacobimatrix im Newtonverfahren durch eine geeignete Approximation Bi ≈ F 0 (x(i) ). Die Approximation Bi+1 wird dann durch eine geeignete Update-Formel aus Bi berechnet. Man erhält so ein superlinear konvergentes Verfahren. • Globalisierung des Newtonverfahrens: Ein Hauptnachteil des Newtonverfahrens ist die nur lokale Konvergenz. Es gibt jedoch Ansätze, die eine Konvergenz des Newtonverfahrens für beliebige Startwerte garantieren (natürlich nur unter bestimmten Voraussetzungen). Die Globalisierung des Newtonverfahrens besteht in der Einführung einer Schrittweite αi > 0 bei der Berechnung von x(i+1) gemäß x(i+1) = x(i) + αi d(i) . Die Schrittweite αi wird hierbei so bestimmt, dass das Residuum kF k beim Übergang von x(i) zu x(i+1) hinreichend stark abnimmt. Diese Technik ist aus der Theorie der nichtlinearen Optimierung bekannt und kann z.B. durch eindimensionale Liniensuche realisert werden. Details finden sich z.B. in Geiger und Kanzow [GK99] und Kosmol [Kos93]. Globalisierte Newtonverfahren sind in der Regel so konstruiert, 79 dass die Schrittweite αi = 1 akzeptiert wird, sobald x(i) den lokalen Konvergenzeinzugsbereich des lokalen Newtonverfahrens erreicht. Insofern erben die globalisierten Varianten die lokalen Konvergenzeigenschaften des Newtonverfahrens. 4.4 Sekantenverfahren Wiederum im eindimensionalen Fall n = 1 verwendet man häufig auch das Sekantenverfahren. Sind zwei Näherungen x(i−1) 6= x(i) gegeben, approximiert man die Funktion F lokal durch die Sekante S(x), die durch die Punkte (x(i−1) , F (x(i−1) )) und (x(i) , F (x(i) )) verläuft, vgl. Abbildung 4.4. Die Sekante ist durch S(x) = F (x(i−1) ) + x − x(i−1) (F (x(i) ) − F (x(i−1) )) x(i) − x(i−1) gegeben. Die nächste Näherung x(i+1) ist durch die Nullstelle der Sekante gegeben: S(x(i+1) ) = 0. Diese lässt sich leicht ausrechnen: x(i) − x(i−1) (i+1) (i−1) x = x − F (x(i−1) ) F (x(i) ) − F (x(i−1) ) x(i−1) F (x(i) ) − x(i) F (x(i−1) ) , i = 1, 2, 3, . . . = F (x(i) ) − F (x(i−1) ) F x(i+1) x(i) x(i−1) x Abbildung 4.4: Sekantenverfahren: Ersetzung der Tangente durch Sekante. Zum Start des Verfahrens werden zwei Anfangswerte x(0) und x(1) mit F (x(0) ) 6= F (x(1) ) benötigt. Ein Problem entsteht, wenn während der Iteration der Fall F (x(i) ) ≈ F (x(i−1) ) auftritt. Deshalb sollte das Verfahren abgebrochen werden, wenn |F (x(i) ) − F (x(i−1) )| ≤ ε|F (x(i) )| mit einer gegebenen Toleranz ε gilt. 80 Satz 4.4.1 Sei x̂ ∈ (a, b) eine Nullstelle einer auf [a, b] zweimal stetig differenzierbar Funktion F : R → R mit F 0 (x̂) 6= 0. Dann konvergiert das Sekantenverfahren lokal gegen x̂. Die √ Konvergenzordnung ist p = (1 + 5)/2 ≈ 1.618. 2 Beweis: siehe Lempio [Lem98, S.140] Gemäß Definition 4.0.3 ist das Sekantenverfahren also lokal superlinear konvergent. Vorteile: • Es werden nur Funktionsauswertungen benötigt (sehr einfach zu implementieren). √ • Die Konvergenzordnung ist p = (1 + 5)/2 ≈ 1.618 > 1. Die Konvergenzordnung ist formal zwar kleiner als die des Newtonverfahrens, berücksichtigt man jedoch den Aufwand zur Durchführung eines Iterationsschrittes, zeigt sich, dass das Sekantenverfahren nur eine Funktionsauswertung benötigt, während das Newtonverfahren eine Funktionsauswertung und eine Ableitung benötigt, die im Allgemeinen teurer als eine Funktionsauswertung ist. Um den gleichen Aufwand zu erhalten, können zwei Schritte des Sekantenverfahrens bei nur einem Newtonschritt durchgeführt werden. Damit erhielte man die Konvergenzordnung p2 ≈ 2.618. Nachteile: • Nur für n = 1 anwendbar. • Das Verfahren konvergiert nur lokal aber nicht global, z.B. nicht für F (x) = arctan(x), • Es ist wegen Auslöschung instabil, falls F (x(i) ) ≈ F (x(i−1) ) gilt. Beispiel 4.4.2 √ Wir betrachten erneut das Problem, 2 als (positive) Nullstelle von F (x) = x2 − 2 zu berechnen. Mit dem Sekantenverfahren und den Startpunkten x(0) = 1, x(1) = 2 erhalten wir Iteration x(i+1) Fehler 1 1.3333333333 0.0808802290 2 1.4000000000 0.0142135624 3 1.4146341463 0.0004205840 4 1.4142114385 0.0000021239 5 1.4142135621 0.0000000003 6 1.4142135624 0.0000000000 81 Kapitel 5 Interpolation Interpolationsaufgaben beschäftigen sich mit der Konstruktion einer Funktion, die an endlich vielen Stellen mit einer im Allgemeinen unbekannten oder schwer zu berechnenden Funktion übereinstimmt. Eine typische Aufgabenstellung ist die Messung eines physikalischen oder technischen Vorgangs, welche n + 1 ∈ N verschiedene reellwertige Stützwerte x : x0 x 1 · · · f (x) : f0 f1 · · · xn fn (5.1) ergebe. Die Funktion f , die dem Vorgang zu Grunde liegt, ist häufig nicht explizit bekannt und liegt somit nur an den Stützstellen x0 , . . . , xn mit Funktionswerten fi = f (xi ), i = 0, . . . , n, vor. Um trotzdem mit Funktionswerten zwischen den Stützstellen arbeiten zu können und nicht jedesmal ein mitunter aufwändiges Experiment durchführren zu müssen, wird eine Approximation an die unbekannte Funktion f konstruiert, die zumindest die folgende Interpolationsaufgabe löst und ggf. weitere sinnvolle Eigenschaften besitzt (z.B. hinreichend hohe Glattheit). Problem 5.0.1 (Interpolationsaufgabe) Bestimme eine geeignete Funktion g : R → R, die die Interpolationsbedingungen g(xi ) = fi , i = 0, . . . , n, (5.2) erfüllt, vgl. auch Abbildung 5.1. f0 f2 f3 f1 g(x) x0 x1 x2 x3 x Abbildung 5.1: Interpolation: Stützwerte und interpolierende Funktion g. 82 Eine Funktion g, die die Interpolationsbedingungen (5.2) erfüllt, heißt interpolierende Funktion. Man sagt auch, g interpoliert die Stützwerte (5.1). Interpolationsaufgaben treten auch in der Computergrafik (zeichnen einer geeigneten Kurve durch vorgegebene Punkte), in technischen Geräten (z.B. Scanner oder Digitalkameras: echte Auflösung 600 × 1200 dpi, interpolierte Auflösung 9600 × 9600 dpi), in der Bildverarbeitung (Kompressionsalgorithmen) und in der Signalverarbeitung (Approximation von Signalen) auf. Interpolierende Funktionen können auch zur Approximation von Funktionen dienen, deren Auswertung sehr kostenintensiv ist. In diesem Fall stellt sich die Frage, wie gut die Approximation durch die interpolierende Funktion ist, d.h. wie groß der Approximationsfehler ist. Beispiel 5.0.2 Die Funktion f (x) = exp(x) soll mittels Interpolation von Funktionswerten approximiert werden. Insbesondere soll f an der Stelle x = 0.4 approximiert werden. (a) Wir verwenden die Funktionswerte 1 = f (0) = exp(0) und e = f (1) = exp(1) und berechnen aus den Interpolationsbedingungen 1 = p(0) = a0 , exp(1) = p(1) = a0 + a1 , das interpolierende Polynom ersten Grades p(x) = a0 + a1 x = 1 + (exp(1) − 1)x. Damit ergibt sich an der Stelle x = 0.4 der Approximationsfehler |p(0.4) − f (0.4)| ≈ |1.6873127 − 1.4918247| = 0.1954880. Dieser Fehler ist verhältnismäßig groß. (b) Wir verwenden nun die Funktionswerte 1 = f (0) = exp(0), 1.6487213 ≈ f (0.5) = exp(0.5) und e = f (1) = exp(1) und berechnen aus den Interpolationsbedingungen 1 = p(0) = a0 , a1 a2 + , 2 4 exp(1) = p(1) = a0 + a1 + a2 , exp(0.5) = p(0.5) = a0 + das interpolierende Polynom p(x) = a0 + a1 x + a2 x2 zweiten Grades mit a0 = 1, a1 = 4 exp(0.5) − exp(1) − 3, a2 = 2 exp(1) + 2 − 4 exp(0.5). Damit ergibt sich an der Stelle x = 0.4 der Approximationsfehler |p(0.4) − f (0.4)| ≈ |1.4853099 − 1.4918247| = 0.0065148. Dieser Fehler ist deutlich kleiner als in (a), vgl. nachstehende Abbildung. 83 4.5 4 f(x)=exp(x) interpolierendes Polynom 1. Grades interpolierendes Polynom 2. Grades 3.5 3 2.5 2 1.5 1 0.5 0 -0.5 5.1 0 0.5 1 1.5 Polynominterpolation Bei der Polynominterpolation werden Polynome als interpolierende Funktion g gewählt. Ein Polynom p vom Höchstgrad r kann allgemein als 2 r p(x) = a0 + a1 x + a2 x + . . . + ar x = r X ai x i (5.3) i=0 mit den r + 1 Koeffizienten a0 , a1 , . . . , ar dargestellt werden. r bezeichnet des Grad des Polynoms. Unser Ziel ist es, ein interpolierendes Polynom p mit p(xi ) = fi , i = 0, . . . , n, zu finden. Folgende Fragestellungen sind dabei relevant: • Existenz und Eindeutigkeit eines interpolierenden Polynoms. • Wie kann das (bzw. ein) Interpolationspolynom berechnet werden? • Wie groß ist der Interpolationsfehler zwischen den Stützstellen? Wir weisen die Existenz eines interpolierenden Polynoms konstruktiv nach, indem wir explizit ein interpolierendes Polynom konstruieren. Dabei setzen wir voraus, dass die Knoten xi , i = 0, . . . , n, in (5.1) paarweise verschieden sind. Definition 5.1.1 (Lagrange Polynom) Für paarweise verschiedene Knoten xi , i = 0, . . . , n, ist das k-te Lagrange Polynom definiert durch n Y x − xi Lk (x) := , k = 0, . . . , n. x k − xi i=0,i6=k Für paarweise verschiedene Knoten xi , i = 0, . . . , n, sind die Lagrange Polynome wohldefiniert, da der Nenner niemals Null wird. Man sieht ebenfalls sofort, dass Lk , k = 0, . . . , n, 84 Polynome vom Grad n sind, da jedes Lk sich aus n Produkten von Polynomen ersten Grades zusammensetzt. Desweiteren gilt für k = 0, . . . , n, ( 1 falls j = k, Lk (xj ) = 0 falls j 6= k. Hieraus folgt sofort, dass p(x) = f0 L0 (x) + f1 L1 (x) + . . . + fn Ln (x) = n X fj Lj (x). (5.4) j=0 ein interpolierendes Polynom ist. Mithin haben wir konstruktiv die Existenz einer Lösung nachgewiesen. Weiter gilt Satz 5.1.2 Sind die Stützstellen xi , i = 0, . . . , n, paarweise verschieden, so gibt es genau ein interpolierendes Polynom vom Höchstgrad n. Beweis: Die Existenz eines interpolierenden Polynoms vom Höchstgrad n wurde bereits durch das Polynom p in (5.4) konstruktiv nachgewiesen. Es verbleibt, die Eindeutigkeit nachzuweisen. Es seien also p1 und p2 interpolierende Polynome vom Höchstgrad n. Dann ist das Polynom p = p1 −p2 ebenfalls ein Polynom mit Höchstgrad n und es hat mindestens n + 1 verschiedene Nullstellen p(xi ) = p1 (xi ) − p2 (xi ) = fi − fi = 0, i = 0, . . . , n. Da ein Polynom mit Höchstgrad n, welches verschieden vom Nullpolynom ist, maximal n reelle Nullstellen haben kann (Fundamentalsatz der Algebra), muss p = p1 − p2 das Nullpolynom sein. Also gilt p1 = p2 und die Eindeutigkeit ist gezeigt. 2 Zumindest theoretisch kann das interpolierende Polynom vom Höchstgrad n gemäß (5.4) berechnet werden. Es zeigt sich jedoch, dass die Berechnung gemäß (5.4) für gegebenes x sehr rundungsfehleranfällig und somit numerisch instabil ist. Deshalb sollte man den Wert eines Interpolationspolynoms an der Stelle x in der Praxis nicht über (5.4) berechnen, sondern auf die Newton’schen dividierten Differenzen des nächsten Abschnitts zurückgreifen. Ein weiterer Nachteil von (5.4) ist, dass das Interpolationspolynom komplett neu berechnet werden muss, wenn nachträglich weitere Stützwerte hinzugefügt werden. Bemerkung 5.1.3 (Interpolation und Vandermonde-Matrix) Wertet man die Interpolationsbedingungen (5.2) mit dem Polynom (5.3) aus, folgt 1 x0 x20 · · · xr0 a0 f0 1 x1 x21 · · · xr1 a1 f1 . . . = . (5.5) .. .. . . . . . . . . . . 1 xn x2n · · · xrn ar fn 85 ein lineares Gleichungssystem für die Koeffizienten a0 , . . . , ar . Die Matrix ist die VandermondeMatrix. Die Frage der Lösbarkeit der Interpolationsaufgabe ist also gleichbedeutend mit der Frage nach der Lösbarkeit des linearen Gleichungssystems. Falls der Grad r des Polynoms kleiner als n ist, besitzt das Gleichungssystem im Allgemeinen keine Lösung. Für r = n und verschiedene Stützstellen ist die Vandermonde-Matrix regulär, da wir die Existenz und Eindeutigkeit eines interpolierenden Polynoms bereits bewiesen haben. Es empfiehlt sich nicht, das lineare Gleichungsystem zu lösen, da die Vandermonde-Matrix extrem schlecht konditioniert ist, wie wir bereits in Beispiel 2.4.1 gesehen haben. Für r > n ist die Interpolationsaufgabe ebenfalls stets lösbar, allerdings nicht mehr eindeutig. 5.1.1 Newton’sche dividierte Differenzen Wir leiten ein effizientes Verfahren zur Bestimmung des interpolierenden Polynoms her und stellen das Interpolationspolynom p hierfür in der Form p(x) = α0 + α1 (x − x0 ) + α2 (x − x0 )(x − x1 ) + . . . + αn (x − x0 ) · · · (x − xn−1 ) (5.6) dar. Beachte, dass die αi , i = 0, . . . , n, im Allgemeinen verschieden von den ai , i = 0, . . . , n, in der Darstellung (5.3) sind. Auswertung der Interpolationsbedingungen (5.2) für p führt auf das lineare Gleichungssystem f0 = α0 f1 = α0 + α1 (x1 − x0 ) f2 = α0 + α1 (x2 − x0 ) + α2 (x2 − x0 )(x2 − x1 ) .. .. .. . . . fn = α0 + α1 (xn − x0 ) + . . . + αn (xn − x0 ) · · · (xn − xn−1 ). Dieses Gleichungssystem hat Dreiecksgestalt und besitzt, falls die Stützstellen verschieden sind, eine eindeutige Lösung für α0 , . . . , αn . Die Lösung kann mit Hilfe der Newton’schen dividierten Differenzen f [xi , . . . , xi+k ] angegeben werden, welche wir im Folgenden herleiten möchten. Es bezeichne pxi xi+1 ···xi+k das Polynom vom Höchstgrad k, welches die Stützwerte (xj , fj ), j = i, . . . , i + k, interpoliert. Nach Satz 5.1.2 ist dieses eindeutig bestimmt. Der folgende Satz gibt eine Rekursionsformel zur Berechnung an. Satz 5.1.4 (Rekursionsformel von Neville) Für die Interpolationspolynome pxi xi+1 ···xi+k gilt pxi (x) = fi für alle i = 0, . . . , n und für k ≥ 1 mit i + k ≤ n gilt die Rekursionsformel pxi xi+1 ···xi+k (x) = (x − xi )pxi+1 ···xi+k (x) − (x − xi+k )pxi xi+1 ···xi+k−1 (x) . xi+k − xi 86 Beweis: Der Beweis wird induktiv nach dem Grad k der interpolierenden Polynome geführt. Für k = 0 ist pxi (x) = fi für alle i = 0, . . . , n nach Definition das interpolierende Polynom. Sei die Behauptung richtig für Polynome vom Höchstgrad k − 1 mit k ≥ 1. Definiere das Polynom q(x) := (x − xi )pxi+1 ···xi+k (x) − (x − xi+k )pxi xi+1 ···xi+k−1 (x) , xi+k − xi welches offenbar den Höchstgrad k hat. Nach Induktionsvoraussetzung gilt: q(xi ) = pxi xi+1 ···xi+k−1 (xi ) = fi , (xj − xi )pxi+1 ···xi+k (xj ) − (xj − xi+k )pxi xi+1 ···xi+k−1 (xj ) q(xj ) = xi+k − xi (xj − xi )fj − (xj − xi+k )fj = = fj , j = i + 1, . . . , i + k − 1, xi+k − xi q(xi+k ) = pxi+1 ···xi+k (xi+k ) = fi+k . Also interpoliert q genauso wie pxi xi+1 ···xi+k die Stützstellen (xj , fj ), j = i, . . . , i + k. Da das Interpolationspolynom nach Satz 5.1.2 eindeutig ist, gilt q ≡ pxi xi+1 ···xi+k . 2 Definition 5.1.5 (Newton’sche dividierte Differenz) Der Koeffizient vor xk in pxi xi+1 ···xi+k heißt Newton’sche dividierte Differenz und wird mit f [xi , . . . , xi+k ] bezeichnet. Durch Koeffizientenvergleich erhalten wir aus Satz 5.1.4 eine Rekursionsformel für die Newton’schen dividierten Differenzen. Satz 5.1.6 (Rekursionsformel für Newton’sche dividierte Differenzen) Die Newtonschen dividierten Differenzen sind rekursiv durch f [xi ] := fi , f [xi+1 , . . . , xi+k ] − f [xi , . . . , xi+k−1 ] f [xi , xi+1 , . . . , xi+k ] := xi+k − xi 87 gegeben und können nach folgendem Schema spaltenweise berechnet werden: x0 f [x0 ] = f0 f [x0 , x1 ] x1 f [x1 ] = f1 f [x0 , x1 , x2 ] .. f [x1 , x2 ] x2 xn−2 f [x2 ] = f2 .. . f [xn−2 ] = fn−2 .. . .. . xn−1 f [xn−1 ] = fn−1 ··· .. f [xn−1 , xn−2 ] . f [x0 , . . . , xn ] . f [xn−2 , xn−1 , xn ] f [xn−1 , xn ] xn f [xn ] = fn Der Aufwand des Schemas beträgt 12 n(n + 1) Divisionen und n(n + 1) Subtraktionen. Ein wesentlicher Vorteil der dividierten Differenzen ist, dass zusätzliche Knoten leicht hinzugefügt werden können. Hierfür muss das Schema lediglich um eine weitere Zeile ergänzt werden. Der folgende Satz klärt den Zusammenhang zwischen dem Interpolationspolynom in (5.6) und den Newton’schen dividierten Differenzen, der oberste Eintrag einer jeden Spalte in obigem Schema gibt nämlich den entsprechenden Koeffizienten im Interpolationspolynoms der Form (5.6) an. Satz 5.1.7 Das Interpolationspolynom p zu den Stützstellen (xi , fi ), i = 0, . . . , n, ist durch (5.6) mit αi = f [x0 , . . . , xi ], i = 0, 1, . . . , n gegeben. Beweis: Es gilt die Darstellung p = px0 ···xn = px0 + (px0 x1 − px0 ) + (px0 x1 x2 − px0 x1 ) + . . . + (px0 x1 ···xn − px0 x1 ···xn−1 ). Das Polynom px0 x1 ···xk+1 −px0 x1 ···xk hat den Höchstgrad k+1, die k+1 Nullstellen x0 , . . . , xk , sowie per Definition den Koeffizienten f [x0 , . . . , xk+1 ] vor xk+1 . Da px0 x1 ···xk ein Polynom vom Grad k ist, gilt px0 x1 ···xk+1 (x) − px0 x1 ···xk (x) = f [x0 , . . . , xk+1 ](x − x0 )(x − x1 ) · · · (x − xk ), woraus die Behauptung folgt. Beispiel 5.1.8 Wir betrachten zwei Interpolationsaufgaben mit 5 bzw. 13 Punkten: 88 2 (i) xi fi -2 0 -1 0 0 1 2 10 0 0 (ii) xi fi -6 0 -5 0 -4 0 -3 0 -2 0 -1 0 1 2 3 4 5 6 0 10 0 0 0 0 0 0 Für den Fall (i) ergibt das Schema der dividierten Differenzen: −2 0 0−0 −1−(−2) −1 10−0 0−(−2) 0 10−0 0−(−1) 0 −10−5 1−(−2) 10 = −5 5−(−5) 2−(−2) = −10 5−(−10) 2−(−1) = −10 0−(−10) 2−0 0 0−0 2−1 2 =5 = 10 −10−10 1−(−1) 0−10 1−0 1 =0 = 5 2 =5 =5 =0 0 Das Interpolationspolynom lautet also p(x) = 0 + 0 · (x − (−2)) + 5(x − (−2))(x − (−1)) − 5(x − (−2))(x − (−1))(x − 0) 5 + (x − (−2))(x − (−1))(x − 0)(x − 1) 2 5 = 5(x + 2)(x + 1) − 5x(x + 2)(x + 1) + x(x + 2)(x + 1)(x − 1) 2 5 = (x + 2)(x + 1)(x − 1)(x − 2). 2 Abbildung 5.2 zeigt dieses und das zu (ii) gehörige Interpolationspolynom. Abbildung 5.2: Interpolationspolynome für die Daten (i) (links) und (ii) (rechts). 12 150 10 100 8 50 6 4 0 2 -50 0 -100 -2 -150 -4 -6 -200 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 -6 -4 -2 0 2 4 6 Beispiel (ii) verdeutlicht eine grundsätzliche Problematik bei der Interpolation mit Polynomen: Mit wachsender Anzahl der Stützwerte wächst auch der Polynomgrad und das Interpolationspolynom weist starke Oszillationen auf, die zum Rand hin immer stärker werden. 89 In der Praxis interpoliert man die gegebenen Stützwerte daher häufig abschnittsweise durch Polynome mit nicht zu hohem Grad und heftet“ die abschnittsweise definierten ” Polynome aneinander. Schreibt man das Polynom p in (5.6) in geschachtelter Form p(x) = α0 + (x − x0 ) {α1 + (x − x1 ) [α2 + (. . . + (αn−1 + (x − xn−1 )αn ) . . .)]} , (5.7) so kann es effizient an einer Stelle x ausgewertet werden: Algorithmus 5.1.9 (Horner Algorithmus) Input: x, αi , i = 0, . . . , n. Output: p = p(x). Setze p = αn . for i = (n − 1) : −1 : 0 p := αi + p(x − xi ) end Ist man nur am Wert des Interpolationspolynoms an einer festen Stelle x interessiert so kommt das Neville-Schema zum Einsatz. Im Gegensatz zu den Newton’schen dividierten Differenzen werden nicht die Koeffizienten des Interpolationspolynoms berechnet, sondern nur der Funktionswert an der Stelle x. Es wird wiederum vorausgesetzt, dass die Stützstellen xi , i = 0, . . . , n, verschieden sind. Das Neville-Schema verwendet die Darstellung des Interpolationspolynoms aus Satz 5.1.4 (x − xi )pxi+1 ···xi+k (x) − (x − xi+k )pxi xi+1 ···xi+k−1 (x) xi+k − xi px ···x (x) − pxi xi+1 ···xi+k−1 (x) = pxi+1 ···xi+k (x) + i+1 i+k (x − xi+k ) xi+k − xi pxi xi+1 ···xi+k (x) = und berechnet den gesuchten Funktionswert p(x) = px0 ,...,xn (x) nach folgendem Schema: x0 px0 (x) = f0 px0 x1 (x) x1 px1 (x) = f1 px0 ,x1 ,x2 (x) .. px1 x2 (x) x2 px2 (x) = f2 .. .. . . xn−2 pxn−2 (x) = fn−2 .. . .. . xn−1 pxn−1 (x) = fn−1 pxn−2 ,xn−1 ,xn (x) pxn−1 ,xn (x) xn ··· .. pxn−1 xn−2 pxn (x) = fn 90 . . px0 ,...,xn (x) Beispiel 5.1.10 (Neville-Schema) Gesucht ist der Funktionswert an der Stelle x = 2 des Interpolationspolynoms zu den Daten xi 0 1 3 fi 1 3 2 Das Neville-Schema liefert: 0 1 3+ 3−1 (2 1−0 − 1) = 5 5 2 1 3 2+ 2−3 (2 3−1 − 3) = + 5/2−5 (2 3−0 − 3) = 10 3 5 2 3 2 Es gilt also p(2) = 10 . 3 5.1.2 Fehlerdarstellung interpolierender Polynome Der folgende Satz gibt Auskunft über die Approximationsgüte des interpolierenden Polynoms vom Höchstgrad n an eine Funktion f zwischen den Stützstellen in (5.1). Satz 5.1.11 Seien f : [a, b] → R (n + 1)-mal stetig differenzierbar und die Stützstellen xi ∈ [a, b], i = 0, . . . , n, verschieden. Für das Interpolationspolynom p vom Höchstgrad n mit p(xi ) = f (xi ), i = 0, . . . , n, gilt f (x) − p(x) = n Y 1 f (n+1) (ξ) · (x − xi ) (n + 1)! i=0 (5.8) mit einer von x abhängigen Stelle ξ mit min{x0 , . . . , xn , x} < ξ < max{x0 , . . . , xn , x}. Beweis: Für x ∈ {x0 , . . . , xn } ist die Behauptung offenbar wahr. Sei nun x 6∈ {x0 , . . . , xn } fest gewählt. Definiere n Y g(t) := f (t) − p(t) − α (t − xi ), i=0 wobei α so bestimmt sei, dass g(x) = 0 gilt. Dies ist wegen x 6∈ {x0 , . . . , xn } möglich. Dann besitzt g auf dem Intervall I := [min{x0 , . . . , xn , x}, max{x0 , . . . , xn , x}] 91 die n + 2 Nullstellen x0 , x1 , . . . , xn , x. Desweiteren ist g (n + 1)-mal stetig differenzierbar. Nach dem Satz von Rolle1 hat g 0 in I noch mindestens n + 1 Nullstellen, g 00 mindestens n Nullstellen, usw. Schließlich hat g (n+1) noch mindestens eine Nullstelle ξ ∈ I. Als Polynom vom Höchstgrad n gilt p(n+1) ≡ 0, also 0 = g (n+1) (ξ) = f (n+1) (ξ) − α(n + 1)! und somit α= f (n+1) (ξ) . (n + 1)! Damit ergibt sich aus g(x) = 0 die Behauptung gemäß f (x) − p(x) = α n n Y f (n+1) (ξ) Y · (x − xi ). (x − xi ) = (n + 1)! i=0 i=0 2 Die Fehlerdarstellung in (5.8) hängt von der Funktion f und dem Polynom n Y ω(x) = (x − xi ) i=0 ab, wobei letzteres nur von den Stützstellen xi , i = 0, . . . , n, abhängt. Insbesondere folgt aus ihr die Abschätzung 1 kf (n+1) k∞ · |ω(x)| (n + 1)! 1 ≤ kf (n+1) k∞ · kωk∞ (n + 1)! |f (x) − p(x)| ≤ (5.9) für alle x ∈ I := [min{x0 , . . . , xn }, max{x0 , . . . , xn }]. Hierin bezeichnet k · k∞ die Supremumsnorm, die für eine stetige Funktion g : [a, b] → R definiert ist durch kgk∞ := max |g(x)|. a≤x≤b Falls die Stützstellen xi , i = 0, . . . , n, frei wählbar sind, kann der Approximationsfehler in (5.9) für gegebenes f minimiert werden, indem die Stützstellen xi , i = 0, . . . , n, optimal gewählt werden. Optimal bedeutet in diesem Zusammenhang, dass kωk∞ = max min{x0 ,...,xn }≤x≤max{x0 ,...,xn } |ω(x)| durch geeignete Wahl von paarweise verschiedenen Stütztellen xi , i = 0, . . . , n, minimal wird. Es gilt also, paarweise verschiedene Stützstellen xi , i = 0, . . . , n, zu bestimmen, die ω 1 Sei g in [a, b] stetig, in (a, b) differenzierbar und es gelte g(a) = g(b). Dann gibt es ξ ∈ (a, b) mit g (ξ) = 0. 0 92 bzgl. der Supremumsnorm minimieren. Es zeigt sich, dass die Nullstellen der sogenannten Tschebyscheff-Polynome Tn (x) := cos(n · arccos(x)), n = 0, 1, 2, . . . , genau dies leisten (wenn man das Problem zuvor auf das Intervall [−1, 1] transformiert). Die Tschebyscheff-Polynome bis zum Grad 5 lauten wie folgt: 1.5 T0 (x) = 1 1 x 2*x**2-1 4*x**3-3*x 8*x**4-8*x**2+1 16*x**5-20*x**3+5*x 1 T1 (x) = x 0.5 2 T2 (x) = 2x − 1 0 T3 (x) = 4x3 − 3x -0.5 T4 (x) = 8x4 − 8x2 + 1 T5 (x) = 16x5 − 20x3 + 5x -1 -1 -0.5 0 0.5 1 Beispiel 5.1.12 Betrachte 1 . + 0.01 Die folgende Abbildung zeigt das interpolierende Polynom peq zu den äquidistanten Stützstellen xi = −1 + i · h, i = 0, . . . , 10, mit h = 2/10, sowie das interpolierende Polynom zu den Tschebyscheff-Stützstellen 2i + 1 xi = cos π , i = 0, . . . , 10. 2n + 2 f (x) = x2 500 equidistant tschebyscheff f 400 300 200 100 0 -100 -1 -0.8 -0.6 -0.4 -0.2 0 93 0.2 0.4 0.6 0.8 1 Es ist deutlich zu sehen, dass das auf Tschebyscheff-Stützstellen basierende Interpolationspolynom die Funktion f besser (im Sinne der Supremumsnorm) approximiert als das auf äquidistanten Stützstellen basierende Interpolationspolynom. Insbesondere am Rand oszilliert letzteres stark. 5.2 Splineinterpolation Im Computer Aided Design (CAD) oder in anderen Interpolationsaufgaben sind vor allem visuell ansprechende“ interpolierende Funktionen erwünscht. Wir haben in Beispiel 5.1.8 ” gesehen, dass Polynominterpolation hier nur bedingt geeignet ist, da mit wachsendem Polynomgrad starke Oszillationen auftreten können. Polynome mit niedrigem Polynomgrad hingegen sehen angenehm aus. Dies motiviert die Verwendung von sogenannten Splines. Ein Spline besteht dabei aus abschnittsweise definierten Polynomen, die an den Stützstellen unter Beachtung von Glattheitsanforderungen zusammengesetzt werden. Definition 5.2.1 (Spline vom Grad k) Eine Funktion s heißt Spline vom Grad k auf [a, b], k ∈ N, falls s folgende Bedingungen erfüllt: • s ist (k − 1)-mal stetig differenzierbar in [a, b]. • Es gibt eine Partition a = x0 < x1 < . . . < xn = b von [a, b], so dass s in jedem Teilintervall [xi , xi+1 ], 0 ≤ i ≤ n − 1, ein Polynom vom Höchstgrad k ist. Für G = {x0 < x1 < . . . < xn } bezeichnet Sk (G) := {s : [x0 , xn ] → R | s ist Spline vom Grad k auf [x0 , xn ]} den Raum der Splines vom Grad k mit Stützstellen in G. Ein interpolierender Spline vom Grad k für die Daten (xi , fi ), i = 0, . . . , n, ist ein Spline s vom Grad k, der die Interpolationsbedingungen s(xi ) = fi , i = 0, . . . , n, erfüllt. Beispiel 5.2.2 (Spline vom Grad 1) Eine stetige, stückweise lineare Funktion auf [a, b] ist ein Spline vom Grad 1. Offenbar ist der interpolierende Spline s vom Grad 1 für (xi , fi ), i = 0, . . . , n, eindeutig bestimmt durch s0 (x), x 0 ≤ x ≤ x1 , s1 (x), x 1 ≤ x ≤ x2 , s(x) = . .. sn−1 (x), xn−1 ≤ x ≤ xn 94 mit x − xi si (x) = s[xi ,xi+1 ] (x) = fi + (fi+1 − fi ), xi+1 − xi a = x0 x1 x2 x3 x4 x5 i = 0, . . . , n − 1. x6 x 7 = b Im Folgenden werden wir uns speziell mit Splines vom Grad 3, den sogenannten kubischen Splines, beschäftigen. Kubische Splines werden insbesondere in der Computergrafik verwendet, um Funktionsverläufe zu erzeugen, die für das menschliche Auge glatt“ erscheinen. Der Hintergrund ist, dass das menschliche Auge Unstetigkeiten in der ” Krümmung (also im Wesentlichen in der zweiten Ableitung) noch wahrnehmen kann, während Funktionen mit stetiger Krümmung (also stetiger zweiter Ableitung) als glatt aufgefasst werden. Zur Motivation betrachten wir folgendes Problem: Ein Balken soll an den Stellen (5.1) durch Lager so eingespannt werden, dass dort nur Kräfte senkrecht zur Biegelinie wirken. Er wird dann eine Lage annehmen, die durch eine minimale Biegeenergie E(s) charakterisiert ist, wobei die Funktion s die Biegekurve bezeichnet. Die Biegeenergie ist durch Z xn s00 (x)2 E(s) = c p 5 dx x0 1 + s0 (x)2 mit einer Konstanten c gegeben. Die Funktion s00 (x) κ(x) = p 3 1 + s0 (x)2 beschreibt gerade die Krümmung der Funktion s. Für |s0 (x)| << 1 gilt näherungsweise Z xn E(s) ≈ c s00 (x)2 dx. (5.10) x0 Die Aufgabe besteht nun darin, eine Funktion s zu finden, die das Integral in (5.10) und damit auch (näherungsweise) die Biegeenergie unter den Nebenbedingungen s(xi ) = 95 1 0 0 1 0 01 1 0 1 0 1 0 1 11 00 0 1 00 11 00 11 1 0 01 1 0 0 1 Abbildung 5.3: Spline vom Grad 3 (kubischer Spline): Kurve mit minimaler Krümmung durch gegebene Punkte. fi , i = 0, . . . , n minimiert. Dies leisten Splines, wie der folgende Satz aussagt, den wir hier nicht beweisen wollen. Satz 5.2.3 (Minimaleigenschaft interpolierender kubischer Splines) Seien (xi , fi ), i = 0, . . . , n gegeben und s ein interpolierender kubischer Spline, sowie g eine beliebige zweimal stetig differenzierbaren Funktion mit g(xi ) = fi für alle i = 0, . . . , n. Ferner gelte eine der folgenden Bedingungen: (a) s00 (x0 ) = s00 (xn ) = 0 (natürlicher Spline); (b) f0 = fn , s0 (x0 ) = s0 (xn ), s00 (x0 ) = s00 (xn ), g 0 (x0 ) = g 0 (xn ) (periodischer Spline); (c) f00 = g 0 (x0 ) = s0 (x0 ), fn0 = g 0 (xn ) = s0 (xn ) für gegebenes f00 und fn0 (Spline mit Hermitescher Datenvorgabe). Dann ist der interpolierende kubische Spline eindeutig. Ferner gilt: Z xn Z xn 00 2 s (x) dx ≤ g 00 (x)2 dx. x0 x0 Konstruktion interpolierender kubischer Splines Ziel ist es nun, einen die Stützstellen (5.1) interpolierenden kubischen Spline zu konstruieren. Nach Definition ist ein interpolierender kubischer Spline s für i = 0, . . . , n − 1 im Intervall [xi , xi+1 ] ein Polynom si vom Höchstgrad 3, etwa s0 (x), x 0 ≤ x ≤ x1 , s1 (x), x 1 ≤ x ≤ x2 , s(x) = (5.11) . .. sn−1 (x), xn−1 ≤ x ≤ xn 96 mit 1 1 si (x) = s[xi ,xi+1 ] (x) = αi + βi (x − xi ) + γi (x − xi )2 + δi (x − xi )3 . 2 6 (5.12) Weiterhin müssen die folgenden Bedingungen erfüllt sein: s0 (x0 ) si (xi+1 ) = si+1 (xi+1 ) sn−1 (xn ) s0i (xi+1 ) s00i (xi+1 ) = = = = = f0 , fi+1 , fn , s0i+1 (xi+1 ), s00i+1 (xi+1 ), i = 0, 1, . . . , n − 2, i = 0, 1, . . . , n − 2, i = 0, 1, . . . , n − 2. Aus der Darstellung (5.12) folgt γi = s00i (xi ) = s00 (xi ), i = 0, 1, . . . , n − 1. Zusätzlich definieren wir γn := s00 (xn ) = s00n−1 (xn ) und hi+1 := xi+1 − xi , i = 0, 1, . . . , n − 1. Zunächst wird gezeigt, wie si (und damit auch s) berechnet werden kann, wenn die sogenannten Momente γi , i = 0, . . . , n bekannt sind. • Aus s00i (xi+1 ) = s00i+1 (xi+1 ), i = 0, . . . , n − 1, folgt s00i (xi+1 ) = γi + δi hi+1 = γi+1 = s00i+1 (xi+1 ). Auflösen nach δi liefert die Darstellung δi = γi+1 − γi , hi+1 i = 0, . . . , n − 1. (5.13) • Aus si (xi ) = fi , i = 0, . . . , n − 1, folgt sofort i = 0, . . . , n − 1. αi = f i , Aus si (xi+1 ) = fi+1 , i = 0, . . . , n − 1, folgt si (xi+1 ) = fi + βi hi+1 + γi 2 γi+1 − γi 3 hi+1 + h = fi+1 . 2 6hi+1 i+1 Auflösen nach βi liefert die Darstellung βi = fi+1 − fi 2γi + γi+1 − hi+1 , hi+1 6 97 i = 0, . . . , n − 1. (5.14) Zusammenfassend sind die Koeffizienten αi , βi und δi in (5.12) nach (5.13) und (5.14) gegeben durch hi+1 = xi+1 − xi , i = 0, . . . , n − 1, i = 0, . . . , n − 1, αi = fi , fi+1 − fi 2γi + γi+1 − hi+1 , i = 0, . . . , n − 1, hi+1 6 γi+1 − γi , i = 0, . . . , n − 1, = hi+1 βi = δi falls die Momente γi , i = 0, . . . , n, bekannt sind. Es verbleibt, die Momente γi , i = 0, . . . , n, zu berechnen. Hierzu werden die Bedingungen s0i (xi ) = s0i−1 (xi ), i = 1, . . . , n − 1, ausgenutzt. Diese führen auf 1 βi = βi−1 + γi−1 hi + δi−1 h2i , 2 i = 1, . . . , n − 1. Einsetzen von βi , βi−1 und δi−1 aus (5.13) und (5.14) liefert fi+1 − fi 2γi + γi+1 fi − fi−1 2γi−1 + γi γi − γi−1 2 − hi+1 = − hi + γi−1 hi + hi , hi+1 6 hi 6 2hi und dies ist äquivalent zu hi γi−1 + 2(hi + hi+1 )γi + hi+1 γi+1 = 6 fi+1 − fi fi − fi−1 − hi+1 hi (5.15) für i = 1, 2, . . . , n − 1. (5.15) definiert insgesamt n − 1 lineare Gleichungen für die n + 1 Momente γi , i = 0, . . . , n. Wir benötigen also noch zwei weitere Bedingungen, um eine eindeutige Lösung zu ermöglichen. Hier kommt eine der Bedingungen (a), (b) oder (c) aus Satz 5.2.3 ins Spiel und je nachdem, welche Bedingung gewählt wird, entsteht ein natürlicher Spline, ein periodischer Spline oder ein Spline mit hermitescher Datenvorgabe. In jedem Fall führt die Auswertung von (5.15) zusammen mit einer der Bedingungen (a), (b) oder (c) aus Satz 5.2.3 jeweils auf ein lineares Gleichungsystem der Form A · γ = d, γ := (γ0 , . . . , γn )> , d = (d0 , . . . , dn )> (5.16) für die Momente γi , i = 0, . . . , n. Wir diskutieren die einzelnen Fälle und die daraus entstehenden linearen Gleichungssysteme. Natürlicher Spline Es wird s00 (x0 ) = s00 (xn ) = 0 gefordert. Per Definition der Momente führt dies auf die beiden Gleichungen γ0 = 0 und γn = 0. Zusammen mit (5.15) ist das resultierende 98 Gleichungssystem (5.16) dann gegeben durch 1 h 2(h + h ) h2 1 1 2 h2 2(h2 + h3 ) h3 .. .. .. A = . . . hn−2 2(hn−2 + hn−1 ) hn−1 hn−1 2(hn−1 + hn ) hn 1 0, für i = 0 und i = n, di = fi+1 − fi fi − fi−1 6 , für i = 1, . . . , n − 1. − hi+1 hi , Periodischer Spline Für einen periodischen Spline muss zunächst f0 = fn gelten. Weiter wird gefordert, dass s0 (x0 ) = s0 (xn ) und s00 (x0 ) = s00 (xn ) gelten. Dies führt auf 0 = s00 (x0 ) − s00 (xn ) = γ0 − γn , 1 0 = s0 (x0 ) − s0 (xn ) = β0 − (βn−1 + γn−1 hn + δn−1 h2n ). 2 Einsetzen von β0 , βn−1 und δn−1 aus (5.13) und (5.14) liefert f1 − f0 2γ0 + γ1 fn − fn−1 2γn−1 + γn γn − γn−1 2 − h1 − + hn − γn−1 hn − hn h1 6 hn 6 2hn und damit f1 − f0 fn − fn−1 − 2h1 γ0 + h1 γ1 + hn γn−1 + 2hn γn = 6 . h1 hn 0= Zusammen mit (5.15) ist das resultierende Gleichungssystem (5.16) dann gegeben durch 1 −1 h 2(h + h ) h2 1 1 2 h 2(h h 2 2 + h3 ) 3 .. .. .. A = , . . . h 2(h + h ) h n−2 n−2 n−1 n−1 hn−1 2(hn−1 + hn ) hn 2h1 di = h1 hn 0, für i = 0, fi+1 − fi fi − fi−1 − , für i = 1, . . . , n − 1, hi hi+1 f1 − f0 fn − fn−1 6 − , für i = n. h1 hn 6 99 2hn Spline mit Hermitescher Datenvorgabe In x0 und xn werden neben den Funktionswerten s(x0 ) = f0 und s(xn ) = fn auch Ableitungen s0 (x0 ) = f00 und s0 (xn ) = fn0 vorgegeben. Mit den obigen Bezeichnungen führt dies auf die beiden Gleichungen f1 − f0 2γ0 + γ1 − h1 , h1 6 1 fn0 = s0 (xn ) = βn−1 + γn−1 hn + δn−1 h2n 2 fn − fn−1 2γn−1 + γn γn − γn−1 2 = − hn + γn−1 hn + hn . hn 6 2hn f00 = s0 (x0 ) = β0 = und schließlich zu 2h1 γ0 + h1 γ1 hn γn−1 + 2hn γn f1 − f0 0 = 6 − f0 h1 fn − fn−1 0 . = 6 fn − hn Zusammen mit (5.15) ist das resultierende Gleichungssystem (5.16) dann gegeben durch 2h1 h1 h 2(h + h ) h2 1 1 2 h2 2(h2 + h3 ) h3 .. .. .. A = . . . hn−2 2(hn−2 + hn−1 ) hn−1 hn−1 2(hn−1 + hn ) hn hn 2hn f1 − f0 6 − f00 , für i = 0, h1 f − f fi − fi−1 i+1 i di = 6 − , für i = 1, . . . , n − 1, hi+1 hi fn − fn−1 0 6 fn − , für i = n. hn , Bemerkung 5.2.4 (Eindeutigkeit der kubischen Splines) Die Matrizen A im Fall natürlicher Splines oder solcher mit Hermitscher Datenvorgabe sind strikt diagonaldominant und daher regulär. Im Fall periodischer Splines kann man, wegen γ0 = γn eine dieser beiden Variablen ersetzen und erhält dann ebenfalls eine strikt diagonaldominante Matrix. Daher ist der interpolierende kubische Spline unter jeder der Zusatzbedingungen eindeutig. 100 5.2.1 Fehlerdarstellung interpolierender kubischer Splines Wie bei der Polynominterpolation interessieren wir uns für Approximationseigenschaften von kubischen Splines bezüglich einer gegebenen Funktion f . Zusätzlich zur Supremumsnorm k · k∞ definieren wir die L2 -Norm Z kgk2 := b 1/2 |g(t)| dt 2 a für eine quadratintegrierbare Funktion g : [a, b] → R. Es gilt folgender Satz, der hier nicht bewiesen werden soll: Satz 5.2.5 Gegeben seien a = x0 < x1 < . . . < xn = b und eine zweimal stetig differenzierbare Funktion f : [a, b] → R. Der interpolierende kubische Spline s mit s(xi ) = f (xi ), i = 0, . . . , n, und einer der Bedingungen (a), (b) oder (c) aus Satz 5.2.3 erfüllt die Fehlerabschätzungen 1 00 kf k2 h2max , 4 1 00 3/2 ≤ kf k2 hmax , 2 kf − sk2 ≤ kf − sk∞ wobei hmax = max (xi+1 − xi ) sei. i=0,...,n−1 101 Kapitel 6 Numerische Integration Viele praktisch relevante Anwendungen, wie z.B. die Lösung einer elliptischen partiellen Differentialgleichung mit der Finite-Element-Methode oder die Fourieranalyse, erfordern die Auswertung von bestimmten Integralen der Form Zb a, b ∈ R, a < b. f (x)dx, I[f ] := (6.1) a Da diese Integrale oftmals nicht exakt berechnet werden können, z.B. nicht bei Zπ/2 1− x 2 r 2 1/2 sin θ dθ, 0 werden numerische Methoden zur approximativen Berechnung von (6.1) benötigt. Hierzu verwendet man in der Regel den Ansatz Zb f (x)dx ≈ N X wi f (xi ) i=0 a mit N ∈ N0 , Gewichten wi und Stützstellen a ≤ x0 ≤ x1 ≤ . . . ≤ xN ≤ b. Die Aufgabe besteht nun darin, geeignete Werte für wi und xi zu finden, so dass die sogenannte Quadraturformel Q[f ] := N X wi f (xi ) (6.2) i=0 den Wert des Integrals möglichst gut approximiert. Beispiel 6.0.1 (Trapezregel) Die Trapezregel basiert auf der linearen Interpolation der Stützstellen (a, f (a)) und (b, f (b)), vgl. Abbildung 6.1. 102 f(b) f(a) f(x) a b Abbildung 6.1: Trapezregel Der Integrand f wird dann durch das lineare Interpolationspolynom p1 (x) = f (a) + x−a (f (b) − f (a)) b−a (6.3) ersetzt. Das resultierende Integral über p1 kann berechnet werden und lautet Z b Z b b−a p1 (x)dx = f (x)dx ≈ (f (a) + f (b)) =: T1 [f ]. 2 a a Geometrisch ist dies gerade der Flächeninhalt des Trapezes, welches durch die Punkte (a, 0), (a, f (a)), (b, f (b)), (b, 0) gegeben ist, vgl. Abbildung 6.1. Die Trapezregel T1 [f ] ist also eine spezielle Quadraturformel (6.2) mit N = 1, w0 = w1 = b−a , x0 = a, x1 = b. 2 Ein Maß für die Güte der Approximation, welches später noch häufiger auftreten wird, ist unter anderem der Exaktheitsgrad einer Quadraturformel. Definition 6.0.2 (Exaktheitsgrad einer Quadraturformel) Eine Quadraturformel Q[f ] hat Exaktheitsgrad p, falls alle Polynome bis zum Grad p exakt integriert werden. Durch einfaches Nachrechnen zeigt sich, dass die Trapezregel den Exaktheitsgrad eins hat. Polynome ersten Grades werden also durch die Trapezregel exakt integriert. Wie Abbildung 6.1 bereits erahnen lässt, ist die Approximation von f durch p1 im Falle der Trapezregel im Allgemeinen nicht sehr genau. Häufig verwendet man daher sogenannte 103 zusammengesetzte Quadraturformeln. Dabei wird das Intervall [a, b] (in der Regel) äquidistant mit Schrittweite h = (b − a)/n für n ∈ N unterteilt und in jedem Teilintervall [xi , xi+1 ] mit xi = a + ih, i = 0, . . . , n − 1, wird dann eine Quadraturformel der Form (6.2) auf Z xi+1 Ii [f ] = f (x)dx, i = 0, . . . , n − 1, xi angewendet. Eine Approximation von I[f ] ergibt sich dann durch Summation über alle Teilintervalle: n−1 Z xi+1 n−1 X X I[f ] = f (x)dx ≈ Q[f |[xi ,xi+1 ] ] = Qn [f ]. i=0 xi i=0 Beispiel 6.0.3 (zusammengesetzte Trapezregel) Das Intervall [a, b] wird durch Einführung von äquidistanten Stützstellen xi = a + ih, i = 0, . . . , n, h = b−a n (6.4) in Teilintervalle [xi , xi+1 ] zerlegt. Die Trapezregel wird dann auf jedes Teilintervall angewendet, vgl. Abbildung 6.2. f6 f5 f0 f4 f1 f2 a = x 0 x1 x2 f3 x3 x4 x5 b = x6 Abbildung 6.2: Zusammengesetzte Trapezregel für n = 6. Dies liefert die zusammengesetzte Trapezregel oder iterierte Trapezregel n−1 X h h f (xi ) + f (b). Tn [f ] := f (a) + h 2 2 i=1 Die zusammengesetzte Trapezregel entspricht der Approximation von f durch eine stetige, stückweise lineare Funktion. Sie ist ebenfalls eine Quadraturformel der Form (6.2) mit N = n, w0 = wn = h , wi = h, i = 1, . . . , n − 1, xi = a + ih, i = 0, . . . , n. 2 104 Bei zusammengesetzten Quadraturformeln spielt die Konvergenzordnung bzgl. der Schrittweite h eine große Rolle. Definition 6.0.4 (Konvergenzordnung zusammengesetzte Quadraturformel) Eine zusammengesetzte Quadraturformel Qn [f ] hat die Konvergenzordnung p (ist von Ordnung p), falls |Qn [f ] − I[f ]| = O(hp ) mit h = (b − a)/n gilt. 6.1 Interpolatorische Quadraturformeln Die Idee der interpolatorischen Quadraturformeln besteht darin, die Funktion f im Integranden von I[f ] durch ein interpolierendes Polynom mit Höchstgrad N zu ersetzen, welches f an den Stützstellen a ≤ x0 < x1 < . . . < xN −1 < xN ≤ b interpoliert. Für dieses interpolierende Polynom kann das Integral dann leicht ausgewertet werden. Als Spezialfall haben wir die Trapezregel bereits kennen gelernt. Allgemeiner wird nun ein Polynom pN höchstens N -ten Grades zur Interpolation der N +1 Stützpunkte (xi , f (xi )), i = 0, . . . , N , verwendet. Dieses Interpolationspolynom ist nach Satz 5.1.2 eindeutig bestimmt und besitzt nach (5.4) die Darstellung pN (x) = N X f (xk )Lk (x), k=0 wobei Lk (x) das k-te Lagrange Polynom bezeichnet. Ersetzen von f durch pN in I[f ] in (5.1) liefert die Approximation Z b Z f (x)dx ≈ I[f ] = a b pN (x)dx, = a N X Z f (xk ) k=0 mit den Koeffizienten a | b Lk (x)dx = {z } N X k=0 =:wk b Z wk = Lk (x)dx. a Insgesamt erhalten wir daraus die interpolatorischen Quadraturformeln. Definition 6.1.1 (interpolatorische Quadraturformeln) Seien a < b, N ∈ N0 und Stützstellen xi , i = 0, . . . , N , mit a ≤ x0 < x1 < . . . < xN −1 < xN ≤ b 105 wk f (xk ) gegeben. Die Quadraturformel Q[f ] = N X Z wk f (xk ), b Lk (x)dx wk = (6.5) a k=0 heißt interpolatorische Quadraturformel. Im Fall äquidistanter Stützstellen b−a N heißt eine interpolatorische Quadraturformel geschlossene Newton-Cotes-Formel, wobei hierbei (a, f (a)) und (b, f (b)) interpoliert werden. Im Fall äquidistanter Stützstellen xi = a + ih, i = 0, . . . , N, h= b−a N +2 heißt eine interpolatorische Quadraturformel offene Newton-Cotes-Formel, wobei hierbei (a, f (a)) und (b, f (b)) nicht interpoliert werden. xi = a + (i + 1)h, i = 0, . . . , N, h= Gemäß Definition 6.1.1 ist die Trapezregel also eine geschlossene Newton-Cotes-Formel. Ein Beispiel für eine offene Newton-Cotes-Formel ist die folgende Mittelpunktsregel. Beispiel 6.1.2 (Mittelpunktsregel) Sei [a, b], a < b, gegeben. Wähle N = 0 und x0 = (a + b)/2 (Mittelpunkt des Intervalls). Das interpolierende Polynom 0-ten Grades zum Stützwert (x0 , f (x0 )) ist gegeben durch p0 (x) = f (x0 ). Damit wird das Integral I[f ] in (6.1) approximiert durch die sogenannte Mittelpunktsregel Z b a+b I[f ] ≈ M0 [f ] = p0 (x)dx = (b − a)f , 2 a vgl. Abbildung 6.3. f(x) a (a+b)/2 b Abbildung 6.3: Mittelpunkstregel 106 Dies ist eine offenee Newton-Cotes-Formel der Form (6.2) mit N = 0, w0 = b − a, x0 = (a + b)/2. Die daraus resultierende zusammengesetze Quadraturformel lautet n−1 n−1 X X xi + xi+1 xi + xi+1 =h f Mn [f ] = (xi+1 − xi )f 2 2 i=0 i=0 mit h = (b − a)/n, vgl. Abbildung 6.4. a = x0 x1 x2 x3 x4 x5 b = x6 Abbildung 6.4: Iterierte Mittelpunktsregel für n = 6. Wir interessieren uns nun für den Approximationsfehler interpolatorischer Quadraturformeln, wofür wir auf Satz 5.1.11 zurück greifen. Satz 6.1.3 Seien f : [a, b] → R (N + 1)-mal stetig differenzierbar und a ≤ x0 < x1 < . . . < xN ≤ b gegeben. Dann gilt 1 max f (N +1) (x) · |I[f ] − Q[f ]| ≤ (N + 1)! x∈[a,b] Z bY N |x − xi |dx a i=0 für die interpolatorische Quadraturformel Q[f ] in (6.5). Insbesondere besitzt die interpolatorische Quadraturformel Q[f ] (mindestens) den Exaktheitsgrad N . Beweis: Es gilt Z b Z b f (x) − pN (x)dx ≤ |f (x) − pN (x)|dx. |I[f ] − Q[f ]| = a a 107 Die Fehlerdarstellung des interpolierendes Polynoms pN aus Satz 5.1.11 liefert |f (x) − pN (x)| ≤ N Y 1 max |f (N +1) (x)| · |x − xi | (N + 1)! x∈[a,b] i=0 für alle x ∈ [a, b]. Integration liefert die erste Behauptung. Da f (N +1) ≡ 0 für alle Polynome mit Höchstgrad N gilt, folgt aus der Fehlerdarstellung |I[f ] − Q[f ]| = 0, so dass Polynome bis zum Höchstgrad N exakt integriert werden. 2 Der maximale Exaktheitsgrad von Q[f ] kann größer als N ausfallen. So besitzt die Simpsonregel (geschlossene Newton-Cotes-Formel mit N = 2) sogar den Exaktheitsgrad N + 1 = 3. Tabelle 6.1 fasst einige gebräuchliche geschlossene Newton-Cotes-Formeln und deren Exaktheitsgrad zusammen. Hierbei gilt xi = a + ih, i = 0, . . . , N , h = (b − a)/N . Tabelle 6.1: Einige gebräuchliche geschlossene Newton-Cotes-Formeln N Gewichte Name Fehler E.-Grad 1 2 3 4 3 w0 = w1 = h/2 Trapezregel − h12 f (2) (ξ) w0 = w2 = h/3, 5 Simpsonregel − h90 f (4) (ξ) w1 = 4h/3 w0 = w3 = 3h/8, 5 f (4) (ξ) 3/8-Regel − 3h 80 w1 = w2 = 9h/8 w0 = w4 = 28h/90, 7 w1 = w3 = 128h/90, Milneregel − 8h f (6) (ξ) 945 w2 = 48h/90 1 3 3 5 Tabelle 6.2 fasst einige gebräuchliche offene Newton-Cotes-Formeln und deren Exaktheitsgrad zusammen. Hierbei gilt xi = a + (i + 1)h, i = 0, . . . , N , h = (b − a)/(N + 2). Tabelle 6.2: Einige gebräuchliche offene Newton-Cotes-Formeln N Gewichte Name Fehler E.-Grad 0 1 2 3 w0 = 2h w0 = w1 = 3h/2, w0 = w2 = 8h/3, w1 = −4h/3 w0 = w3 = 55h/24, w1 = w2 = 5h/24 Mittelpunktsregel 108 h3 (2) f (ξ) 3 3h3 (2) f (ξ) 4 1 1 28h5 (4) f (ξ) 90 3 95h5 (4) f (ξ) 144 3 Für jede Regel gilt jeweils ξ ∈ [a, b] mit eigenem ξ. 6.2 Extrapolationsverfahren Sei f : [a, b] → R Riemann integrierbar und (2m+2) Mal stetig differenzierbar. Betrachten wir die zusammengesetzte Trapezregel als Funktion der Schrittweite h := b−a n n−1 X h h f (a + ih) + f (b). T (h) := f (a) + h 2 2 i=1 Da f Riemann-integrierbar ist, gilt Z f (x) dx. lim T (h) = I[f ] = h→0 b (6.6) a Ferner kann man unter der Voraussetzung der (2m+2)-fachen stetigen Differenzierbarkeit die asymptotische Fehlerdarstellung der zusammengesetzten Trapezregel T (h) = I[f ] + γ2 h2 + γ4 h4 + . . . + γ2m h2m + O(h2m+2 ) mit von h unabhängigen Konstanten γ2i , i = 1, . . . , m zeigen. Beispiel 6.2.1 Für die zusammengesetzte Trapezsumme zu verschiedenen Schrittweiten erhält man h0 := h h h1 := 2 T (h) = I[f ] + γ2 h2 + γ4 h4 + O(h6 ) γ2 γ4 T h2 = I[f ] + h2 + h4 + O(h6 ). 4 16 Subtrahiert man nun vom Vierfachen der zweiten die erste Gleichung, so folgt 3 h − T (h) = 3I[f ] − γ4 h4 + O(h6 ). 4T 2 4 Damit hat man also durch geeignete Zusammensetzung vergenzordnung 2 mit 1 h T̂ (h) := 4T − T (h) = I[f ] − 3 2 zweier Trapezsummen der Kon- 1 γ4 h4 + O(h6 ). 4 eine Approximation des Integrals I[f ] mit der verbesserten Konvergenzordnung 4 erhalten. , ni ∈ N Die Idee aus Beispiel 6.2.1 lässt sich verallgemeinern. Zu der Folge hi = b−a ni wollen wir die Trapezsummen T (hi ), i ∈ N berechnen. Sind T (h0 ), T (h1 ), . . . , T (hk ) bereits bekannt, so bestimme das Polynom T0,1,...,k (h2 ) mit T0,1,...,k (h2j ) = T (hj ) ∀j = 0, 1, . . . , k. 109 Wir extrapolieren das Polynom in den Nullpunkt und wählen T0,1,...,k (0) als neue Näherung für I[f ], vgl. (6.6). Zur Berechnung: Bezeichne Ti,i+1,...,i+k (h2 ) das Polynom mit Ti,i+1,...,i+k (h2j ) = T (hj ) ∀j = i, i + 1, . . . , i + k. Setze Ti,i+k := Ti,i+1,...,i+k (0). T (h0 ) T (h3 ) T0,1,...,3 (h2 ) T (h2 ) T (h ) 1 T0,3 0 h3 h2 h1 h0 h Abbildung 6.5: Prinzip der Extrapolation für k = 3. Dann liefert das Neville-Schema aus Abschnitt 5.1.1 (mit xi = h2i , x = 0) Ti,i = T (hi ), (Ti+1,i+k − Ti,i+k−1 ) (0 − h2i+k ) h2i+k − h2i (Ti+1,i+k − Ti,i+k−1 ) = Ti+1,i+k + . 2 hi − 1 hi+k Ti,i+k = Ti+1,i+k + Speziell für die in der Praxis häufig verwendete Romberg-Folge hi+1 := wegen hi+k = 21k hi , Ti,i+k = Ti+1,i+k + hi 2 hat man, (Ti+1,i+k − Ti,i+k−1 ) . 4k − 1 Ein wesentlicher Vorteil der Romberg-Folge liegt darin, dass man T (hi+1 ) unter Verwen110 dung von T (hi ) berechnen kann. Mit hi = T (hi+1 ) = hi+1 2 hi = 4 hi = 4 b−a ni und ni+1 = 2ni bzw. hi = 2hi+1 folgt ni+1 −1 f (a) + 2 X ! f (a + jhi+1 ) + f (b) j=1 f (a) + 2 2n i −1 X f j=1 f (a) + 2 1 hi = T (hi ) + 2 2 n i −1 X j=1 n i X ! + f (b) ! f (a + jhi ) + f (b) f j=1 hi a+j 2 ni hi X 2j − 1 + f a+ hi 2 j=1 2 2j − 1 a+ hi . 2 Damit erhalten wir Algorithmus 6.2.2 (Extrapolationsverfahren mit Romberg-Folge) n = 1; h = b − a; k = 0; p(0) = 0.5 ∗ h(f (a) + f (b)); repeat k = k + 1; s = 0; h = 0.5 ∗ h; n = 2 ∗ n; q = 1; for i = 1 : 2 : (n − 1) s = s + f (a + i ∗ h); end p(k) = 0.5 ∗ p(k − 1) + s ∗ h; for i = (k − 1) : −1 : 0 q = 4 ∗ q; p(i) = p(i + 1) + (p(i + 1) − p(i))/(q − 1); end until Abbruch Integral = p(0); Als Abbruchkriterium verwendet man beispielsweise für ε = 10−6 |T0,k − T0,k−1 | ≤ ε. |T0,k | 111 Beispiel 6.2.3 (Romberg-Verfahren) Betrachte Z 1 f (x)dx = [− ln(cos(x))]10 = − ln(cos(1)) ≈ 0.6156264704. f (x) = tan(x), 0 Das Romberg-Verfahren liefert folgende Approximation 0.7787038623 0.6625031761 0.6237696140 0.6279861833 0.6164805191 0.6159945794 0.6187665645 0.6156933582 0.6156408808 0.6156352666 0.6164148747 0.6156309781 0.6156268194 0.6156265962 0.6156265622 und folgenden Fehler 1.63e-01 4.69e-02 8.14e-03 1.24e-03 8.54e-04 3.68e-04 3.14e-04 6.69e-05 1.44e-05 8.80e-06 7.88e-05 4.51e-06 3.49e-07 1.26e-07 9.19e-08 Eine genauere Betrachtung der Fehler in Spalte 1 zeigt, dass sich der Fehler bei jeder Schritthalbierung viertelt, was auf die Approximationsordnung 2 der iterierten Trapezregel schließen lässt. In der zweiten Spalte erkennt man (mit etwas gutem Willen), dass sich der Fehler bei Schrittweitenhalbierung etwa um den Faktor 1/16 verringert, was auf die Ordnung 4 hindeutet. 6.3 Gauß’sche Quadraturformeln Ziel ist es nun, Quadraturformeln der Form (6.2) mit maximalem Exaktheitsgrad herzuleiten, wobei neben den Gewichten wi auch die Stützstellen xi frei wählbar sind. Wir betrachten hier etwas allgemeinere Integrale der Form Z b Iα [f ] := α(x)f (x)dx a die durch die Quadraturformel QN [f ] = N X wk f (xk ) k=0 approximiert werden sollen, wobei die Gewichtsfunktion α : [a, b] → R+ genau dann zulässig ist, wenn sie stetig und positiv ist. Ist a = −∞ oder b = +∞, so verlangen Rb wir zusätzlich noch, dass die Momente a α(x)xk dx existieren und endlich sind. Wir haben bei den interpolatorischen Quadraturformeln (Satz 6.1.3) gesehen, dass bei fest 112 gewählten Knoten xk , k = 0, . . . , N die N + 1 Freiheitsgrade wk , k = 0, . . . , N ausreichen um Polynome vom Höchstgrad N exakt zu integrieren (evtl. ist aber auch mehr möglich). Da wir nun auch die Stützstellen frei wählen können, haben wir 2N +2 Freiheitsgrade und vermuten daher, dass Polynome mit Höchstgrad 2N + 1 exakt integriert werden können. Mehr geht nicht, wie das folgende Lemma beweist. Lemma 6.3.1 PN Zu jeder Quadraturformel QN [f ] = k=0 wk f (xk ) existiert ein Polynom p vom Grad 2N + 2 mit Iα [p] 6= QN [p]. Q 2 Beweis: Sei p(x) := N k=0 (x − xk ) . Dann ist p ein Polynom vom Grad 2N + 2 und es gilt p(x) ≥ 0 für alle x ∈ [a, b]. Es folgt b Z Iα [p] = α(x)p(x) dx > 0. a Andererseits ist QN [p] = N X wk p(xk ) = k=0 N X k=0 wk N Y (xk − xk )2 = 0, k=0 2 und damit gilt Iα [p] 6= QN [p]. Also beträgt der maximale Exaktheitsgrad von QN [f ] höchstens 2N + 1. Ziel ist es zu zeigen, dass er exakt 2N + 1 ist. Dann interessieren wir uns natürlich für die Wahl der Stützstellen. Wir gehen konstruktiv vor und bestimmen diese zuerst. Sei dazu für alle stetigen Funktionen f, g : [a, b] → R das Skalarprodukt Z hf, giα := b α(x)f (x)g(x) dx a definiert. Wir suchen nun Polynome qk , k = 0, 1, . . . vom Grad k mit Höchstkoeffizient 1, die paarweise orthogonal bezüglich dieses Skalarprodukts sind, d.h. die hqk , q` iα = 0 für k 6= ` erfüllen. Nach Definition ist q0 ≡ 1. Um die weiteren Polynome zu bestimmen, verwenden man das Gram-Schmidtschen Orthogonalisierungsverfahrens und die Polynome 1, x, x2 , . . . , die eine Basis des Raums aller Polynome bilden. Man verwendet also die Vorschrift k−1 X hqj , xk iα q k = xk − qj , k = 1, 2, . . . , hq , q i j j α j=0 und erhält dann als Ergebnis. Satz 6.3.2 (Orthogonale Polynome) 113 Die durch die Rekursion q0 (x) = 1, q1 (x) = x − a1 , qk (x) = (x − ak )qk−1 (x) − bk qk−2 (x), k = 2, 3, . . . , mit hxqk−1 , qk−1 iα , hqk−1 , qk−1 iα hxqk−1 , qk−2 iα = hqk−2 , qk−2 iα ak = bk definierten Polynome sind bezüglich des Skalarprodukts h·, ·iα orthogonal. Beweis: Der Beweis wird per vollständiger Induktion geführt: Aus der Rekursion ergibt sich induktiv, dass qk ein Polynom vom Grad k und somit nicht das Nullpolynom ist. Daher sind die Nenner in ak und bk ungleich Null. Wir zeigen induktiv, dass hqk , qi iα = 0 für alle i < k gilt. Für k = 0 ist nichts zu zeigen. Für k = 1 gilt mit der Definition von a1 hq1 , q0 iα = hx − a1 , q0 iα = hxq0 , q0 iα − a1 hq0 , q0 iα = 0. Sei die Behauptung gezeigt für ein k ≥ 1. Dann gilt mit den Definitionen für ak+1 und bk+1 : hqk+1 , qk iα = hxqk , qk iα − ak+1 hqk , qk iα − bk+1 hqk−1 , qk iα = 0, | {z } =0 hqk+1 , qk−1 iα = hxqk , qk−1 iα − ak+1 hqk , qk−1 iα −bk+1 hqk−1 , qk−1 iα = 0. | {z } =0 Für i = 1, . . . , k − 2 folgt hqk+1 , qi iα = hxqk , qi iα − ak+1 hqk , qi iα −bk+1 hqk−1 , qi iα | {z } | {z } =0 =0 = hqk , xqi iα = hqk , qi+1 + ai+1 qi + bi+1 qi−1 iα = 0, wobei die Rekursionsgleichung ausgenutzt wurde, um xqi darzustellen. Schließlich gilt noch hqk+1 , q0 iα = hxqk , q0 iα − ak+1 hqk , q0 iα − bk+1 hqk−1 , q0 iα = hqk , xq0 iα = hqk , q1 + a1 q0 iα = 0, 114 womit die Behauptung gezeigt ist. Für klassische Wahlen der Gewichtsfunktion ergeben sich die folgenden Polynome. [a, b] [−1, 1] [−1, 1] α(x) 1 √ 1 1−x2 [0, ∞) exp(−x) (−∞, ∞) exp(−x2 ) 2 Name erste Vertreter Legendre-Polynome 1, x, x2 − 13 , x3 − 53 x Tschebyscheff-Polynome 1, x, x2 − 12 , x3 − 43 x Laguerre-Polynome Hermite-Polynome 1, x − 1, x2 − 4x + 2, x3 − 9x2 + 18x − 6 1, x, x2 − 12 , x3 − 23 x Die Nullstellen des Polynoms qN +1 sollen nun die Stützstellen für QN [f ] werden. Dies macht nur Sinn, wenn diese Nullstellen einfach sind und im Intervall [a, b] liegen. Lemma 6.3.3 Das gemäß Satz 6.3.2 definierte orthogonale Polynom qN +1 hat lauter einfache Nullstellen in (a, b). Beweis: Seien a < x0 < x1 < . . . < xk < b die Nullstellen von qN +1 in (a, b), an denen Q qN +1 sein Vorzeichen wechselt. Angenommen es gilt k < N. Dann hat p(x) = ki=0 (x − xi ) den Grad k + 1 < N + 1. Es gilt also p ∈ span{q0 , q1 , . . . , qk+1 }, da die qi orthogonal und damit linear unabhängig sind. Da qN +1 orthogonal zu allen qi ist, folgt hqN +1 , piα = 0. Ferner ändert die Funktion α(x)qN +1 (x)p(x) auf [a, b] ihr Vorzeichen nicht, da bei jedem Vorzeichenwechsel von qN +1 (x) nach Definition auch einer von p(x) auftritt und α(x) > 0 für alle x ∈ [a, b] gilt. Damit folgt aber Z b α(x)qN +1 (x)p(x) dx 6= 0, hqN +1 , piα = a also ein Widerspruch. Die Annahme war falsch und es gilt k = N. 2 Wir definieren unsere Stützstellen x0 , x1 , . . . , xN also als die Nullstellen von qN +1 , und PN bestimmen nun die Gewichte wk , k = 0, 1, . . . , N in QN [f ] := k=0 wk f (xk ) so, dass zumindest alle Polynome vom Grad N exakt integriert werden können. Betrachten wir hierzu die Lagrange-Polynome Lk (x) = N Y x − xj . x k − xj j=0,j6=k Unsere Forderung besagt dann, dass für alle k = 0, 1, . . . , N Z b ! α(x)Lk (x) dx = Iα [Lk ] = QN [Lk ] = a N X j=0 115 wj Lj (xk ) = wk | {z } =δkj gelten muss. Damit sind auch die Gewichte eindeutig bestimmt und wir haben die sogenannte Gauß-Quadraturformel hergeleitet. Es gilt nun Satz 6.3.4 (Gauß-Quadratur) Seien α eine zulässige Gewichtungsfunktion und q ein nichttriviales Polynom vom Grad N + 1 mit Zb xk α(x)q(x)dx = 0 (0 ≤ k ≤ N ), (6.7) a d.h. q ist α-orthogonal zu allen Polynomen vom Höchstgrad N . Seien x0 , x1 , . . . , xN die Nullstellen von q. Dann ist die Gauß-Quadraturformel QN [f ] = N X Zb wk f (xk ) mit α(x)Lk (x)dx wk = k=0 (6.8) a exakt für alle Polynome vom Höchstgrad 2N + 1. Beweis: Per Konstruktion ist QN [f ] = Iα [f ] für alle Polynome f vom Höchstgrad N. Sei nun f ein Polynom vom Höchstgrad 2N + 1. Division mit Rest von f durch q liefert f = qp+r mit Polynomen p und r vom Höchstgrad N . Da alle xk , k = 0, . . . , N nach Definition Nullstellen von q sind, gilt f (xk ) = r(xk ). Nun liefert die folgende Gleichungskette die Behauptung: QN [f ] = N X wk f (xk ) = k=0 = N X N X k=0 wk [q(xk ) p(xk ) + r(xk )] | {z } =0 wk r(xk ) = QN [r] k=0 Z b = Iα [r] = α(x)r(x) dx a Z b Z b = α(x)q(x)p(x) dx + α(x)r(x) dx a |a {z } =0 Z b = α(x)f (x) dx = Iα [f ]. a 2 Tabelle 6.3 fasst Knoten und Gewichte von einigen Gauß-Quadraturformeln zusammen. Zu beachten ist, dass die angegebenen Stützstellen und Koeffizienten ausschließlich für das Intervall [−1, 1] und die Gewichtsfunktion α ≡ 1 gültig sind. 116 Tabelle 6.3: Einige Gauß -Quadraturformeln für α ≡ 1 und [a, b] = [−1, 1]. N Knoten xi Gewichte Ai 0 1 x0 x0 x1 x0 x1 x2 w0 w0 w1 w0 w1 w2 2 3 x0 x1 x2 x3 =0 p = − 1/3 p = + 1/3 p = − 3/5 =0 p = + 3/5 r q = − 37 + 71 24 5 r q = − 37 − 71 24 5 r q = + 37 − 17 24 5 r q = + 37 + 17 24 5 =2 =1 =1 = 5/9 = 8/9 = 5/9 w0 = √ 18− 30 36 w1 = √ 18+ 30 36 w2 = √ 18+ 30 36 w3 = √ 18− 30 36 Beispiel 6.3.5 Approxiere das Integral Z 2 I= exp(−t2 )dt. −2 Da das Intervall auf [a, b] := [−2, 2] und nicht auf [−1, 1] definiert ist, muss es zunächst auf [−1, 1] transformiert werden. Dies wird erreicht durch die lineare Transformation t(x) = a(1 − x) + b(1 + x) −2(1 − x) + 2(1 + x) = = 2x 2 2 ⇒ dt = 2dx. Transformation des Integrals: Z 2 Z 2 1 Z 1 2 exp(−4x2 )dx. exp(−(2x) ) · 2dx = exp(−t )dt = −2 2 −1 −1 Anwendung der Formeln: Z 2 2 exp(−t )dt = −2 Z 1 2 exp(−4x2 ) dx ≈ {z } −1 | =f (x) 117 N X k=0 wk f (xk ). Für N = 1 folgt Z 2 2 Z exp(−t )dt = −2 1 2 exp(−4x2 ) dx {z } −1 | =f (x) ≈ w0 f (x0 ) + w1 f (x1 ) p p = 1 · f (− 1/3) + 1 · f ( 1/3) = 2 exp(−4/3) + 2 exp(−4/3) = 4 exp(−4/3) ≈ 1.054388553. Für N = 2 erhalten wir den Wert 1.979373230 und für N = 3 1.714546068. Der exakte Wert des Integrals beträgt ungefähr I ≈ 1.764162782. 6.4 Adaptive Verfahren am Beispiel der Simpson-Regel In den vorangehenden Abschnitten wurde stets eine feste Unterteilung des Intervalls [a, b] gewählt. In diesem Abschnitt wird diese Unterteilung nicht a priori festgelegt, sondern adaptiv in Abhängigkeit vom Approximationsfehler gewählt. Der Grund hierfür ist, dass in einigen Abschnitten von [a, b] mitunter (z.B. wenn f dort stark variiert) eine sehr feine Unterteilung notwendig ist, um eine gute Approximation zu erreichen, während in anderen Bereichen (z.B. weil f sich dort kaum ändert) eine grobe Unterteilung ausreicht, vgl. Abbildung 6.6. Abbildung 6.6: Idee eines adaptiven Verfahrens: Die Schrittweite wird lokal adaptiv an die Eigenschaften von f angepasst. Zur Demonstration wählen wir die 118 elementare Simpson-Regel: Z a b b−a a+b f (x)dx ≈ S(a, b) := f (a) + 4f + f (b) . 6 2 Falls f in [a, b] 4-mal stetig differenzierbar ist, dann besitzt die Simpson-Regel folgende Fehlerdarstellung: Fehler: 1 I[f ] − S(a, b) = − 90 b−a 2 5 f (4) (ξ) (ξ ∈ (a, b)). Die adaptive Simpson-Regel basiert auf einer Methode zur numerischen Schätzung des Interpolationsfehlers: Adaptive Simpson-Regel: (i) Seien das Intervall [a, b] und eine Toleranz ε > 0 gegeben. (ii) Wende die Simpson-Regel auf [a, b] an und berechne S (1) := S(a, b). (iii) Wende die elementare Simpsonregel auf [a, c] und [c, b] mit c = (a + b)/2 an und berechne S (2) := S(a, c) + S(c, b) (zusammengesetzte Simpson-Regel). (iv) Falls die Ungleichung |S (2) − S (1) | < 15ε (6.9) erfüllt ist, STOP (der Approximationsfehler ist hinreichend klein) mit Z b f (x)dx ≈ S (2) + a 1 S (2) − S (1) . 15 (6.10) Falls die Ungleichung nicht erfüllt ist: – Teile [a, b] in zwei Teilintervalle [a, c] und [c, b] auf, – Wende diesen Algorithmus rekursiv auf [a, c] und [c, b] mit ε ersetzt durch ε/2 an. Hintergrund der Fehlerschätzung (6.9): Seien h = b − a und C := f (4) (c)1 . Mit Hilfe der Fehlerformel der elementaren Simpson1 f (4) Falls f 5-mal stetig differenzierbar ist, liefert Taylorentwicklung f (4) (ξ) = f (4) (c) + f (5) (ζ)(ξ − c) = (c) + O(h) für jedes ξ ∈ (a, b) und ζ ∈ (c, ξ). 119 Regel folgt E (1) := I[f ] − S (1) E (2) := I[f ] − S (2) 5 h C + O(h6 ) 2 5 5 1 h/2 1 h/2 = − C− C + O(h6 ) 90 2 90 2 5 ! 1 h 1 − C + O(h6 )). = 16 90 2 1 = − 90 (6.11) (6.12) Ein Vergleich der rechten Seiten liefert I[f ] − S (1) ≈ 16E (2) , I[f ] − S (2) = E (2) . Subtraktion der zweiten von der ersten Gleichung liefert S (2) − S (1) ≈ 15E (2) , und daher I[f ] − S (2) = E (2) ≈ Deshalb kann der Wert 1 S (2) − S (1) . 15 1 S (2) − S (1) 15 (6.13) als Approximation des Fehlers I[f ] − S (2) verwendet werden. Beachte, dass E (2) nicht bekannt ist, aber wir können (6.13) berechnen. Gleichung (6.9) sichert (approximativ), dass E (2) kleiner als ε ist. Die Approximation in (6.10) hat den Fehler O(h6 ). Dies folgt aus den Fehlerdarstellungen (6.11) und (6.12) (mulitpliziere die erste Gleichung mit −1/15, die zweite mit 16/15 und addiere beide Gleichungen). Die Division von ε durch zwei in Schritt (iv) soll garantieren, dass der Fehler in jedem der zwei Teilintervalle kleiner als ε/2 verbleibt, so dass der Fehler des gesamten Intervalls kleiner als ε/2 + ε/2 = ε ist. 6.5 Mehrdimensionale Integrale Ziel dieses Abschnitts ist es, n-dimensionale Integrale der Form Z b1 Z bn I[f ] = ··· f (x1 , . . . , xn )dxn · · · dx1 a1 an zu approximieren. Zur Vereinfachung sei im folgenden n = 2 gewählt, d.h. Z bZ d I[f ] = f (x, y)dxdy. a c 120 (6.14) Man kann leicht Quadraturformeln entwickeln, indem die Quadraturformeln für eindimensionale Integrale verwendet werden. Dies soll am Beispiel der Trapezregel verdeutlicht werden. Zunächst wird das innere Integral Z d f (x, y)dx c für festes y durch die zusammengesetzte Trapezregel approximiert: Zd f (x, y)dx ≈ h1 c n1 −1 f (x0 , y) X f (xn1 , y) + f (xi , y) + 2 2 i=1 ! =: n1 X Ai f (xi , y) i=0 wobei h1 = (d − c)/n1 und xi = c + ih1 , i = 0, 1, . . . , n1 . Dann wird diese Approximation in das ursprüngliche Integral eingesetzt und die iterierte Trapezregel wird angewendet auf das äußere Integral mit h2 = (b − a)/n2 , yj = a + jh2 , j = 0, 1, . . . , n2 : Zb Zd f (x, y)dxdy ≈ a c Zb X n1 a = ≈ Ai f (xi , y)dy i=0 n1 X Ai f (xi , y)dy i=0 a n1 X n2 X Ai i=0 = Zb n1 X n2 X ! Bj f (xi , yj ) j=0 Ai Bj f (xi , yj ). i=0 j=0 Entsprechend kann für n > 2 verfahren werden. Ebenso können auch Newton-CotesFormeln oder zusammengesetzte Verfahren verwendet werden. Es gibt auch Verfahren, die auf der Gauß-Quadratur basieren. 121 Kapitel 7 Eigenwert-Probleme Definition 7.0.1 (Eigenwert) Sei A ∈ Rn×n eine gegebene Matrix. Eine (eventuell komplexe) Zahl λ heißt Eigenwert von A, wenn es einen (eventuell komplexen) Vektor x 6= 0 gibt mit Ax = λx. Jedes solche x 6= 0 heißt Eigenvektor von A zum Eigenwert λ. Aus der linearen Algebra sind folgende Charakterisierungen eines Eigenwertes bekannt: λ ist ein Eigenwert von A ⇔ ∃x 6= 0 : Ax = λx ⇔ ∃x 6= 0 : (A − λI)x = 0 ⇔ (A − λI) ist singulär ⇔ det(A − λI) = 0 ⇔ λ ist Nullstelle des charakteristischen Polynoms pA (t) := det(tI − A). Aus der letzten Charakterisierung folgt, dass ähnliche Matrizen dieselben Eigenwerte besitzen: Sind nämlich A, B ∈ Rn×n ähnlich, also B = X −1 AX für eine reguläre Matrix X ∈ Rn×n , so folgt aus bekannten Eigenschaften der Determinante pB (t) = det(tI − B) = det(X −1 (tI − A)X) = det(X −1 ) det(tI − A) det(X) 1 det(tI − A) det(X) = det(X) = pA (t), d.h. die charakteristischen Polynome von A und B stimmen überein. Damit sind auch die Eigenwerte der Matrizen gleich. Nach dem Fundamentalsatz der Algebra hat jedes Polynom vom Grad n, genau n komplexe Nullstellen. Damit hat jede Matrix in Rn×n genau n komplexe Eigenwerte. Es sei jedoch bemerkt, dass eine reelle Matrix keine reellen Eigenwerte haben muss. 122 Beispiel 7.0.2 ! 0 1 Die Matrix A = ∈ R2×2 hat das charakteristische Polynom pA (t) = t2 + 1 und −1 0 damit die beiden rein imaginären Eigenwerte +i und −i. Betrachten wir eine Anwendung eines Eigenwertproblems. Beispiel 7.0.3 Zwei Armeen A und B liefern sich ein Gefecht. Es bezeichne A(t) bzw. B(t) die Stärke (etwa die Anzahl der Soldaten) der Armeen A bzw. B zum Zeitpunkt t ≥ 0. In t = 0 mögen die Anfangsstärken A0 bzw. B0 sein. Im Laufe des Gefechts dezimieren sich die Armeen gegenseitig proportional zur Stärke der jeweils verfeindeten Armee, so dass wir das Differentialgleichungssystem Ȧ(t) = −βB(t), A(0) = A0 , Ḃ(t) = −αA(t), B(0) = B0 erhalten, wobei α > 0 bzw. β > 0 ein (von der Ausrüstung der Armee abhängiges) Maß der Zerstörungskraft der Armee A bzw. B bezeichne. Wir erhalten das System ! ! ! ! ! Ȧ(t) 0 −β A(t) A(0) A0 = , = . Ḃ(t) −α 0 B(t) B(0) B0 Zur Lösung dieses Problems benötigt man bekanntlich die Eigenwerte der Matrix 0 −β −α 0 die sich zu λ1 = p αβ, p λ2 = − αβ berechnen. Damit erhält man die Lösung p p √ 1 √ A(t) = √ [ αA0 − βB0 ]eλ1 t + [ αA0 + βB0 ]eλ2 t , 2 α p p √ 1 √ λ1 t λ2 t B(t) = √ −[ αA0 − βB0 ]e + [ αA0 + βB0 ]e 2 β für alle t ≥ 0. Wegen λ2 < 0 verschwindet sowohl für A(t) als auch für B(t) der zweite Term für t → ∞. Die Stärken beider Armeen A und B hängen mit zunehmender Zeit √ √ daher fast nur vom ersten Term ab. Gilt αA0 − βB0 > 0, so bleibt A(t) positiv, √ √ während B(t) nach endlicher Zeit ausgelöscht ist. Im Fall αA0 − βB0 = 0 tritt eine Pattsituation ein, es werden zwar beide Armeen dezimiert, die Armeestärken bleiben aber √ √ positiv. Schließlich wird A(t) nach endlicher Zeit ausgelöscht, falls αA0 − βB0 < 0 ist, während B(t) dann stets positiv bleibt. In diesem Beispiel konnten wir die Eigenwerte der Matrix explizit bestimmen. Im Allgemeinen ist dies jedoch nicht möglich. Mit Hilfe der Algebra kann man beweisen, dass es 123 ! , keine explizite Lösungsformel für die Berechnung von Nullstellen von allgemeinen Polynomen vom Grad größer gleich 5 gibt. Daher kann es auch keinen endlichen Algorithmus zur Berechnung der Eigenwerte einer beliebigen Matrix geben, denn andernfalls könnte man diesen auf die Matrix 0 . . . 0 −a0 . .. 1 . . ... . A= ∈ Rn×n ... 0 −an−2 0 1 −an−1 anwenden, und somit die Nullstellen von deren charakteristischen Polynom det(tI − A), (immer nach letzter Zeile entwickeln) pA (t) = det(tI − A) = tn + an−1 tn−1 + . . . + a1 t + a0 bestimmen. Somit hätte man eine Algorithmus zur Bestimmung der Nullstellen eines beliebigen Polynoms vom Grad n mit führendem Koeffizienten 1. Daher sind alle Verfahren zur Eigenwertberechnung iterative Verfahren, die sukzessive, bessere Näherungen der Eigenwerte berechnen. Da die Eigenwerte einer Matrix die Nullstellen des charakteristischen Polynoms sind, hat man mit jedem Verfahren zur Lösung einer nichtlinearen Gleichung auch ein Verfahren zur Eigenwertberechnung. Im Folgenden werden wir aber effizientere Verfahren konstruieren. 7.1 Matrixzerlegungen Bei manchen Matrizen kann man die Eigenwerte sehr leicht ablesen. So stehen bei Diagonalmatrizen oder auch bei oberen oder unteren Dreiecksmatrizen die Eigenwerte auf der Diagonalen. Wir haben bereits gesehen, dass ähnliche Matrizen dieselben Eigenwerte haben. Daher ist es naheliegend zu versuchen durch Ähnlichkeitstransformationen eine Matrix auf Diagonalgestalt zu bringen. Definition 7.1.1 (Diagonalisierbarkeit) Eine Matrix A ∈ Rn×n heißt diagonalisierbar, wenn es eine (eventuell komplexe) reguläre Matrix X ∈ Rn×n und eine (eventuell ebenfalls komlexe) Diagonalmatrix D ∈ Rn×n gibt, mit A = X D X −1 . Bezeichnen wir die Spalten von X mit xi und die Diagonalelemente von D mit λi , so folgt für eine diagonalisierbare Matrix A AX = XD ⇔ Axi = λi xi 124 ∀i = 1, . . . , n. Wir erkennen, dass die Diagonalelemente gerade die Eigenwerte von A und die Spalten von X die dazugehörigen Eigenvektoren sind. Da X regulär ist, müssen die Eigenvektoren also linear unabhängig sein. Damit sieht man leicht, dass nicht jede Matrix diagonalisierbar ! 0 1 ist. Die Matrix A = ist nicht diagonalisierbar, da λ = 0 zwar doppelter Eigenwert 0 0 ist, der zugehörige Eigenraum aber eindimensional ist. Es bleibt also die Frage, ob wir zumindest auf eine Dreiechsgestalt mittels Ähnlichkeitstransformationen kommen können. Satz 7.1.2 (Schur-Zerlegung) Sei A ∈ Rn×n eine Matrix mit lauter reellen Eigenwerten. Dann existiert eine orthogonale Matrix Q ∈ Rn×n derart, dass T = Q> AQ eine obere Dreiecksmatrix ist. Beweis: Wir beweisen den Satz mittels vollständiger Induktion. Für n = 1 ist die Behauptung mit Q = I richtig. Sei daher n ≥ 2 und die Aussage für alle Matrizen mit lauter reellen Eigenwerten der Dimension n−1 richtig. Sei A ∈ Rn×n eine gegebene Matrix mit lauter reellen Eigenwerten, und sei λ1 ∈ R einer dieser Eigenwerte mit zugehörigem Eigenvektor x1 ∈ Rn . Ohne Einschränkung können wir diesen normieren, so dass kx1 k2 = 1 gilt. Wir erweitern x1 mit den Vektoren x2 , . . . , xn zu einer Orthonormalbasis des Rn . Ferner seien λ2 , . . . , λn die übrigen Eigenwerte von A. Definiere die Matrix V ∈ Rn×(n−1) durch | | V = x2 . . . xn . | | Dann ist die Matrix (x1 |V ) ∈ Rn×n eine orthogonale Matrix mit A(x1 |V ) = (Ax1 |AV ) = (λ1 x1 |AV ). > Wegen x> 1 x1 = 1 und V x1 = 0 folgt (x1 |V )> A(x1 |V ) = λ1 x> 1 AV 0 V > AV ! , und die Eigenwerte von V > AV sind gerade λ2 , . . . , λn , da (x1 |V )> A(x1 |V ) dieselben Eigenwerte wie A besitzt (orthogonale Ähnlichkeitstransformation). Nach Induktionsvoraussetzung existiert daher eine orthogonale Matrix U ∈ R(n−1)×(n−1) derart, dass R = U > V > AV U eine obere Dreiecksmatrix ist. Definiert man ! 1 0 Q := (x1 V ) ∈ Rn×n , 0 U 125 so ist T := Q> AQ = = λ1 0 = λ1 0 ! λ 1 x> 1 AV 0 V > AV ! x> AV U 1 > > U V AV U ! x> AV U 1 R 1 0 0 U> ! 1 0 0 U ! eine obere Dreiecksmatrix und der Induktionsschluss bewiesen. 2 Eine wichtige Konsequenz aus diesem Satz ist der sogenannte Spektralsatz. Satz 7.1.3 (Spektralsatz) Sei A ∈ Rn×n eine symmetrische Matrix. Dann existiert eine orthogonale Matrix Q ∈ Rn×n derart, dass D = Q> AQ eine Diagonalmatrix ist. Beweis: Eine reelle symmetrische Matrix (A> = Ā) hat nur reelle Eigenwerte, denn für einen Eigenwert λ ∈ C mit zugehörigem normierten Eigenvektor x ∈ Cn gilt: λ = λx> x̄ = (Ax)> x̄ = x> A> x̄ = x> Āx̄ = x> (λ̄x̄) = λ̄, also λ = λ̄ und damit λ ∈ R. Damit lässt sich Satz 7.1.2 anwenden, und es existiert eine orthogonale Matrix Q ∈ Rn×n derart, dass T = Q> AQ eine obere Dreiecksmatrix ist. Aus der Symmetrie von A folgt T > = (Q> AQ)> = Q> A> Q = Q> AQ = T, 2 also muss T eine Diagonalmatrix sein. Bezeichnen wir die Spalten der Matrix Q mit v1 , . . . , vn und die Diagonalelemente von D mit λ1 , . . . , λn , so folgt aus dem Spektralsatz A = QDQ> und daraus die Spektraldarstellung n X A= λi vi vi> . i=1 Aus dem Spektralsatz ergibt sich das folgende Korollar. Korollar 7.1.4 Sei A ∈ Rn×n eine symmetrische Matrix mit größtem Eigenwert λmax und kleinstem Eigenwert λmin . Dann gilt für alle x ∈ Rn : λmin x> x ≤ x> Ax ≤ λmax x> x. 126 Beweis: Nach dem Spektralsatz hat man D = Q> AQ mit einer orthogonalen Matrix Q, und die Eigenwerte von A stehen auf der Diagonalen von D. Es folgt mit y := Q> x x> Ax = x> QDQ> x = y > Dy = n X λi yi2 . i=1 Damit folgt die Behauptung aus λmin x> x = λmin x> Q> Qx = λmin y > y ≤ n X λi yi2 ≤ λmax y > y = λmax x> x. i=1 2 Hieraus erhält man eine Abschätzung für den größten und den kleinsten Eigenwert einer Matrix. Da man jeweils einen entsprechenden Eigenvektor wählen kann, gilt Gleichheit λmin = min x6=0 x> Ax x> x und λmax = max x6=0 x> Ax x> x Um die Lage der Eigenwerte abschätzen zu können eignet sich der folgende Satz. Satz 7.1.5 (Satz von Gerschgorin) Seien A ∈ Rn×n und λ ∈ C ein beliebiger Eigenwert von A. Dann ist λ∈ n [ Ki , i=1 mit den Gerschgorin-Kreisen Ki := {z ∈ C | |z − aii | ≤ n X |aij |} (i = 1, . . . , n) j=1,j6=i Beweis: Nach Voraussetzung existiert ein Eigenvektor x 6= 0 mit Ax = λx. Sei i0 ∈ {1, . . . , n} ein Index mit |xi0 | = max1≤i≤n |xi |. Insbesondere ist dann |xi0 | = 6 0. Es folgt λxi0 = [Ax]i0 = n X ai 0 j x j , j=1 und somit nach Division durch xi0 n X xj |λ − ai0 i0 | = ai 0 j ≤ xi 0 j=1,j6=i0 127 n X j=1,j6=i0 |ai0 j |. Also ist λ ∈ Ki0 ⊆ Sn i=1 2 Ki . Beispiel 7.1.6 1 0.1 −0.1 Betrachten wir die Matrix A = 0 2 0.4 . Wir erhalten die Gerschgorin-Kreise −0.2 0 3 K1 = {z ∈ C | |z − 1| ≤ 0.2}, K2 = {z ∈ C | |z − 2| ≤ 0.4}, K3 = {z ∈ C | |z − 3| ≤ 0.2}. Die drei Eigenwerte von A liegen also in dem durch die Kreise gekennzeichneten Gebiet, jedoch folgt aus dem Satz nicht, dass in jedem der Kreise genau ein Eigenwert liegen muss. Man kann den Satz auch noch auf die transponierte Matrix A> anwenden, und erhält die Gerschgorin-Kreise K1> = {z ∈ C | |z − 1| ≤ 0.2}, K2> = {z ∈ C | |z − 2| ≤ 0.1}, K3> = {z ∈ C | |z − 3| ≤ 0.5}. Man kann also den Kreis um 2 verkleinern, d.h. eine bessere Abschätzung erhalten. 7.2 Potenz-Methoden Sei A ∈ Rn×n symmetrisch und v ∈ Rn eine Näherung für einen Eigenvektor von A. Dann liefert die Lösung des linearen Ausgleichsproblems min kAv − λvk2 λ eine zugehörige Näherung für den Eigenwert. Lemma 7.2.1 (Rayleigh-Quotient) Sei A ∈ Rn×n symmetrisch und v ∈ Rn \ {0}. Dann ist der Rayleigh-Quotient λ := λ(v) = 128 v > Av v>v (7.1) die Lösung von (7.1). Beweis: (7.1) ist äquivalent zu der Gauß’schen Normalengleichung λv > v = v > Av. > Also löst λ = vv>Av nach Satz 3.1.3 das lineare Ausgleichsproblem (7.1). 2 v Wie erhält man eine Näherung an einen Eigenvektor von A? Sei A= QDQ> dienach dem Spektralsatz existierende Zerlegung mit orthogonaler Matrix | | Q = q1 . . . qn und Diagonalmatrix D = diag(λ1 , . . . , λn ). O.B.d.A. sei |λ1 | ≥ |λ2 | ≥ | | . . . ≥ |λn |. Ferner sei der Eigenwert λ1 dominant, d.h. |λ1 | > |λ2 | ≥ . . . ≥ |λn |. P Sei v (0) ∈ Rn \ {0}. Dann gibt es Koeffizienten αi ∈ R mit v (0) = ni=1 αi qi . Dabei habe v (0) einen nichttrivialen Anteil in Richtung q1 , d.h. α1 6= 0. Setzt man v (k+1) := Av (k) für alle k = 0, 1, . . . , so folgt induktiv v (k) = Ak v (0) und weiter für alle k = 0, 1, . . . v (k) = Ak v (0) = n X αi Ak qi i=1 = n X αi (λi )k qi i=1 k n X αi λi qi = α1 (λ1 ) α1 λ1 i=1 k k !> α λ α λn 2 2 n = α1 (λ1 )k Q 1, ,..., α 1 λ1 α 1 λ1 {z } | k =:x(k) k (k) = α1 (λ1 ) Qx . Dabei gilt nach Voraussetzung x(k) → e1 := (1, 0, . . . , 0)> und somit {Qx(k) } → q1 . Also konvergiert v (k) gegen ein Vielfaches von q1 . Der Name Potenz-Methode erklärt sich, da man Potenzen der Matrix A berechnen muss. Durch Normierung der Länge des Vektors v (k) erhält man den folgenden Algorithmus: Algorithmus 7.2.2 (Potenz-Methode, Vektor-Iteration) (i) Wähle v (0) ∈ Rn mit kv (0) k2 = 1, und setze k := 0. (ii) Setze w(k+1) := Av (k) . (iii) Setze v (k+1) := w(k+1) . kw(k+1) k2 (iv) Setze λ(k+1) := (v (k+1) )> Av (k+1) . 129 (iv) Setze k := k + 1 und gehe zu (ii). Als Abbruchkriterium kann man bei der numerischen Umsetzung beispielsweise kAv (k) − λ(k) v (k) k2 ≤ ε (7.2) wählen. Es gilt der folgende Konvergenzsatz Satz 7.2.3 Sei A ∈ Rn×n symmetrisch mit dominantem Eigenwert λ1 . Der Startvektor v (0) habe einen nichttrivialen Anteil in Richtung q1 . Dann gilt: (a) Jeder Häufungspunkt einer durch die Potenz-Methode erzeugten Folge {v (k) } liegt in span{q1 }. (b) Die durch die Potenz-Methode erzeugte Folge {λ(k) } konvergiert gegen λ1 Beweis: Wegen kv (k) k2 = 1 ist die Menge {v (k) } beschränkt und hat somit einen Häufungspunkt. Nach der Herleitung gilt {v (k) } = {Ak v (0) } ∈ span{q1 }, also (a). Wegen kv (k) k2 = 1 gilt für jeden Häufungspunkt v̂ von {v (k) } kv̂k2 = 1 = kq1 k, also v̂ = ±q1 . Damit folgt λ(k) = (v (k) )> Av (k) → q1> Aq1 = λ1 q1> q1 = λ1 , 2 und somit (b). Bemerkung 7.2.4 • In der Herleitung der Potenz-Methode erkennt man, dass man auch einen mehrfachen betragsgrößten Eigenwert zulassen kann, d.h. λ1 = λ2 = . . . = λr für r ≤ n. • Was allerdings nicht funktioniert ist, wenn der betragsgrößte Eigenwert einmal po k sitiv und einmal negativ auftritt, d.h. λ1 = −λ2 da dann die Quotienten λλ21 = (−1)k oszillieren, und somit keine Konvergenz von v (k) erzielt werden kann. • Die Matrix A bleibt bei der Potenz-Methode unverändert. • Man benötigt nur eine Matrix-Vektor-Multiplikation pro Iteration. • Die Normierung verhindert einen overflow. 130 • Die Potenz-Methode konvergiert gegen den betragsgrößten Eigenwert und jeder Häufungspunkt der Folge {v (k) } ist ein zugehöriger Eigenvektor. • Die Konvergenzgeschwindigkeit hängt von der Größe |λ1 | deutlich größer als die übrigen Eigenwerte ist. |λ2 | |λ1 | ab, sie ist also gut, wenn Beispiel 7.2.5 ! −2 0 Wir wollen am Beispiel der Matrix A = zeigen, dass man die Konvergenz 0 1 des Eigenvektors nicht erzwingen kann. Der Eigenwert λ1 = −2 ist dominant und v (0) = (1, 0)> ist ein Eigenvektor. Wendet man die Potenzmethode auf diesen Vektor an, so erhält man v (k) = (−1)k v (0) . Man erhält also zwei Häufungspunkte, aber keinen Grenzwert. Dahingehend lässt sich Satz 7.2.3 also nicht verschärfen. Nun wollen wir noch ein numerisches Beispiel betrachten. Beispiel 7.2.6 (Potenz-Methode) 5 4 3 4 6 0 Wir wollen die Eigenwerte der Matrix A = 3 0 7 2 4 6 1 3 5 1 > Methode mit dem Startvektor √5 (1, 1, 1, 1, 1) und dem 10−6 ergibt den in der Tabelle dargestellten Verlauf. k kAv (k) − λ(k) v (k) k2 1 1.42e+00 2 4.61e-01 3 1.52e-01 4 5.07e-02 5 1.69e-02 6 5.67e-03 7 1.90e-03 8 6.37e-04 9 2.14e-04 10 7.17e-05 11 2.40e-05 12 8.06e-06 13 2.70e-06 14 9.06e-07 131 2 1 4 3 6 5 berechnen. Die Potenz 8 7 7 9 Abbruchkriterium (7.2) mit ε = λ(k) 22.2732784755 22.3927301087 22.4053275309 22.4067033079 22.4068560629 22.4068731480 22.4068750649 22.4068752803 22.4068753045 22.4068753072 22.4068753075 22.4068753075 22.4068753075 22.4068753075 Man erkennt lineare Konvergenz gegen den Eigenwert λ1 = 22.4068753075. Als zugehörigen Eigenvektor erhalten wir v ≈ (0.2458779720, 0.3023960827, 0.4532145027, 0.5771771472, 0.5563845677)> . Die Potenzmethode liefert nur den betragsgrößten Eigenwert von A. Wendet man den Algorithmus auf die Inverse A−1 an, so erhält man den betragsgrößten Eigenwert von A−1 , der gleich 1 durch den betragskleinsten Eigenwert von A ist. Wählt man einen Shift-Parameter µ, der dichter an λi liegt als an allen anderen Eigenwerten, so ist λi1−µ der betragsgrößte Eigenwert der Matrix (A − µI)−1 , und man kann wieder die PotenzMethode anwenden um den Eigenwert λi von A zu berechnen. Dies führt auf Algorithmus 7.2.7 (Inverse Iteration von Wieland) (i) Wähle v (0) ∈ Rn mit kv (0) k2 = 1, sowie µ ∈ R und setze k := 0. (ii) Bestimme w(k+1) als Lösung des linearen Gleichungssystems (A − µI)w = v (k) . (iii) Setze v (k+1) := w(k+1) . kw(k+1) k2 (iv) Setze λ(k+1) := (v (k+1) )> Av (k+1) . (iv) Setze k := k + 1 und gehe zu (ii). Bemerkung 7.2.8 • Schritt (ii) ist nur durchführbar, wenn µ noch kein Eigenwert von A ist. • Das Verfahren ist teurer als die Potenz-Methode, da in jedem Schritt ein lineares Gleichungssytem gelöst werden muss. Die Matrix (A − µI) wird jedoch nicht verändert. Daher genügt es, diese einmal zu Faktorisieren, so dass das Gleichungssystem relativ günstig zu lösen ist. • Bei vorhandener Näherung µ an einen Eigenwert von A eignet sich die Inverse Iteration besonders zur Berechnung eines zugehörigen Eigenvektors. Beispiel 7.2.9 (Inverse Iteration von Wielandt) Wir wenden die Inverse Iteration von Wieland auf die Matrix aus dem Beispiel 7.2.6 mit dem Startvektor v (0) = √15 (1, 1, 1, 1, 1)> , der Eigenwertnäherung µ = 5 und dem Abbruchkriterium (7.2) mit ε = 10−6 an. Wir erhalten den in der Tabelle dargestellte Iterationsverlauf. 132 k kAv (k) − λ(k) v (k) k2 1 1.19e+00 2 2.22e-02 3 1.18e-03 4 7.07e-05 5 4.24e-06 6 2.55e-07 Wir haben also einen anderen Eigenwert als in erhalten wir λ(k) 4.959971665915 4.849095349254 4.848950631563 4.848950122172 4.848950120323 4.848950120316 Beispiel 7.2.6 berechnet. Als Eigenvektor v ≈ (0.5471727958, −0.3125699200, 0.6181120763, −0.1156065936, −0.4554937467)> . Wählt man µ variabel und setzt hierbei insbesondere µ(k) gleich dem Rayleigh-Quotienten, so erhält man: Algorithmus 7.2.10 (Rayleigh-Quotienten-Iteration) (i) Wähle v (0) ∈ Rn mit kv (0) k2 = 1, und setze k := 0. (ii) Setze µ(k) := (v (k) )> Av (k) . (iii) Bestimme w(k+1) als Lösung des linearen Gleichungssystems (A − µ(k) I)w = v (k) . (iv) Setze v (k+1) := w(k+1) . kw(k+1) k2 (v) Setze λ(k+1) := (v (k+1) )> Av (k+1) . (vi) Setze k := k + 1 und gehe zu (ii). Der Aufwand der Rayleigh-Quotienten-Iteration ist deutlich größer als bei der Inversen Iteration von Wieland, da sich die Matrix im linearen Gleichungssystem in jedem Schritt ändert. Dafür lässt sich eine sehr schnelle Konvergenz des Verfahrens beweisen. Wir wollen nur das Resultat angeben und auf den Beweis verzichten. Satz 7.2.11 Sei A ∈ Rn×n symmetrisch. Die durch die Rayleigh-Quotienten-Iteration erzeugte Folge {v (k) } konvergiere gegen einen Eigenvektor von A. Dann konvergiert die Folge {µ(k) } lokal kubisch gegen den zugehörigen Eigenwert, d.h. es existiert ein c > 0 mit |µ(k+1) − λ| ≤ c|µ(k) − λ|3 133 für alle k ∈ N hinreichend groß. Beispiel 7.2.12 (Rayleigh-Quotienten-Iteration) Wenden wir die Rayleigh-Quotienten-Iteration auf die Matrix aus dem Beispiel 7.2.6 mit dem Startvektor v (0) = √15 (1, 1, 1, 1, 1)> , und dem Abbruchkriterium (7.2) mit ε = 10−6 an, so beobachten wir sehr schnelle Konvergenz in nur drei Iterationen. k kAv (k) − λ(k) v (k) k2 1 4.81e-01 2 4.54e-04 3 3.96e-13 Als Eigenvektor erhalten wir λ(k) 22.392146487612 22.406875294193 22.406875307580 v ≈ (0.2458779385, 0.3023960396, 0.4532145234, 0.5771771523, 0.5563845840)> . 7.3 QR-Algorithmus Der am häufigsten verwendete Algorithmus zur Eigenwertberechnung ist der sogenannte QR-Algorithmus, den wir in diesem Abschnitt beschreiben wollen. In seiner Grundform lautet er wie folgt: Algorithmus 7.3.1 (Grundform QR-Algorithmus) Setze A(0) := A ∈ Rn×n , V (0) := I. for k = 0, 1, . . . Bestimme eine QR-Zerlegung von A(k) , also A(k) = Q(k) R(k) . Berechne A(k+1) = R(k) Q(k) . Setze V (k+1) = V (k) Q(k) . end Bemerkung 7.3.2 Die QR-Zerlegung kann z.B. nach Householder oder mit Givens-Rotationen erfolgen. Die Matrix V (k+1) muss nur berechnet werden, wenn auch eine Eigenvektornäherung zu bestimmen ist. Der Beweis eines Konvergenzsatzes für den QR-Algorithmus ist etwas aufwendiger, so dass hier auf ihn verzichtet wird. Um ihn formulieren zu können, benötigen wir folgende Bezeichnung: Definition 7.3.3 (Vorzeichenmatrix) Eine Matrix S ∈ Rn×n heißt Vorzeichenmatrix, wenn S eine Diagonalmatrix ist, deren 134 Einträge sii ∈ {−1, +1} für alle i = 1, . . . , n sind. Die Vorzeichenmatrizen sind notwendig, da die QR-Zerlegung einer regulären Matrix A nur bis auf eine Vorzeichenmatrix eindeutig bestimmt ist: Für zwei QR-Zerlegungen A = Q1 R1 und A = Q2 R2 gibt es eine Vorzeichenmatrix S mit Q1 = Q2 S und R2 = SR1 . Satz 7.3.4 (Konvergenz des QR-Algorithmus) Sei A ∈ Rn×n eine reguläre diagonalisierbare Matrix mit betragsmäßig einfachen Eigenwerten λ1 , . . . , λn ∈ R, die ohne Einschränkung betragsmäßig fallend angeordnet seien, d.h. |λ1 | > |λ2 | > . . . |λn | > 0. | | X = v1 . . . vn bezeichne eine Matrix, deren Spalten vi Eigenvektoren zu den Eigen| | werten λi von A sind. Besitzt X −1 eine LR-Zerlegung ohne Pivotisierung, so existieren Vorzeichenmatrizen Sk ∈ Rn×n mit kA(k) − Sk U Sk k2 → 0 für k → ∞, wobei U ∈ Rn×n eine von k unabhängige obere λ1 ∗ λ2 U = Dreiecksmatrix der Form ··· ∗ . .. . .. .. . ∗ λn ist. Insbesondere approximieren die Diagonalelemente von A(k) die betragsmäßig fallend sortierten Eigenwerte von A. Bemerkung 7.3.5 • Die Forderung der Existenz der LR-Zerlegung ohne Pivotisierung ist nur notwendig damit die Eigenwerte betragsmäßig geordnet erscheinen. Für die Konvergenz des Verfahrens benötigt man sie nicht. • Es wurde explizit gefordert, dass die Eigenwerte reell sind. Eine reelle, diagonalisierbare Matrix kann auch komplexe Eigenwerte haben. Dann ist aber immer auch das komplex konjugierte ein Eigenwert und somit hätte man zwei Eigenwerte mit dem gleichen Betrag, was den Voraussetzungen im Satz widerspricht. • Hat man zwei betragsgleiche Eigenwerte λp , λp+1 so konvergiert die Folge der Matrizen A(k) nur noch in den Spalten, die nicht zu den betragsgleichen Eigenwerten 135 gehören. In den anderen Spalten kann sie divergieren, jedoch konvergieren die Eigenwerte des auf der entsprechenden Diagonalen befindlichen 2x2-Blocks noch gegen λp und λp+1 . • Im Spezialfall symmetrischer Matrizen folgt aus der Symmetrie aller A(k) und der Bedingung kA(k) − Sk U Sk k2 → 0 für k → ∞, dass U eine Diagonalmatrix sein muss und somit ist Sk U Sk = U. Der Aufwand des QR-Algorithmus in seiner Grundform ist zu hoch, da in jedem Schritt eine QR-Zerlegung bestimmt werden muss. Diesen kann man verringern, indem man die Matrix A ∈ Rn×n zunächst in eine obere Hessenberg-Matrix (alle Elemente unterhalb der ersten unteren Nebendiagonalen sind 0) mit gleichen Eigenwerten transformiert, so dass die QR-Zerlegung sehr effizient (mittels n − 1 Givens-Rotationen) durchführbar ist. Ferner lässt sich die relativ langsame Konvergenzgeschwindigkeit verbessern, indem man Shift-Parameter einführt. Algorithmus 7.3.6 (QR-Algorithmus) Bestimme eine orthogonale Matrix U ∈ Rn×n derart, dass A(0) := U > AU eine obere Hessenberg-Matrix ist. Setze V (0) := I. for k = 0, 1, . . . Wähle einen Shift-Parameter µ(k) . Bestimme eine QR-Zerlegung von A(k) − µ(k) I, also A(k) − µ(k) I = Q(k) R(k) . Berechne A(k+1) = R(k) Q(k) + µ(k) I. Setze V (k+1) = V (k) Q(k) . end Um den Algorithmus durchführen zu können benötigt man eine Methode um A auf obere Hessenberg-Matrixform zu transformieren. Für symmetrische Matrizen bedeutet dies sogar Tridiagonalgestalt. Die Idee hierzu besteht darin geeignete Householder-Spiegelungen Qk ∈ Rn×n zu berechnen, die unterhalb der ersten unteren Nebendiagonalen die notwendigen Nullen erzeugen. Dann setzt man U := Q1 Q2 · · · Qn−2 . Schematisch sieht das ganze bei einer 5x5 Matrix wie folgt aus, wobei wir mit + die Elemente kennzeichnen, die sich beim jeweiligen Schritt verändern, und * die (im Allge136 meinen) von Null verschiedenen Elemente kennzeichnet: ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ A1 := A = ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ + + + + + + ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ + ∗ ∗ ∗ ∗ > → Q> A = → A := Q A Q = 0 ∗ ∗ ∗ ∗ 0 ∗ ∗ ∗ ∗ 1 2 1 1 1 1 0 ∗ ∗ ∗ ∗ 0 ∗ ∗ ∗ ∗ 0 ∗ ∗ ∗ ∗ 0 ∗ ∗ ∗ ∗ + + ∗ ∗ ∗ + + + + + + + ∗ ∗ ∗ + + + + + > 0 + ∗ ∗ ∗ → A := Q A Q = → Q> ∗ ∗ ∗ ∗ 3 2 2 2 2 A2 = 0 0 0 ∗ ∗ ∗ 0 0 ∗ ∗ ∗ 0 0 ∗ ∗ ∗ 0 0 ∗ ∗ ∗ + + + ∗ ∗ + + + + + + + + ∗ ∗ + + + + + > → Q> 3 A3 = 0 + + + + → A4 := Q3 A3 Q3 == 0 + + ∗ ∗ 0 0 + ∗ ∗ 0 0 ∗ ∗ ∗ 0 0 0 ∗ ∗ 0 0 0 ∗ ∗ Würde man, was zunächst naheliegt, die Matrizen Qk so wählen, dass alle Elemente (k) unterhalb der Diagonalen bei der Berechnung von Q> Null werden, so würde man die kA (k) erzeugten Nullen im nächsten Schritt bei der Berechnung von Q> k A Qk im Allgemeinen wieder zerstören, würde also am Ende keine Diagonalgestalt bekommen. Ein detaillierter Algorithmus läst sich mit Hilfe der Householder-Spiegelungen wie im Kapitel über QRZerlegungen herleiten. Ferner benötigt der QR-Algorithmus einen Shift-Parameter. Bekannte Wahlmöglichkeiten für diesen sind: (k) • Rayleigh-Quotienten-Shift: µ(k) = ann ; (k) • Wilkinson-Shift: µ(k) ist der am dichtesten an ann gelegene Eigenwert von ! (k) (k) an−1,n−1 an−1,n . (k) (k) an,n−1 an,n In beiden Fällen wird versucht den kleinsten Eigenwert von A möglichst gut zu approximieren. 137 Bemerkung 7.3.7 • Man kann (mit Hilfe der Givens-Rotationen) beweisen, dass alle im Algorithmus erzeugten Matrizen A(k) obere Hessenberg-Matrizen sind. • Statt mit U, V (k) kann man die Eigenvektoren günstiger sukzessive mit der Inversen Iteration von Wielandt mit µ = λ1 , . . . , µ = λn berechnen. Bemerkung 7.3.8 (Deflation) Hat eine Matrix eine Blockgestalt der Form M= A C 0 B ! mit quadratischen Matrizen A und B, so gilt ! ! ! A C I 0 A C M= = . 0 B 0 B 0 I Damit ist ersichtlich, dass die Menge der Eigenwerte der Matrix M genau die Menge der Eigenwerte der Matrix A vereinigt mit der Menge der Eigenwerte der Matrix B ist. Man kann also das Problem der Eigenwertberechnung von M auf zwei kleinere Eigenwert-Probleme zurückführen. Sind bei einer Matrix A(k) im QR-Algorithmus also untere Nebendiagonalelemente Null (betragsmäßig sehr klein), so kann man die Matrix in Blöcke zerlegen und auf jeden einzelenen den QR-Algorithmus anwenden, wobei die obere Hessenberg-Matrixform in jedem Block vorliegt. Dieser Vorgang heißt Deflation. Beispiel 7.3.9 (QR-Algorithmus) Gesucht sind die Eigenwerte der Tridiagonalmatrix 12 1 0 0 1 9 1 0 A= 0 1 6 1 0 0 1 3 0 0 0 1 0 0 0 . 1 0 Wir verwenden den QR-Algorithmus mit Wilkinson-Shift und ε = 10−16 als Parameter für die Durchführung der Deflation. Man erhält: (k) (k) k a54 a55 µ(k) 0 1.00000000000000 0.00000000000000 -0.30277563773199 1 0.00454544295112 -0.31686978239085 -0.31687587422668 2 0.00000000010677 -0.31687595261688 -0.31687595261688 3 0.00000000000000 -0.31687595261688 138 (k) Nach drei Iterationen ist das untere Nebendiagonalelement a54 < ε. Wir nutzen die Deflation und rechnen mit dem oberen 4x4 Block weiter. (k) (k) k a43 a44 µ(k) 3 0.14372385063470 2.99069135876582 2.98389967723150 4 0.00000171156219 2.98386369683915 2.98386369683818 5 0.00000000000000 2.98386369683818 Wir haben erneut Deflation und rechnen mit dem oberen 3x3 Block weiter. (k) (k) k a32 a33 µ(k) 5 0.07800880529022 6.00201597257969 6.00000324472688 6 0.00000008388556 6.00000000000000 6.00000000000000 7 0.00000000000000 6.00000000000000 Nach Deflation bleibt eine 2x2 Matrix, deren Eigenwerte sich analytisch bestimmen lassen. Nach insgesamt 7 Iterationen haben wir die fünf Eigenwerte: λ1 λ2 λ3 λ4 λ5 ≈ 12.31687595261687, ≈ 9.01613630316183, ≈ 6.00000000000000, ≈ 2.98386369683818, ≈ −0.31687595261688. Wendet man den QR-Algorithmus in seiner Grundform ohne Shift und ohne Deflation auf dieses Beispiel an, so benötigt man 43 Iterationen um mit 12.316875952 0.000001559 0 0 0 0.000001559 9.016136303 0.000000024 0 0 0 0.000000024 6.000000000 0 0 0 0 0 2.983863696 0 0 0 0 0 −0.316875952 eine Matrix zu erhalten, deren Diagonalelemente weitgehend mit den Eigenwerten übereinstimmen, wenngleich einige Nebendiagonalelemente noch deutlich von Null verschieden sind. 139 Literaturverzeichnis [DH91] Deuflhard, P. und Hohmann, A. Numerische Mathematik. de Gruyter, Berlin, 1991. [GK99] Geiger, C. und Kanzow, C. Numerische Verfahren zur Lösung unrestringierter Optimierungsaufgaben. Springer, Berlin-Heidelberg-New York, 1999. [Kel95] Kelley, C. T. Iterative Methods for Solving Linear and Nonlinear Equations, volume 16 of Frontiers in Applied Mathematics. SIAM, Philadelphia, 1995. [Kos93] Kosmol, P. Methoden zur numerischen Behandlung nichtlinearer Gleichungen und Optimierungsaufgaben. Teubner, Stuttgart, 2nd edition, 1993. [Lem98] Lempio, F. Numerische Mathematik II – Methoden der Analysis. volume 55 of Bayreuther Mathematische Schriften. Bayreuth, 1998. 140
© Copyright 2025 ExpyDoc