Ing Mathematik: Numerisches Lösen nichtlinearer Gleichungen

Nullstellenprobleme

Nullstellen (oder Wurzeln) einer Funktion sind diejenigen Punkte, welche die x-Achse berühren oder schneiden (y=0). Wichtig sind insbesondere die Fragen

• existieren Nullstellen,
• in welchem Bereich liegen sie und
• wieviele Nullstellen gibt es?

Zwischenwertsatz

Bernard Bolzano (böhmischer Mathematiker, 1781-1848)

Eine stetige Funktion $f:[a,b]\to ℝ$ nimmt jeden Wert $s$ zwischen $f(a)$ und $f(b)$ an. Mit $s=0$ ist der Zwischenwertsatz für die Lösung von Nullstellenproblemen wichtig. Man nennt ihn dann auch Nullstellensatz (von Bolzano). Er besagt dann: Haben $f(a)$ und $f(b)$ verschiedene Vorzeichen, so existiert mindestens eine Nullstelle von $f$ im Intervall $(a,b)$.

Siehe zu diesem Thema insbesondere Vorlage:W.

Bisektion

Das Bisektions- oder Intervallhalbierungsverfahren kann benutzt werden, um Nullstellen grob zu bestimmen. Es ist bei Weitem nicht das schnellste Verfahren um Nullstellen genau zu bestimmen. Allerdings konvergieren die schnelleren Verfahren nicht immer und das Bisektionsverfahren kann hier ggf. Hilfestellung leisten. Das Bisektionsverfahren läuft wie folgt ab:

Sei $f:[a,b]\to ℝ$ mit $f(a)\cdot f(b)<0$ stetig.

$[a_0,b_0]=[a,b]$

$[a_{i+1},b_{i+1}]=\{\begin{array}{cc}\left[\begin{array}{c}{a}_i,\frac{a_i+b_i}{2}\end{array}\right]& \text{wenn}f\left(\begin{array}{c}\frac{a_i+b_i}{2}\end{array}\right)\cdot f(a_i)\le 0\\ \left[\begin{array}{c}\frac{a_i+b_i}{2},b_i\end{array}\right]& \text{sonst}\end{array}$

In jedem der erzeugten Intervalle gibt es eine Nullstelle von $f$.

Natürlich kann man auch das Bisektionsverfahren durch ein Python-Programm realisieren:

import numpy as np

def funktion(x):
    y = x**3-x+0.3
    return y

def bisektion(a, b, maxit):
    ai = a
    bi = b

    for i in np.arange(0, maxit):
        if funktion((ai+bi)/2.) * funktion(ai) <= 0.:
            bi = (ai+bi) / 2.
        else:
            ai = (ai+bi) / 2.
        print(ai, bi)

bisektion(0., 0.5, 10)

Ausgabe:

0.25 0.5
0.25 0.375
0.3125 0.375
0.3125 0.34375
...
0.3388671875 0.33935546875

Die Zahlenwerte für das Polynom in funktion und das Startintervall wurden aus einem Beispiel in Knorrenschild: Seite 27 entnommen. Die anderen Nullstellen lassen sich genauso bestimmen. Es muss nur das Startintervall entsprechend gesetzt werden.

Auch SciPy bietet eine vorgefertigte Funktion namens scipy.optimize.bisect:

from scipy.optimize import bisect

def funktion(x):
    y = x**3-x+0.3
    return y

print (bisect(funktion, 0., 0.5))

Ausgabe:

0.33893624159645697

Die Nullstellen des Beispiels und den Graphen allgemein zeigt folgender Plot:

Fixpunktiteration

Gleichungen der Form $x=f(x)$ (Fixpunktgleichungen) sollen gelöst werden. In diese Form muss die Gleichung umgeformt werden. Anschließend wird eine Startnäherung $x_0$ gewählt und $x_1=f(x_0)$ berechnet und weiter iterativ verfeinert $x_2=f(x_1)$ etc. Unter bestimmten Voraussetzungen nähert sich die Folge $x_0,x_1,x_2\dots$ einer Lösung des ursprünglichen Problems an.

Für weitere Details siehe z.B. Vorlage:W.

Beispiel:

$\log x+x-10=0$

Fixpunktgleichung:

$x_{n+1}=10-\log x_n$

Startwert:

$x_0=1$

Iterationen:

$x_1=10$

$x_2=9$

$x_3=9.046$

$x_4=9.0436$

$x_5=9.0437$

$\dots$

Sehr einfach lässt sich das auch mittels Python-Programm rechnen:

import numpy as np

x_i = 1. 

for i in np.arange(0, 10):
    x_i = 10 - np.log10(x_i)
    print(x_i)

Ausgabe:

10.0
9.0
9.045757490560675
...
9.043655967157635

Banachscher Fixpunktsatz

Stefan Banach (polnischer Mathematiker, 1892-1945)

Mit dem banachschen Fixpunktsatz lässt sich die Konvergenz von iterativen Verfahren (wie z.B. dem Newton-Verfahren) zeigen. S. Banach veröffentlichte diesen Satz im Jahr 1922.

Eine Abbildung $f$ heißt kontrahierend, wenn $\Vert f(x_1)-f(x_2)\Vert \le K\Vert x_1-x_2\Vert$ mit $0<K<1$. $K$ heißt auch Kontraktionszahl von $f$. Daraus folgt: $f$ hat genau einen Fixpunkt. Die Iterationsfolge $(x_n)$ konvergiert gegen diesen Fixpunkt. Der Startwert sei $x_0\in ℝ$.

Es gelten die

• a-priori-Fehlerabschätzung: $\Vert x_n-x\Vert \le \frac{K^n}{1-K}\Vert x_1-x_0\Vert$

und die

• a-posteriori-Fehlerabschätzung: $\Vert x_n-x\Vert \le \frac{K}{1-K}\Vert x_n-x_{n-1}\Vert$

mit $n=0,1,\dots$

Siehe auch Vorlage:W. Für die Beweise siehe z.B. Meister: Seite 30ff oder auch Burg, Haf, Wille, Meister: Bd. I, Seite 72ff, Burg, Haf, Wille, Meister: Bd. VI, Seite 20f.

Newton-Verfahren

Isaac Newton (englischer Mathematiker und Naturwissenschaftler, 1643-1727)

• 
• 

Steigung der Tangente: $\tan {\alpha }_0=f^{\prime }(x_0)$

Mit der Tangente wird ein Dreieck aufgespannt:

$\tan {\alpha }_0=f^{\prime }(x_0)=\frac{f(x_0)}{x_0-x_1}$

$x_0-x_1=\frac{f(x_0)}{f^{\prime }(x_0)}$

$x_1=x_0-\frac{f(x_0)}{f^{\prime }(x_0)}$

oder allgemein

$x_{n+1}=x_n-\frac{f(x_n)}{f^{\prime }(x_n)};n=0,1,\dots$

Obiges Beispiel mit der von SciPy bereitgestellten Funktion newton liefert wieder das gewünschte Ergebnis. Der Solver verwendet je nach gelieferten Informationen das Newton-Raphson-Verfahren (kurz: Newton-Verfahren) oder das später vorgestellte Sekantenverfahren oder das Halley-Verfahren. Details siehe auch [1]. Zum Halley-Verfahren siehe Vorlage:W. Es benötigt auch die zweite Ableitung. Wird nur die erste Ableitung angegeben, so kommt das Newton-Verfahren zum Einsatz. Wird keine Ableitung angegeben, so wird das Sekantenverfahren verwendet.

from scipy.optimize import newton

def funktion(x):
    y = x**3 - x + 0.3
    return y

print (newton(funktion, 0.))

Ausgabe:

0.3389362415949407

Siehe auch Vorlage:W.

Heron-Verfahren

Zur Berechnung von $\sqrt{a}$ kann das Heron-Verfahren Anwendung finden.

$f(x)=x^2-a$

$f^{\prime }(x)=2x$

Newton: $x_{i+1}=x_i-\frac{x_i^2-a}{2x_i}$

Iteration für $a=5$ und $x_0=1$:

$x_1=1-\frac{1^2-5}{2\cdot 1}=3.0$

$x_2=3-\frac{3^2-5}{2\cdot 3}=2.3333$

$x_3=2.3333-\frac{2.3333^2-5}{2\cdot 2.3333}=2.2381$

$\dots$

Python-Code: Es soll $\sqrt{5}$ berechnet werden, $x_0=1$

import numpy as np

def funktion(x, a):
    y = x**2 - a
    return y

def funktion_diff(x, a):
    y = 2*x
    return y

def heron(xstart, a, maxit):
    x_i = xstart

    for i in np.arange(0, maxit):
        x_i = x_i - funktion(x_i, a) / funktion_diff(x_i, a)
        print(x_i)

heron(1.0, 5.0, 10)

Ausgabe:

3.0
2.3333333333333335
2.238095238095238
...
2.23606797749979
2.23606797749979

Zum Prinzip siehe auch Vorlage:W.

Sekantenverfahren

Iterationsvorschrift: $x_{n+1}=x_n-\underset{\approx 1/f^{\prime }(x_n)}{\underbrace{\frac{x_n-x_{n-1}}{f(x_n)-f(x_{n-1})}}}\cdot f(x_n)$

Startwerte: $x_0,x_1$.

Zur automatischen Berechnung kann einerseits wiederum die SciPy-Funktion newton verwendet werden. Oder man programmiert das Verfahren von Grund auf selbst. Dies sei hier aber nicht gezeigt. Es ist sehr ähnlich zu den Ausführungen beim Newton- oder Heron-Verfahren.

Siehe auch Vorlage:W.

Regula falsi

Die Regula falsi kombiniert das Sekanten- mit dem Bisektionsverfahren. Näheres siehe Vorlage:W.

Konvergenz

Iterationsverfahren konvergieren i.A. nicht mit jedem Startwert gegen die gesuchte Lösung $\tilde{x}$. Man unterscheidet zwischen anziehenden ($F^{\prime }(\tilde{x})<1$) und abstoßenden Fixpunkten ($F^{\prime }(\tilde{x})>1$, konvergiert nicht). Vereinfacht gesagt: Wenn der Startwert nahe genug bei der Lösung liegt, so ist der Fixpunkt anziehend.

Die Konvergenzgeschwindigkeit oder Konvergenzordnung $q\ge 1$: $\Vert x_{n+1}-\tilde{x}\Vert \le c\cdot \Vert x_n-\tilde{x}{\Vert }^q$ für $c\in {ℝ}^{+}$ und $n\in \mathbb{N}$. Für $q=1$ muss $c<1$ sein.

Ist

• $q=1$, so liegt lineare Konvergenz vor
• $q=2$, es liegt quadratische Konvergenz vor.

Das Newton-Verfahren konvergiert quadratisch, ausgenommen es liegt eine mehrfache Nullstelle vor. Dann konvergiert es nur noch linear. Das Sekantenverfahren: konvergiert mit $q=1.618$ (Quelle: Knorrenschild: Seite 37).

Nichtlineare Gleichungssysteme - Newton für Systeme

Es soll folgendes Gleichungssystem gelöst werden:

$f(x)=f(x_1,\dots ,x_n)=\left(\begin{array}{c}{f}_1(x_1,\dots ,x_n)\\ ⋮\\ f_n(x_1,\dots ,x_n)\end{array}\right)=\left(\begin{array}{c}0\\ ⋮\\ 0\end{array}\right)=0$

Newton: $x^{(n+1)}=x^{(n)}-(f^{\prime }(x^{(n)}){)}^{-1}f(x^{(n)})$

mit der Jacobi-Matrix

$f^{\prime }(x)={\left(\begin{array}{c}\frac{\partial f_i}{\partial x_j}\end{array}\right)}_{ij}\in {ℝ}^{𝕟\mathbb{\times }𝕟}=\left(\begin{array}{cccc}\frac{\partial f_1}{\partial x_1}& \frac{\partial f_1}{\partial x_2}& \dots & \frac{\partial f_1}{\partial x_n}\\ \frac{\partial f_2}{\partial x_1}& \frac{\partial f_2}{\partial x_2}& \dots & \frac{\partial f_2}{\partial x_n}\\ ⋮& ⋮& \ddots & ⋮\\ \frac{\partial f_n}{\partial x_1}& \frac{\partial f_n}{\partial x_2}& \dots & \frac{\partial f_n}{\partial x_n}\end{array}\right).$

Die Jacobi-Matrix muss invertierbar sein. Bei der konkreten Anwendung wird aber nicht die Inverse der Jacobi-Matrix berechnet, sondern die Lösung in Form eines linearen GS berechnet.

$f{\text{'}}^{(n)}{\delta }^{(n)}=-f(x^{(n)})$

mit

$x^{(n+1)}=x^{(n)}+{\delta }^{(n)}$

Mit dem gedämpften Newton-Verfahren kann eine bessere Konvergenz erzielt werden (Schrittweitendämpfung).

Mit dem vereinfachtes Newton-Verfahren erspart man sich die laufende Auswertung der Jacobi-Matrix: $f{\text{'}}^{(0)}{\delta }^{(n)}=-f(x^{(n)})$. Die Konvergenzgeschwindigkeit ist aber schlechter (nur noch linear anstatt quadratisch).

Nachfolgend ein Python-Beispiel, das von den Funktionen her an Hanke-Bourgeois: Seite 174f angelehnt ist. Es wird der Einfachheit halber die SciPy-Funktion fsolve verwendet.

import numpy as np

from scipy.optimize import fsolve

def func(x):
    return [x[0] - 0.25*np.cos(x[0]) + 0.25*np.sin(x[1]),
            x[1] - 0.25*np.cos(x[0]) + 0.5*np.sin(x[1])]

root = fsolve(func,  [0, 0])
print(root)

Ausgabe:

[0.20412903 0.16344858]

Gedruckte Werke (auszugsweise)

• Burg, Haf, Wille, Meister: Höhere Mathematik für Ingenieure, Band I: Analysis. 9. Auflage, Springer Vieweg, 2011, ISBN 978-3-8348-1218-6
• Burg, Haf, Wille, Meister: Partielle Differentialgleichungen und funktionalanalytische Grundlagen, Höhere Mathematik für Ingenieure, Naturwissenschaftler und Mathematiker (Band VI). 5. Auflage, Vieweg+Teubner, 2010, ISBN 978-3-8348-1294-0
• Hanke-Bourgeois: Grundlagen der Numerischen Mathematik und des Wissenschaftlichen Rechnens. 3. Aufl., Vieweg+Teubner, 2009, ISBN 978-3-8348-0708-3
• Knorrenschild: Numerische Mathematik, Eine beispielorientierte Einführung. 6. Aufl., Hanser, 2017, ISBN 978-3-446-45161-2
• Meister: Numerik linearer Gleichungssysteme. 4. Aufl., Vieweg+Teubner, 2011, ISBN 978-3-8348-1550-7