In [20]:
plot(Loes_DGL.rhs,xlim=(0,2),ylim=(0,6),ylabel="y(t)");
Dieses Jupyter Notebook basiert auf den Materialien der Vorlesung "Einführung in die Programmierung für Studierende der Physik (Introduction to Programming for Physicists)" (siehe Vorlesung 8 ).
Zunächst wird das Python Modul "sympy" eingebunden, das ein Computer-Algebra-System für Python bereitstellt und eine Vielzahl an symbolischen Berechnungen im Bereich der Mathematik und Physik relativ einfach möglich macht (weiteres siehe Vorlesung 6 Jupyter Notebooks und das Rechnen mit symbolischen Ausdrücken). Falls Sie das "sympy" Modul das erste Mal verwenden, müssen Sie es zunächst in Ihrer Python 3 Umgebung installieren (z.B. in einem Linux Terminal mit "pip3 install sympy").
from sympy import *
init_printing()
Wir betrachten in diesem Jupyter Notebook das numerische Lösen einer Differentialgleichung (DGL) erster Ordnung der Form
$$ \begin{equation} \dot{y}(t) = \frac{d y(t)}{dt} = f(t,y(t)) \quad, \, \hbox{mit:} \,\, a \leq t \leq b \, , \,\, y(a)=\alpha \quad. \end{equation} $$Die Funktion $f(t,y(t))$ bestimmt die DGL und somit das Verhalten der gesuchten Funktion $y(t)$. Es wird hierbei vorausgesetzt, dass $f(t,y(t))$ auf einer Teilmenge ${\cal D}=\{ (t,y) | a \leq t \leq b \, , \,\, -\infty \leq y \leq \infty \}$ kontinuierlich definiert ist. Weiter wird angenommen, dass so definierte Anfangswertproblem "well-posed" ist und eine eindeutige Lösung $y(t)$ existiert ("well-posed" bedeutet hier, dass die Differentialgleichung eine Struktur hat, bei der kleine Störungen im Anfangszustand nicht exponentiell anwachsen).
Wir betrachten im Folgenden die DGL:
$$ \begin{equation} \frac{d y(t)}{dt} = f(t,y(t)) = y(t)-t^2 +1 \quad, \, \hbox{mit:} \,\, 0 \leq t \leq 2 \, , \,\, y(0)=\alpha = 0.5 \quad. \end{equation} $$Zunächst definieren wir uns die DGL als sympy-Gleichung
t= symbols('t',real=True)
y = Function('y')(t)
DGL = Eq(y.diff(t), y-t**2 +1)
DGL
und versuchen mittels des Befehls "dsolve()" eine allgemeine analytische Lösung zu finden:
dsolve(DGL)
Die Konstante $C_1$ ist durch die oben angegebene Anfangsbedingung ($y(0)= \alpha$) bestimmt. Einarbeiten dieser Anfangsbedingung:
alpha= symbols('alpha')
y_0 = alpha
Loes_DGL = dsolve(DGL,ics={y.subs(t,0):y_0})
Loes_DGL
Diese allgemeine Lösung können wir später zum Vergleich in unserem C++ Programm verwenden. Wir schreiben sie somit als C++-Ausdruck:
print(ccode(Loes_DGL.rhs))
pow(t, 2) + 2*t + (alpha - 1)*exp(t) + 1
bzw.
from sympy.codegen.ast import Assignment
wert= symbols('wert')
print(ccode(Assignment(wert, Loes_DGL.rhs)))
wert = pow(t, 2) + 2*t + (alpha - 1)*exp(t) + 1;
Es wird somit eine allgemeine analytische Lösung gefunden, sodass wir die spezielle Anfangsbedingung ($y(0)= 0.5$, $\alpha=0.5$) einarbeiten können, um den noch unbestimmten Parameter $C_1$ bzw. $\alpha$ festzulegen.
y_0 = 0.5
Loes_DGL = dsolve(DGL,ics={y.subs(t,0):y_0})
Loes_DGL
Wir stellen uns die berechnete analytische Lösung der DGL dar:
plot(Loes_DGL.rhs,xlim=(0,2),ylim=(0,6),ylabel="y(t)");
Nicht jede DGL lässt sich analytisch lösen und falls der Befehl "dsolve()" keine sinnvollen Resultate liefert, muss man die zeitliche Entwicklung der Funktion $y(t)$ numerisch berechnen. Um zu verdeutlichen, wie eine solche numerische Lösung erzeugt wird, benutzen wir die DGL aus dem vorigen Beispiel (obwohl man hier eine analytische Lösung findet).
Wir benutzen hierfür das Python-Modul SciPy, welches eine breite Kollektion von mathematischen Algorithmen und Funktionen bereitstellt. Im Speziellen werden wir die Funktion 'solve_ivp(...)' verwenden, die sich im Untermodul 'scipy.integrate' befindet, welches Funktionen zum Lösen von gewöhnlichen Differentialgleichungen bereitstellt. Zusätzlich berechnen wir die numerische Lösung auch mittels der älteren Methode "integrate.odeint()", die ebenfalls in dem Python-Modul "scipy" definiert ist.
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
from scipy import integrate
Wir definieren uns die DGL als eine Funktion
def DGL_num(t, y):
dxdt = y-t**2 +1
return dxdt
Zusätzlich auch für die ältere Methode "integrate.odeint()"
def DGL_num_odeint(y,t):
dydt = y-t**2 +1
return dydt
und lösen das Anfangswertproblem im Bereich $t \in [0,2]$ unter Verwendung von $N=1000$ zeitlichen Gitterpunkten (mesh points).
N = 1000
tval = np.linspace(0, 2, N)
Loes_DGL_num = integrate.solve_ivp(DGL_num, [0, 2], [y_0], t_eval=tval)
Zusätzlich auch für die Methode "integrate.odeint()"
Loes_DGL_num_odeint = integrate.odeint(DGL_num_odeint, y_0, tval)
Wir stellen uns die numerische Lösung grafisch dar:
params = {
'figure.figsize' : [8,5],
# 'text.usetex' : True,
'axes.titlesize' : 14,
'axes.labelsize' : 16,
'xtick.labelsize' : 14 ,
'ytick.labelsize' : 14
}
matplotlib.rcParams.update(params)
plt.xlabel(r"$\rm t$")
plt.ylabel(r"$\rm y(t)$")
plt.plot(Loes_DGL_num.t, Loes_DGL_num.y[0],c="blue", label=r'$\rm solve-ivp() \, Loesung$')
plt.plot(tval, Loes_DGL_num_odeint,c="green", label=r'$\rm odeint() \, Loesung$')
plt.legend(frameon=True, loc="lower right",fontsize=15);
Wir stellen uns den Fehler $\Delta y$ der simulierten Datenpunkte zur analytischen Lösung dar ($\Delta y = y_{analytisch}-y$):
Loes_DGL_lam = lambdify(t, Loes_DGL.rhs)
plt.xlabel(r"$\rm t$")
plt.ylabel(r"$\rm \Delta y(t)$")
plt.plot(Loes_DGL_num.t, Loes_DGL_lam(Loes_DGL_num.t) - Loes_DGL_num.y[0],c="blue", label=r'$\rm solve-ivp() \, Loesung$')
plt.legend(frameon=True, loc="lower left",fontsize=15);
und zusätzlich auch für die Methode "integrate.odeint()"
plt.xlabel(r"$\rm t$")
plt.ylabel(r"$\rm \Delta y(t)$")
plt.plot(tval, Loes_DGL_lam(tval) - Loes_DGL_num_odeint[:,0],c="green", label=r'$\rm odeint() \, Loesung$')
plt.legend(frameon=True, loc="lower right",fontsize=15);
Wir vergleichen diese numerische Lösung mit den Resultaten des C++ Programms DGL_1.cpp. Wobei wir bei der Erzeugung der Daten zunächst nur 10 Gitterpunkte verwendeten und die Ergebnisse in eine Datei "DGL_1.dat" schrieben.
data = np.genfromtxt("./DGL_1.dat")
import matplotlib.gridspec as gridspec
# Bildabmessungen usw.
params = {
'figure.figsize' : [14,10],
'text.usetex' : True,
'axes.titlesize' : 14,
'axes.labelsize' : 16,
'xtick.labelsize' : 14 ,
'ytick.labelsize' : 14
}
matplotlib.rcParams.update(params)
fig = plt.figure() # Hauptbild
gs = gridspec.GridSpec(2, 2, width_ratios=[1,1], wspace=0.3, hspace=0.3) # Anordnung der vier Unterbilder
ax1 = plt.subplot(gs[0])
ax2 = plt.subplot(gs[1])
ax3 = plt.subplot(gs[2])
ax4 = plt.subplot(gs[3])
ax1.set_title(r'Euler Methode') # Titel der Abbildung in ax1
ax1.set_xlabel(r"$\rm t$")
ax1.set_ylabel(r"$\rm y(t)$")
ax2.set_title(r'Mittelpunkt Methode') # Titel der Abbildung in ax2
ax2.set_xlabel(r"$\rm t$")
ax2.set_ylabel(r"$\rm y(t)$")
ax3.set_title(r'Modifizierte Euler Methode') # Titel der Abbildung in ax3
ax3.set_xlabel(r"$\rm t$")
ax3.set_ylabel(r"$\rm y(t)$")
ax4.set_title(r'Runge-Kutta Ordnung vier') # Titel der Abbildung in ax4
ax4.set_xlabel(r"$\rm t$")
ax4.set_ylabel(r"$\rm y(t)$")
l_width = 0.8 # Festlegung der Plot-Liniendicke
alp = 0.7 # Festlegung der Transparenz der Kurven
ax1.plot(data[:,1],data[:,2], color="blue", linewidth=l_width, linestyle='-', alpha=alp, label=r'$\rm Euler-Loesung$')
ax1.plot(data[:,1],data[:,6], color="black", linewidth=l_width, linestyle='-', alpha=alp, label=r'$\rm Analytische \, Loesung$')
ax2.plot(data[:,1],data[:,3], color="blue", linewidth=l_width, linestyle='-', alpha=alp, label=r'$\rm Mittelpunkt-Loesung$')
ax2.plot(data[:,1],data[:,6], color="black", linewidth=l_width, linestyle='-', alpha=alp, label=r'$\rm Analytische \, Loesung$')
ax3.plot(data[:,1],data[:,4], color="blue", linewidth=l_width, linestyle='-', alpha=alp, label=r'$\rm Modifizierte \,Euler-Loesung$')
ax3.plot(data[:,1],data[:,6], color="black", linewidth=l_width, linestyle='-', alpha=alp, label=r'$\rm Analytische \, Loesung$')
ax4.plot(data[:,1],data[:,5], color="blue", linewidth=l_width, linestyle='-', alpha=alp, label=r'$\rm Runge-Kutta \, Ordnung \, vier \, Loesung$')
ax4.plot(data[:,1],data[:,6], color="black", linewidth=l_width, linestyle='-', alpha=alp, label=r'$\rm Analytische \, Loesung$')
ax1.legend(frameon=True, loc="lower right",fontsize=10) # Anordnung der Legende auf ax1
ax2.legend(frameon=True, loc="lower right",fontsize=10) # Anordnung der Legende auf ax2
ax3.legend(frameon=True, loc="lower right",fontsize=10) # Anordnung der Legende auf ax3
ax4.legend(frameon=True, loc="lower right",fontsize=10); # Anordnung der Legende auf ax4
Wir stellen uns den Fehler $\Delta y$ der simulierten Datenpunkte zur analytischen Lösung dar ($\Delta y = y_{analytisch}-y$):
fig = plt.figure() # Hauptbild
gs = gridspec.GridSpec(2, 2, width_ratios=[1,1], wspace=0.3, hspace=0.3) # Anordnung der vier Unterbilder
ax1 = plt.subplot(gs[0])
ax2 = plt.subplot(gs[1])
ax3 = plt.subplot(gs[2])
ax4 = plt.subplot(gs[3])
ax1.set_title(r'Euler Methode') # Titel der Abbildung in ax1
ax1.set_xlabel(r"$\rm t$")
ax1.set_ylabel(r"$\rm \Delta y(t)$")
ax2.set_title(r'Mittelpunkt Methode') # Titel der Abbildung in ax2
ax2.set_xlabel(r"$\rm t$")
ax2.set_ylabel(r"$\rm \Delta y(t)$")
ax3.set_title(r'Modifizierte Euler Methode') # Titel der Abbildung in ax3
ax3.set_xlabel(r"$\rm t$")
ax3.set_ylabel(r"$\rm \Delta y(t)$")
ax4.set_title(r'Runge-Kutta Ordnung vier') # Titel der Abbildung in ax4
ax4.set_xlabel(r"$\rm t$")
ax4.set_ylabel(r"$\rm \Delta y(t)$")
l_width=0.8 # Festlegung der Plot-Liniendicke
alp=0.7 # Festlegung der Transparenz der Kurven
ax1.plot(data[:,1],data[:,6]-data[:,2], color="blue", linewidth=l_width, linestyle='-', alpha=alp, label=r'$\rm Euler-Loesung$')
ax2.plot(data[:,1],data[:,6]-data[:,3], color="blue", linewidth=l_width, linestyle='-', alpha=alp, label=r'$\rm Mittelpunkt-Loesung$')
ax3.plot(data[:,1],data[:,6]-data[:,4], color="blue", linewidth=l_width, linestyle='-', alpha=alp, label=r'$\rm Modifizierte Euler-Loesung$')
ax4.plot(data[:,1],data[:,6]-data[:,5], color="blue", linewidth=l_width, linestyle='-', alpha=alp, label=r'$\rm Runge-Kutta \, Ordnung \, vier \, Loesung$')
ax1.legend(frameon=True, loc="lower right",fontsize=10) # Anordnung der Legende auf ax1
ax2.legend(frameon=True, loc="lower right",fontsize=10) # Anordnung der Legende auf ax2
ax3.legend(frameon=True, loc="lower right",fontsize=10) # Anordnung der Legende auf ax3
ax4.legend(frameon=True, loc="lower right",fontsize=10); # Anordnung der Legende auf ax4
Wir erhöhen nun die Anzahl der Gitterpunkte ($N=10000$) und betrachten uns wieder den Fehler $\Delta y$.
data = np.genfromtxt("./DGL_1a.dat")
fig = plt.figure() # Hauptbild
gs = gridspec.GridSpec(2, 2, width_ratios=[1,1], wspace=0.3, hspace=0.3) # Anordnung der vier Unterbilder
ax1 = plt.subplot(gs[0])
ax2 = plt.subplot(gs[1])
ax3 = plt.subplot(gs[2])
ax4 = plt.subplot(gs[3])
ax1.set_title(r'Euler Methode') # Titel der Abbildung in ax1
ax1.set_xlabel(r"$\rm t$")
ax1.set_ylabel(r"$\rm \Delta y$")
ax2.set_title(r'Mittelpunkt Methode') # Titel der Abbildung in ax2
ax2.set_xlabel(r"$\rm t$")
ax2.set_ylabel(r"$\rm \Delta y$")
ax3.set_title(r'Modifizierte Euler Methode') # Titel der Abbildung in ax3
ax3.set_xlabel(r"$\rm t$")
ax3.set_ylabel(r"$\rm \Delta y$")
ax4.set_title(r'Runge-Kutta Ordnung vier') # Titel der Abbildung in ax4
ax4.set_xlabel(r"$\rm t$")
ax4.set_ylabel(r"$\rm \Delta y$")
l_width=0.8 # Festlegung der Plot-Liniendicke
alp=0.7 # Festlegung der Transparenz der Kurven
ax1.plot(data[:,1],data[:,6]-data[:,2], color="blue", linewidth=l_width, linestyle='-', alpha=alp, label=r'$\rm Euler-Loesung$')
ax2.plot(data[:,1],data[:,6]-data[:,3], color="blue", linewidth=l_width, linestyle='-', alpha=alp, label=r'$\rm Mittelpunkt-Loesung$')
ax3.plot(data[:,1],data[:,6]-data[:,4], color="blue", linewidth=l_width, linestyle='-', alpha=alp, label=r'$\rm Modifizierte Euler-Loesung$')
ax4.plot(data[:,1],data[:,6]-data[:,5], color="blue", linewidth=l_width, linestyle='-', alpha=alp, label=r'$\rm Runge-Kutta Ordnung vier-Loesung$')
ax1.legend(frameon=True, loc="lower right",fontsize=10) # Anordnung der Legende auf ax1
ax2.legend(frameon=True, loc="lower right",fontsize=10) # Anordnung der Legende auf ax2
ax3.legend(frameon=True, loc="lower right",fontsize=10) # Anordnung der Legende auf ax3
ax4.legend(frameon=True, loc="lower right",fontsize=10); # Anordnung der Legende auf ax4
Die in den oberen Abbildungen dargestellten Fehler $\Delta y$ hängen stark von dem zugrundeliegenden nummerischen Verfahren ab. Die Runge-Kutta Ordnung vier Lösung liefert hier mit Abstand die besten Ergebnisse und die mit ihr berechneten Werte sind 14 Nachkommastellen genau.
Bei dem numerischen Lösen in Python (mittels solve_ivp(...) bzw. odeint(...)) kann der Benutzer die relativen und absoluten Fehler-Toleranzen der berechneten nummerischen Werte ebenfalls festlegen, wober der relative Fehler mit der Zusatzoption 'rtol' und der absolute Fehler mit der Zusatzoption 'atol' kontrolliert wird (näheres siehe scipy.integrate.solve_ivp).
In Folgenden setzen wir die relativen und absoluten Fehler-Toleranzen (rtol und atol) auf $10^{-10}$, $10^{-12}$ und $10^{-13}$ (noch kleinere Werte konnten nicht ausgeführt werden) und stellen uns den Fehler $\Delta y$ der Simulationen wieder grafisch dar.
fehler = 10**(-10)
N = 10000
tval = np.linspace(0, 2, N)
Loes_DGL_num_F10 = integrate.solve_ivp(DGL_num, [0, 2], [y_0], t_eval=tval, rtol=fehler, atol=fehler)
Loes_DGL_num_odeint_F10 = integrate.odeint(DGL_num_odeint, y_0, tval, rtol=fehler, atol=fehler)
fehler = 10**(-12)
Loes_DGL_num_F12 = integrate.solve_ivp(DGL_num, [0, 2], [y_0], t_eval=tval, rtol=fehler, atol=fehler)
Loes_DGL_num_odeint_F12 = integrate.odeint(DGL_num_odeint, y_0, tval, rtol=fehler, atol=fehler)
fehler = 10**(-13)
Loes_DGL_num_F13 = integrate.solve_ivp(DGL_num, [0, 2], [y_0], t_eval=tval, rtol=fehler, atol=fehler)
Loes_DGL_num_odeint_F13 = integrate.odeint(DGL_num_odeint, y_0, tval, rtol=fehler, atol=fehler)
fig = plt.figure() # Hauptbild
gs = gridspec.GridSpec(2, 2, width_ratios=[1,1], wspace=0.3, hspace=0.3) # Anordnung der vier Unterbilder
ax1 = plt.subplot(gs[0])
ax2 = plt.subplot(gs[1])
ax3 = plt.subplot(gs[2])
ax4 = plt.subplot(gs[3])
ax1.set_title(r'Fehler-Toleranzen rtol=atol=$10^{-10}$') # Titel der Abbildung in ax1
ax1.set_xlabel(r"$\rm t$")
ax1.set_ylabel(r"$\rm \Delta y$")
ax2.set_title(r'Fehler-Toleranzen rtol=atol=$10^{-12}$') # Titel der Abbildung in ax2
ax2.set_xlabel(r"$\rm t$")
ax2.set_ylabel(r"$\rm \Delta y$")
ax3.set_title(r'Fehler-Toleranzen rtol=atol=$10^{-13}$') # Titel der Abbildung in ax3
ax3.set_xlabel(r"$\rm t$")
ax3.set_ylabel(r"$\rm \Delta y$")
ax4.set_title(r'Runge-Kutta Ordnung vier') # Titel der Abbildung in ax4
ax4.set_xlabel(r"$\rm t$")
ax4.set_ylabel(r"$\rm \Delta y$")
l_width = 0.8 # Festlegung der Plot-Liniendicke
alp = 0.7 # Festlegung der Transparenz der Kurven
ax1.plot(tval, Loes_DGL_lam(tval) - Loes_DGL_num_odeint_F10[:,0],c="green", label=r'$\rm odeint() \, Loesung$')
ax1.plot(Loes_DGL_num_F10.t, Loes_DGL_lam(Loes_DGL_num_F10.t) - Loes_DGL_num_F10.y[0],c="blue", label=r'$\rm solve-ivp() \, Loesung$')
ax2.plot(tval, Loes_DGL_lam(tval) - Loes_DGL_num_odeint_F12[:,0],c="green", label=r'$\rm odeint() \, Loesung$')
ax2.plot(Loes_DGL_num_F12.t, Loes_DGL_lam(Loes_DGL_num_F12.t) - Loes_DGL_num_F12.y[0],c="blue", label=r'$\rm solve-ivp() \, Loesung$')
ax3.plot(tval, Loes_DGL_lam(tval) - Loes_DGL_num_odeint_F13[:,0],c="green", label=r'$\rm odeint() \, Loesung$')
ax3.plot(Loes_DGL_num_F13.t, Loes_DGL_lam(Loes_DGL_num_F13.t) - Loes_DGL_num_F13.y[0],c="blue", label=r'$\rm solve-ivp() \, Loesung$')
ax3.plot(data[:,1],data[:,6]-data[:,5], color="black", linewidth=l_width, linestyle='-', alpha=alp, label=r'$\rm Runge-Kutta \, Ordnung \, vier-Loesung$')
ax4.plot(data[:,1],data[:,6]-data[:,5], color="black", linewidth=l_width, linestyle='-', alpha=alp, label=r'$\rm Runge-Kutta \, Ordnung \, vier-Loesung$')
ax1.legend(frameon=True, loc="upper left",fontsize=10) # Anordnung der Legende auf ax1
ax2.legend(frameon=True, loc="upper left",fontsize=10) # Anordnung der Legende auf ax2
ax3.legend(frameon=True, loc="upper left",fontsize=10) # Anordnung der Legende auf ax3
ax4.legend(frameon=True, loc="lower left",fontsize=10); # Anordnung der Legende auf ax4
