Musterlösung: Aufgabe 2, Übungsblatt 7¶
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¶
Frankfurt am Main 17.05.2026¶
Aufgabe 2
Diese Aufgabe ist angelehnt an das Kapitel 23 "Der gedämpfte harmonische Oszillator" des Buches von Prof. Walter Greiner, Mechanik (Teil 1) [5. Auflage, 1989, siehe Seite 226- 237]. Siehe auch Vorlesungsskript von Prof. Rischke auf Seite 117- 126 http://itp.uni-frankfurt.de/~drischke/Skript_MI_WiSe2022-2023.pdf ). Wir betrachten im Folgenden den gedämpften harmonischen Oszillator am Beispiel eines reibungsfrei gelagerten Wagens (Masse=$M$) auf den eine Rückstellkraft einwirkt (die proportional zu seiner Auslenkung $x$ ist (Proportionalitätskonstante $k$)), wobei zusätzlich eine geschwindigkeitsabhängige Reibungskraft auf den Wagen einwirkt (z.B. verursacht durch den auf den Wagen einwirkende Luftwiderstand, Stokesscher Ansatz: Proportionalitätskonstante $\alpha$). Aufgrund der Rückstellkraft, besitzt das zugrundeliegende Potenzial $V(x)$ die Form einer Parabel $V(x)=\frac{k \, x^2}{2}$. Die zeitliche Entwicklung des linearen harmonischen Oszillators mit Dämpfung wird mittels der folgenden Differenzialgleichung zweiter Ordnung beschrieben (wir setzen $\omega_0^2=\frac{k}{M}$ und $\beta = \frac{\alpha}{2M}$):
$$ \begin{equation} \ddot{x}(t) = - \omega_0^2 \, x(t) - 2 \beta \, \dot{x}(t) \end{equation} $$
Die Anfangsbedingungen seien zunächst noch allgemein gehalten: $x(0) = x_0 \,\, , \,\, \dot{x}(0) = v_0$. Bestimmen Sie die allgemeine Lösung der Differenzialgleichung mittels eines eigenen Jupyter Notebooks. Geben Sie dann die spezielle Lösung der Differenzialgleichung bei festgelegten Parameterwerten ($\omega_0^2=3$ und $\beta = 0.25$) und Anfangsbedingungen ($x_0 = 0$ und $v_0 = 40$) an und visualisieren Sie diese in einem x-t Diagramm. An welchem Ort befindet sich der Wagen zur Zeit $t=10$ ( $x(10)$ )? Erstellen Sie zusätzlich eine Animation, indem Sie einen der zuvor festgelegten Parameter (z.B. $x_0$, $v_0$, $\omega_0$ oder $\beta$) in einem gewissen Wertebereich verändern.
Berechnen Sie nun die allgemeine Lösung für den aperiodischen Grenzfall (kritische Dämpfung, $\beta = \omega_0$) und erstellen Sie eine weitere Animation, welche die Abbildung 2.23 im Vorlesungsskript von Prof. Rischke (Seite 124) verdeutlicht.
Wir berechnen zunächst die allgemeine Lösung der Differenzialgleichung (DGL):
from sympy import *
init_printing()
Definition der DGL:
t = symbols('t',real=True)
omega_0, beta = symbols('omega_0, beta', positive = True, real = True)
x = Function('x')(t)
DGL = Eq(x.diff(t).diff(t), -omega_0**2*x - 2*beta*x.diff(t))
DGL
Lösen der DGL (allgemein):
dsolve(DGL)
Die Konstanten $C_1$ und $C_2$ werden durch die oben angegebene Anfangsbedingung ($x(0) = x_0 \,\, , \,\, \dot{x}(0) = v_0$) festgelegt. Die allgemeine Lösung der DGL mit diesen Anfangsbedingung lautet dann:
x0, v0 = symbols(r'x_0, v_0', real = True)
Loes_allg = dsolve(DGL,ics={x.subs(t,0):x0,x.diff(t).subs(t, 0): v0})
Loes_allg
bzw. mittels 'simplify()'
Loes_allg.simplify()
Festlegung der Parameter und Anfangsbedingung auf die oben angegebenen Werte:
Loes_allg.subs({(omega_0,sqrt(3)),(beta,0.25),(x0,0),(v0,40)}).simplify()
Wir sind nur an dem reelwertigen Teil der Lösung interessiert:
Loes_allg_FestlPara = re(Loes_allg.rhs.subs({(omega_0,sqrt(3)),(beta,0.25),(x0,0),(v0,40)})).simplify()
Eq(Loes_allg.lhs,Loes_allg_FestlPara)
Wir stellen uns die berechnete analytische Lösung der DGL dar:
plot(Loes_allg_FestlPara, (t, 0, 15), xlabel = 't', ylabel = 'x(t)');
bzw. Darstellung mittels Matplotlib:
import matplotlib.pyplot as plt
import numpy as np
tvals = np.linspace(0, 15, 300)
func = lambdify(t, simplify(Loes_allg_FestlPara))
plt.ylabel('x(t)')
plt.xlabel('t')
plt.plot(tvals,func(tvals));
Zur Zeit $t=10$ befindet sich der Wagen bei $x(10)=$
Loes_allg_FestlPara.subs(t,10)
Wir visualisieren nun die Auswirkungen einer unterschiedlichen Dämpfung, indem wir den Parameter $\beta$ im Bereich von $\beta \in [0,1]$ verändern. Wir wählen dafür zunächst die spb-Python-Bibliothek (das SymPy Plotting Backends-Modul) welche interaktive Animationen von SymPy-Gleichungen direkt erstellen kann (siehe SPB Modules Reference).
Wir legen wie in der vorigen speziellen Lösung die Parameter und Anfangsbedingung auf die oben angegebenen Werte fest, lassen jedoch den Parameter $\beta$ dabei noch offen.
Loes_allg_beta = re(Loes_allg.rhs.subs({(omega_0,sqrt(3)),(x0,0),(v0,40)}))
Eq(Loes_allg.lhs,Loes_allg_beta)
Diesen analytischen Ausdruck visualisieren wir nun in einer interaktiven Animation, wobei $\beta \in [0,1]$.
from spb import plot as spb_plot
from IPython.display import IFrame
anim = spb_plot( Loes_allg_beta, (t, 0, 10),
params={beta: (0, 1)},
animation={"fps": 20, "time": 5},
title = (r"$\beta$ = {:.2f}", beta),
xlabel = ('t'), ylabel= ('x(t)'),
ylim = (-25, 25),
imodule = "panel", show = False
)
#anim.show().save("Oszi.html", embed=True)
IFrame(src="Oszi.html", width="100%", height="600px")
Alternativ kann man auch in Matplotlib die animierte Grafik erzeugen, speichern und im Jupyter Notebook darstellen (näheres siehe matplotlib.animation und Beispiele von Animationen).
Möchte man die Lösung mittels "matplotlib" visualisieren, muss man zunächst den analytischen Ausdruck mittels "lambdify(...)" in eine numerische Funktion umwandelt. Da die analytische Lösung $x(t)$ komplexwertige Wurzeln, während der internen Berechnung erzeugen kann, benutzen wir eine spezielle NumPy-Wurzelfunktion numpy.lib.scimath.sqrt, die komplexe Zahlen explizit zulässt.
func_beta = lambdify((t,beta), Loes_allg_beta, modules=[{'sqrt': np.lib.scimath.sqrt}, 'numpy'])
from matplotlib.animation import FuncAnimation
from IPython.display import HTML
fig, ax = plt.subplots(figsize=(7, 5))
frames_N = 20
def animate(i):
ax.cla()
beta_s = i/frames_N
ax.set_title(r"$\beta$ = {:.2f}".format(beta_s))
ax.set(xlabel='t', ylabel='x(t)', ylim=(-25, 25))
ax.plot(tvals, func_beta(tvals, beta_s))
return fig,
ani = FuncAnimation(fig, animate, frames=frames_N, interval=200)
ani.save("./Oszi.gif", writer='pillow')
plt.close(fig)
HTML(ani.to_html5_video())
Wir betrachten nun den aperiodischen Grenzfall (kritische Dämpfung, $\beta = \omega_0$) und berechnen zunächst den analytischen Ausdruck der allgemeinen Lösung der Differenzialgleichung.
Im aperiodischen Grenzfall lautet die Differenzialgleichung
DGL.subs(omega_0,beta)
mit der allgemeinen Lösung:
Loes_allg_ap = dsolve(DGL.subs(omega_0,beta),ics={x.subs(t,0):x0,x.diff(t).subs(t, 0): v0})
Loes_allg_ap
Dieser Ausdruck entspricht genau der Gleichung (2.123) im Vorlesungsskript von Prof. Rischke (Seite 124). Die Nullstelle der Lösung ($x(t_N) = 0$) berechnet sich allgemein zu:
t_N = solve(Eq(Loes_allg_ap.rhs,0),t)[0]
t_N
Für $t_N > 0$ besitzt die Lösung somit nur dann einen Nulldurchgang, falls $v_0 < − \beta \, x_0 < 0$.
Um die Abbildung 2.23 im Vorlesungsskript von Prof. Rischke (Seite 124, http://itp.uni-frankfurt.de/~drischke/Skript_MI_WiSe2022-2023.pdf) in einer Animation zu verdeutlichen, setzen wir die Anfangsauslenkung $x_0=1$ und den Wert der Dämpfung auf $\beta = 0.25$, sodass die Lösung jetzt lediglich von der Anfangsgeschwindigkeit $v_0$ abhängt.
Loes_allg_ap_FestlPara = Loes_allg_ap.subs({(x0,1),(beta,0.25)}).rhs
Eq(Loes_allg.lhs,Loes_allg_ap_FestlPara)
Für $t_N > 0$ besitzt die Lösung somit nur dann einen Nulldurchgang, falls $v_0 < − 0.25$.
t_N_FestlPara = t_N.subs({(x0,1),(beta,0.25)})
t_N_FestlPara
In der folgenden Animation visualisieren wir den aperiodischen Grenzfall in einer interaktiven Animation, wobei wir den Animationsparameter, die Anfangsgeschwindigkeit $v_0$, im Bereich von $v_0 \in [0,-0.7]$ variieren.
anim = spb_plot( Loes_allg_ap_FestlPara, (t, 0, 30),
params={v0: (0, -0.7)},
animation={"fps": 20, "time": 5},
title = (r"$v_0$ = {:.2f}, $t_N$ = {:.2f}", v0, t_N_FestlPara),
xlabel = ('t'), ylabel= ('x(t)'),
ylim = (-0.4, 1.1),
imodule = "panel", show = False
)
#anim.show().save("Oszi_ap.html", embed=True)
IFrame(src="Oszi_ap.html", width="100%", height="600px")
Wieder visualisieren wir auch die Animation mittels "matplotlib".
func_v0 = lambdify((t,v0), Loes_allg_ap_FestlPara)
frames_N = 40
tvals = np.linspace(0, 30, 300)
fig, ax = plt.subplots(figsize=(7, 5))
def animate(i):
ax.cla()
v0_s = -0.7*i/frames_N
t_N_s = t_N_FestlPara.subs(v0,v0_s)
ax.set_title(r"$v_0$ = {:.2f}, $t_N$ = {:.2f}".format(v0_s,t_N_s))
ax.set(xlabel='t', ylabel='x(t)', xlim=(0, 30), ylim=(-0.4, 1.1))
ax.plot(tvals, func_v0(tvals, v0_s))
ax.scatter(t_N_s,0, marker='o', color="red", s=25);
ax.plot([0,tvals[-1]], [0,0], color="black", linestyle=":")
return fig,
ani = FuncAnimation(fig, animate, frames=frames_N, interval=200)
ani.save("./Oszi_ap.gif", writer='pillow')
plt.close(fig)
HTML(ani.to_html5_video())