4.3.2. Gewöhnliche DGL#
Show code cell source
import numpy as np
np.set_printoptions(precision=2, linewidth=65)
import matplotlib.pyplot as plt
plt.rc('figure', dpi=150)
import seaborn as sns
sns.set()
sns.set_style('ticks')
sns.set_context("notebook", font_scale=1.2, rc={"lines.linewidth": 1.2})
import scipy.integrate
Bei gewöhnlichen Differentialgleichungen werden zwei Fälle unterschieden: Anfangswertproblem und Randwertproblem. Für beide gilt, dass die Lösung auf einem Intervall gesucht wird, jedoch wird im ersten Fall nur der Anfangswert vorgegeben und im zweiten Fall die Lösungswerte an den (im eindimensionalen Fall zwei) Randpunkten des Intervalls.
In diesem Kapitel wird nur auf das Anfangswertproblem eingegangen. Als Beispiel dafür und dessen Lösung werden im Folgenden zwei Gleichungen vorgestellt und numerisch gelöst.
Hinweis: Alle physikalischen Konstanten sind in diesem Kapitel vernachlässigt, um den Fokus auf die numerischen Lösungsmethoden zu setzten.
Temperaturentwicklung#
Im ersten Beispiel wird die zeitliche Temperaturentwicklung \(\sf T(t)\) mittels der Gleichung
vorgeschrieben und ein Anfangswert, d.h. \(\sf T(t=0) = T_0\), gegeben. Dabei ist \(\sf T_\infty\) eine Konstante und \(\sf Q(t)\) eine Funktion der Zeit. Die Gleichung beschreibt die Änderung der Temperatur aufgrund des Unterschieds zur Umgebung (\(\sf T_\infty\)) und eines vorgegebenen Wärmestroms (\(\sf Q(t)\)).
Die Diskretisierung in der Zeit mit einem frei gewählten \(\sf \Delta t\) führt zu folgender Näherung der Ableitung
Hierbei wird der Index, welcher der Zeitpunkt angibt, oben aufgeführt. Dies ist eine Konvention, welche später bei partiellen Differentialgleichungen die Indizierung lesbarer macht.
Die Terme auf der rechnten Seite der Gleichung sind aber noch nicht diskretisiert. Hierbei gibt es zwei Möglichkeiten: die Auswertung zum Zeitpunkt \(\sf n\) oder \(\sf n+1\). Wird die erste gewählt, so führt dies zum expliziten Euler-Verfahren
Die zweite Möglichkeit, d.h. die Auswertung der Terme zum Zeitpunkt \(\sf n+1\) führt zum impliziten Euler-Verfahren, auf welches hier nicht weiter eingegangen wird.
Die Umformung der Gleichung, welche nach dem expliziten Euler-Verfahren erzeugt wurde, nach dem nächsten Wert der Lösung, d.h. \(\sf T^{n+1}\), ergibt
Damit kann mittels der Informationen zum Zeitpunkt \(\sf n\) die Lösung zum Zeitpunkt \(\sf n+1\) direkt (explizit) berechnet werden. Mit dem vorgegebenen Anfangswert, \(\sf T^0 = T_0\) kann der nächste Lösungswert \(\sf T^1\), daraus dann \(\sf T^2\) und so weiter berechnet werden.
Als erstes Lösungsbeispiel wird
gewählt und die Lösung soll bis \(\sf t=10\) mit \(\sf \Delta t=0.5\) gesucht werden.
# Setzten der Konstanten
T0 = 400
Tinfty = 300
dt = 0.5
# Anlegen des Zeit- und Lösungsarrays
zeit = np.arange(0, 10, dt)
T = np.zeros_like(zeit)
# Setzten des Anfangswerts
T[0] = T0
# Lösung nach dem expliziten Euler-Verfahren
for n in range(0, len(zeit)-1):
T[n+1] = T[n] - dt * (T[n] - Tinfty)
# Visualisierung der Lösung
plt.plot(zeit, T, marker='.')
plt.xlabel("t")
plt.ylabel("T");
Im zweiten Lösungsbeispiel wird nun die Funktion \(\sf Q(t)\) als Stufenfunktion hinzugenommen, mit
und die Lösung soll bis \(\sf t=25\) berechnet werden.
# Anlegen des Zeit- und Lösungsarrays
zeit = np.arange(0, 25, dt)
T = np.zeros_like(zeit)
# Setzten des Anfangswerts
T[0] = T0
# Lösung nach dem expliziten Euler-Verfahren
for n in range(0, len(zeit)-1):
t = zeit[n]
# Quellterm Q
Q = 0
if t >= 5 and t <= 15:
Q = 100
# Zeititeration
T[n+1] = T[n] - dt * (T[n] - Tinfty) + Q
# Visualisierung der Lösung
plt.plot(zeit, T, marker='.')
plt.xlabel("t")
plt.ylabel("T");
Schwingung#
Genauso wie das vorangegangene Beispiel mit der Temperaturentwicklung, kann auch die Gleichung für eine gedämpfte harmonische Schwingung mit dem expliziten Euler-Verfahren gelöst werden. Die gewöhnliche Differentialgleichung zweiter Ordnung für die Auslenkung \(\sf u(t)\) lautet
Oder als System von zwei Gleichungen erster Ordnung
Die diskretisierte Form lautet
Die Lösung soll mit dem Anfangswerten \(\sf u(t=0) = u_0 = 1\) und \(\sf v(t=0) = v_0 = 0\) bis zur Zeit \(\sf t=25\) mit \(\sf \Delta t = 0.25\) berechnet werden.
# Setzten der Konstanten
u0 = 1
v0 = 0
dt = 0.25
# Anlegen des Zeit- und Lösungsarrays
zeit = np.arange(0, 25, dt)
u = np.zeros_like(zeit)
v = np.zeros_like(zeit)
# Setzten der Anfangswerte
u[0] = u0
v[0] = v0
# Lösung nach dem expliziten Euler-Verfahren
for n in range(0, len(zeit)-1):
u[n+1] = u[n] + dt * v[n]
v[n+1] = v[n] - dt * (v[n] + u[n])
# Visualisierung der Lösung
plt.plot(zeit, u, marker='.')
plt.xlabel("Zeit t")
plt.ylabel("Amplitude u");
Lösung mit dem scipy-Modul#
Das oben vorgestellte Lösungsverfahren mit der expliziten Euler-Methode ist eines der einfachsten Verfahren. Es existieren viele, vor allem bessere, Verfahren zur Lösung von gewöhnlichen Differentialgleichungen, welche hier aber nicht weiter vorgestellt werden. Einige der häufig eingesetzten Löser, wie z.B. das Runge-Kutta-Verfahren sind bereits im scipy-Modul implementiert. Anhand der obigen Anfangswertprobleme, wird deren Verwendung demonstriert.
Die Form, in welcher die zu lösenden Differentialgleichungen erster Ordnung vorliegen müssen ist
mit den Anfangsbedingungen für \(\sf y(t=0) = y_0\). Dabei wird für \(\sf f\) oft der Begriff ‘Rechte-Hand-Seite’, englisch right hand side (RHS), verwendet. Wird ein System von Gleichungen gelöst, so sind \(\sf y\) und \(\sf f\) entsprechend Vektoren.
Mit dieser Darstellung lautet die Funktion \(\sf f\) für die Temperaturgleichung
Die Funktion scipy.integrate.solve_ivp benötigt folgende Informationen für die Lösung der Gleichung
die Funktion \(\sf f\)
das Intervall auf dem die Gleichung gelöst werden soll, d.h. hier \(\sf t \in [0, t_{ende}]\)
den Anfangswert \(\sf T_0\)
Die Funktion \(\sf f\) muss dabei als eine Python-Funktion implementiert werden, welche zu gegebenen \(\sf t\) und \(\sf T\) eine entsprechende Rückgabe liefert. Dabei muss das erste Argument die unabhängige Variable sein, hier die Zeit \(\sf t\).
# Definition der Funktion f, d.h. der RHS
def f_T(t, T):
Tinfty = 300
Q = 0
if t >= 5 and t <= 15:
Q = 100
return -(T-Tinfty) + Q
# Definition der Konstanten
t_ende = 25
T0 = 400
Mit diesen Informationen kann nun der Löser aufgerufen werden. Das dabei verwendetet Lösungsverfahren kann mit dem optionalen Argument method ausgewählt werden, as Standardwert wird ein Runge-Kutta-Verfahren der fünften Ordnung (RK45) verwendet.
# Aufruf des Lösers
res = scipy.integrate.solve_ivp(f_T, [0, t_ende], [T0])
print(res)
message: The solver successfully reached the end of the integration interval.
success: True
status: 0
t: [ 0.000e+00 1.320e-01 ... 2.335e+01 2.500e+01]
y: [[ 4.000e+02 3.876e+02 ... 3.001e+02 3.000e+02]]
sol: None
t_events: None
y_events: None
nfev: 128
njev: 0
nlu: 0
Die Rückgabe enthält neben dem Lösungsarray (y) auch die unabhängige Variable (t) als auch weitere Informationen, wie z.B. die Anzahl der Funktionsauswertungen (nfev) der Funktion \(\sf f\).
Die graphische Ausgabe ist über die Variablen des Rückgabewerts direkt möglich. Dabei ist zu beachten, dass res.y in der allgemeinen Form eines mehrdimiensionalen Arrays vorliegt. In diesem Beispiel, in welchem nur eine Gleichung gelöst wird, muss trotzdem die erste und einzige Komponente der Lösung addressiert werden.
# Visualisierung der Lösung
plt.plot(res.t, res.y[0], marker='.')
plt.xlabel("t")
plt.ylabel("T");
Da es sich im Beispiel der Schwingungsgleichung von oben um ein Gleichungssystem handelt, ist die Funktion \(\sf \vec{f}\) folgende Vektorform an
Die Implementierung sieht entsprechend aus:
def f_osz(t, z):
u = z[0]
v = z[1]
# u'
up = v
# v'
vp = - v - u
return [up, vp]
Mit den gleichen Konstanten wie oben, kann die scipy-Funktion aufgerufen werden.
u0 = 1
v0 = 0
z0 = [u0, v0]
t_end = 25
res = scipy.integrate.solve_ivp(f_osz, [0, t_end], z0)
print(res)
message: The solver successfully reached the end of the integration interval.
success: True
status: 0
t: [ 0.000e+00 9.990e-04 ... 2.382e+01 2.500e+01]
y: [[ 1.000e+00 1.000e+00 ... 3.020e-06 -2.385e-06]
[ 0.000e+00 -9.985e-04 ... -7.475e-06 -1.761e-06]]
sol: None
t_events: None
y_events: None
nfev: 170
njev: 0
nlu: 0
Die Visualisierung kann genauso wie oben direkt durchgeführt werden.
# Visualisierung der Lösung
plt.plot(res.t, res.y[0], marker='.')
plt.xlabel("Zeit t")
plt.ylabel("Amplitude u");
