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

Zurück zu Numerische Berechnung bestimmter Integrale |

Hoch zu Gesamtinhaltsverzeichnis |

Vor zu FFT und Wavelets

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'\approx {\frac {\Delta y}{\Delta x}}={\frac {y(x+h)-y(x)}{(x+h)-x}}={\frac {y_{(k+1)}-y_{(k)}}{h}}$
Rückwärts-Differenzenquotient: $y'\approx {\frac {\Delta y}{\Delta x}}={\frac {y(x)-y(x-h)}{x-(x-h)}}={\frac {y_{(k)}-y_{(k-1)}}{h}}$
Zentraler Differenzenquotient: $y'\approx {\frac {\Delta y}{\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 Numerische Differentiation, Differenzenquotient.

Satz von Taylor

Taylorpolynom: $T_{n}f(x;a)=\sum _{k=0}^{n}{f^{(k)}(a) \over 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 Taylor-Formel, 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'(x_{0})h+{\frac {f''(\xi )}{2}}h^{2}$

$\underbrace {\frac {f(x_{0}+h)-f(x_{0})}{h}} _{D_{1}f(x_{0},h)}=f'(x_{0})+{\frac {f''(\xi )}{2}}h$

Diskretisierungsfehler: $|D_{1}f(x_{0},h)-f'(x_{0})|$

(Fehler)Ordnung $k$ : $|D_{1}f(x_{0},h)-f'(x_{0})|\leq Kh^{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 Lokaler Diskretisierungsfehler.

Gewöhnliche Differentialgleichungen - Allgemeines

$F\left(x,y(x),y'(x),\dots ,y^{(n)}(x)\right)=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 Differentialgleichung, Gewöhnliche Differentialgleichung 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'=f(t,y)$ mit $\quad y(0)=y_{0},\quad 0\leq t\leq 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)

$\|f(t,y)-f(t,z)\|_{2}\leq L\|y-z\|_{2}$ heißt lipschitzstetig, wenn eine Konstante (die Lipschitz-Konstante) $L\in \mathbb {R} _{0}$ existiert, die für alle $x_{1},x_{2}\in \mathbb {R}$ gilt.

Siehe auch Lipschitz-stetige Funktion, Burg, Haf, Wille, Meister: Seite 21, Hanke-Bourgeois: Seite 553

Picard-Lindelöf

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 ( Satz von Picard-Lindelöf, Autoren siehe dort):

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

{\overline {B}}(y_{0},R):=\{z\in E\ |\ \|z-y_{0}\|\leq R\}

die abgeschlossene Kugel um $y_{0}$ mit Radius $R$ . Ist

M:=\max\{\|f(x,y)\|\ |\ (x,y)\in [a,b]\times {\overline {B}}(y_{0},R)\}

und

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

dann existiert genau eine Lösung des Anfangswertproblems

y'=f(x,y),\quad y(a)=y_{0}

auf dem Intervall $[a,a+\alpha ]$ ; sie hat Werte in ${\overline {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}+\int _{x_{0}}^{x}f(t,y(t))dt$ ersetzt.

$y_{0}(x)=y_{0}$

$y_{i}(x)=y_{0}+\int _{x_{0}}^{x}f(t,y_{n-1}(t))dt;\quad 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})=\int _{t_{i}}^{t_{i+m}}f(t,y(t))dt$ . Das Integral kann durch numerische Quadraturformeln gelöst werden. Mit der Rechteckregel ergibt sich:

$\int _{t_{i}}^{t_{i+m}}f(t,y(t))dt\approx \int _{t_{i}}^{t_{i+m}}r(t)dt=\underbrace {(t_{i+m}-t_{i})} _{m\Delta t}\underbrace {r(t_{i})} _{f(t_{i},y(t_{i}))}=m\Delta tf(t_{i},y(t_{i}))$

Differenziationsmethoden

$y'(t_{i})\approx {\frac {y(t_{i+1})-y(t_{i})}{\Delta t}}$

$\underbrace {y'(t_{i})} _{f(t_{i},y(t_{i}))}\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.

Einschrittverfahren

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.

Mehrschrittverfahren

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 Steifes Anfangswertproblem. Im Englischen werden derartige DGL als stiff bezeichnet.

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

Euler

Explizit

$y(t_{i+m})-y(t_{i})=\int _{t_{i}}^{t_{i+m}}f(t,y(t))dt$

$\int _{t_{i}}^{t_{i+m}}f(t,y(t))dt\approx \int _{t_{i}}^{t_{i+m}}r(t)dt=\underbrace {(t_{i+m}-t_{i})} _{m\Delta t}\underbrace {r(t_{i})} _{f(t_{i},y(t_{i}))}=m\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}'={\frac {\Delta y}{\Delta t}}={\frac {y_{1}-y_{0}}{\Delta t}}$

$y_{1}=y_{0}+\Delta tf(t_{0},y_{0})$

oder allgemein

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

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

Siehe auch Explizites Euler-Verfahren, Burg, Haf, Wille, Meister: Seite 43ff, Knorrenschild: Seite 147ff oder Hanke-Bourgeois: Seite 557ff

Implizit

$y_{i+1}=y_{i}+\Delta tf(t_{i+1},y_{i+1});\quad 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 Implizites Euler-Verfahren, Burg, Haf, Wille, Meister: Seite 45ff oder Hanke-Bourgeois: Seite 560ff

Runge

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

$y_{i+1}=y_{i}+\Delta tf\left(t_{i}+{\frac {\Delta t}{2}},y_{i}+{\frac {\Delta t}{2}}f(t_{i},y_{i})\right);\quad i=0,\dots ,n-1$

Herleitung:

$y'_{i+1}={\frac {\Delta y}{\Delta t}}={\frac {y_{i+1}-y_{i}}{\Delta t}}=f(t_{i+1},y_{i+1})=f\left(t_{i}+{\frac {\Delta t}{2}},y_{i}+{\frac {\Delta t}{2}}f(t_{i},y_{i})\right)$

$y_{i+1}=y_{i}+\Delta tf\left(t_{i}+{\frac {\Delta t}{2}},y_{i}+{\frac {\Delta t}{2}}f(t_{i},y_{i})\right)$

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

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'=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'=t^{2}+0.1\cdot y,\quad y(-1.5)=y_{0}=0,\quad 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 Runge-Kutta-Verfahren, 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

Zurück zu Numerische Berechnung bestimmter Integrale |

Hoch zu Gesamtinhaltsverzeichnis |

Vor zu FFT und Wavelets