Ing Mathematik: Näherungsweises Lösen gewöhnlicher Dgl.

Numerisches Differenzieren

Die findiff-Funktion

Python-Code:

import findiff
import numpy as np

x = np.linspace(0, 1, 100)

f = np.sin(x)  
d_dx = findiff.FinDiff(0, x, 1)
df_dx = d_dx(f) 

print(df_dx[-1])

Ausgabe:

0.5403208977697744

Differenzenquotient

• Vorwärts-Differenzenquotient: $y^{\prime }\approx \frac{\mathrm{\Delta }y}{\mathrm{\Delta }x}=\frac{y(x+h)-y(x)}{(x+h)-x}=\frac{y_{(k+1)}-y_{(k)}}{h}$
• Rückwärts-Differenzenquotient: $y^{\prime }\approx \frac{\mathrm{\Delta }y}{\mathrm{\Delta }x}=\frac{y(x)-y(x-h)}{x-(x-h)}=\frac{y_{(k)}-y_{(k-1)}}{h}$
• Zentraler Differenzenquotient: $y^{\prime }\approx \frac{\mathrm{\Delta }y}{\mathrm{\Delta }x}=\frac{y(x+h)-y(x-h)}{(x+h)-(x-h)}=\frac{y_{(k+1)}-y_{(k-1)}}{2\cdot h}$

Siehe auch Vorlage:W, Vorlage:W.

Satz von Taylor

Brook Taylor (englischer Mathematiker, 1685-1731)

Joseph-Louis Lagrange (italienisch-französischer Mathematiker, 1736-1813)

Taylorpolynom: $T_{n}f(x;a)=\sum _{k=0}^n\frac{f^{(k)}(a)}{k!}(x-a{)}^k$

Lagrangesche Restgliedformel: $R_{n}f(x;a)=\frac{f^{(n+1)}(\xi )}{(n+1)!}(x-a{)}^{n+1}$

$f(x)=T_{n}f(x;a)+R_{n}f(x;a)$

Nachfolgend seien Entwicklungen der ${\text{e}}^x$-Funktion als Taylorpolynome gezeigt. Zuerst die Linearisierung um den Arbeitspunkt und eine Parabel. Anschließend Polynome höheren Grades.

Zugehöriger Python-Code:

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import approximate_taylor_polynomial

x = np.arange(-10.0, 10.0, 0.01)
degree1 = 5
degree2 = 11

taylor1 = approximate_taylor_polynomial(np.exp, 0, degree1, 1, degree1+2)
taylor2 = approximate_taylor_polynomial(np.exp, 0, degree2, 1, degree2+2)

plt.plot(x, np.exp(x), label = r"$e^x$")
plt.plot(x, taylor1(x), label="Taylor, %d" % degree1)
plt.plot(x, taylor2(x), label="Taylor, %d" % degree2)
plt.plot(0, 1, "o", label = "Entwicklungspunkt")
plt.legend()
plt.tight_layout()
plt.grid()
plt.axis([-10, 10, -5, 5])
plt.show()

Siehe auch Vorlage:W, Knorrenschild: Seite 112ff

Diskretisierungsfehler

Problem: Wählt man $h$ zu klein, so ist ggf. der Rundungsfehler zu groß. Ist $h$ zu groß, dann ist offensischtlich der Diskretisierungsfehler hoch.

Herleitung: Mit der Taylorformel erhält man für $x=x_0+h;n=1$

$f(x_0+h)=f(x_0)+f^{\prime }(x_0)h+\frac{f^{″}(\xi )}{2}{h}^2$

$\underset{D_{1}f(x_0,h)}{\underbrace{\frac{f(x_0+h)-f(x_0)}{h}}}=f^{\prime }(x_0)+\frac{f^{″}(\xi )}{2}h$

Diskretisierungsfehler: $|D_{1}f(x_0,h)-f^{\prime }(x_0)|$

(Fehler)Ordnung $k$: $|D_{1}f(x_0,h)-f^{\prime }(x_0)|\le K{h}^k;$ mit $K>0$.

Je höher die Fehlerordnung, desto genauer ist i. A. die Formel. Die Schrittweite h ist optimal, wenn der Gesamtfehler aus Diskretisierungs- und Rundungsfehler minimal wird. Lt. Korrenschild: Seite 115 ist das Optimum bei $h_{opt}\approx 2\sqrt{\frac{|f(x_0)|}{|f^{″}(x_0)|}\cdot \text{eps}}$ erreicht. eps sei die Maschinengenauigkeit.

Siehe auch Vorlage:W.

Gewöhnliche Differentialgleichungen - Allgemeines

$F(x,y(x),y^{\prime }(x),\dots ,y^{(n)}(x))=0$ ist eine gewöhnliche Differenzialgleichung n-ter Ordnung. Es sind sämtliche Lösungen dieser Gleichung zu bestimmen. Häufig gelingt das nicht mehr analytisch, sondern die Gleichung muss numerisch bearbeitet werden.

Siehe auch Vorlage:W, Vorlage:W und Burg, Haf, Wille, Meister

Anfangswertprobleme

Allgemeines

Die Lösung eines Anfangswertproblems ist die Lösung der Differentialgleichung plus Berücksichtigung von vorgegebenen Anfangswerten.

$y^{\prime }=f(t,y)$ mit $y(0)=y_0,0\le t\le T$.

Wie nachfolgendes Bild zeigt, gibt es auch hier eine Vielzahl an numerischen Methoden (das Bild beinhaltet nicht alle verfügbaren Methoden, Quelle: Tabellen in Burg, Haf, Wille, Meister: Seite 75f).

Lipschitz-Bedingung

Rudolf Lipschitz (deutscher Mathematiker, 1832-1903)

$\Vert f(t,y)-f(t,z){\Vert }_2\le L\Vert y-z{\Vert }_2$ heißt lipschitzstetig, wenn eine Konstante (die Lipschitz-Konstante) $L\in {ℝ}_0$ existiert, die für alle $x_1,x_2\in ℝ$ gilt.

Siehe auch Vorlage:W, Burg, Haf, Wille, Meister: Seite 21, Hanke-Bourgeois: Seite 553

Picard-Lindelöf

Charles Émile Picard (französischer Mathematiker, 1856-1941

Ernst Leonard Lindelöf (finnischer Mathematiker, 1870-1946)

• Existenzproblem: Finden wir stets eine Lösung des Anfangswertproblems?
• Eindeutigkeit: Ist dies die einzige Lösung?
• Ist diese Lösung lokal oder lässt sie sich auf einen größeren Bereich fortsetzen?

Die dritte Frage lässt sich so beantworten: Existenzaussagen sind i.A. nur lokal gültig. Die ersten beiden beantwortet der Existenz- und Eindeutigkeitssatz von Picard-Lindelöf.

Die nachfolgende Variante des Satzes von Picard-Lindelöf stammt aus Wikipedia (Vorlage:W, Autoren siehe dort):

Sei $E$ ein Banachraum, $G\subset ℝ\times E$, $y_0\in E,R>0$ mit $[a,b]\times \bar{B}(y_0,R)\subset G$ und $f:G\to E$ stetig und lokal Lipschitz-stetig in der zweiten Variablen. Hierin bezeichnet

$\bar{B}(y_0,R):=\{z\in E|\Vert z-y_0\Vert \le R\}$

die abgeschlossene Kugel um $y_0$ mit Radius $R$. Ist

$M:=\max \{\Vert f(x,y)\Vert |(x,y)\in [a,b]\times \bar{B}(y_0,R)\}$

und

$\alpha :=\min \{b-a,\frac{R}{M}\},$

dann existiert genau eine Lösung des Anfangswertproblems

$y^{\prime }=f(x,y),y(a)=y_0$

auf dem Intervall $[a,a+\alpha ]$; sie hat Werte in $\bar{B}(y_0,R)$.

Eine etwas einfacher verständliche Variante findet sich in Burg, Haf, Wille, Meister: Seite 15ff. Auch in Hanke-Bourgeois: Seite 553ff findet sich dieser Satz.

Rekursionsformel nach Picard-Lindelöf:

Das Anfangswertproblem wird durch die Integralgleichung $y(x)=y_0+\underset{x_0}{\overset{x}{\int }}f(t,y(t))dt$ ersetzt.

$y_0(x)=y_0$

$y_i(x)=y_0+\underset{x_0}{\overset{x}{\int }}f(t,y_{n-1}(t))dt;i\in \mathbb{N}$.

Dabei muss gelten: $f,\frac{\partial f}{\partial y}$ stetig.

Mit dem Verfahren von Picard-Lindelöf ist zwar schon eine Methode zur numerischen Lösung von Anfangswertproblemen gegeben. In Burg, Haf, Wille, Meister: Seite 22ff ist ein einfaches Beispiel gezeigt. Allerdings lassen sich oft die Integrale nicht einfach explizit bestimmen. Es gibt daher eine Vielzahl anderer Methoden, die dieses Problem vermeiden.

Integrationsmethoden

$y(t_{i+m})-y(t_i)=\underset{t_i}{\overset{t_{i+m}}{\int }}f(t,y(t))dt$. Das Integral kann durch numerische Quadraturformeln gelöst werden. Mit der Rechteckregel ergibt sich:

$\underset{t_i}{\overset{t_{i+m}}{\int }}f(t,y(t))dt\approx \underset{t_i}{\overset{t_{i+m}}{\int }}r(t)dt=\underset{m\mathrm{\Delta }t}{\underbrace{(t_{i+m}-t_i)}}\underset{f(t_i,y(t_i))}{\underbrace{r(t_i)}}=m\mathrm{\Delta }tf(t_i,y(t_i))$

Differenziationsmethoden

$y^{\prime }(t_i)\approx \frac{y(t_{i+1})-y(t_i)}{\mathrm{\Delta }t}$

$\underset{f(t_i,y(t_i))}{\underbrace{y^{\prime }(t_i)}}\mathrm{\Delta }t\approx y(t_{i+1})-y(t_i)$

Einschrittverfahren

Zur Berechnung der Größe $y_{i+1}$ ist stets die Verfügbarkeit von $y_i$ ausreichend.

Mehrschrittverfahren

Zur Berechnung der Lösung $y_{i+m}$ sind neben $y_{i+m-1}$ mit $y_{i+m-2},\dots ,y_i$ noch weitere Stützpunkte zu berücksichtigen.

Explizite Verfahren

Weist die Verfahrensfunktion keine Abhängigkeit von $y_{i+m}$ auf, so heißt das Verfahren explizit.

Implizite Verfahren

Weist die Verfahrensfunktion eine Abhängigkeit von $y_{i+m}$ auf, so heißt das Verfahren implizit.

Steife Differenzialgleichungen

Bei diesen versagen explizite Einschrittverfahren (Quelle: Hanke-Bourgeois: Seite 587). Die Lösungskomponenten setzen sich aus verschieden stark exponentiell abklingenden Anteilen zusammen (Quelle: Bronstein, Semendjajew, Musil, Mühlig: Taschenbuch der Mathematik. 7. Aufl., Harri Deutsch, 2008, Seite 978). Die Lipschitzkonstante nimmt große Werte $L\gg 1$ an.

Siehe auch Vorlage:W. Im Englischen werden derartige DGL als stiff bezeichnet.

Ein PDF-File zum Thema Steife Differentialgeichungen findet sich auch auf [1].

Euler

Leonhard Euler (schweizer Mathematiker und Naturwissenschaftler, 1707-1783)

Explizit

Augustin-Louis Cauchy (französischer Mathematiker, 1789-1857)

$y(t_{i+m})-y(t_i)=\underset{t_i}{\overset{t_{i+m}}{\int }}f(t,y(t))dt$

$\underset{t_i}{\overset{t_{i+m}}{\int }}f(t,y(t))dt\approx \underset{t_i}{\overset{t_{i+m}}{\int }}r(t)dt=\underset{m\mathrm{\Delta }t}{\underbrace{(t_{i+m}-t_i)}}\underset{f(t_i,y(t_i))}{\underbrace{r(t_i)}}=m\mathrm{\Delta }tf(t_i,y(t_i))$

Mit $m=1$ kann dieses Verfahren aus obigen Formeln abgeleitet werden. Oder noch einfacher via nachfolgendem Bildchen.

$f(t_0,y_0)={y_0}^{\prime }=\frac{\mathrm{\Delta }y}{\mathrm{\Delta }t}=\frac{y_1-y_0}{\mathrm{\Delta }t}$

$y_1=y_0+\mathrm{\Delta }tf(t_0,y_0)$

oder allgemein

$y_{i+1}=y_i+\mathrm{\Delta }tf(t_i,y_i);i=0,\dots ,n-1$

Dieses Verfahren wird auch Euler-Cauchy-, Vorwärts-Euler-Verfahren oder eulersches Polygonzugverfahren genannt.

Siehe auch Vorlage:W, Burg, Haf, Wille, Meister: Seite 43ff, Knorrenschild: Seite 147ff oder Hanke-Bourgeois: Seite 557ff

Implizit

$y_{i+1}=y_i+\mathrm{\Delta }tf(t_{i+1},y_{i+1});i=0,\dots ,n-1$

Auch dieses Verfahren lässt sich aus dem obigen Bild herleiten. Diesmal wird aber nicht der Punkt $i$ betrachtet, sondern $i+1$. Es wird auch Rückwärts-Euler-Verfahren genannt.

Siehe auch Vorlage:W, Burg, Haf, Wille, Meister: Seite 45ff oder Hanke-Bourgeois: Seite 560ff

Runge

Carl David Tolmé Runge (deutscher Mathematiker, 1856-1927)

Es ist ein verbessertes Polygonzugverfahren. Es lässt sich auch mit der Mittelpunktregel herleiten.

$y_{i+1}=y_i+\mathrm{\Delta }tf(t_i+\frac{\mathrm{\Delta }t}{2},y_i+\frac{\mathrm{\Delta }t}{2}f(t_i,y_i));i=0,\dots ,n-1$

Herleitung:

$y{\text{'}}_{i+1}=\frac{\mathrm{\Delta }y}{\mathrm{\Delta }t}=\frac{y_{i+1}-y_i}{\mathrm{\Delta }t}=f(t_{i+1},y_{i+1})=f(t_i+\frac{\mathrm{\Delta }t}{2},y_i+\frac{\mathrm{\Delta }t}{2}f(t_i,y_i))$

$y_{i+1}=y_i+\mathrm{\Delta }tf(t_i+\frac{\mathrm{\Delta }t}{2},y_i+\frac{\mathrm{\Delta }t}{2}f(t_i,y_i))$

Der folgende Python-Code dient nur der Visualisierung und zu Testzwecken und ist nicht für den Produktiveinsatz gedacht. Dazu siehe das Kapitel Runge-Kutta und die vorgefertigte SciPy-Funktion solve_ivp.

Grafikausgabe des nachfolgenden Python-Codes:

import numpy as np
import matplotlib.pyplot as plt

def f(t,y):
    return(y)

def runge(dt, n, y0):
    t = np.zeros(n+2)
    y = np.zeros(n+2)
    t[0] = 0.
    y[0] = y0

    for i in np.arange(0, n+1):
        y[i+1] = y[i] + dt * f(t[i]+dt/2, y[i]+dt/2*f(t[i],y[i]))
        t[i+1] = t[i] + dt

    plt.plot(t,y, label="Runge-Näherung")
    return y,t

y, t = runge(1., 9, 1.)

print("Runge-Näherung = ", y[-1])

x = np.arange(0., 10.1, .1)
y = np.exp(x)

plt.plot(x,y, label="exaktes Ergebnis")
plt.grid()
plt.legend(loc="best")
plt.show()

print("exaktes Ergebnis = ", y[-1])

Konsolenausgabe:

Runge-Näherung =  9536.7431640625
exaktes Ergebnis =  22026.465794806718

Python-Code:

import numpy as np

def f(t,y):
    return(y)

def runge(dt, n, y0):
    t = np.zeros(n+2)
    y = np.zeros(n+2)
    t[0] = 0.
    y[0] = y0

    for i in np.arange(0, n+1):
        y[i+1] = y[i] + dt * f(t[i]+dt/2, y[i]+dt/2*f(t[i],y[i]))
        t[i+1] = t[i] + dt
    return y

y = runge(.01, 1000, 1.)
print(y[-1])

Ausgabe:

22244.151807018065

Die exakte Lösung wäre $y=e^x=e^{10}=22026.47$

Siehe auch Burg, Haf, Wille, Meister: Seite 53f.

Runge-Kutta

Martin Wilhelm Kutta (deutscher Mathematiker, 1867-1944)

Das obige Bild zeigt auch, dass das Runge-Kutta-Verfahren schon sehr genaue Werte liefert (verglichen mit den anderen dargestellten Verfahren).

Python-Programm (es soll die DGL $y^{\prime }=y$ mit dem Anfangswert $y(0)=1$ gelöst werden. Siehe dazu auch Burg, Haf, Wille, Meister: Seite22f. Dort wird auch die exakte Lösung mit vollständiger Induktion hergeleitet. Dies kann aber auch mit grundlegenden Mathematikkenntnissen selbst durchgeführt werden -> Separation der Variablen):

import scipy
import numpy as np

def dy_dx(x, y):
    return y

y0 = [1]
xi = [0, 10]

z = scipy.integrate.solve_ivp(dy_dx, xi, y0, method="RK45")
y = z.y 

print(y[-1, -1])

Ausgabe:

22016.12509092586

Die exakte Lösung wäre: $y=e^x=e^{10}=22026.47$

Knorrenschild gibt auf Seite 157 seines Buches eine Tabelle mit Zahlenwerten für verschiedene Verfahren an. Es soll die DGL $y^{\prime }=t^2+0.1\cdot y,y(-1.5)=y_0=0,t=[-1.5,-0.9,-0.3,0.3,0.9,1.5]$ u.a. mit dem Runge-Kutta-Verfahren gelöst werden. Nachfolgendes Python-Programm zeigt die Vorgehensweise. Die gelieferten Werte liegen sehr nahe an der Knorrenschild-Tabelle (und auch an den exakten Werten).

Python-Code:

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

def dy_dx(x, y):
    return x*x + 0.1*y

y0 = [0]
x = [-1.5, 1.5]

z = scipy.integrate.solve_ivp(dy_dx, x, y0, method="RK45" ,
                              t_eval=[-1.5,-0.9,-0.3,0.3,0.9,1.5])
y = z.y

print(y)

Ausgabe:

0.         0.91345018 1.21332982 1.30691808 1.62660618 2.63179452

Siehe auch Vorlage:W, Burg, Haf, Wille, Meister: Seite 55ff, Knorrenschild: Seite 157ff, Hanke-Bourgeois: Seite 565ff und [2].

Gedruckte Bücher (auszugsweise)

• Burg, Haf, Wille, Meister: Höhere Mathematik für Ingenieure, Band III: Gewöhnliche Differentialgleichungen, Distributionen, Integraltransformationen. 5. Auflage, Vieweg+Teubner, 2009, ISBN 978-3-8348-0565-2
• 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

Vorlage:Navigation zurückhochvor buch