Die Newton-Iteration
xk+1 =xn− f(xk)
f′(xk), k ∈N0 (151)
zur Approximation einer Nullstelle der differenzierbaren reellen Funktionf wurde von Isaac Newton um 1670 herum am Beispiel einer kubischen Gleichung entwickelt. Joseph Raphson beschrieb das Verfahren allgemein in einer Arbeit über das Lösen von Gleichungen, und die abstrakte Form (151) erhielt sie gegen 1735 von Thomas Simpson88.
Das Verfahren beruht auf der Linearisierung der Gleichung f(x) = 0
bei der Approximationxk:
f(x)≈f(xk) +f′(xk)(x−xk) und der ersatzweisen Lösung von
f(xk) +f′(xk)(x−xk) = 0 (152) durch xk+1, was schließlich auf (151) führt.
86Wir setzen hier einmal globale Stetigkeit voraus, um nicht zu viele Fälle berücksichtigen zu müssen.
87Mit demSatz von Arzela-Ascoli: Sei(R, d(·,·))ein kompakter metrischer Raum undFeine Familie stetiger Funktionen von R in einen Banachraum (V,∥ · ∥), die im Raum der stetigen Funktionen von R nach V abgeschlossen ist. Dann istF genau kompakt, wennF beschränkt und gleichmäßig gleichgradig stetig ist.
Lezteres bedeutet, dass es zu jedemε >0einδ >0unabhängig vonx∈R(gleichmäßig!) und unabhängig vonf ∈ F (gleichgradig!) gibt, so dass
∥f(y)−f(x)∥ ≤εfür alled(y, x)< δ und für allef ∈ F.
88Bekannt durch die Simpsonsche Integrationsformel, die wiederum eigentlich auf Kepler (1615) zurück-geht.
−2.5 −2 −1.5 −1 −0.5 0 0.5 1 1.5 Newton−Iteration zur Bestimmung einer Nullstelle von f(x) = cos(x)+x/2
x1
x3
x2 x
0
f(x)=cos(x)+x/2
Abbildung 22: Newton-Iteration für f :R7−→R Die Linearisierung (152) verwendet man auch für Gleichungssysteme
F(x) = 0,
wennF eine - mindestens - differenzierbare Funktion ist mit
F :Rn −→Rm (153)
oder auch
F : (Ba,∥ · ∥a)−→(Bb,∥ · ∥b, (154) worin Ba und Bb Banachräume mit zugehörigen Normen sind.
In diesem Fall lautet die Linearisierung
F(x)≈F(xk) +F′(xk)(x−xk), (155) worin F′(xk), die Frechet-Ableitung89 oder einfach die Ableitung von F bei xk ist.
Im Fall von (153) und
F(x) =
f1(x1, . . . , xn) f2(x1, . . . , xn)
...
fm(x1, . . . , xn)
(156)
ist die lineare Abbildung F′(x) durch die Jacobi-Matrix
F′(x) = J(x) :=
∂f1
∂x1(x) ∂f∂x1
2(x), . . . , ∂x∂f1
n(x)
∂f2
∂x1(x) ∂f∂x2
2(x), . . . , ∂x∂f2
n(x) ... ... . . . ...
∂fm
∂x1(x) ∂f∂xm
2(x), . . . , ∂f∂xm
n(x)
(157)
89Eine Funktion (154) ist in einem inneren Punktxihres DefinitionsbereichesD Frechet-ableitbar oder auch einfach nur ableitbar, wenn es eine lineare AbbildungA: (B1,∥ · ∥a)−→(B2,∥ · ∥2)gibt, mit der
F(x+h) =F(x) +Ah+R(h)für alleh∈B1 mitx+h∈D gilt, wobei die RestfunktionR die Bedingung
lim
x+h∈D,∥h∥1→0
∥R(h)∥2
∥h∥1
= 0
erfüllt.
Die Zuordnung der spezifischen linearen AbbildungAzum Differentiationspunktxmacht man durch die BezeichnungF′(x)kenntlich.
′ ′
gegeben90. Für den Fall
F(x) = 0, F ∈C1(Rn,Rn),
also bei einem System mit endlich vielen Unbekannten und ebensovielen Gleichungen haben wir das Newton-Verfahren sowie zwei approximative Varianten hiervon schon kurz auf der Seite 63 kennengelernt und mit dem Satz von Ostrowski sogar schon einer ersten Konvergenzanalyse unterzogen.
Solchen Aufgaben wird zunächst weiter unser Hauptaugenmerk gelten.
Wir weisen aber schon hier darauf hin, dass die Linearisierung (155) und ihre (gegebenen-falls näherungsweise Lösung) auch Grundlage von Verfahren ist, wenn
(i) die Anzahl der Gleichungen die der Unbekannten überwiegt oder (ii) es weniger Gleichungen als Unbekannte gibt.
Im ersten Fall wird man die Linearisierung als Ausgleichsaufgabe ansehen und dadurch zum sogenannten „Gauss-Newton-Verfahren“ geführt. Im zweiten Fall hat das Gleichungssystem (im Falle maximalen Zeilenranges) eine lineare Mannigfaltigkeit als Lösungsgesamtheit und man wählt durch geeignete Zusatzbedingungen unter diesen Lösungen die Schrittrichtung aus.
Zunächst arbeiten wir aber weiter am Fall n=m.
3.2.1 Lokale Konvergenz des Newton-Verfahrens
Mit dem Satz von Ostrowski wurde oben gezeigt, dass das Newton-Verfahren für Funktio-nen F : Rn ⊃ D −→ Rn, die in einer Umgebung einer regulären Nullstelle91 x∗ zweifach stetig differenzierbar sind, lokal quadratische gegen die Nullstelle konvergiert.
Diese Voraussetzungen sind etwas zu stark. Üblicherweise setzt man voraus die Differen-zierbarkeit in einer Umgebung92U der regulären Nullstellex∗ sowie die Annahme, dass sich die Ableitung F′(x)in dieser Umgebung mit x „nicht zu schnell ändert“. Eine klassische93 Realisierung dieser Voraussetzung ist die Lipschitzbeschränktheit der Ableitung in dieser Umgebung
∥F′(x)−F′(y)∥ ≤γ∥x−y∥für alle x, y ∈U. (158) Dann können wir wie folgt versuchen, das lokale Fehlerverhalten94 in den Griff zu bekom-men. Aus
xk+1 :=xk−F′(xk)−1F(xk)
90Achtung: Man beachte, dass die Existenz aller partiellen Ableitungen ∂x∂fi
j(x)als Elemente von J(x) für die Differenzierbarkeit von F in xund alsoF′(x) =J(x) nicht reicht. Eine hinreichende Bedingung dafür ist aber schon die Stetigkeit der partiellen Ableitungen inx.
91Das hieß, es warF(x∗) = 0 unddetF′(x∗)̸= 0
92Zur Übung wiederholen wir hier, dass man eine MengeU in der Mathematik gewöhnlich als Umgebung eine Punktes bezeichnet, wenn neben dem Punkt selbst ein ganze Normkugel um den Punkt zuU gehört.
93Die Klassik liegt in der Numerischen Mathematik noch gar nicht so weit zurück, nämlich in den Tegen der ersten (relativ-) Großrechner in der Mitte des letzten Jahrhunderts.
94Das bedeutet das Verhalten in einem Bereich, in dem Newton wie gleich gesehen wird, von allen Startpunkten aus gegenx∗ konvergiert.
schließen wir für die Fehler xk+1−x∗ und xk−x∗ das Folgende:
xk+1−x∗ = xk−F′(xk)−1F(xk)−x∗ =xk−x∗−F′(xk)−1(F(xk)−F(x∗))
= F′(xk)−1(
F(x∗)−[F(xk) +F′(xk)(x∗−xk)])
| {z }
T
. (159)
Wenn wir uns im Eindimensionalen befänden und wenn F zweimal stetig differenzierbar wäre, so wären wir schon fertig, denn der gekennzeichnete Ausdruck T ist nichts anders als der Fehler der Auswertung der Linearisierung von F bei xk an der Stelle x∗, und wir könnten schließen
T =F(x∗)−[F(xk) +F′(xk)(x∗−xk)] = F′′(ζ)
2 (x∗−xk)2.
Hier gehen wir etwas anders vor, gruppieren die Elemente inT etwas um und finden - wobei wir U als konvex voraussetzen - mit dem Mittelwertsatz von Seite 122 die Darstellung
xk+1−x∗ = F′(xk)−1 (
[F(x∗)−F(xk)]−F′(xk)(x∗−xk)])
= F′(xk)−1 ∫1
0
(F′(xk+t(x∗−xk))−F′(xk))
(x∗−xk)dt. (160) Indem wir zu Normen übergehen, (158) verwenden und - gegebenenfalls unter Verkleine-rung der UmgebungU - die Norm der Inversen F′(xk) über
∥F′(xk)−1∥ ≤2· ∥F′(x∗)−1∥:=β abschätzen95, gelangen wir zu
∥xk+1−x∗∥ ≤ ∥F′(xk)−1∥∫1
0 F′(xk+t(x∗−xk))−F′(xk)· ∥x∗ −xk∥dt
≤ ∥F′(xk)−1∥ ∫1
0 γt∥x∗−xk∥dt∥x∗−xk∥ ≤βγ∥x∗−xk∥2. (161) Aus dieser Abschätzung gewinnen wir erstens einmal überhaupt lokale Konvergenz, indem wir x0 so nahe bei x∗ annehmen, dass βγ· ∥x0−x∗∥ < 1 ist, so dass die Folge {xk}k∈N0
im Kreis mit Radiusr:=∥x0−x∗∥bleibt und dort bei jedem Iterationsschritt mindestens um den Faktor r anx∗ heranrückt.
Zweitens hat man natürlich auch die quadratische Ordnung der dann als konvergent er-kannten Folge.
Fassen wir zusammen:
Satz 3.29 (Lokale quadratische Konvergenz des Newton-Verfahrens)
Sei D⊂Rn offen und konvex. Sei F ∈C1,1(D,Rn) mit F(x∗) = 0 für ein x∗ ∈D. Sei
∥F′(x∗)−1∥ ≤ β 2 und
∥F′(x)−F′(y)∥ ≤γ∥x−y∥ für alle x, y ∈D. (162) Dann gibt es ein r > 0, so das für alle x0 in der Kr(x0) := {x ∈ Rn | ∥x−x∗∥ ≤ r} die Folge
xk+1 :=xk−F′(xk)−1F(xk) inKr wohldefiniert ist, und es ist
∥xk+1−x∗∥ ≤βγ ∥xk−x∗∥2, k ≥0.
3.2.2 Affin-Invarianz
Wenn wir das nichtlineare System
F(x) = 0, F ∈C1,1(R,Rn) (163) mit einer regulären (n, n)-MatrixA multiplizieren, ist das sich ergebende Problem
G(x) = 0, G(x) = AF(x) (164)
sicherlich äquvalent zu (163) in dem Sinne, dass beide dieselben Lösungen haben.
Wendet man das Newton-Verfahren auf (164) an, so ergibt sich überdies nach
xk+1 :=xk−G′(xk)−1G(xk) =xk−(AF′(xk))−1AF(xk) = xk−F′(xk)−1F(xk), dass es ganz unbeeindruckt vom affinen Übergang von (163) zu (164) ist und dieselben Iterationsvektoren liefert.
Man sagt, dass die Lösungsmenge von (163) und das Newton-Verfahren invariant unter der affinen Transformation (164) seien oder auch einfach, sie seien affin-invariant.
Die Affininvarianz besagt zum Beispiel, dass wir auf Gleichungssysteme (163) lineare Trans-formationen anwenden können - etwa solche, wie wir sie bei der Gauss-Elimination wendeten, Vielfache einer Gleichung von einer anderen subtrahieren oder Gleichungen ver-tauschen - ohne dass sich die Lösungen änderten und ohne dass das die Newton-Iteration beeinflusste96.
Fallbeispiel:
Als erstes einfaches Beispiel für die Nützlichkeit der Affininvarianz der Lösungs-menge eines Gleichungssystemes betrachten wir das nichtlineare System
(x+ 1)2+y2−2 = 0,
(x−1)2+y2−2 = 0. (165) Hierin hat die erste Gleichung als Lösungsgesamtheit einen Kreis mit Radius
√2 um den Punkt (xa, ya) = (−1,0) und die zweite Gleichung den Kreis um (xb, yb) = (1,0) mit demselben Radius.
Die Lösungen des Gleichungssystems sind offenbar (vgl. Abbildung 23) die Punkte (x1, y1) := (0,1) und (x2, y2) = (0,−1).
−3 −2 −1 0 1 2 3
−2
−1.5
−1
−0.5 0 0.5 1 1.5 2
(x+1)2+y2=2
(x−1)2+y2=2
Abbildung 23: Gleichungssystem „Zwei Kreise“
96Dass das Letzte durchaus nicht so selbstverständlich ist, sieht man etwa bei den einfachsten Iterati-onsverfahren für lineare Systeme. Gesamt- und Einzelschrittverfahren reagieren extrem empfindlich auf Änderungen des Systems; schon das Vertauschen von Zeilen oder Spalten kann das Verhalten vollständig ändern.
Führen wir nun einen „Gauss-Eliminationsschritt“ aus, indem wir die erste Glei-chung von der zweiten abziehen97, so entsteht das neue System
(x+ 1)2+y2−2 = 0,
−4x = 0. (166)
Hier ist die zweite nichtlineare Gleichung ersetzt worden durch die lineare Glei-chung
−4x= 0.
−3 −2 −1 0 1 2 3
−2
−1.5
−1
−0.5 0 0.5 1 1.5
2 (x−1)2+y2=2
−4x=0
Abbildung 24: Äquivalent zu „Zwei Kreise“
Aus dem Wissen, dass die Newton-Iteration für beide Systeme (165) und (166) dieselbe Iterationsfolge liefert und mit der Erinnerung, dass Newton lokal die Linearisierung des Systems löst, folgert man, dass jeder Newtonschritt - egal woher - sofort auf der Gerade x = 0 landen muss und die Iterationsfolge an-schließend dort verbleibt.
Da jede nichttriviale Linearkombination der zwei Kreisgleichung (165) wieder eine Kreisgleichung ergibt98 und diese Gleichung mit den beiden Grundglei-chungen durch die Punkte (0,1) und (0,−1) gelöst wird, führt jede reguläre Transformation (164) auf die Gleichung zweier verschiedener Kreise des Kreis-büschels durch die beiden Punkte (0,1) und (0,−1) (vgl. Abbildung 25).
−3 −2 −1 0 1 2 3
−2
−1.5
−1
−0.5 0 0.5 1 1.5 2
Einige Mitglieder des Kreisbüschels durch (0,1) und (0,−1)
Abbildung 25: Kreisbüschel Gleichungssystem: Zwei Kreise
97Das entspricht der Multiplikation des Systems mit der regulären Matrix
( 1 0
−1 1 )
.
98Oder die Gleichung der Geradex= 0, eines Kreises mit Radius∞.
Die Punkte eines jeden der Kreise sind die Lösung einer speziellen zugehörigen Linearkombination
0 = (a
b )T
F(x, y) = a·(
(x+ 1)2+y2−2) +b·(
(x−1)2+y2−2) .
Das heißt, dass F(x, y) überall orthogonal zu (a
b )
sein muss. Im vorliegenden einfachen Fall F : R2 −→ R2 muss äquivalent F(x, y) überall auf dem Kreis eine Vielfaches von
(−b a
) sein:
F(x, y) =µ (−b
a )
, µ∈R (167)
Ist(x0, y0)ein Punkt dieses Kreises mitF(x0, y0)̸= 0, so gilt für den durch ihn verlaufenden Kreis sicher
F(x0, y0) = ˆµ (−b
a )
mit einem reellenµˆ ̸= 0.
Dies zeigt, dass der Vektor (−b
a )
ein nichttriviales Vielfaches vonF(x0, y0)ist.
Daher können wir das parameterabhängige Gleichungssystem (167) alternativ aufschreiben als
F(x, y) = (1−λ)F(x0, y0), λ∈R. (168) Wir werden diese parameterabhängige Gleichung etwas später sehr intensiv weiter untersuchen.
Während die Gleichung (163) und das Newton-Verfahren unempfindlich gegenüber dem Übergang nach (164) sind, ist es die Fehleranalyse aus dem letzten Konvergenzlemma nicht.
Ersetzt man in (161)F durch AF, so erhält man dort - wenn man die Lipschitzbedingung (158) beibehalten will (und ansonsten naiv weiter abschätzt) das Ergebnis
∥xk+1−x∗∥ ≤ ∥(AF′(xk))−1∥∫1
0 AF′(xk+t(x∗ −xk))−AF′(xk)· ∥x∗−xk∥dt
≤ ∥F′(xk)−1A−1∥ ∫1
0 ∥A∥γt∥x∗−xk∥dt ∥x∗−xk∥
≤ ∥A−1∥ · ∥A∥ ·βγ∥x∗−xk∥2.
Dies Ergebnis hat sich um den Faktor cond∥·∥(A) = ∥A−1∥ · ∥A∥verschlechtert.
Wenn wir die zu (161) führende Fehlerdarstellung (160) näher anschauen xk+1−x∗ = F′(xk)−1 (
[F(x∗)−F(xk)]−F′(xk)(x∗−xk)])
= ∫1
0 F′(xk)−1(
F′(xk+t(x∗−xk))−F′(xk))
(x∗−xk)dt, (169) so sehen wir, dass wir in (161) letztlich eine Abschätzung
∥F′(xk)−1(
F′(xk+t(x∗−xk))−F′(xk))
(x∗−xk)∥ ≤ωt∥x∗−xk∥2 (170) gesichert haben. Diese Abschätzung ist offenbar affin-invariant, denn jeder reguläre Faktor A von F hebt sich automatisch wieder heraus.
Indem wir (169) direkt mit (170) abschätzen, erhalten wir :
Lemma 3.30 (Affin-invariante lokale Konvergenz)
Es sei D ∈ Rn offen und konvex. F : D −→ R sei in D stetig differenzierbar. Es seien x∗ ∈ D mit F(x∗) = 0. Es sei F′(x) regulär in ganz D und erfülle die affin-invariante Lipschitzbedingung
∥F′(x)−1(F′(x+t(y−x))−F′(x)) (y−x)∥ ≤ωt∥y−x∥2, für alle x, y ∈D und t∈[0,1].
(171) Dann gibt es eine Kugel Kr(x∗) :={x∈D| ∥x−x∗∥ ≤r} ⊂D, so dass
Φ(x) :=x−F′(x)−1F(x) Kr(x∗)in sich abbildet, und es gilt
∥x∗−Φ(x)∥ ≤ ω
2∥x∗ −x∥2 für alle x∈Kr(x∗).
Anmerkungen 3.31
Dem Leser mag die Bemerkung auf der Zunge liegen, dass die Lipschitzbedingung (171) aber um einige komplizierter aussähe als die Bedingung (158). Dem ist zu entgegegnen, dass die Lipschitzkonstanten in beiden Fällen praktisch nie genau gefunden werden können, da hier ja letztendlich Maxima über alle unendlich vielen Elemente in Dzu bilden sind.
Zugegebenerweise sieht der Ausdruck der in (171) beschränkt werden muss, komplizierter aus als der in(158). Man wird in der Praxis aber für beide Werte nur Approximationen durch die Auswertung an endlich vielen Stellenx, y bilden können, und es stellt sich heraus, dass die komplizierter wirkende Lipschitzbedingung (171) im Verlauf der Iteration nicht schwieriger zu bilden ist als (158).
3.2.3 Globales Verhalten des Newtonverfahrens
Mit den letzten Konvergenzsätzen haben wir die schöne Gewissheit, dass das Newton-Verfahren jede reguläre Nullstelle approximiert, wenn man nur nahe genug bei ihr startet und dass - wenn das Newton-Verfahren gegen eine reguläre Nullstelle konvergiert - dies auch mit befriedigender zunehmender Geschwindigkeit geschieht.
Ob und wann die Newton-Iteration in die Region quadratischer Konvergenz eintritt, ist aber leider für weiter von einer Nullstelle entfernte Startwerte ganz unklar.
In der Abbildung 26 ist z.B. dargestellt, wie sich das Newton-Verfahren für die Berechnung von Nullstellen der Cosinus-Funktion verhalten kann. Man kann sich vorstellen, dass man durch geeignete Wahl des Anfangswertes beliebig lange Strecken durchmessen kann, bevor sich die Iteration bei einer der unendlich vielen Nullstellen schließlich fängt.
−8 −6 −4 −2 0 2 4 6 8 10 12
−1
−0.8
−0.6
−0.4
−0.2 0 0.2 0.4 0.6 0.8 1
7 Newton−Steps to solve cos(x)=0 starting at x
1=0.15303.
x2
x3
x4
x5
x6
x1
x7
x8
Abbildung 26: Konvergenzprobleme bei Newton
Das das auch im Mehrdimensionalen ähnlich ist, zeigt die Abbildung 27 für das Itera-tionsverhalten von Newton in der komplexen Ebene zur Lösung der kubische Gleichung z3−1 = 0.
Abbildung 27: Julia-Menge der Newton-Iteration in C für f(z) = z3−1.
Die Punkte, von denen aus das Newton-Verfahren in die Lösung z1 := 1 konvergiert, sind blau eingefärbt, die zu z2 = exp(i2π/3) führenden Startpunkte sind rot, und die zu z3 = exp(i4π/3)leitenden sind schließlich grün.
Man sieht - konsistent mit den Sätzen über die lokale quadratische Konvergenz - dass es um diese Nullstellen herum Gebiete Kreise gibt, die vollständig die entsprechenden Farben tragen. Diese werden also alle in die nahe gelegene Nullstelle geworfen.
Man sieht aber auch, dass die Konvergenzgebiete anders als man vielleicht erwartet, nicht einfach durch die Geraden
{rexp(iπ
3)|r ∈R+}, {−r |r∈R+}und {rexp(−iπ
3)|r ∈R+}
getrennt werden. Vielmehr sind die Trennungsmengen komplizierte sogenannte „Julia-Mengen“ (vgl. [JUL]). Neben den oft zitierten Selbstähnlichkeiten99 hat die Menge die bemerkenswerte Eigenschaft, dass jeder Punkt, von dem aus selbst keine Konvergenz zu einer der drei Nullstellen stattfindet, Häufungspunkt sowohl der roten, als auch der blau-en als ebanfalls auch der grünblau-en Mblau-enge ist. Wblau-enn man nahe eines solchblau-en Punktes dblau-en Startwert beliebig wenig stört, kann man in eine jede der drei Nullstellen konvergieren. Die Literatur über dies Phänomen ist riesig.
Nun mögen diese Sachverhalte vom Standpunkt der reinen Mathematik grundsätzlich in-teressant sein100und auch gut in Form hübscher Bilder in allen möglichen Zeitschriften ver-öffentlichbar sein101, der professionelle Anwender der Newton-Verfahrens selbst ist durch diese Hüpferei aber eher gestört und erwartet von mathematischen Iterations-Algorithmen, dass sie möglichst fix in die Nullstelle laufen, die er haben will.
Das ist verständlich, aber für den die Befriedigung der Anwenderwünsche anstrebenden Angewandten Mathematiker macht diese Formulierung nun auch Probleme. Meistens ha-ben nichtlineare Gleichungssysteme nicht nur eine sondern mehrere Lösungen. Wie soll man nun die gewünschte Lösung aussortieren, vor allem dann, wenn man sie alle selbst noch nicht kennt?
Eine Methode, den Anwender bei der Auswahl mithelfen zu lassen, ist es, ihm eine Start-näherung abzuverlangen102. Mit dieser Näherung an der Hand sollte der Numeriker nun versuchen, die nächstgelegene Näherung anzusteuern.
Wie sich herausstellen wird, ist dabei schon gleich völlig unklar, was denn unter der nächst-gelegenen Nullstelle zu verstehen ist,
Aber dazu kommen wir später.
3.2.4 Globalisierungsstrategien: Dämpfung, Homotopie, Trust-Region
Ziel der sogenannten „Globalisierung“ ist es, Iterationen103, für die - wie in den letzten beiden Lemmata für das Newton-Verfahren - nachgewiesen werden kann, dass sie für hin-reichend nahe bei der Lösung gelegene Startnäherungen gegen die Lösung konvergieren, Variationen zu schaffen, die auch für weiter von der Lösung entfernt liegende Startnähe-rungen einsetzbar sind, indem sie sich einen Weg zu den Konvergenzgebieten erarbeiten und die dort dann in die Originaliterationen übergehen.
Die Methoden zur Globalisierung des Newton-Verfahrens sind im Wesentlichen durch die drei Schlagwörter „Dämpfung“, „Homotopie“ und „Trustregion-Globalisierung“ angerissen.
99Wenn man in diese Bilder hineinzoomt, stößt man immer wieder auf sich wiederholende Strukturen.
Vgl. [SAE].
100Die Analyse dieser Mengen spielt eine Rolle in der Chaos-Theorie, die wiederum schon eine Menge praktischer Anwendungen gefunden hat.
101Mit dem zweifelhaften Nebenerfolg, die Allgemeinheit weiter in ihrem Vorurteil zu bestärken, dass Mathematik nur Spielerei sei, mit der man letztlich nicht viel anfangen kann.
102Jeder ernstzunehmende Mathematikanwender, der ein nichtlineares Gleichungssysten lösen will, kennt eine vernünftige Näherung seiner Lösung. Ist dies nicht der Fall, hat er sich die Aufgabe nicht richtig über-legt. Den Ingenieur, der nach einem „robusten Iterationsverfahren verlangt“, das unabhängig von weiterer Information immer eine Lösung liefert, würde ich bei keiner Konstruktion sicherheitsrelevanter techni-scher Erzeugnisse wie bei der von Kernkraftwerken, Hochbrücken, Flugzeugen oder Küchen-Mixgeräten
Wir werden alle drei Methodenbereiche zunächst kurz skizzieren, bevor wir die ersten beiden davon etwas tiefer analysieren.
A. Dämpfung Unter Dämpfung der 1D-Newton-Methode xk+1 =xk−f′(xk)−1f(xk)
versteht man eine (möglicherweise) wiederholte Verkürzung der Schrittweite, wenn der durch das Newton-Verfahren vorgeschlagene Schritt keine verbesserte Näherung erbringt, wobei man die Güte durch den Betrag des Funktionswertes misst: Je kleiner desto besser.
−6 −4 −2 0 2 4 6 8 10 12
Newton für f(x)=x*cos(x), x0 =−3.34
x1 f(x0) =3.27
f(x1)= 5.84
x0
Zu groß!
Abbildung 28: Zu großer Newton-Schritt
Die einfachste Strategie, dies zu erzwingen, stammt aus dem Jahre 1966 von einem Herrn L. Armijo (vgl. [ARM]), heißt darum auch Armijo-Strategie und schreibt vor, dass im Falle eines Anwachsens der (absoluten) Größe der Funktion beim vorgeschlagenen Newton-Schritt die aktuelle Newton-Schrittweite halbiert wird. Ist also
|f(x0+N(x0))|>|f(x0)|,
so werden sukzessive|f(x0+12N(x0))|,|f(x0+14N(x0))|,|f(x0+18N(x0))|, .... getestet, bis der erste Wert|f(x0)| unterschreitet. Dass diese Rechnung nach endlich viel Halbierungen zum Ziel kommt104, liegt daran, dass die Ableitung von|f(x))|2 beix0 in Newton-Richtung
[d
dtf(x0+tN(x0))2]
|t=0 = 2f(x0)f′(x0)N(x0)
= −2f(x0)f′(x0)f′(x0)−1f(x0) =−2f(x0)2 <0 (172) negativ ist. |f| fällt also nahe x0 in Newton-Schritt-Richtung, und es gibt deshalb ein δ0, so dass |f(x0+sN)|<|f(x0)| für alle s∈(0, δ).
Diese Art der Schrittweitenreduktion zielt einerseits darauf, den Newton-Schritt zu ver-wenden, wenn er möglich ist; denn man möchte natürlich seine lokale quadratische Kon-vergenzrate auskosten.
Dennis und Schnabel halten das Primat des Newton-Schrittes für so wichtig, dass sie den Leser ihres Buch „Numerical methods for unconstrained optimization and nonlinear equa-tions“ [DS] mit
„Try Newton first“
104Jedenfalls bei theoretisch exakter Rechnung; bei der konkreten Rechnung mit beschränkter Stellenzahl muss man die Schleife abbrechen, wenn die Schrittweite unter die Rundungsgenauigkeit des Rechners fällt.
darauf einzuschwören versuchen.
Andererseits sorgt die sukzessive Verkleinerung dafür, dass ein Mindestanteil des mit Ab-stieg möglichen Weges auch gegangen wird105. Dies stellt sicher, dass dass Verfahren nicht einfach stehenbleibt, weil die Schrittweiten zu schnell gegen Null gefahren werden.
1 2 3 4 5 6 7 8 9 10
Immer noch zu groß zu groß
Newton−Schritt halber Schritt
Viertel−Schritt GUT
Abbildung 29: Armijo klappt Allerdings kann dies auch
−6 −4 −2 0 2 4 6 8 10 12
Newton für f(x)=x*cos(x), x 0 =−3.34
x1
f(x0) =3.27
f(x1)= 5.84
x0
(x0+x
1)/2
Blow up
Abbildung 30: Armijo steuert Nullstelle in Zoom-Viereck an zu einer Konvergenz gegen eine entferntere Nullstelle führen:
105Viele Implementierungen wählen deshalb bei der Schrittweitenverkleinerung eine Verkleinerungsfaktor, der zwar kleiner als 1 aber doch größer als12 ist. Zugleich kann man natürlich bei mehrfacher Verkleinerung die berechneten Funktionswerte verwenden, um – z.B. durch Interpolation – ein genaueres Bild über das Verhalten der Funktion in Newton-Richtung zu erhalten, um damit gegebenenfalls einen besseren Schritt zu finden. Vgl. dazu die „backtracking“-Algorithmen in [DS].
0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6
Blow−Up
(x0+x1)/2
Konvergenz, aber gegen die falsche Nullstelle
Abbildung 31: Zoom: Lokale Konvergenz
Um diesem unerwünschten Effekt entgegenzuwirken, versucht man, in dem an den Wertx0 anschließenden „ersten Tal“ der Betragsfunktion zu f(x) zu bleiben. Damit die Funktion, mit der man dies „testet“ nett glatt wird, wählt man anstelle von|f(x)|die „Testfunktion“
T(x) :=f(x)2.
Für die Analyse der Algorithmen definiert man erstens die Level- oder Niveau-Menge L(f, x0) oder kürzer L(x0) als die Menge all der Punkte im Definitionsbereich von D, an denen die Testfunktion nicht größer ist alsT(x0),
L(f, x0) :={x∈D|T(x)≤T(x0)}, (173) und zweitens die Zusammenhangskomponente Z(f, x0) von x0 in L(f, x0) über
Z(f, x0) :={x∈D|[x, x0]⊂L(f, x0) oder[x0, x]⊂L(f, x0)}. (174) Dies ist die Menge all der Punkte in D, die mit x0 durch einen stetigen Weg in L(f, x0) verbunden werden können106. In Abbildung (32) findet man Skizzen zu diesen Begriffen:
0 0.5 1 1.5 2 2.5
B(s):= f(0)−2f(0)2 x + C* x2
T(x):= f(x)2 x0
f(x)
Levelmenge L(x 0) Zusammenhangskomponente von L(x
0), die x
0 enthält
Abbildung 32: L(f, x0), Z(f, x0)und obere Schranke für Testfunktion f(x)2
Steht eine obere Schranke C für die zweite Ableitung von T(x0 + sN) in Z(f, x0) zur Verfügung, so findet man mit T(0) = f(x0)2 und dtdT′(x0 +tN) = −2f(x0)2 nach (172) eine obere Schranke
T(x0+sN)≤B(s) :=f(0)2−2f(0)2s+C 2s2
106Diese Sprechweise zielt auf die Definition von Zusammenhangskomponenten in allgemeineren Räumen alsR. InRsind solche Wege schlicht Intervalle.
die in N-Richtung erst einmal absteigt. Wählt man sˆ= f(0)C2, so dass B minimiert wird, so steigt man mit dem Schritt
x0 −→x1 =x0+ ˆsN(x0) sicher ab, bleibt aber auch wegen
Tˆ(x0+sN)≤B(s)≤T(x0) für alle s ∈[0,s]ˆ
sicher in der ZusammenhangskomponenteZ(f, x0). Zusätzlich zum Abstieg kann man zei-gen, dass dieser Schritt für den Nachweis der Konvergenz gegen ein Minimum von f(x)2 ausreichend groß ist. Leider ist hierbei (und auch bei anderen Dämpfungsstrategien ) nicht sicherzustellen, dass ein solches angesteuertes Minimum einer Nullstelle vonf ist.
Außerdem zeigt die nächste Abbildung 33 für das schon einmal weiter oben behandelte Beispiel f(x) =xcos(x), dass diese Strategie auch nicht garantiert, wirklich die nächstge-legene Nullstelle zu approximieren107. In Abbildung 33 liegen drei Nullstellen in Z(f, x0).
−4 −2 0 2 4 6 8 10
−4 −2 0 2 4 6 8 10
x0
Levelmenge L(x
0)
Zusammenhangskomponente von x 0
Drei Nullstellen in Zusammenhangskomponente
f(x)
[f(x)]
2Abbildung 33: Zusammenhang ist nicht genug
Wir werden im Abschnitt „Halbglobale Konvergenz des Newton Verfahrens“ die Konver-genzaussagen tatsächlich nur ein wenig erweitern können, müssen dafür aber im Wesentli-chen voraussetzen, dass man einer Nullstelle schon so nahe ist, dass die die x0 enthaltende Zusammenhangskomponente der Levelmenge als kritischen Punkt nur die Nullstelle ent-hält.
Zuvor kümmern wir uns allerdings erst einmal um zwei andere Ansätze der Globalisierung.
107Zu dieser Unsicherheit kommt hinzu, dass eine Schranke C selten wirklich zur Verfügung steht. Als
B.1. Homotopie Homotopie-Methoden108 benutzt man in der Numerik zur Globalisie-rung der Newton-Iteration, indem man die zu lösende Aufgabe
F(x) = 0, F :D−→Rn, D⊂Rn offen (175) mit einer durch ein x0 ∈D schon gelösten Aufgabe
G(x) = 0, G:D−→Rn (176)
homotop verbindet, wie wir das ja schon bei der Anwendung des Brouwerschen Fixpunkt-grades kennengelernt haben.
Es wird mit einer auf D×[0,1]mindestens stetigen Funktion und bezüglichx differenzier-baren Funktion H :D×[0,1]−→Rn das Kontinuum von Aufgaben
H(x, t) = 0, t∈[0,1] (177)
betrachtet, wobei mit
H(x,0) =G(x) und H(x,1) =F(x)
die „fortzusetzende“ Start-Gleichung (176) und die zu erreichende Zielgleichung (175) in das Kontinuum von Gleichungen (177) „eingebettet“ sind.
Gibt es109 für die parameterabhängige Gleichung (177) ein stetiges x: [0,1]−→Rn, mit H(x(t), t) = 0, x(0) =x0
so kann man versuchen, durch sukzessives Lösen mehrerer Aufgaben H(x, tk) = 0, 0 =t0 < t1 < . . . < tm−1 < tm = 1 am Ende auch die Aufgabe (175) zu lösen.
Im Falle eines Newton-Zuganges macht man sich die lokal (quadratische) Konvergenz zu-nutze, indem man den Schritt vontknachtk+1 jeweils so klein wählt, dass die Lösungx(tk) des Systems H(x, tk) = 0 in den Konvergenzbereich der Iteration
xj+1(tk+1) =xj(tk+1)−Hx(xj(tk+1), tk+1)−1H(xj(tk+1), tk+1), j ≥0, x0(tk+1) :=x(tk) fällt.
108auch unter den Namen Fortsetzungs- oder Einbettungsmethoden „gehandelt“
109Zu sichern etwa mit dem Satz über implizite Funktionen.
0 0.2 0.4 0.6 0.8 1 1.2
x(t)
Fortsetzung x(t
2) zu Startpunkt der Iteration bei (t
3)
Newton−Iterationen
G(x(0))=0
F(x(1))=0
x(t2) x(t
3) x(t
4)
Abbildung 34: Fortsetzungsmethode
Wie diese Fortsetzung effizient gestaltet werden kann, wie z.B. die Schrittlängen zwischen tk und tk+1 adaptiv gewählt werden, um ein möglichst rechenzeitsparendes Verfahren zu erhalten, werden wir später im Abschnitt 3.3 über die Behandlung parameterabhängiger Gleichungen ohnehin noch genauer besprechen. Der/die daran interessierte Leser/in wird hierfür auf diesen Vorlesungsabschnitt verwiesen.
Um von einer GleichungF(x) = 0und einem Startwertx0 zu einer Homotopie H(x, t)mit einer StartaufgabeG(x) =H(x,0) = 0zu kommen, dieG(x0) = 0erfüllt, gibt es sehr viele verschiedene Vorgehensweise.
Die einfachsten der verwendeten Funktionen G sind bei gegebenem Startwert x0 unter vielen weiteren die Funktionen
G1(x) := x−x0,
G2(x) := F′(x0)−1(x−x0), G3(x) := F(x)−F(x0).
(178)
Als Homotopie verwendet man meistens die Konvexkombination
H(x, t) = t·F(x) + (1−t)·G(x), t∈[0,1], (179) die man äquivalent auch wie folgt schreibt:
H(x, t) = G(x) +t(F(x)−G(x)), t∈[0,1]. (180) Für den FallG(x) = G3(x) = F(x)−F(x0)bekommt die Homotopie (180) die Form
Dies ist die Homotopie, die wir weiterhin betrachten werden. Sie ist unter verschiedenen Namen bekannt. Schwetlick nennt sie „defektreduzierende Einbettung“ (vgl. [SCHW]), weil Sie den Defekt ∥F(x(t))∥= (1−t)∥F(x(0))∥in jeder Norm linear in treduziert. Wir nen-nen Sie aus noch darzulegendem Grund „Newton-Homotopie".
Anmerkungen 3.32
(i) Es sei an dieser Stelle darauf hingewiesen, dass wir einen Anwendungsfall der Homo-topie (181) schon kennen. Im auf der Seite 83 beginnenden Fallbeispiel hatten wir ermittelt, dass die Kreise des zur Aufgabe (165) gehörenden Kreisbüschels genau die Lösungen dieser Homotopie sind (vgl. (167)).
(ii) Anders als bei dem Beispiel (165) führen die Lösungskurvenx(t)von (181) durchaus nicht bei allen Aufgaben zu einer Lösung des Problems F(x) = 0.
Dies ist von vornherein klar, wenn das Problem F(x) = 0 überhaupt keine Lösung hat wie zum Beispiel die eindimensionale Aufgabe
f(x) =x2+ 1 = 0.
Mit der Startnäherung x0 = 1 wird die Homotopie (181) zu
H(x, t) = x2−1 + 2t. (182)
Für t= 0 ist x0 = 1 offenbar wie gefordert eine Lösung von H(x,0) =x2−1 = 0.
Die Fortsetzung in Richtung größer werdender t-Werte versagt aber bei t∗ = 12 und zugehörigemx∗ = 0.
Da der x-Anteil der Jacobi-„Matrix“
H′(x∗, t∗) = (2x,2)|(x,t)=(0,1
2) = (0,2)
verschwindet, ist eine Newton-Iteration zur Anpassung von xin diesem Punkt nicht möglich.
Tatsächlich ist ja fürt = 12 +δ mit δ > 0die Gleichung 0 =H(x, t) = x2−1 + 2·
(1 2+δ
)
=x2+ 2δ nicht reell nachx auflösbar.
DaH′(x∗, t∗) = (0,2)aber vollen Rang hat, sagt der Satz über implizite Funktionen, dass der Lösungsast im (x, t)-Raum anstandslos fortsetzbar ist. Nur muss jetzt eine andere Richtung als die t-Richtung als Parameter gewählt werden.
Das wäre natürlich sofort aus der Gleichung (182) zu sehen gewesen. Während diese nur für t-Werte kleiner oder gleich 12 lokal nach x auflösbar ist, ist sie offenbar für alle x-Werte einfach nach t auflösbar.
t= 1−x2 2 .
−0.2 0 0.2 0.4 0.6 0.8 1
x0
t Startlinie t=0
Keine x−Lösung rechts dieser Linie Zilinielinie t=1
Umkehrpunkt (x*,t*) =(0,1/2) Versuch der
Fortsetzung bis t=1
Abbildung 35: Umkehrpunkt
Sehen wir t als abhängige Variable an, wird das Problem ganz einfach. Es handelt sich um eine einfache Parabel110. Vom Systemparameter x aus gesehen kommt es dagegen zu einem Zusammenbruch der einfachen Fortsetzung durch Erhöhung der t-Komponente.
Der hier beobachtete Effekt, dass der Lösungsast in t-Richtung „aufhört“, wobei er tatsächlich bezüglich der t-Richtung umkehrt111, wird häufig bei Homotopien beob-achtet.
Wenn der Parameter t nicht nur ein künstlicher numerischer Parameter ist, son-dern eine mit dem behandelten System verbundene Bedeutung hat (vgl. das Bratu-Problem oben), so sagt uns die Detektion eines solchen Umkehrpunktes, dass es rechts vont∗ (lokal) keine Lösung der Problems vorhanden ist. Wenn das Gleichungssystem z.B. den Zustand eines Reaktors beschreibt, sollte man sich vielleicht besser nicht in seiner Nähe befinden, wenn die Mathematik angibt, dass er nicht mehr existiert. Aus diesem Grund wird der Parameterwert t∗ gern auch kritischer Parameter genannt.
B.1. Kontinuierliche Fortsetzungsmethode Anstatt Werte der impliziten Funktion x(t) aus (181) also
F(x(t)) = (1−t)F(x0), x(0) =x0 (183) an endlich vielen t-Stellen zu bestimmen und so eine von vornherein t-diskrete Approxi-mation fürx(t)zu suchen, kann man durch Differentiation von (183) nacht die Homotopie (181) in eine Differentialgleichung für x(t) überführen.
Einfaches Differenzieren ergibt nämlich
F′(x(t))x′(t) = −F(x0).
Setzen wir voraus, dass F′(x(t))regulär ist112, so können wir schreiben
x′ =−F′(x(t))−1F(x0). (184) Diese „Davidenko-Differentialgleichung“ genannte Gleichung bildet mit der Anfangsbedin-gung
x(0) =x0
eine Anfangswertaufgabe, welche man mit Integrationsmethoden der numerischen Mathe-matik113 angehen kann, um sie vont = 0 bist = 1 numerisch zu integrieren.
110Beachten Sie bitte, dass die Auflösung nach t (allein) nicht möglich ist, wenn x ein Vektor (hoher Dimension) ist. Das Beispiel wird nur zu Demonstrationszwecken so einfach gehalten.
111um - im vorliegenden Fall - zu einer zweiten Lösung des Ausgangsproblems beit= 0zurück zu laufen