Einführung in die Programmierung für Studierende der Physik¶

(Introduction to Programming for Physicists)¶

Vorlesung gehalten an der J.W.Goethe-Universität in Frankfurt am Main¶

(Sommersemester 2026)¶

von Dr.phil.nat. Dr.rer.pol. Matthias Hanauske¶

Erster Vorlesungsteil:¶

Numerisches Lösen von Differentialgleichungen (das Anfangswertproblem)¶

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 ).

Allgemeine Betrachtungen¶

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").

In [1]:
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).

Beispiel: Analytische Lösung¶

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

In [2]:
t= symbols('t',real=True)
y = Function('y')(t)
DGL = Eq(y.diff(t), y-t**2 +1)
DGL
Out[2]:
$\displaystyle \frac{d}{d t} y{\left(t \right)} = - t^{2} + y{\left(t \right)} + 1$

und versuchen mittels des Befehls "dsolve()" eine allgemeine analytische Lösung zu finden:

In [3]:
dsolve(DGL)
Out[3]:
$\displaystyle y{\left(t \right)} = C_{1} e^{t} + t^{2} + 2 t + 1$

Die Konstante $C_1$ ist durch die oben angegebene Anfangsbedingung ($y(0)= \alpha$) bestimmt. Einarbeiten dieser Anfangsbedingung:

In [4]:
alpha= symbols('alpha')
y_0 = alpha
Loes_DGL = dsolve(DGL,ics={y.subs(t,0):y_0})
Loes_DGL
Out[4]:
$\displaystyle y{\left(t \right)} = t^{2} + 2 t + \left(\alpha - 1\right) e^{t} + 1$

Diese allgemeine Lösung können wir später zum Vergleich in unserem C++ Programm verwenden. Wir schreiben sie somit als C++-Ausdruck:

In [5]:
print(ccode(Loes_DGL.rhs))

pow(t, 2) + 2*t + (alpha - 1)*exp(t) + 1

bzw.

In [6]:
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.

In [7]:
y_0 = 0.5
Loes_DGL = dsolve(DGL,ics={y.subs(t,0):y_0})
Loes_DGL
Out[7]:
$\displaystyle y{\left(t \right)} = t^{2} + 2 t - 0.5 e^{t} + 1$

Wir stellen uns die berechnete analytische Lösung der DGL dar:

In [8]:
plot(Loes_DGL.rhs,xlim=(0,2),ylim=(0,6),ylabel="y(t)");
No description has been provided for this image

Beispiel: Numerische Lösung¶

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.

In [9]:
import numpy as np
import matplotlib.pyplot as plt 
from matplotlib import rcParams
from scipy import integrate

Wir definieren uns die DGL als eine Funktion

In [10]:
def DGL_num(t, y):
    dxdt = y-t**2 +1
    return dxdt

Zusätzlich auch für die ältere Methode "integrate.odeint()"

In [11]:
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).

In [12]:
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()"

In [13]:
Loes_DGL_num_odeint = integrate.odeint(DGL_num_odeint, y_0, tval)

Wir stellen uns die numerische Lösung grafisch dar:

In [14]:
rcParams.update({
    'figure.figsize'    : [8,5],
    'axes.titlesize' : 14,
    'axes.labelsize' : 16,  
    'xtick.labelsize' : 14 ,
    'ytick.labelsize' : 14 
})
In [15]:
plt.xlabel("t")
plt.ylabel("y(t)")
plt.plot(Loes_DGL_num.t, Loes_DGL_num.y[0],c="blue", label='solve_ivp()-Lösung')
plt.plot(tval, Loes_DGL_num_odeint,c="green", label='odeint()-Lösung')
plt.legend(frameon=True, loc="lower right",fontsize=15);
No description has been provided for this image

Wir stellen uns den Fehler $\Delta y$ der simulierten Datenpunkte zur analytischen Lösung dar ($\Delta y = y_{analytisch}-y$):

In [16]:
Loes_DGL_lam = lambdify(t, Loes_DGL.rhs)
In [17]:
plt.xlabel("t")
plt.ylabel(r"$\Delta y$")
delta_y = Loes_DGL_lam(Loes_DGL_num.t) - Loes_DGL_num.y[0]
plt.plot(Loes_DGL_num.t, delta_y, c="blue", label='solve_ivp()-Lösung')
plt.legend(frameon=True, loc="lower left",fontsize=15);
No description has been provided for this image

und zusätzlich auch für die Methode "integrate.odeint()"

In [18]:
plt.xlabel("t")
plt.ylabel(r"$\Delta y$")
delta_y = Loes_DGL_lam(tval) - Loes_DGL_num_odeint[:,0]
plt.plot(tval, delta_y, c="green", label=r'odeint()-Lösung')
plt.legend(frameon=True, loc="lower right",fontsize=15);
No description has been provided for this image

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.

In [19]:
data = np.genfromtxt("./DGL_1.dat")
In [20]:
import matplotlib.gridspec as gridspec
# Bildabmessungen usw.
rcParams.update({
    'figure.figsize'    : [14,10],
    'axes.titlesize' : 14,
    'axes.labelsize' : 16,  
    'xtick.labelsize' : 14 ,
    'ytick.labelsize' : 14 
})
In [21]:
fig = plt.figure()
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])

for ax in [ax1, ax2, ax3, ax4]:
    ax.set_xlabel("t")
    ax.set_ylabel("y(t)")

# Titel der Unterbilder
ax1.set_title("Euler-Verfahren")
ax2.set_title("Mittelpunkt-Verfahren")
ax3.set_title("Modifiziertes Euler-Verfahren")
ax4.set_title("RK4-Verfahren") # Runge-Kutta-Verfahren 4. Ordnung

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="Euler-Lösung")
ax1.plot(data[:,1], data[:,6], color="black", linewidth=l_width, linestyle='-', alpha=alp, label="Analytische Lösung")

ax2.plot(data[:,1], data[:,3], color="blue", linewidth=l_width, linestyle='-', alpha=alp, label="Mittelpunkt-Lösung")
ax2.plot(data[:,1], data[:,6], color="black", linewidth=l_width, linestyle='-', alpha=alp, label="Analytische Lösung")

ax3.plot(data[:,1], data[:,4], color="blue", linewidth=l_width, linestyle='-', alpha=alp, label="Modifizierte Euler-Lösung")
ax3.plot(data[:,1], data[:,6], color="black", linewidth=l_width, linestyle='-', alpha=alp, label="Analytische Lösung")

ax4.plot(data[:,1], data[:,5], color="blue", linewidth=l_width, linestyle='-', alpha=alp, label="RK4-Lösung")
ax4.plot(data[:,1], data[:,6], color="black", linewidth=l_width, linestyle='-', alpha=alp, label="Analytische Lösung")

# Anordnung der Legenden
ax1.legend(frameon=True, loc="lower right", fontsize=12)
ax2.legend(frameon=True, loc="lower right", fontsize=12)
ax3.legend(frameon=True, loc="lower right", fontsize=12)
ax4.legend(frameon=True, loc="lower right", fontsize=12);
No description has been provided for this image

Wir stellen uns den Fehler $\Delta y$ der simulierten Datenpunkte zur analytischen Lösung dar ($\Delta y = y_{analytisch}-y$):

In [22]:
fig = plt.figure()
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])

for ax in [ax1, ax2, ax3, ax4]:
    ax.set_xlabel("t")
    ax.set_ylabel(r"$\Delta y$")

# Titel der Unterbilder
ax1.set_title("Euler-Verfahren")
ax2.set_title("Mittelpunkt-Verfahren")
ax3.set_title("Modifiziertes Euler-Verfahren")
ax4.set_title("RK4-Verfahren") # Runge-Kutta-Verfahren 4. Ordnung

ax1.plot(data[:,1], data[:,6]-data[:,2], color="blue", linewidth=l_width, linestyle='-', alpha=alp, label="Euler-Lösung")
ax2.plot(data[:,1], data[:,6]-data[:,3], color="blue", linewidth=l_width, linestyle='-', alpha=alp, label="Mittelpunkt-Lösung")
ax3.plot(data[:,1], data[:,6]-data[:,4], color="blue", linewidth=l_width, linestyle='-', alpha=alp, label="Modifizierte Euler-Lösung")
ax4.plot(data[:,1], data[:,6]-data[:,5], color="blue", linewidth=l_width, linestyle='-', alpha=alp, label="RK4-Lösung")

# Anordnung der Legenden
ax1.legend(frameon=True, loc="lower right", fontsize=12)
ax2.legend(frameon=True, loc="lower right", fontsize=12)
ax3.legend(frameon=True, loc="lower right", fontsize=12)
ax4.legend(frameon=True, loc="lower right", fontsize=12);
No description has been provided for this image

Wir erhöhen nun die Anzahl der Gitterpunkte ($N=10000$) und betrachten uns wieder den Fehler $\Delta y$.

In [23]:
data = np.genfromtxt("./DGL_1a.dat")
In [24]:
fig = plt.figure()
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])

for ax in [ax1, ax2, ax3, ax4]:
    ax.set_xlabel("t")
    ax.set_ylabel(r"$\Delta y$")

# Titel der Unterbilder
ax1.set_title("Euler-Verfahren")
ax2.set_title("Mittelpunkt-Verfahren")
ax3.set_title("Modifiziertes Euler-Verfahren")
ax4.set_title("RK4-Verfahren") # Runge-Kutta-Verfahren 4. Ordnung

ax1.plot(data[:,1], data[:,6]-data[:,2], color="blue", linewidth=l_width, linestyle='-', alpha=alp, label="Euler-Lösung")
ax2.plot(data[:,1], data[:,6]-data[:,3], color="blue", linewidth=l_width, linestyle='-', alpha=alp, label="Mittelpunkt-Lösung")
ax3.plot(data[:,1], data[:,6]-data[:,4], color="blue", linewidth=l_width, linestyle='-', alpha=alp, label="Modifizierte Euler-Lösung")
ax4.plot(data[:,1], data[:,6]-data[:,5], color="blue", linewidth=l_width, linestyle='-', alpha=alp, label="RK4-Lösung")

# Anordnung der Legenden
ax1.legend(frameon=True, loc="lower right", fontsize=12)
ax2.legend(frameon=True, loc="lower right", fontsize=12)
ax3.legend(frameon=True, loc="lower right", fontsize=12)
ax4.legend(frameon=True, loc="lower left", fontsize=12);
No description has been provided for this image

Die in den oberen Abbildungen dargestellten Fehler $\Delta y$ hängen stark von dem zugrundeliegenden nummerischen Verfahren ab. Das Runge-Kutta-Verfahren (4. Ordnung) 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.

In [25]:
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)
In [26]:
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)
In [27]:
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)
In [28]:
fig = plt.figure()
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])

for ax in [ax1, ax2, ax3, ax4]:
    ax.set_xlabel("t")
    ax.set_ylabel(r"$\Delta y$")

# Titel der Unterbilder
ax1.set_title(r"Fehler-Toleranzen rtol=atol=$10^{-10}$")
ax2.set_title(r"Fehler-Toleranzen rtol=atol=$10^{-12}$")
ax3.set_title(r"Fehler-Toleranzen rtol=atol=$10^{-13}$")
ax4.set_title("RK4-Verfahren") # Runge-Kutta-Verfahren 4. Ordnung

ax1.plot(tval, Loes_DGL_lam(tval) - Loes_DGL_num_odeint_F10[:,0],c="green", label='odeint()-Lösung')
ax1.plot(Loes_DGL_num_F10.t, Loes_DGL_lam(Loes_DGL_num_F10.t) - Loes_DGL_num_F10.y[0],c="blue", label='solve_ivp()-Lösung')
ax2.plot(tval, Loes_DGL_lam(tval) - Loes_DGL_num_odeint_F12[:,0],c="green", label='odeint()-Lösung')
ax2.plot(Loes_DGL_num_F12.t, Loes_DGL_lam(Loes_DGL_num_F12.t) - Loes_DGL_num_F12.y[0],c="blue", label='solve_ivp()-Lösung')
ax3.plot(tval, Loes_DGL_lam(tval) - Loes_DGL_num_odeint_F13[:,0],c="green", label='odeint()-Lösung')
ax3.plot(Loes_DGL_num_F13.t, Loes_DGL_lam(Loes_DGL_num_F13.t) - Loes_DGL_num_F13.y[0],c="blue", label='solve_ivp()-Lösung')
ax3.plot(data[:,1],data[:,6]-data[:,5], color="black", linewidth=l_width, linestyle='-', alpha=alp, label="RK4-Lösung")
ax4.plot(data[:,1],data[:,6]-data[:,5], color="black", linewidth=l_width, linestyle='-', alpha=alp, label="RK4-Lösung")

# Anordnung der Legenden
ax1.legend(frameon=True, loc="upper left", fontsize=12)
ax2.legend(frameon=True, loc="upper left", fontsize=12)
ax3.legend(frameon=True, loc="upper left", fontsize=12)
ax4.legend(frameon=True, loc="lower left", fontsize=12);
No description has been provided for this image