Das Projekt Achterbahn ist ein Anwendungsfall aus der klassischen Mechanik. Das System besteht aus zwei Achterbahn-Wagen unterschiedlicher Massen $m_i \, , \,\, i \in [0,1]$, die reibungsfrei auf einer Metallschiene gleiten können. Die Achterbahn hat eine Länge von $l \, [m]$ und besitzt keine Kurven, sodass man den Verlauf der Metallschiene der Achterbahn als eine eindimensionale Funktion $h(x) \, , \,\, x \in [-l/2, l/2]$ beschreiben kann, wobei $h(x)$ die Höhe der Metallschiene der Achterbahn an der Position $x$ in der Einheit Metern beschreibt. Die Koordinaten des Wagens seien $\vec{r} = \left( x(t), \, y(t) \right) = \left( x(t), \, h(x(t)) \right)$ und seine Masse sei durch den Parameter $m$ gekennzeichnet. Die Herleitung der Bewegungsgleichungen erfolgt am elegantesten mittels der Euler-Lagrange Gleichungen, bzw. mittels der Hamilton Theorie. Das System besitzt die folgende Lagrangefunktion: $$ \begin{eqnarray} & L = T - V \quad \hbox{mit:} &\\ & T = \frac{1}{2} m \left[ \dot{x}(t)^2 + \dot{y}(t)^2 \right] \quad , \quad V(t) = m \, g \, y(t) \,\, , & \\ & \hbox{mit:} \quad \dot{y}(t) = \frac{d y}{dt} = \frac{d h (x(t))}{dt} = \frac{d x(t)}{dt} \cdot \frac{d h (x)}{dx} = \dot{x}(t) \cdot \frac{d h (x)}{dx}\quad , & \end{eqnarray} $$ wobei $g$ der Wert der Erdbeschleunigung ist. Nun kann man für eine beliebige Form des Verlaufs der Achterbahn (beschrieben durch die Funktion $h(x)$) die Bewegungsgleichung des Wagens mittels der Lagrange-Gleichungen $$ \begin{eqnarray} &\frac{d }{dt} \left( \frac{\partial L}{\partial \dot{x}} \right) - \frac{\partial L}{\partial x} \, =\, 0 \quad \hbox{mit:} & \\ &L = \frac{1}{2} m \left[ \dot{x}(t)^2 + \left( \frac{d h(x(t))}{dt} \right)^2 \right] - m \, g \, h(t) = \frac{1}{2} m \left[ \dot{x}^2 + \left( \dot{x} \cdot \frac{d h(x)}{dx} \right)^2 \right] - m \, g \, h(x) = \frac{1}{2} m \, \dot{x}^2 \cdot\left[ 1 + \left( \frac{d h(x)}{dx} \right)^2 \right] - m \, g \, h(x)& \end{eqnarray} $$ herleiten.
Wir berechnen zunächst auf analytischem Weg, unter Verwendung der SymPy Bibliothek, die Euler-Lagrange Gleichungen der Achterbahn und schreiben diese dann in ein System von erster Ordnung Differentialgleichungen um. Dabei betrachten wir zunächst den Fall einer allgemeinen Funktion $h(x)$, damit eine, in späteren Teilaufgaben möglicherweise abgeänderte Form der Achterbahn, einfach in den Algorithmus eingebaut werden kann.
from sympy import *
init_printing()
t= symbols('t',real=True)
m, g = symbols('m, g', positive = True, real = True)
x = Function('x')(t)
h = Function('h')(x)
T = 1/2 * m * diff(x,t)**2 * ( 1 + diff(h,x)**2)
V = m * g * h
L = T - V
E = T + V
expand(L)
simplify(L)
Mittels der so definierten Lagrangefunktion $L$ und den Euler-Lagrange Gleichungen bestimmen wir die Bewegungsgleichungen des Systems. Die Euler-Lagrange Gleichung lautet wie folgt:
ELG_1 = diff( diff(L, diff(x,t) ) ,t) - diff(L, x)
Eq(0,ELG_1)
Dies ist eine Differentialgleichung zweiter Ordnung in der Zeit. Wir lösen die Gleichung nach $\frac{d^2 x}{dt^2}$ auf und erhalten die folgende DGL.
f=solve(ELG_1, diff(x,t,t))[0]
Eq(diff(x,t,t),f)
Wir betrachten nun eine spezielle Form der Achterbahn und legen die Funktion $h(x)$ fest.
h = exp(sqrt(x**2/10)) * (sin(x**2/10))**2
T = 1/2 * m * diff(x,t)**2 * ( 1 + diff(h,x)**2)
V = m * g * h
L = T - V
E = T + V
expand(L)
simplify(L)
Die Euler-Lagrangegleichungen lauten:
ELG = diff( diff(L, diff(x,t) ) ,t) - diff(L, x)
Eq(0,ELG)
Wir lösen die Gleichung wieder nach $\frac{d^2 x}{dt^2}$ auf und erhalten die folgende DGL.
f=solve(ELG, diff(x,t,t))[0]
Eq(diff(x,t,t),f)
Um die entstandene DGL zweiter Ordnung numerisch lösen zu können, schreiben wir sie in ein System von zwei DGLs erster Ordnung um. Wir betrachten nun die numerische Lösung und machen die folgenden Variablenumbenennung: $u_1(t)=x(t)$, und $u_2(t)=\dot{x}(t)$ und definieren das System von DGLs wie folgt:
$$ \begin{eqnarray} \dot{u_1}(t) &=& \frac{d u_1}{dt} = \frac{d x}{dt} = u_2(t) \\ \dot{u_2}(t) &=& \frac{d u_2}{dt} = \frac{d \dot{x}}{dt} = \ddot{x}(t) := f(u_1,u_2) \end{eqnarray} $$Wir bringen zunächst die DGL der Achterbahn in eine Form $\ddot{x}(t) = \frac{d^2 x}{dt^2} = f(u_1,u_2)$. Die DGL bestimmende Funktion $f(u_1,u_2)$ lautet:
u_1, u_2 = symbols('u_1, u_2', real = True)
f_c = f.subs(diff(x,t),u_2).subs(x,u_1)
f_c
Diesen Ausdruck exportieren wir in die Sprache C++ um ihn dann später im C++ Programm verwenden zu können.
print(ccode(f_c))
0.20000000000000001*u_1*(-200.0*g*pow(u_1, 2)*exp(0.316227766016838*fabs(u_1))*cos(0.10000000000000001*pow(u_1, 2)) - 158.11388300841898*g*exp(0.316227766016838*fabs(u_1))*sin(0.10000000000000001*pow(u_1, 2))*fabs(u_1) + 32.0*pow(u_1, 4)*pow(u_2, 2)*exp(0.63245553203367599*fabs(u_1))*pow(sin(0.10000000000000001*pow(u_1, 2)), 2)*cos(0.10000000000000001*pow(u_1, 2)) - 16.0*pow(u_1, 4)*pow(u_2, 2)*exp(0.63245553203367599*fabs(u_1))*cos(0.10000000000000001*pow(u_1, 2)) + 75.894663844041105*pow(u_1, 2)*pow(u_2, 2)*exp(0.63245553203367599*fabs(u_1))*pow(sin(0.10000000000000001*pow(u_1, 2)), 3)*fabs(u_1) + 80.0*pow(u_1, 2)*pow(u_2, 2)*exp(0.63245553203367599*fabs(u_1))*pow(sin(0.10000000000000001*pow(u_1, 2)), 3) - 60.0*pow(u_1, 2)*pow(u_2, 2)*exp(0.63245553203367599*fabs(u_1))*pow(sin(0.10000000000000001*pow(u_1, 2)), 2)*cos(0.10000000000000001*pow(u_1, 2)) - 63.245553203367592*pow(u_1, 2)*pow(u_2, 2)*exp(0.63245553203367599*fabs(u_1))*sin(0.10000000000000001*pow(u_1, 2))*fabs(u_1) - 80.0*pow(u_1, 2)*pow(u_2, 2)*exp(0.63245553203367599*fabs(u_1))*sin(0.10000000000000001*pow(u_1, 2)) - 15.811388300841898*pow(u_2, 2)*exp(0.63245553203367599*fabs(u_1))*pow(sin(0.10000000000000001*pow(u_1, 2)), 3)*fabs(u_1) - 63.245553203367592*pow(u_2, 2)*exp(0.63245553203367599*fabs(u_1))*pow(sin(0.10000000000000001*pow(u_1, 2)), 2)*cos(0.10000000000000001*pow(u_1, 2))*fabs(u_1))*sin(0.10000000000000001*pow(u_1, 2))/(100.0*pow(u_1, 2) + 16.0*pow(pow(u_1, 2)*cos(0.10000000000000001*pow(u_1, 2)) + 0.79056941504209488*sin(0.10000000000000001*pow(u_1, 2))*pow(fabs(u_1), 1.0), 2)*exp(0.63245553203367599*fabs(u_1))*pow(sin(0.10000000000000001*pow(u_1, 2)), 2))
Zum Lösen dieses gekoppelten Systems von zwei Differenzialgleichungen erster Ordnung benutzen wir das Python-Modul SciPy, welches eine breite Kollektion von mathematischen Algorithmen und Funktionen bereitstellt. Im Speziellen werden wir in diesem Semester die Funktion 'solve_ivp(...)' verwenden, die sich im Untermodul 'scipy.integrate' befindet, welches Funktionen zum Lösen von gewöhnlichen Differenzialgleichungen bereitstellt.
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
from scipy.integrate import solve_ivp
Um die mittels sympy gewonnenen analytische Ausdrücke direkt beim numerischen Lösen zu verwenden, schreiben wir die Funktion $f(u_1,u_2)$ in eine numerische Lambda-Funktion um. Zusätzlich definieren wir auch eine Lambda-Funktion der kinetischen Energie, potentiellen Energie, Gesamtenergie und Höhe $h(x)$. Auch legen wir ab jetzt den Wert der Erdbeschleunigung fest $g = 9.81$.
set_g = 9.81
set_m = 1
DGL_L = lambdify((u_1,u_2), f_c.subs(g,set_g))
E_L = lambdify((u_1,u_2), E.subs(g,set_g).subs(m,set_m).subs(diff(x,t),u_2).subs(x,u_1))
T_L = lambdify((u_1,u_2), T.subs(g,set_g).subs(m,set_m).subs(diff(x,t),u_2).subs(x,u_1))
V_L = lambdify((u_1,u_2), V.subs(g,set_g).subs(m,set_m).subs(diff(x,t),u_2).subs(x,u_1))
h_L = lambdify((x), h)
Wir definieren uns das System der DGLs als eine Funktion
def DGLs(t,u_vec):
u1, u2 = u_vec
du1_dt = u2
du2_dt = DGL_L(u1,u2)
return np.array([du1_dt,du2_dt])
und lösen das Anfangswertproblem im Bereich $t \in [0,15]$ unter Verwendung von $N=10000$ zeitlichen Gitterpunkten (mesh points). Wir wählen als Anfangsbedingungen die folgenden Werte: $x(t=0)=0.01$ [m] , $v_x(0) = 8.5$ [m/s].
x_0 = 0.01
v_0 = 8.5
initialval = np.array([x_0,v_0])
N=10000
t_end=15
t_val = np.linspace(0, t_end, N+1)
Bei dem numerischen Lösen von Differenzialgleichungen in Python (mittels solve_ivp(...)) kann der Benutzer zusätzlich die relativen und absoluten Fehler-Toleranzen der berechneten numerischen Werte festlegen, wobei 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^{-13}$.
fehler = 10**(-13)
Loes = solve_ivp(DGLs, [0, t_end], initialval, t_eval=t_val, rtol=fehler, atol=fehler)
Wir stellen uns die numerische Lösung grafisch dar:
params = {
'figure.figsize' : [10,5],
# 'text.usetex' : True,
'axes.titlesize' : 14,
'axes.labelsize' : 16,
'xtick.labelsize' : 14 ,
'ytick.labelsize' : 14
}
matplotlib.rcParams.update(params)
Die folgende Abbildung zeigt die zeitliche Entwicklung der x-Koordinate $x(t)$ des Wagens, der Geschwindigkeit in x-Richtung ($\dot{x} = v_x = \frac{d x}{dt}$) und seiner Höhe $h(x)$.
plt.xlabel(r"$\rm t$")
plt.ylabel(r"$\rm x \,\, \dot{x} \,\, h$")
plt.plot(Loes.t, Loes.y[0],c="blue", label=r"$\rm x(t)$");
plt.plot(Loes.t, Loes.y[1],c="red", label=r"$\rm \dot{x}(t)$");
plt.plot(Loes.t, h_L(Loes.y[0]),c="green", label=r"$\rm h(x)$");
plt.legend(loc='lower right',fontsize=16);
In der unteren Abbildung erkennt man, dass sich die potentielle und kinetische Energie periodisch austauschen und die Gesamtenergie des Systems erhalten bleibt.
plt.xlabel(r"$\rm t$")
plt.ylabel(r"$\rm Energie \, E \, , \, \, T \, , \, \, V$")
plt.plot(Loes.t, E_L(Loes.y[0],Loes.y[1]),c="black", label=r"$\rm Gesamtenergie \, E$");
plt.plot(Loes.t, T_L(Loes.y[0],Loes.y[1]),c="green", label=r"$\rm kinetische \, Energie \, T$");
plt.plot(Loes.t, V_L(Loes.y[0],Loes.y[1]),c="orange", label=r"$\rm potentielle \, Energie \, V$");
plt.legend(loc='upper right',fontsize=16);
Wir entwerfen nun eine Visualisierung der Achterbahn.
import matplotlib.animation as animation
from IPython.display import HTML
fig = plt.figure()
ax = fig.gca()
set_frames=250
stepi=int(N/set_frames)
plot_max_y=np.max(h_L(Loes.y[0])) + 0.5
plot_max_x=np.max(Loes.y[0]) + 0.5
def init():
ax.scatter(Loes.y[0][0], h_L(Loes.y[0][0]), s=400, marker='o', c="red")
return fig,
def animate(i):
ax.cla()
ax.scatter(Loes.y[0][stepi*i], h_L(Loes.y[0][stepi*i]), s=400, marker='s', c="red")
ax.plot(Loes.y[0],h_L(Loes.y[0]), c="black")
ax.set_xlim(-plot_max_x,plot_max_x)
ax.set_ylim(-0.1,plot_max_y)
ax.set_xlabel(r"$\rm x $")
ax.set_ylabel(r"$\rm y $")
ax.set_title('Die Achterbahn')
return fig,
ani = animation.FuncAnimation(fig,animate,init_func=init,frames=set_frames,interval=70)
plt.close(ani._fig)
HTML(ani.to_html5_video())
Nun erhöhen wir die Anfangsgeschwindigkeit auf $v_x(0) = 13.2$ [m/s].
x_0 = 0.01
v_0 = 13.2
initialval = np.array([x_0,v_0])
fehler = 10**(-13)
Loes1 = solve_ivp(DGLs, [0, t_end], initialval, t_eval=t_val, rtol=fehler, atol=fehler)
Die folgende Abbildung zeigt wieder die zeitliche Entwicklung der x-Koordinate $x(t)$ des Wagens, der Geschwindigkeit in x-Richtung ($\dot{x} = v_x = \frac{d x}{dt}$) und seiner Höhe $h(x)$.
plt.xlabel(r"$\rm t$")
plt.ylabel(r"$\rm x \,\, \dot{x} \,\, h$")
plt.plot(Loes1.t, Loes1.y[0],c="blue", label=r"$\rm x(t)$");
plt.plot(Loes1.t, Loes1.y[1],c="red", label=r"$\rm \dot{x}(t)$");
plt.plot(Loes1.t, h_L(Loes1.y[0]),c="green", label=r"$\rm h(x)$");
plt.legend(loc='lower right',fontsize=16);
In der unteren Abbildung erkennt man wieder, dass sich die potentielle und kinetische Energie periodisch austauschen und die Gesamtenergie des Systems erhalten bleibt.
plt.xlabel(r"$\rm t$")
plt.ylabel(r"$\rm Energie \, E \, , \, \, T \, , \, \, V$")
plt.plot(Loes1.t, E_L(Loes1.y[0],Loes1.y[1]),c="black", label=r"$\rm Gesamtenergie \, E$");
plt.plot(Loes1.t, T_L(Loes1.y[0],Loes1.y[1]),c="green", label=r"$\rm kinetische \, Energie \, T$");
plt.plot(Loes1.t, V_L(Loes1.y[0],Loes1.y[1]),c="orange", label=r"$\rm potentielle \, Energie \, V$");
plt.legend(loc='upper right',fontsize=16);
Wir stellen uns beide numerischen Lösungen in einer Animation dar:
fig = plt.figure()
ax = fig.gca()
set_frames=250
stepi=int(N/set_frames)
plot_max_y=np.max(h_L(Loes1.y[0])) + 0.5
plot_max_x=np.max(Loes1.y[0]) + 0.5
plot_min_x=np.min(Loes1.y[0]) - 0.5
x_L = np.linspace(plot_min_x, plot_max_x, 4000)
def init():
ax.scatter(Loes1.y[0][0], h_L(Loes1.y[0][0]), s=400, marker='o', c="green")
ax.scatter(Loes.y[0][0], h_L(Loes.y[0][0]), s=400, marker='o', c="red")
return fig,
def animate(i):
ax.cla()
ax.scatter(Loes1.y[0][stepi*i], h_L(Loes1.y[0][stepi*i]), s=400, marker='s', c="green")
ax.scatter(Loes.y[0][stepi*i], h_L(Loes.y[0][stepi*i]), s=400, marker='s', c="red")
ax.plot(x_L,h_L(x_L), c="black")
ax.set_xlim(-plot_max_x,plot_max_x)
ax.set_ylim(-0.1,plot_max_y)
ax.set_xlabel(r"$\rm x $")
ax.set_ylabel(r"$\rm y $")
ax.set_title('Die Achterbahn')
return fig,
ani = animation.FuncAnimation(fig,animate,init_func=init,frames=set_frames,interval=70)
plt.close(ani._fig)
HTML(ani.to_html5_video())
Am Beispiel des grünen Wagens vergleichen wir nun die numerischen Berechnungen mit den Ergebnissen des C++ Programmes (siehe Achterbahn.cpp). Bei der folgenden Simulation verwendeten wir $N=10000$ zeitliche Gitterpunkte.
data = np.genfromtxt("./Achterbahn_low.dat")
plt.xlabel(r"$\rm t$")
plt.ylabel(r"$\rm x \,\, \dot{x}$")
plt.plot(Loes1.t, Loes1.y[0],c="blue", label=r"$\rm x(t)$");
plt.plot(data[:,1], data[:,2],c="green", linestyle=":", label=r"$\rm x (C++)$");
plt.plot(Loes1.t, Loes1.y[1],c="red", label=r"$\rm \dot{x}(t)$");
plt.plot(data[:,1], data[:,3],c="black", linestyle=":", label=r"$\rm \dot{x} (C++)$");
plt.legend(loc='lower right',fontsize=16);
Die obere Abbildung zeigt, dass sich die Ergebnisse des C++ Programms ab $t \approx 2$ merklich von der Python Simulation unterscheidet. Der Grund liegt dabei an einer zu geringen Anzahl von Gitterpunkten im C++ Programm. Bei der folgenden Simulation verwenden wir eine größere Anzahl an zeitlichen Gitterpunkten ($N=100000$).
data = np.genfromtxt("./Achterbahn_med.dat")
plt.xlabel(r"$\rm t$")
plt.ylabel(r"$\rm x \,\, \dot{x}$")
plt.plot(Loes1.t, Loes1.y[0],c="blue", label=r"$\rm x(t)$");
plt.plot(data[:,1], data[:,2],c="green", linestyle=":", label=r"$\rm x (C++)$");
plt.plot(Loes1.t, Loes1.y[1],c="red", label=r"$\rm \dot{x}(t)$");
plt.plot(data[:,1], data[:,3],c="black", linestyle=":", label=r"$\rm \dot{x} (C++)$");
plt.legend(loc='lower right',fontsize=16);
Die obere Abbildung zeigt, dass die Ergebnisse des C++ Programms gut mit der Python Simulation übereinstimmen. Um dies genauer zu untersuchen, stellen wir in der folgenden Abbildung die absoluten Fehler der zeitlichen Verläufe der berechneten Größen dar ($\Delta x$ und $\Delta \dot{x}$).
N=100000
t_val = np.linspace(0, t_end, N+1)
Loes1 = solve_ivp(DGLs, [0, t_end], initialval, t_eval=t_val, rtol=fehler, atol=fehler)
plt.xlabel(r"$\rm t$")
plt.ylabel(r"$\rm \Delta x \,\, \Delta \dot{x}$")
plt.plot(Loes1.t, Loes1.y[0]-data[:,2],c="blue", label=r"$\rm \Delta x(t)$");
plt.plot(Loes1.t, Loes1.y[1]-data[:,3],c="red", label=r"$\rm \Delta \dot{x}(t)$");
plt.legend(loc='lower right',fontsize=16);
Die obere Abbildung zeigt, dass sich in einigen Bereichen noch große Abweichungen zwischen den Ergebnissen des C++ Programms und der Python Simulation ergeben. Wir Erhöhen deshalb nochmals die Anzahl an zeitlichen Gitterpunkten ($N=500000$).
N=500000
t_val = np.linspace(0, t_end, N+1)
Loes1 = solve_ivp(DGLs, [0, t_end], initialval, t_eval=t_val, rtol=fehler, atol=fehler)
data = np.genfromtxt("./Achterbahn_high.dat")
plt.xlabel(r"$\rm t$")
plt.ylabel(r"$\rm x \,\, \dot{x}$")
plt.plot(Loes1.t, Loes1.y[0],c="blue", label=r"$\rm x(t)$");
plt.plot(data[:,1], data[:,2],c="green", linestyle=":", label=r"$\rm x (C++)$");
plt.plot(Loes1.t, Loes1.y[1],c="red", label=r"$\rm \dot{x}(t)$");
plt.plot(data[:,1], data[:,3],c="black", linestyle=":", label=r"$\rm \dot{x} (C++)$");
plt.legend(loc='lower right',fontsize=16);
plt.xlabel(r"$\rm t$")
plt.ylabel(r"$\rm x \,\, \dot{x} \,\, h$")
plt.plot(Loes1.t, Loes1.y[0]-data[:,2],c="blue", label=r"$\rm x(t)$");
plt.plot(Loes1.t, Loes1.y[1]-data[:,3],c="red", label=r"$\rm \dot{x}(t)$");
plt.legend(loc='lower right',fontsize=16);
