4 Nichtlineare Gleichungen

Beschäftigen uns mit dem numerischen Lösen nichtlinearer Gleichungen

wobei der Definitionsbereich $[a,b]$ von $f:[a,b]\to\bbR$ so eingeschränkt sei, dass es genau eine Lösung gibt. Nichtlineare Gleichungen können immer auf diese Form gebracht werden, sodass das Lösen nichtlinearer Gleichungen äquivalent zum Finden von Nullstellen ist.

4.1Kondition¶

Formale Untersuchungen zur Kondition des Problems sind mit unseren Mittel nicht möglich. Die Eingabe ist hier eine Funktion. Die Ausgabe ist eine Nullstelle der Funktion. Wir müssten uns also mit Abständen in Funktionenräume beschäftigen um überhaupt über Eingabefehler reden zu können. Grundsätzlich gilt aber, dass die Kondition des Problems vor allem an der konkreten Gestalt der Funktion $f$ festgelegt wird.

Für quadratische Gleichungen haben wir die Kondition des Nullstellenfindens genauer untersucht. Allerdings haben wir dort nicht die gesamte quadratische Funktion als Eingabe betrachtet, sondern nur die in der Funktionsdefinition auftretenden Parameter $p$ und $q$. Somit konnten wir Eingabefehler mit unseren Mittel ausdrücken.

Eine andere Frage ist die nach der Kondition beim Auswerten der Funktion $f$ an einer Stelle $x$. Diese Frage können wir mit unseren Mitteln beantworten, sofern eine konkrete Funktion betrachtet wird. Für eine allgemeine Funktion $f$ sind keine Aussagen möglich. Da alle numerische Lösungsverfahren die Funktion $f$ an mehreren Stellen auswerten müssen, ist die Frage nach der Kondition dieser Auswertungen sehr relevant für die Stabilität von Lösungsalgorithmen.

4.2Iterative Verfahren¶

4.2.1Idee¶

Alle Lösungsverfahren für nichtlineare Gleichungen arbeiten iterativ, d.h. sie arbeiten nach folgendem Muster:

1. Wähle einen Startwert $x_0\in[a,b]$.

2. Für $k=1,2\ldots$ wiederhole:

   1. Berechne eine neue Näherungslösung $x_k$ aus $x_{k-1}$.

Im Prinzip kann $x_k$ nicht nur von $x_{k-1}$ sondern auch von mehreren Vorgängern abhängen, ist aber hier nicht Thema.

Die Folge $x_0,x_1,x_2,\ldots$ liefert immer bessere Näherungen an die tatsächliche Lösung, sofern das Verfahren sinnvoll konstruiert ist.

Um aus der obigen Verfahrensidee konkrete Algorithmen zu erhalten, müssen wir folgende Fragen beantworten:

• Wie wählen wir $x_0$?

• Wie berechnen wir $x_k$ aus $x_{k-1}$?

• Wann brechen wir die Iteration ab?

4.2.2Abbruchkriterium¶

Die Frage nach dem Abbruch der Iteration ist relevant, weil ein Computer nur endlich viele Rechenschritte durchführen kann. Die Antwort ist eine Abwägung zwischen verfügbarer Rechenzeit und akzeptiertem Lösungsfehler. Folgende Varianten sind üblich:

• Falls vor allem die Rechenzeit relevant ist, bricht man nach einer vorgegebenen Anzahl an Iterationen ab.

• Falls vor allem die Genauigkeit der Lösung relevant ist, bricht man erst ab, wenn die gewünschte Genauigkeit erreicht ist. Da man die exakte Lösung nicht kennt, ist die Machbarkeit dieses Ansatzes vom konkreten Verfahren abhängig.

• Falls vor allem die Genauigkeit der Lösung relevant ist, aber dise bei dem gewählten Verfahren nicht zuverlässig bestimmt werden kann, bricht man ab, wenn sich die Näherungslösungen im Rahmen der gewünschten Genauigkeit nicht mehr ändern.

• Praktisch wird man abbrechen, wenn die Genauigkeit hoch genug ist, und zusätzlich eine (hohe) maximale Iterationsanzahl vorgeben um auch in ungünstigen Fällen (langsame Konvergenz) die Rechenzeit in akzeptablem Rahmen zu halten.

4.2.3Konvergenzordnung¶

Aus der Beziehung $x^\ast=\lim\limits_{k\to\infty}x_k$ lassen sich keine Erkenntnisse zur Konvergenzgeschwindigkeit gewinnen. Deshalb ist man an konkreten Abschätzungen des Lösungsfehlers interessiert:

Aus beiden Abschätzungen folgt insbesondere $x^\ast=\lim\limits_{k\to\infty}x_k$, also Konvergenz.

Je höher die Konvergenzordnung, desto schneller die Konvergenz. Entsprechend sind bei hoher Konvergenzordnung weniger Iterationen nötig um einen vorgegebenen maximalen Lösungsfehler zu erreichen als bei niedriger Konvergenzordnung.

Abschätzungen der obigen Form sind jedoch selten direkt zu bekommen. Meist fügt man einen Zwischenschritt folgender Form ein:

wobei $p\geq 1$ und $C>0$. Daraus folgt dann Konvergenz mit Konvergenzordnung $p$, sofern folgende Bedingungen an die Konstante $C$ und den Startfehler $|x_0-x^\ast|$ erfüllt sind (IDVID 420):

• $p=1$: Gilt $C<1$, so folgt lineare Konvergenz mit $c=|x_0-x^\ast|$ und $\gamma=C$.

• $p>1$: Gilt

  $$|x_0-x^\ast|<\begin{cases} C^{-\frac{1}{p}},&\text{falls }C<1,\\ C^{-\frac{1}{p-1}},&\text{falls }C\geq 1, \end{cases}$$

  (6)

  so folgt Konvergenz mit Konvergenzordnung $p$, wobei $c=1$ und

  $$\gamma=\begin{cases} C^{\frac{1}{p}}\,|x_0-x^\ast|,&\text{falls }C<1,\\ C^{\frac{1}{p-1}}\,|x_0-x^\ast|,&\text{falls }C\geq 1,\\ \end{cases}$$

  (7)

  gilt.

Beispiel: Gilt für ein festes $k$, dass $|x_{k-1}-x^\ast|\approx 0.1$, so sinkt der Fehler in den nächsten Schritten bei linear Konvergenz ($p=1$) stets nur um den Faktor $C<1$. Bei quadratischer Konvergenz ($p=2$) hingegen beträgt der Fehler in den nächsten Schritte etwa 0.01, 0.0001, 0.00000001.

4.3Bisektionsverfahren¶

4.3.1Verfahren¶

Das Bisektionsverfahren (auch: Intervallschachtelung) ist ein einfaches Verfahren zur Nullstellenbestimmung.

Die Bedinung der Form $f(a)\,f(b)\leq 0$ sichert, dass entweder im Intervall $(a,b)$ eine Nullstelle liegt (bei strenger Ungleichheit) oder ein Intervallende eine Nullstelle ist (bei Gleichheit). Liegen mehrere Nullstellen im Startintervall $[a_0,b_0]$, so liefert das Verfahren nur eine dieser Nullstellen, wobei keine allgemeine Aussagen darüber möglich sind, welche Nullstelle geliefert wird.

4.3.2Konvergenz und Abbruchkriterium¶

Ist der gewünschte maximale (absolute) Lösungsfehler $\Delta>0$ vorgegeben, so kann wegen

die Iteration gestoppt werden, sobald $b_k-a_k\leq 2\,\Delta$. Wegen

ist die Iterationsanzahl dann

insbesondere ist sie im Vorfeld exakt bekannt. Ein zusätzliche Begrenzung der Iterationsanzahl als weiteres Abbruchkriterium ist nicht nötig.

4.3.3Algorithmus¶

Ein konkreter Algorithmus zur Umsetzung des Bisektionsverfahrens kann wie folgt aussehen (ohne explizites Erwähnen des Rundens nach jeder Rechenoperation):

Eine vollständige Stabilitätsanalyse ist trotz des relativ einfachen Algorithmus recht mühsam. Das Verfahren gilt als sehr stabil. Wenn die Funktion $f$ also fehlerfrei auswertbar ist (fehlerfreie Eingabe), dann liegt der Ausgabefehler in der Größenordnung des durch übliche fehlerbehaftete Eingaben auch bei exakter Rechnung nicht zu vermeidenden Ausgabefehlers (Kondition!), sofern $\Delta$ im Abbruchkriterium entsprechend klein gewählt wird. Einige Hinweise zur Einschätzung der Stabilität:

• Die Summen $a+b$ beim Berechnen der Intervallmitten sind anfällig für Auslöschung falls $0\in(a,b)$. Der relative Fehler der Intervallmitten kann dann recht groß sein. Die Werte $x_k$ werden also nicht exakt in der Mitte liegen, was aber aufgrund der Arbeitsweise des Bisektionsverfahrens keinen Einfluss auf das Endergebnis hat.

• Rundungsfehler bei der Berechnung der Intervallmitten haben ebenfalls keinen Einfluss auf das Endergebnis. Insbesondere werden sie nicht von Schritt zu Schritt verstärkt, da stets eine neue (im Rahmen der Stabilitätsanalyse exakte) Auswertung von $f$ stattfindet.

• Auslöschungseffekte bei $b-a$ im Abbruchkriterium können zu leichten Abweichungen bei der theoretisch ermittelten Iterationsanzahl führen. Da die Fehlerschranke $\Delta$ in der Praxis ohnehin nur eine grobe Vorgabe/Schätzung ist, spielt dieses Problem in der Praxis keine Rolle.

• Formuliert man das Abbruchkriterium bzgl. des relativen Fehlers, also z.B.

  $$b-a\leq\Delta_{\mathrm{rel}}\,\frac{|a+b|}{2},$$

  (13)

  so können zusätzlich Auslöschungseffekte bei $a+b$ auftreten, wenn $0\in(a,b)$. Die damit eventuell verbundenen Abweichungen bei den Iterationsanzahlen sind wiederum praktisch irrelevant.

4.4Newton-Verfahren¶

4.4.1Verfahren¶

Das Newton-Verfahren nutzt zusätzlich zu den Funktionswerten von $f$ auch die erste Ableitung $f'$ um Informationen über die Lage einer Nullstelle zu erhalten.

Dieses Verfahren wählt $x_k$ also Nullstelle der Tangent an $f$ an der Stelle $x_{k-1}$ (IDVID 450).

4.4.2Konvergenz und Abbruchkriterium¶

Aus dieser Abschätzung folgt quadratische Konvergenz $|x_k-x^\ast|\leq c\,\gamma^{(2^k)}$ mit $c=1$ und

wobei

gelten muss. Damit das Newton-Verfahren konvergiert müssen also insbesondere zwei Bedingungen erfüllt sein:

• Der Startwert $x_0$ liegt nah genug bei der (unbekannten) Lösung $x^\ast$.

• Die Iterierten $x_k$ müssen in der Nähe von $x^\ast$ bleiben.

Sind diese Bedingungen erfüllt, so konvergiert das Newton-Verfahren quadratisch. Man spricht auch von lokaler quadratischer Konvergenz. Im Gegensatz dazu ist das Bijektionsverfahren global konvergent, allerdings nur linear.

Da die Konstanten $c_1$ und $c_2$ praktisch nicht zugänglich sind, kann aus der Abschätzung für die Konvergenzordnung kein in der Praxis einsetzbares Abbruchkriterium gewonnen werden. Man behilft sich mit der Tatsache, dass bei Konvergenz der Abstand $|x_k-x_{k-1}|$ zwischen zwei Iterierten gegen Null geht. Abgebrochen wird also, wenn dieser Abstand unter eine vorgegebene Schranke fällt.

4.4.3Algorithmus¶

Ein konkrete Umsetzung des Newton-Verfahrens kann wie folgt aussehen (ohne explizites Erwähnen des Rundens nach jeder Rechenoperation):

Dieser Algorithmus ist stabil. Insbesondere werden Rundungsfehler im Laufe der Iterationen nicht verstärkt. Jede neue Iteration kann als Neustart des Algorithmus mit anderem Startwert betrachtet werden. Da die Wahl des Startwertes (innerhalb gewisser Grenzen) keinen Einfluss auf die Lösung hat, spielen Rundungsfehler des Startwertes keine Rolle.

4.4.4Beispiele¶

Um das Problem der nicht globalen Konvergenz zu verdeutlichen führen wir das Newton-Verfahren für viele verschiedene Startwerte $x_0$ aus und stellen die Konvergenz bzw. Nichtkonvergenz grafisch dar. Dabei sehen wir eine Näherungsfolge $x_0,x_1,x_2,\ldots$ als konvergent an, wenn innerhalb von 1000 Iterationen der Abstand $|x_k-x_{k-1}|$ unter 10^-10 sinkt. Zusätzlich stellen wir dar, gegen welche von ggf. mehreren Nullstellen die Folge konvergiert.

In den folgenden Diagrammen sind neben der betrachteten Funktion $f$ die Iterationsanzahlen bis zum erstmaligen Unterschreiten der Fehlerschranke und (farblich) die angenäherten Nullstellen ausgehend von verschiedenen Startwerten $x_0$ dargestellt.

import numpy as np
import matplotlib.pyplot as plt

def visualize_newton(f, f1, zeros, a, b, n=1000):

    # grid of different initial values
    x0s = np.linspace(a, b, n)
    
    # Newton method for each initial value
    delta = 1e-10
    max_iter = 1000
    limits = []
    iters = []
    for x0 in x0s:
        x_old = x0
        for k in range(max_iter):
            x = x_old - f(x_old) / f1(x_old)
            if np.abs(x - x_old) < delta:
                limits.append(round(x))
                iters.append(k)
                break
            x_old = x
        else:
            limits.append(np.nan)
            iters.append(max_iter)
    
    # convert limits to colors
    colors = ['#ff0000', '#40c000', '#0000ff']
    no_conv_color = '#808080'  # no convergence
    if len(zeros) > len(colors):
        print('More zeros than colors! Ignoring some zeros.')
    zero_to_color = dict(zip(zeros, colors))
    limits_color = [zero_to_color.get(limit, no_conv_color) for limit in limits]
    
    # plot
    fig, ax1 = plt.subplots(figsize=(10, 5))
    ax2 = ax1.twinx()
    ax1.bar(x0s, iters, color=limits_color, width=(b - a) / n)
    ax2.plot(x0s, f(x0s), '-k')
    for zero in zeros:
        ax2.plot(zero, f(zero), 'o', color=zero_to_color[int(zero)], markersize=10)
    ax1.set_xlim(a, b)
    ax1.set_ylim(0, max([i for i in iters if i < max_iter]))
    ax1.grid(axis='x')
    ax2.grid(axis='y')
    ax1.set_xlabel('$x$')
    ax1.set_ylabel('Iterationen')
    ax2.set_ylabel('$f$')
    ax1.set_title('Einzugsbereiche der Nullstellen und Konvergenzgeschwindigkeit')
    plt.show()

4.4.4.1Polynom 3. Grades¶

Konvergenzverhalten für das Polynom

visualize_newton(
    lambda x: x ** 3 + 2 * x ** 2 - x - 2,
    lambda x: 3 * x ** 2 + 4 * x - 1,
    [-2, -1, 1],
    -2.5, 1.3
)

An den Extremstellen gilt $f'(x)=0$, sodass das Newton-Verfahren dort nicht anwendbar ist (alle Berechnungen liefern nan, sodass die Abbruchbedingung nie erfüllt ist). Für $x_0$, für die weitere Iterierte auf eine Extremstelle fallen, ist das Newton-Verfahren ebenfalls nicht anwendbar. Ist $x_0$ in der Nähe einer Extremstelle, so führen schon sehr kleine Änderungen an $x_0$ zum Wechsel der angestrebten Nullstelle.

4.4.4.2Sinus¶

Konvergenzverhalten für die Sinusfunktion

visualize_newton(
    lambda x: np.sin(np.pi * x),
    lambda x: np.pi * np.cos(np.pi * x),
    [1, 2, 3],
    0.5, 3.5
)

Liegt $x_0$ in der Nähe einer Extremstelle, so verläuft die Tangente sehr schwach. Entsprechend weit entfernt ist $x_1$, sodass auch Startwerte innerhalb des dargestellten Intervalls zu Näherungen von weit entfernt liegenden Nullstellen führen können. Liegt die angestrebte Nullstelle nicht im dargestellten Intervall, so wird der entsprechende Startwert grau eingefärbt (gleiche Farbe wie bei Nicht-Konvergenz).