Beispiel Projekt: Die Achterbahn
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. Der Einstiegspunkt in die Wagen befindet sich bei $x=0, \, y=0$ und die Achterbahn ist so konstruiert, dass sie die zwei Wagen zu unterschiedlichen Anfangszeitpunkten $\tau_i$ in die Achterbahnschiene katapultiert, wobei die Anfangsgeschwindigkeiten $\vec{v}_i(\tau_i)$ der Wagen (gemessen mit einem Sensor bei $x= 10 \, [cm]$) sich dabei unterscheiden. Jeder der Achterbahn-Wagen besitzt am vorderen und hinterem Bereich eine starke Metallfeder, die zum Einsatz kommt, wenn sich zwei der Wagen treffen. Es wird dabei angenommen, dass es sich bei einem Zusammenprall der Wagen um einen elastischen Stoß zweier Körper handelt.
Betrachten wir uns zunächst den Fall mit nur einem Wagen (Koordinaten: $\vec{r} = \left( x(t), \, y(t) \right) = \left( x(t), \, h(x(t)) \right)$, Masse $m$). Die Herleitung der Bewegungsgleichungen erfolgt am elegantesten mittels der Euler-Lagrange Gleichungen, bzw. mittels der Hamilton Theorie. Die Lagrangefunktion des Systems lautet:
$$
\begin{eqnarray}
& L = T - V \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.
Mögliche Teilaufgaben des Projektes
Mögliche Teilaufgaben des Projektes
Berechnen Sie zunächst auf analytischem Weg, mittels eines Jupyter Notebooks unter Verwendung der SymPy Bibliothek, die Euler-Lagrange Gleichungen der Bewegung eines Wagens auf der Achterbahn, wobei der Verlauf der Metallschiene der Achterbahn durch die Funktion $h(x) = e^{\frac{\sqrt{x^2}}{10}} \cdot \left( \hbox{sin}\left( \frac{x^2}{10} \right) \right)^2$ bestimmt ist. Schreiben Sie dann die Bewegungsgleichung in ein System von erster Ordnung Differentialgleichungen um.
Der Wagen ($m=100$ [Kg]) werde zur Zeit $t=0$ [s], mit einer Anfangsgeschwindigkeit von $\dot{x}(0) = 30.42$ [km/h] bei $x=0.01$ [m] in die Achterbahn katapultiert. Lösen Sie das System von Differentialgleichungen erster Ordnung zunächst mittels Python und visualisieren Sie die Simulation in einem Jupyter Notebook. Erhöhen Sie nun die Anfangsgeschwindigkeit auf $\dot{x}(0) = 60$ [km/h] und stellen beide numerischen Lösungen in einer Animation dar.
Lösen Sie nun das System von Differentialgleichungen erster Ordnung mittels eines C++ Programms und benutzen Sie das Runge-Kutta Ordnung vier Verfahren.
Vergleichen Sie beide Simulationen quantitativ.
Verallgemeinern Sie das C++ Programm so, dass mehrere Wagen auf der Achterbahn fahren können. Simulieren Sie ein System bestehend aus zwei wechselwirkungsfreien Wagen (siehe obige Animation) und benutzen Sie dabei für den ersten Wagen $i=0 \, , \, \tau_0 = 0 $ und $\dot{x}_0(0) = 16.45 \, [m/s]$ (blauer Wagen) und für den zweiten Wagen den ersten Wagen $i=1 \, , \, \tau_1=2 \, [s] $ und $\dot{x}_1(2) = 5.5 \, [m/s]$ (grüner Wagen).
Von den Euler-Lagrange Gleichungen zur Bewegungsgleichung¶
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
expand(L)
simplify(L)
ELG_1 = diff( diff(L, diff(x,t) ) ,t) - diff(L, x)
ELG_1.simplify()
Die DGL bestimmende Funktion lautet somit:
f=solve(ELG_1, diff(x,t,t))[0]
f
Wir exportieren uns für die anderen Teilaufgaben des Projektes die Funktion $f$ als C++ Export.
c_u_1, c_h, c_hs, c_hss = symbols('u_1, c_h, c_hs, c_hss', real = True)
f_c=f.subs(diff(h,x,x),c_hss).subs(diff(h,x),c_hs).subs(diff(h,x),c_hs).subs(diff(x,t),c_u_1)
f_c
from sympy.codegen.ast import Assignment
print(ccode(f_c))
Beispiel: $h(x) = e^{\frac{\sqrt{x^2}}{10}} \cdot \left( \hbox{sin}\left( \frac{x^2}{10} \right) \right)^2$¶
x_= symbols('x',real=True)
wert= symbols('wert')
h_ = exp(sqrt(x_**2/10)) * (sin(x_**2/10))**2
print(ccode(Assignment(wert, h_)))
h_
h_s = diff(h_,x_)
print(ccode(Assignment(wert, h_s)))
h_s
Den oberen Ausdruck für $\frac{d h (x)}{dx}$ werden wir in dem C++ Programm benötigen.
h_ss = diff(h_,x_,x_)
print(ccode(Assignment(wert, h_ss)))
h_ss
Die Diracsche Delta-Funktion kann man sich ja einfach selbst programmieren.
Eine andere Vorgehensweise wäre, die Funktion $h(x)$ schon beim Aufstellen der Bewegungsgleichungen festzulegen:
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
expand(L)
simplify(L)
Die Euler-Lagrangegleichungen lauten:
ELG = diff( diff(L, diff(x,t) ) ,t) - diff(L, x)
ELG.simplify()
Umschreiben in ein System von Differentialgleichungen erster Ordnung¶
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:
f=solve(ELG, diff(x,t,t))[0]
f
Diesen Ausdruck kann man dann wieder nach C++ exportieren (zuvor ist eine Variablenumbennenung nötig).
c_u_2 = symbols('u_2', real = True)
f_c=f.subs(diff(x,t),c_u_2).subs(x,c_u_1).subs(g,9.81)
f_c
print(ccode(Assignment(wert, f_c)))
Die obige Ausgabe verwendeten wir dann schliesslich im C++ Programm Achterbahn_1.cpp und Achterbahn_2.cpp.
Numerisches Lösen und Visualisierung der Simulation mittels Python¶
u_1 = Function('u_1')(t)
u_2 = Function('u_2')(t)
DGL=f.subs(diff(x,t),u_2).subs(x,u_1)
Eq(diff(u_2,t),DGL)
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
from scipy import integrate
DGL_L = lambdify((u_1,u_2,m,g), DGL)
Wir definieren uns das System der DGLs als eine Funktion
set_g=9.81 #Fallbeschleunigung g auf der Erde
set_m=100 # Masse des Wagens in Kg
def DGLs(u_vec,t):
u1, u2 = u_vec
du1_dt = u2
du2_dt = DGL_L(u1,u2,set_m,set_g)
return np.array([du1_dt,du2_dt])
und lösen das Anfangswertproblem im Bereich $t \in [0,15]$ unter Verwendung von $N=100000$ zeitlichen Gitterpunkten (mesh points). Wir wählen als Anfangsbedingungen die oben angegebenen ($m=100$ [Kg] , $v_x(0) = 30.42$ [Km/h] und $x(t=0)=0.01$ [m].
x_0 = 0.01
v_0 = 30.42*1000/(60*60)
print("Anfangsgeschwindigkeit in m/s", v_0)
N=100000
tend=15
tval = np.linspace(0, tend, N+1)
initialval = np.array([x_0,v_0])
Loes = integrate.odeint(DGLs, initialval, 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 x \,\, \dot{x}$")
plt.plot(tval, Loes[:, 0],c="blue", label=r"$\rm x$");
plt.plot(tval, Loes[:, 1],c="red", label=r"$\rm \dot{x}$");
plt.legend(loc='lower right',fontsize=16);
Wir entwerfen nun eine Visualisierung der Achterbahn.
import matplotlib.animation as animation
from IPython.display import HTML
import matplotlib.gridspec as gridspec
params = {
'figure.figsize' : [13,8],
# 'text.usetex' : True,
'axes.titlesize' : 18,
'axes.labelsize' : 15,
'xtick.labelsize' : 15 ,
'ytick.labelsize' : 15
}
matplotlib.rcParams.update(params)
h_L = lambdify((x), h)
h_L(Loes[0:10, 0])
fig = plt.figure()
ax = fig.gca()
set_frames=250
stepi=int(N/set_frames)
plot_max_y=np.max(h_L(Loes[:, 0])) + 0.5
plot_max_x=np.max(Loes[:, 0]) + 0.5
def init():
ax.scatter(Loes[0,0], h_L(Loes[0,0]), s=400, marker='o', c="red")
return fig,
def animate(i):
ax.cla()
ax.scatter(Loes[stepi*i,0], h_L(Loes[stepi*i,0]), s=400, marker='s', c="red")
ax.plot(Loes[:, 0],h_L(Loes[:, 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())
Wir erhöhen nun die Anfangsgeschwindigkeit auf $v_x(0) = 60$ [Km/h] und stellen uns beide numerischen Lösungen in einer Animation dar:
x_0 = 0.01
v_0 = 60*1000/(60*60)
print("Anfangsgeschwindigkeit in m/s", v_0)
N=100000
tend=15
tval = np.linspace(0, tend, N+1)
initialval = np.array([x_0,v_0])
Loes1 = integrate.odeint(DGLs, initialval, tval)
fig = plt.figure()
ax = fig.gca()
set_frames=250
stepi=int(N/set_frames)
plot_max_y = np.max(h_L(Loes1[:, 0])) + 0.5
plot_max_x = np.max(Loes1[:, 0]) + 0.5
plot_min_x = np.min(Loes1[:, 0]) - 0.5
x_L = np.linspace(plot_min_x, plot_max_x, 4000)
def init():
ax.scatter(Loes1[0,0], h_L(Loes1[0,0]), s=400, marker='s', c="green")
ax.scatter(Loes[0,0], h_L(Loes[0,0]), s=400, marker='o', c="red")
return fig,
def animate(i):
ax.cla()
ax.scatter(Loes1[stepi*i,0], h_L(Loes1[stepi*i,0]), s=400, marker='s', c="green")
ax.scatter(Loes[stepi*i,0], h_L(Loes[stepi*i,0]), 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())
Wir vergleichen nun die numerischen Berechnungen mit den Ergebnissen des C++ Programmes (Achterbahn.cpp mit dsolve_Sys.hpp). Wir vergleichen die Bewegung des roten Wagens und verwendeten $N=5000$ Zeit-Gitterpunkte.
data = np.genfromtxt("./Achterbahn_1.dat") # Einlesen der berechneten Daten von Achterbahn.cpp
params = {
'figure.figsize' : [10,7],
# '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 x \,\, \dot{x}$")
plt.plot(tval, Loes[:, 0],c="blue", label=r"$\rm x$");
plt.plot(data[:,1], data[:,2],c="black", label=r"$\rm x (C++)$");
plt.plot(tval, Loes[:, 1],c="red", label=r"$\rm \dot{x}$");
plt.plot(data[:,1], data[:,3],c="black", linestyle=":", label=r"$\rm \dot{x} (C++)$");
plt.legend(loc='lower left',fontsize=16);
Es wurde dann das C++ Programm auf zwei Wagen erweitert:
data = np.genfromtxt("./Achterbahn_2.dat") # Einlesen der berechneten Daten von Achterbahn_2.cpp
plt.xlabel(r"$\rm t$")
plt.ylabel(r"$\rm x \,\, \dot{x}$")
plt.title('Die Achterbahn (Achterbahn_2.cpp)')
plt.plot(data[:,1], data[:,2],c="blue", label=r"$\rm x \, , \,\, Wagen \, 1$");
plt.plot(data[:,1], data[:,3],c="red", label=r"$\rm \dot{x} \, , \,\, Wagen \, 1$");
plt.plot(data[:,1], data[:,4],c="blue", linestyle=":", label=r"$\rm x \, , \,\, Wagen \, 2$");
plt.plot(data[:,1], data[:,5],c="red", linestyle=":", label=r"$\rm \dot{x} \, , \,\, Wagen \, 2$");
plt.legend(loc='lower right',fontsize=14);
import matplotlib.gridspec as gridspec
params = {
'figure.figsize' : [13,8],
# 'text.usetex' : True,
'axes.titlesize' : 18,
'axes.labelsize' : 15,
'xtick.labelsize' : 15 ,
'ytick.labelsize' : 15
}
matplotlib.rcParams.update(params)
fig = plt.figure()
ax = fig.gca()
N=len(data[:,0])
set_frames=700
stepi=int(N/set_frames)
plot_max_y = np.max(h_L(data[:,2])) + 4
plot_max_x = np.max(data[:,2]) + 0.5
plot_min_x = np.min(data[:,2]) - 0.5
x_L = np.linspace(plot_min_x, plot_max_x, 4000)
def init():
for n in range(6):
ax.scatter(data[0,2], h_L(data[0,2]), s=400, marker='o', c="blue")
ax.scatter(data[0,4], h_L(data[0,4]), s=400, marker='o', c="green")
ax.plot(x_L,h_L(x_L), c="black")
def animate(i):
ax.cla()
ax.scatter(data[stepi*i,2], h_L(data[stepi*i,2]), s=400, marker='s', c="blue")
ax.scatter(data[stepi*i,4], h_L(data[stepi*i,4]), s=400, marker='s', c="green")
ax.plot(x_L,h_L(x_L), c="black")
ax.set_xlim(plot_min_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=40)
plt.close(ani._fig)
HTML(ani.to_html5_video())