Nichtlineare Gleichungssysteme - Numerische Mathematik (i282)

Wir werfen einen kurzen Blick auf das numerische Lösen nichtlinearer Gleichungssysteme

F(x)=0,\quad F=(f_1,\ldots,f_d):D\to\bbR^d,

(1)

wobei $D\subseteq\bbR^d$ der Definitionsbereich der vektorwertigen Funktion $F$ ist. Beschränken uns auf den Fall, dass die Anzahl der Unbekannten mit der Anzahl der Gleichungen übereinstimmt.

Gelegentlich treten nichtlineare Gleichungssysteme mit nur kleinem $d$ auf. Meist entstehen die Systeme jedoch als endlich Näherung kontinuierlicher Probleme (Strömungsmechanik usw.) und weisen dann sehr viele Unbekannte auf.

Ziel dieses Kapitels ist nicht die detaillierte behandlung konkreter Verfahren, sondern das Aufzeigen von zusätzlichen Fragen und Problemen, die im eindimensionalen Fall nicht auftraten.

5.1Bisektion vs. Dimension¶

Prinzipiell lässt sich das Bisektionsverfahren auf höhere Dimensionen übertragen, allerdings nicht so einfach wie es zunächst erscheinen mag; siehe z.B. A rapid Generalized Method of Bisection for solving Systems of Non-linear Equations. Hauptproblem ist die mit der Dimension $d$ exponentiell wachsende Anzahl zu untersuchender Punkte. Während im eindimensionalen Fall der Suchbereich durch zwei Punkte begrenzt ist, werden in $d$ Dimensionen $2^d$ Punkte benötigt. Dies entspricht gerade der Anzahl der Ecken eines $d$ -dimensionalen Quaders.

Beispiel: Für ein nichtlineares Gleichungssystem mit nur $d=30$ Unbekannten und Gleichungen wäre das Suchgebiet durch ca. 1 Milliarde Punkte begrenzt!

5.2Newton-Verfahren¶

Das Newton-Verfahren lässt sich leicht auf Gleichungssysteme übertragen. Das Newton-Verfahren und daraus abgeleitete Verfahren (siehe unten) sind das wichtigsten Verfahren für das Lösen auch großer nichtlinearer Gleichungssysteme.

5.2.1Verfahren¶

Wie bei $d=1$ lösen wir eine linearisierte Version des nichtlinearen Gleichungsystems, wobei wir die Linearisierungspunkte sukzessive ausgehend von einem Startpunkt $x_0$ als Lösung des vorhergehenden linearisierten Systems berechnen.

Für gegebenen Punkt $x_{k-1}$ liefert der Satz von Taylor bei zweimal stetig differenzierbarem $F$ die Näherung

F(x)\approx F(x_{k-1})+J_F(x_{k-1})\,(x-x_{k-1}),

(2)

wobei

J_F(x_{k-1}):=\begin{bmatrix} \frac{\partial f_1}{\partial x_1}(x_{k-1})&\cdots&\frac{\partial f_1}{\partial x_d}(x_{k-1})\\ \vdots&&\vdots\\ \frac{\partial f_d}{\partial x_1}(x_{k-1})&\cdots&\frac{\partial f_d}{\partial x_d}(x_{k-1}) \end{bmatrix}

(3)

die Matrix der ersten partiellen Ableitungen der Komponenten $f_1,\ldots,f_d$ von $F$ ist (“Jacobi-Matrix”). Das nichtlineare Gleichungssystem $F(x)=0$ wird also in einer (kleinen) Umgebung von $x_{k-1}$ näherungsweise durch das lineare Gleichungssystem

J_F(x_{k-1})\,x=J_F(x_{k-1})\,x_{k-1}-F(x_{k-1})

(4)

beschrieben. Ist $J_F(x_{k-1})$ invertierbar, so besitzt dieses die eindeutige Lösung

x_k:=x_{k-1}-J_F(x_{k-1})^{-1}\,F(x_{k-1}).

(5)

Dies ist die Iterationsvorschrift für das mehrdimensionale Newton-Verfahren.

5.2.2Probleme¶

Analog zum eindimensionalen Fall (aber technisch aufwendiger) kann man zeigen, dass das mehrdimensionale Newton-Verfahren quadratisch konvergiert, also

|x_k-x^\ast|\leq C\,|x_{k-1}-x^\ast|^2

(6)

gilt, sofern die Funktion $F$ gewisse, recht strenge Bedingungen erfüllt und $x_0$ nah genug an der Lösung $x^\ast$ liegt. Beachte, dass $x^\ast$ und alle $x_k$ hier Vektoren aus $\bbR^d$ sind, in der Ungleichung nun also Längen von Vektoren verglichen werden.

Wie man $x_0$ in der Praxis wählen muss um quadratische Konvergenz oder überhaupt Konvergenz zu erhalten, ist nicht so klar. Der Satz von Kantorovich liefert hier praktisch verwertbare Anhaltspunkte. Letztlich wird es aber meist auf “Versuch und Irrtum” hinauslaufen.

Nebendem auch im eindimensionalen Fall vorhandenen Konvergenzproblem macht zusätzlich die Jacobi-Matrix $J_F(x_{k-1})$ Probleme. Selbst wenn diese als auswertbare Formel vorliegt, müssen in jedem Interationsschritt die $d^2$ Einträge neu berechnet werden (eigentlich nur $\frac{n\,(n+1)}{2}$ , da symmetrisch). Oft ist $F$ so beschaffen, dass man $J_F$ gar nicht explizit in Formeln ausdrücken kann. Dann muss $J_F(x_{k-1})$ numerisch bestimmt werden. Die Probleme um $J_F$ haben vielfältige Lösungen hervorgeracht, die zu unterschiedlichen, sogenannten Newton-ähnlichen Verfahren geführt haben.

5.3Newton-ähnliche Verfahren¶

Bezüglich der Jacobi-Matrix $J_F$ im mehrdimensionalen Newton-Verfahren sind zwei Probleme zu lösen:

Ersetze $J_F$ durch eine Matrix mit möglichst vielen Nullen (weniger Rechenaufwand) aber annäherend gleicher Wirkung.
Bestimme $J_F$ numerisch stabil (!) mit möglichst geringem Aufwand.

Der zweite Punkt ist besonders schwierig, weil numerische Differentiation ein schlecht konditioniertes Problem ist (siehe Veranstaltung zum wissenschaftlichen Rechnen).

Lösungsansätze:

Beim vereinfachten Newton-Verfahren wird $J_F(x_{k-1})$ nicht bei jedem Schritt neu berechnet, sondern über einige Schritte konstant gehalten. Hintergrund ist, dass $J_F$ sich von Schritt zu Schritt ohnehin meist nur wenig ändert.
Quasi-Newton-Verfahren bestimmen eine Nährung für $J_F$ aus den zurückliegenden Iterationsschritten. Sie “sammeln” also im Laufe der Zeit Informationen über die Ableitungen von $F$ . Dies ist deutlich effizienter (aber technisch schwieriger) als numerische Differentiation ohne Nutzung der ohnehin bereits bekannten Funktionswerte.
Bei manchen Anwendungen, z.B. beim Lösen von partiellen Differentialgleichungen (vgl. Veranstaltung zum wissenschaftlichen Rechnen), hat $J_F$ eine gut verstandene Struktur, die den Einsatz schneller Löser für das lineare Gleichungssystem in jedem Iterationsschritt zulässt. Dies ist z.B. möglich, wenn man weiß, dass $J_F$ nur in der Nähe der Diagonale von Null verschieden ist.

5.4Ausflug: Newton-Fraktal¶

Um die Konvergenz bzw. Nichtkonvergenz des Newton-Verfahrens in Abhängigkeit vom Startpunkt zu visualisieren betrachten wir verschiedene Beispiele mit $d=2$ .

Wir generieren ein regelmäßiges Rechteckgitter aus Startpunkten und führen das Newton-Verfahren mit jedem dieser Punkte als Startpunkt aus. Jeder Lösung der Gleichung $F(x,y)=(0,0)$ wird eine Farbe zugeordnet. Die Startpunkte werden in der Farbe der Lösung gefärbt, zu der die Iteration führt. Je langsamer die Konvergenz (hohe Iterationsanzahl), desto dunkler die Farbe.

Startpunkte, die nicht zur Konvergenz führen, werden schwarz erscheinen. Zur leichteren Wahrnehmung dieser Punkte erstellen wir eine zweite Grafik, die diese Startpunkte in weiß auf schwarzem Grund darstellt. Grautöne entsprechen langsamer Konvergenz.

import numpy as np
import matplotlib.pyplot as plt
import sympy

def newton_fractal(f1, f2, x, y, roots):
    # f1, f2, x, y are SymPy objects; roots is a list of coordinates
    
    # canvas size
    xmin, xmax, nx = -2, 2, 1000
    ymin, ymax, ny = -2, 2, 1000

    # colors for zeros
    colors = np.array([[1, 0, 0], [0, 0, 1], [0, 1, 0], [1, 1, 0], [0, 1, 1], [1, 0, 1]])

    # Jacobian
    f1x = sympy.diff(f1, x)
    f1y = sympy.diff(f1, y)
    f2x = sympy.diff(f2, x)
    f2y = sympy.diff(f2, y)
    
    # Newton step
    x_next = x - 1 / (f1x * f2y - f2x * f1y) * (f2y * f1 - f1y * f2)
    y_next = y - 1 / (f1x * f2y - f2x * f1y) * (-f2x * f1 + f1x * f2)
    x_next = sympy.simplify(x_next)
    y_next = sympy.simplify(y_next)
    
    # make Python function from SymPy result
    newton_update = sympy.utilities.lambdify([x, y], (x_next, y_next))

    # grid of initial values for Newton iteration and arrays for storing iteration counts (color-coded)
    grid_x, grid_y = np.meshgrid(np.linspace(xmin, xmax, nx), np.linspace(ymin, ymax, ny))
    canvas = np.zeros((ny, nx, 3))
    iters = np.zeros((ny, nx), dtype=int)

    # Newton interation for each initial value
    delta = 1e-10
    max_iter = 100
    for i in range(ny):
        for j in range(nx):

            # Newton iteration
            x_old, y_old = grid_x[i, j], grid_y[i, j]
            for k in range(max_iter):
                x, y = newton_update(x_old, y_old)
                if (x - x_old) ** 2 + (y - y_old) ** 2 < delta ** 2:
                    iters[i, j] = k
                    break
                x_old, y_old = x, y
            else:
                iters[i, j] = max_iter

            # find closest root (for color selection)
            dists = [(x - root_x) ** 2 + (y - root_y) ** 2 for root_x, root_y in roots]
            idx_root = np.array(dists).argmin()

            # store color-coded iteration count (** 0.5 is for color-scaling)
            canvas[i, j, :] = 1 / iters[i, j] ** 0.5 * colors[idx_root, :]

    # show convergence to roots
    fig, ax = plt.subplots(figsize=(10,10))
    ax.imshow(canvas / canvas.max(), extent=(xmin, xmax, ymin, ymax))
    ax.plot(*zip(*roots), 'ow')
    ax.axis('equal')
    plt.show()

    # show fractal (initial values without convergence)
    fig, ax = plt.subplots(figsize=(10,10))
    ax.imshow(iters, extent=(xmin, xmax, ymin, ymax), cmap='gray')
    ax.plot(*zip(*roots), 'ow')
    ax.axis('equal')
    plt.show()

Schauen uns zunächst die Ergebnisse für die Funktion

F(x,y):=\begin{bmatrix}x^2-y^2-1\\2\,x\,y\end{bmatrix}

(7)

an. Diese Funktion hat zwei Nullstellen bei $(-1,0)$ und $(1, 0)$ (weiße Punkte).

x, y = sympy.symbols('x y', real=True)

#z = x + y * sympy.I
#f = z ** 2 - 1
#f1 = sympy.re(f)
#f2 = sympy.im(f)

f1 = x ** 2 - y ** 2 - 1
f2 = 2 * x * y

roots = [[1, 0], [-1, 0]]

newton_fractal(f1, f2, x, y, roots)

Das Newton-Verfahren konvergiert für alle Startpunkte, die nicht auf der vertikalen Koordinatenachse liegen. Je näher der Startpunkt an einer Nullstelle liegt, desto schneller die Konvergenz.

Erhöhen wir den Polynomgrad auf 3, z.B.

F(x,y):=\begin{bmatrix}x^3-3\,x\,y^2-1\\3\,x^2\,y-y^3\end{bmatrix},

(8)

so wird das Konvergenzverhalten deutlich komplexer. Die Funktion hat 3 Nullstellen (weiße Punkte).

x, y = sympy.symbols('x y', real=True)

#z = x + y * sympy.I
#f = z ** 3 - 1
#f1 = sympy.re(f)
#f2 = sympy.im(f)

f1 = x ** 3 - 3 * x * y ** 2 - 1
f2 = 3 * x ** 2 * y - y ** 3

roots = [[1, 0], [-0.5, 0.5 * np.sqrt(3)], [-0.5, - 0.5 * np.sqrt(3)]]

newton_fractal(f1, f2, x, y, roots)

Bei Start in der Nähe einer Nullstelle konvergiert das Newton-Verfahren zu dieser Nullstellen. Liegt der Startpunkt jedoch etwa gleich weit entfernt von zwei Nullstellen, so führen schon sehr kleine Änderungen am Startwert zu einem Wechsel der angenäherten Nullstelle.

Die Menge der Startpunkte von denen aus das Newton-Verfahren nicht konvergiert heißt Newton-Fraktal. Diese Menge hat eine äußerst komplexe sich sowohl in der Ausdehnung als auch im Detail immer wieder wiederholende Struktur. Das Newton-Fraktal ist weder Linie, noch Fläche, sondern eine Menge mit gebrochener Dimension zwischen 1 und 2 bei ca. 1.42. Siehe z.B. Fractal Basins of Attraction Associated with a Damped Newton’s Method für Details.

Zur besseren Visualisierung des Newton-Fraktals (Zoom usw.) existieren verschiedene Computerprogramme, z.B. XaoS (open-source für alle gängigen Plattformen).

Polynom 6. Grades:

x, y = sympy.symbols('x y', real=True)

#z = x + y * sympy.I
#f = z ** 6 - 1
#f1 = sympy.re(f)
#f2 = sympy.im(f)

f1 = x ** 6 - 15 * x ** 4 * y ** 2 + 15 * x ** 2 * y ** 4 - y ** 6 - 1
f2 = 6 * x ** 5 * y - 20 * x ** 3 * y ** 3 + 6 * x * y ** 5

roots = [
    [1, 0], [-0.5, 0.5 * np.sqrt(3)], [-0.5, - 0.5 * np.sqrt(3)],
    [-1, 0], [0.5, 0.5 * np.sqrt(3)], [0.5, - 0.5 * np.sqrt(3)]
]

newton_fractal(f1, f2, x, y, roots)

Für manche Funktionen gibt es ganze Flächen an Startpunkten, bei denen das Newton-Verfahren nicht konvergiert. Beispiel:

F(x,y):=\begin{bmatrix}x^3-3\,x\,y^2-2\,x+2\\3\,x^2\,y-y^3-2\,y\end{bmatrix}.

(9)

x, y = sympy.symbols('x y', real=True)

#z = x + y * sympy.I
#f = z ** 3 - 2 * z + 2
#f1 = sympy.re(f)
#f2 = sympy.im(f)

f1 = x ** 3 - 3 * x * y ** 2 - 2 * x + 2
f2 = 3 * x ** 2 * y - y ** 3 - 2 * y

roots = [[-1.7693, 0], [0.88465, -0.58974], [0.88465, 0.58974]]

newton_fractal(f1, f2, x, y, roots)

References¶

Vrahatis, M. N., & Iordanidis, K. I. (1986). A rapid Generalized Method of Bisection for solving Systems of Non-linear Equations. Numerische Mathematik, 49(2–3), 123–138. 10.1007/bf01389620

5 Nichtlineare Gleichungssysteme