Kurs:Numerik I/5 Numerische Lösung nichtlinearer Gleichungssysteme

5.4 Das Newton-Verfahren im $\mathbb {R} ^{n}$

4.1 Grundlagen

Es sei $D\subseteq \mathbb {R} ^{n}$ eine offene Menge und $F:D\subseteq \mathbb {R} ^{n}\to \mathbb {R} ^{n}$ eine Funktion mit

$F(x):=(F_{1}(x),\ldots ,F_{n}(x))^{T},\quad x\in D$

und stetigen partiellen Ableitungen $(\partial F_{i}/\partial x_{j})(x)$ $(i,j=1,\ldots ,n)$, es sei also $F\in C^{1}(D,\mathbb {R} ^{n})$. Die zu $F$ gehörige Jacobi-Matrix bezeichnen wir mit

${\mathcal {J}}_{F}(x):=\left({\frac {\partial F_{i}}{\partial x_{j}}}(x)\right)_{i,j=1,\ldots ,n}\in \mathbb {R} ^{n\times n}.$

Mit $\|\cdot \|$ ist im ganzen Abschnitt 5.4 wieder die Euklidische Norm auf dem $\mathbb {R} ^{n}$ bzw. die durch sie induzierte Spektralnorm gemeint.

Schließlich werden wir von dem folgenden Lemma Gebrauch machen.

Lemma 5.17

Sei $v:[a,b]\to \mathbb {R} ^{n}$ eine stetige vektorwertige Funktion, sei $v(t):=(v_{1}(t),\ldots ,v_{n}(t))^{T}$ für $t\in [a,b]$ und sei $u:=\int \limits _{a}^{b}v(t)\,dt$ der Vektor mit den Komponenten $u_{i}:=\int \limits _{a}^{b}v_{i}(t)\,dt$. Dann gilt:

$\left\|\int \limits _{a}^{b}v(t)\,dt\right\|\leq \int \limits _{a}^{b}\|v(t)\|\,dt$

Beweis.

Es sei $K:=\|u\|$. Dann kann man unter Verwendung des Standardskalarprodukts

$\langle x,y\rangle :=\sum _{i=1}^{n}x_{i}y_{i},\quad x,y\in \mathbb {R} ^{n}$

auf dem $\mathbb {R} ^{n}$ und der Cauchy-Schwarz-Ungleichung abschätzen:

$K^{2}=\langle u,u\rangle =\left\langle \int \limits _{a}^{b}v(t)\,dt,u\right\rangle =\sum _{i=1}^{n}u_{i}\int \limits _{a}^{b}v_{i}(t)\,dt=\int \limits _{a}^{b}\sum _{i=1}^{n}u_{i}v_{i}(t)\,dt=\int \limits _{a}^{b}\langle u,v(t)\rangle \,dt\leq \int \limits _{a}^{b}\|u\|\|v(t)\|\,dt$

$=K\int \limits _{a}^{b}\|v(t)\|\,dt.$

q.e.d.

5.4.2 Das Verfahren

Es sei wieder $D\subseteq \mathbb {R} ^{n}$ eine offene Menge und es sei $F\in C^{1}(D,\mathbb {R} ^{n})$ mit

$F(x)=(F_{1}(x),\ldots ,F_{n}(x))^{T},\quad x\in D$

und Jacobi- bzw. Funktionalmatrix ${\mathcal {J}}_{F}(x)$ gegeben. Es soll nun das Newton-Verfahren zur Bestimmung einer Lösung $x^{*}\in D$ des Gleichungssystems

$F(x)=0\Leftrightarrow F_{i}(x_{1},\ldots ,xn)=0\quad (i=1,\ldots ,n)$

vorgestellt und seine Konvergenz untersucht werden. Für $n=1$ hatten wir dies bereits in Abschnitt 5.2.4 getan.

Die Iterationsvorschrift des Newton-Verfahrens lautete im Fall $n=1$

$x_{k+1}:=x_{k}-[f'(x_{k})]^{-1}f(x_{k}),$

wobei sich $x_{k+1}$ als Nullstelle einer linearen Approximation von $f$, der Tangente bei $x_{k}$ an $f$, ergab. Ähnlich kann man für eine Funktion $F\in C^{1}(D,\mathbb {R} ^{n})$ mit beliebigem $n\in \mathbb {N}$ das Newton-Verfahren dadurch motivieren, dass man $x^{k+1}$ als Nullstelle der linearen Approximation von $F$ bei $x^{k}$

$F_{k}(x):=F(x^{k})+{\mathcal {J}}_{F}(x^{k})\left(x-x^{k}\right)$

wählt. Dieses Vorgehen führt zu der allgemeinen Iterationsvorschrift

(5.23) $x^{k+1}:=x^{k}-\left[{\mathcal {J}}_{F}(x^{k})\right]^{-1}F(x^{k}),\quad k\in \mathbb {N} _{0}$

des Newton-Verfahrens. Wir gehen hier implizit davon aus, dass die Jacobi-Matrix ${\mathcal {J}}_{F}(x^{k})$ des Systems für jedes $k\in \mathbb {N} _{0}$ nichtsingulär ist. Da man, wenn immer möglich, die Berechnung der Inversen einer Matrix vermeiden sollte, geht man praktisch bei der Berechnung von $x^{k+1}$ von der zu (5.23) äquivalenten Gleichung

${\mathcal {J}}_{F}(x^{k})(\underbrace {x^{k+1}-x^{k}} _{=:h^{k}})=-F(x^{k})$

aus und bestimmt man die eindeutige Lösung $h^{k}$ des linearen Gleichungssystems

${\mathcal {J}}_{F}(x^{k})h=-F(x^{k}).$

Anschließend setzt man

$x^{k+1}:=x^{k}+h^{k}.$

Das Newton-Verfahren lautet somit wie folgt:

Algorithmus 9 (Newton-Verfahren)

(0) Wähle $x^{0}\in D$ und ein $\varepsilon >0$. Berechne $F(x^{0})$ und setze $k:=0$.

(1) Berechne ${\mathcal {J}}_{F}(x^{k})$ und bestimme die eindeutige Lösung $h^{k}\in \mathbb {R} ^{n}$ von

(5.24) ${\mathcal {J}}_{F}(x^{k})h=-F(x^{k}).$

(2) Setze $x^{k+1}:=x^{k}+h^{k}$ und berechne $F(x^{k+1})$.

(3) Falls $\left\|F(x^{k+1})\right\|_{2}\leq \varepsilon$, stop!

(4) Setze $k:=k+1$ und gehe nach (1).

Der folgende Satz besagt, dass das Newton-Verfahren unter geeigneten Voraussetzungen durchführbar, d. h. für alle $k$ insbesondere $x^{k}\in D$ und ${\mathcal {J}}_{F}(x^{k})$ nichtsingulär ist und dass es superlinear bzw. quadratisch konvergiert.

Satz 5.18

Es sei $D\subset \mathbb {R} ^{n}$ offen und $F\in C^{1}(D,\mathbb {R} ^{n})$. Ferner existiere ein $x^{*}\in D$, für welches $F(x^{*})=0$ und ${\mathcal {J}}_{F}(x^{*})$ nichtsingulär sei. Dann gibt es eine Umgebung ${\mathcal {U}}_{\delta }(x^{*})$ von $x^{*}$ für ein $\delta >0$, so dass das Newton-Verfahren, Algorithmus 9, für jeden Startpunkt $x^{0}\in {\mathcal {U}}_{\delta }(x^{*})$ durchführbar ist und die durch ihn ohne das Abbruchkriterium (3) erzeugte Iteriertenfolge $(x^{k})$ superlinear gegen $x^{*}$ konvergiert. Gilt mit einem $L>0$

(5.25) $\|{\mathcal {J}}_{F}(x)-{\mathcal {J}}_{F}(x^{*})\|\leq L\|x-x^{*}\|,\quad x\in {\mathcal {U}}_{\delta }(x^{*}),$

so konvergiert $(x^{k})$ gegen $x^{*}$ sogar quadratisch.

Beweis.

Wegen der Stetigkeit von ${\mathcal {J}}_{F}(x)$ auf $D$ können wir zunächst $\eta >0$ so klein wählen, dass gilt:

$\|{\mathcal {J}}_{F}(x)-{\mathcal {J}}_{F}(x^{*})\|\leq {\frac {1}{2\left\|[{\mathcal {J}}_{F}(x^{*})]^{-1}\right\|}},\quad x\in {\mathcal {U}}_{\eta }(x^{*}).$

Für $x\in {\mathcal {U}}_{\eta }(x^{*})$ ergibt sich damit und mit $\beta :=\left\|[{\mathcal {J}}_{F}(x^{*})]^{-1}\right\|$ gemäß Korollar 2.21 die Invertierbarkeit der Matrix

${\mathcal {J}}_{F}(x)={\mathcal {J}}_{F}(x^{*})+[{\mathcal {J}}_{F}(x)-{\mathcal {J}}_{F}(x^{*})]$

sowie die Abschätzung

(5.26) $\left\|[{\mathcal {J}}_{F}(x)]^{-1}\right\|\leq {\frac {\left\|[{\mathcal {J}}_{F}(x^{*})]^{-1}\right\|}{1-\left\|[{\mathcal {J}}_{F}(x^{*})]^{-1}\right\|\|{\mathcal {J}}_{F}(x)-{\mathcal {J}}_{F}(x^{*})\|}}\leq 2\beta .$

Sei nun

${\mathcal {N}}(x):=x-[{\mathcal {J}}_{F}(x)]^{-1}F(x),\quad x\in {\mathcal {U}}_{\eta }(x^{*})$

die Iterationsfunktion des lokalen Newton-Verfahrens, die nach dem Gezeigten auf ${\mathcal {U}}_{\eta }(x^{*})$ wohldefiniert ist. Mit $F(x^{*})=0$ und den Identitäten

$\int \limits _{0}^{1}{\mathcal {J}}_{F}(x^{*}+s(x-x^{*}))(x-x^{*})\,ds=F(x^{*}+s(x-x^{*})){\Big |}_{0}^{1}=F(x)-F(x^{*})=F(x)$

schließen wir als nächstes

${\mathcal {N}}(x)-x^{*}=x-x^{*}-[{\mathcal {J}}_{F}(x)]^{-1}[F(x)-F(x^{*})]$

$=x-x^{*}-[{\mathcal {J}}_{F}(x)]^{-1}\left\{{\mathcal {J}}_{F}(x^{*})(x-x^{*})+\int \limits _{0}^{1}[{\mathcal {J}}_{F}(x^{*}+s(x-x^{*}))-{\mathcal {J}}_{F}(x^{*})](x-x^{*})\,ds\right\}$

$=-[{\mathcal {J}}_{F}(x)]^{-1}[{\mathcal {J}}_{F}(x^{*})-{\mathcal {J}}_{F}(x)](x-x^{*})-[{\mathcal {J}}_{F}(x)]^{-1}\int \limits _{0}^{1}[{\mathcal {J}}_{F}(x^{*}+s(x-x^{*}))-{\mathcal {J}}_{F}(x^{*})](x-x^{*})\,ds.$

Für

(5.27) $\varepsilon (x):=2\beta \left\{\|{\mathcal {J}}_{F}(x^{*})-{\mathcal {J}}_{F}(x)\|+\int \limits _{0}^{1}\|{\mathcal {J}}_{F}(x^{*}+s(x-x^{*}))-{\mathcal {J}}_{F}(x^{*})\|\,ds\right\}$

leiten wir daraus unter Anwendung von Lemma 5.17 mit (5.26) die folgende Abschätzung ab:

(5.28) $\|{\mathcal {N}}(x)-x^{*}\|\leq \varepsilon (x)\|x-x^{*}\|.$

Wegen der Stetigkeit von ${\mathcal {J}}_{F}(x)$ auf ${\mathcal {U}}_{\eta }(x^{*})$ existiert ein $\delta \in (0,\eta ]$, so dass $\varepsilon (x)\leq 1/2$ auf ${\mathcal {U}}_{\delta }(x^{*})$ ist und damit gilt:

$\|{\mathcal {N}}(x)-x^{*}\|\leq {\frac {1}{2}}\|x-x^{*}\|,\quad x\in {\mathcal {U}}_{\delta }(x^{*}).$

Beginnend mit $x^{0}\in {\mathcal {U}}_{\delta }(x^{*})$, liegt folglich mit $x^{k}\in {\mathcal {U}}_{\delta }(x^{*})$ auch $x^{k+1}:={\mathcal {N}}(x^{k})$ in ${\mathcal {U}}_{\delta }(x^{*})$ und konvergiert die Folge $(x^{k})$ linear gegen $x^{*}$. Die Konvergenz von $(x^{k})$ impliziert weiter die Konvergenz $\varepsilon (x^{k})\to 0$ $(k\to \infty )$. Da gemäß (5.28)

(5.29) $\left\|x^{k+1}-x^{*}\right\|\leq \varepsilon (x^{k})\left\|x^{k}-x^{*}\right\|$

für alle k gilt, folgt schließlich die superlineare Konvergenz von $(x^{k})$.

Gilt nun (5.25) auf ${\mathcal {U}}_{\delta }(x^{*})$, dann liegt für jedes $k$ mit $x^{*}$ und $x^{k}$ auch $x^{*}+s(x^{k}-x^{*})$ für alle $s\in [0,1]$ in ${\mathcal {U}}_{\delta }(x^{*})$ und folgt somit

$\left\|{\mathcal {J}}_{F}(x^{*}+s(x^{k}-x^{*}))-{\mathcal {J}}_{F}(x^{*})\right\|\leq L\left\|x^{k}-x^{*}\right\|.$

Aus (5.27) gewinnt man damit für alle $k$ die Abschätzung

$\varepsilon (x^{k})\leq 2\beta \left\{L+{\frac {1}{2}}L\right\}\left\|x^{k}-x^{*}\right\|=3\beta L\left\|x^{k}-x^{*}\right\|$

Letzteres zeigt zusammen mit (5.29) die quadratische Konvergenz der Folge $(x^{k})$.

q.e.d.

Beispiel 5.19

Gesucht sei die Lösung $x^{*}:=(x_{1}^{*},x_{2}^{*})^{T}$ der beiden Gleichungen

(5.30) ${\begin{matrix}F_{1}(x_{1},x_{2}):=x_{1}^{2}+x_{2}^{2}+0.6x_{2}-0.16=0,\\F_{2}(x_{1},x_{2}):=x_{1}^{2}-x_{2}^{2}+x_{1}-1.6x_{2}-0.14=0,\end{matrix}}$

für die $x_{1}^{*},x_{2}^{*}>0$ gilt, wobei wir hier keine Abbruchschranke $\varepsilon >0$ angeben. Die Jacobi-Matrix von $F$ lautet

${\mathcal {J}}_{F}(x)={\begin{pmatrix}2x_{1}&2x_{2}+0.6\\2x_{1}+1&-2x_{2}-1.6\end{pmatrix}}.$

Für $(x_{1}^{0},x_{2}^{0})^{T}:=(0.6,0.25)^{T}$ erhält man somit das lineare Gleichungssystem (5.24)

${\begin{pmatrix}1.2&1.1\\2.2&-2.1\end{pmatrix}}{\begin{pmatrix}h_{1}\\h_{2}\end{pmatrix}}=-{\begin{pmatrix}0.4125\\0.3575\end{pmatrix}}.$

Dieses besitzt die Lösung

${\begin{pmatrix}h_{1}^{0}\\h_{2}^{0}\end{pmatrix}}={\begin{pmatrix}-0.254\ 960\\-0.096\ 862\end{pmatrix}},$

so dass sich

${\begin{pmatrix}x_{1}^{1}\\x_{2}^{1}\end{pmatrix}}={\begin{pmatrix}0.60\\0.25\end{pmatrix}}+{\begin{pmatrix}-0.254\ 960\\-0.096\ 862\end{pmatrix}}={\begin{pmatrix}0.345\ 040\\0.153\ 138\end{pmatrix}}$

ergibt mit dem Defekt

${\sqrt {[F_{1}(x_{1}^{1},x_{2}^{1})]^{2}+[F_{2}(x_{1}^{1},x_{2}^{1})]^{2}}}=0.092\ 882\ 7.$

Mit $(x_{1}^{1},x_{2}^{1})^{T}$ verfährt man nun analog usw. Für die ersten vier Iterationen ergibt sich insgesamt die folgende Tabelle:

${\begin{array}{|c|c|c|c|c|c|}\hline k&x_{1}^{k}&x_{2}^{k}&\left\|F(x^{k})\right\|&h_{1}^{k}&h_{2}^{k}\\\hline 0&0.600000+0&0.250000+0&0.545859+0&-0.254960+0&-0.096862+0\\\hline 1&0.345040+0&0.153138+0&0.928827-1&-0.675094-1&-0.306747-1\\\hline 2&0.277531+0&0.122463+0&0.658124-2&-0.564594-2&-0.279860-2\\\hline 3&0.271885+0&0.119664+0&0.464212-4&-0.406023-4&-0.210055-4\\\hline 4&0.271845+0&0.119643+0&0.241346-8&&\\\hline \end{array}}$

Das Newton-Verfahren, Algorithmus 9, ist invariant gegenüber affin-linearen Transformationen (Übung!). Dies bedeutet, wenn $A\in \mathbb {R} ^{n\times n}$ eine beliebige reguläre Matrix und $c\in \mathbb {R} ^{n}$ irgendein Vektor: Ist $\{x^{k}\}$ die durch das lokale Newton-Verfahren für den Startpunkt $x^{0}$ erzeugte Iteriertenfolge zur Bestimmung einer Lösung des Gleichungssystems $F(x)=0$, so erzeugt das Verfahren bei Anwendung auf das System

$G(z):=F(Az+c)=0$

für den Startpunkt $z^{0}:=A^{-1}(x^{0}-c)$ die Iteriertenfolge $\{z^{k}\}$ mit

$z^{k}=A^{-1}\left(x^{k}-c\right)\Leftrightarrow x^{k}=Az^{k}+c.$

Verfahren, die invariant gegenüber affin-linearen Transformationen sind, gelten gegenüber Verfahren, die diese Eigenschaft nicht besitzen, insofern als robuster, als ihre Konvergenzgeschwindigkeit weit weniger von den gerade vorliegenden speziellen Daten abhängt. Anders als z. B. bei dem CG-Verfahren zur Lösung linearer Gleichungssysteme mit symmetrischer, positiv definiter Matrix (s. Kanzow) ändert sich beim lokalen Newton-Verfahren insbesondere durch eine (affin-)lineare Transformation der Variablen die Konvergenzgeschwindigkeit des Verfahrens nicht. Denn ist $\varepsilon >0$ eine vorgegebene Abbruchschranke, so gilt aufgrund der oben beschriebenen Invarianz gegenüber affin-linearen Transformationen und der sich daraus ergebenden Identitäten

$\left\|G(z^{k})\right\|=\left\|F(Az^{k}+c)\right\|=\left\|F(x^{k})\right\|$

die Äquivalenz

$\left\|F(x^{k})\right\|\leq \varepsilon \Leftrightarrow \left\|G(z^{k})\right\|\leq \varepsilon .$

Bei Verfahren, die wie die CG-Verfahren nicht invariant gegenüber affin-linearen Transformationen sind, kann man zwar möglicherweise die Konvergenzgeschwindigkeit durch eine geeignete Wahl der Matrix $A$ erheblich beschleunigen, ist es aber häufig nicht vorhersehbar, ob das Verfahren für die aktuellen Daten langsam konvergiert, oder ist es nicht klar, ob gegebenenfalls eine geeignete Transformation zur Konvergenzbeschleunigung gefunden werden kann. Mit

$p^{k}:=\left[{\mathcal {J}}_{F}(x^{k})\right]^{-1}F(x^{k})$

lautet die Iterationsvorschrift des Newton-Verfahrens

$x^{k+1}:=x^{k}+p^{k},\quad k=0,1,\ldots .$

Die Richtung $p^{k}$ bezeichnet man dabei auch als Newton-Richtung in $x^{k}$.

Es gibt eine große Zahl von Varianten des Newton-Verfahrens, die zum Ziel haben, den Konvergenzbereich des Verfahrens zu vergrößern und/oder seinen numerischen Aufwand zu reduzieren. So kann man das Newton-Verfahren in gewisser Weise globalisieren, indem man eine geeignete Schrittweite $t_{k}>0$ einführt und

$x^{k+1}:=x^{k}+t_{k}p^{k},\quad k=0,1,\ldots$

definiert. Dabei wählt man $t_{k}$ beispielsweise als (Näherungs-)Lösung des Problems

(5.31) $\min _{t\geq 0}\left\|F(x^{k}+tp^{k})\right\|_{2},$

da ja $\left\|F(x^{k+1})\right\|_{2}$ möglichst klein werden sollte. Von dem so modifizierten sog. gedämpften Newton-Verfahren kann man unter relativ schwachen Voraussetzungen für jeden Startpunkt $x^{0}$ einer geeigneten, hinreichend großen Menge Konvergenz zeigen. (Man hat sich dabei zu überlegen, dass ein solches $t_{k}$ existiert und eine positive Zahl ist. Da eine Lösung des globalen Optimierungsproblems (5.31) im Allgemeinen nicht realistisch ist, wählt man die Schrittweite $t_{k}$ häufig aber auf eine andere Weise; siehe z. B. Stoer, wo für eine solche andere Schrittweitenwahl auch ein Konvergenzsatz zu finden ist.)

Eine weitere praktisch relevante Modifikation des Newton-Verfahrens besteht darin, die numerisch aufwendig zu berechnende Jacobi-Matrix ${\mathcal {J}}_{F}(x^{k})$ im Verfahren durch eine geeignete Näherung $H_{k}\in \mathbb {R} ^{n\times n}$ zu ersetzen, wobei $H_{k}$ alleine aus den Daten $x^{k},x^{k-1},F(x^{k})$ und $F(x^{k-1})$ numerisch relativ günstig berechnet werden kann und somit insbesondere keine partiellen Ableitungen benötigt werden. Wir kennen ein solches Vorgehen schon vom Sekantenverfahren her, bei dem der Faktor

${\frac {x_{k}-x_{k-1}}{f(x_{k})-f(x_{k-1})}},$

der eine Näherung für $1/f'(x_{k})$ darstellt, in der Iterationsvorschrift vorkommt. Verfahren dieses Typs werden als Sekanten- oder Quasi-Newton-Verfahren bezeichnet. Das bekannteste ist das Broyden-Verfahren.

Quasi-Newton-Verfahren haben vor allem im Zusammenhang mit der Bestimmung von Extremalpunkten von $F$ und damit der Lösung des Systems $\nabla F(x)=0$ große Bedeutung, da man ja bei Anwendung des Newton-Verfahrens in einem solchen Fall in jeder Iteration die Hesse-Matrix $\nabla ^{2}F(x^{k})$ für $F$, also etwa $n^{2}/2$ partielle Ableitungen zweiter Ordnung zu berechnen hat. Von solchen Quasi-Newton-Verfahren gibt es eine Reihe von Varianten, die sich durch die Wahl der $H_{k}$ unterscheiden, wobei man im Fall der Lösung von Optimierungsaufgaben zusätzlich bestrebt ist, die positive oder negative (Semi-)Definitheit der Hesse-Matrix in einer Umgebung des Extremalpunktes auch für die $H_{k}$ zu erreichen. Die verbreitetste Methode ist das BFGS-Verfahren, das nach seinen Erfindern Broyden, Fletcher, Goldfarb und Shanno benannt wurde, die das Verfahren 1970 unabhängig voneinander vorschlugen. Es hat sich herausgestellt, dass dieses unter allen Quasi-Newton-Verfahren das wohl unempfindlichste gegenüber der Schrittweitenwahl ist. (Es gibt Alternativen zu der numerisch teuren Berechnung der Minimumschrittweite in (5.31).)

Von Quasi-Newton-Verfahren kann man keine quadratische Konvergenz erwarten, da sie ja weniger Information der Funktion als das Newton-Verfahren verwenden. Unter geeigneten Voraussetzungen lässt sich aber für sie im Allgemeinen superlineare Konvergenz nachweisen. Die schlechtere Konvergenzrate gegenüber dem Newton-Verfahren wird jedoch durch den pro Iteration erforderlichen, geringeren numerischen Aufwand kompensiert.

Allgemein kann man sagen, dass das Newton-Verfahren wohl das erfolgreichste und verbreitetste Verfahren der Mathematik ist. Es wurde im Hinblick auf die Lösung zahlloser Probleme, auch solcher in unendlich-dimensionalen Räumen, verallgemeinert und modifiziert.