Problemstellung:
Man wirft einen Ball senkrecht nach oben bzw. lässt diesen aus einer gegebenen Höhe fallen. Die Dämpfung beim Aufprall sei der Einfachheit wegen hierals linear angenommen. In unserem ist unser Dämpfungsfaktor k=0.9.
Gesucht ist nun die Bewegung des Balles vom Zeitpunkt t0=0 bis zu einem gegebenen Zeitpunkt Tende.
Vorbereitungen:
Um ein klein wenig Physik kommt man hier nicht herum:
Numerische Lösung:
Die Lösung des Systems lassen wir einen Algorithmus übernehmen. In unserem Fall ist es das Runge-Kutta-Verfahren (obwohl hier sicherlich auch jedes andere Verfahren funktionieren würde). Auch wenn die Implementation dieses Algorithmus kein großes Problem darstellt, lassen wir das in diesem Fall von
erledigen.Benötigte Pakete:
Benötigt wird Python mit numpy, scipy und matplotlib. Installiert werden diese Paket am einfachsten mit
oder man ladet sich direkt , welches neben den Paketen unter anderen auch Spyder mitbringt. Insgesamt ist Pythonxy sehr zu empfehlen.Matplotlib wird hierbei nicht zwingend benötigt, ist aber wirklich ratsam, wenn man seine Ergebnisse grafisch darstellen möchte.
Das Programm:
Benötigt werden folgende includes:
Code:
from scipy.integrate import ode import matplotlib.pyplot as plt import numpy as np
pyplot aus dem matplotlib Packet wird benötigt, um das Ergebnis später zu visualisieren.
numpy wird hier benötigt, um später ein Array zu erzeugen (kann auch drauf verzichtet werden).
Nun definieren wir uns unser System:
Code:
def f(t, y):
return [y[1], -9.81]
Die Rückgabe ist gleichzusetzen mit Gleichung (3), wobei der Vektor in Python als Liste behandelt wird, und y[1] unser x2 aus Gleichung (3) ist (liegt daran, dass in Python von der 0 beginnend indiziert wird).
Nun zur Hauptroutine, erst ein Teil vom Code, dann die Erklärung:
Code:
solver = ode(f).set_integrator('vode')
solver.set_initial_value([10, 0], 0)
dt = 0.01
end = 10
position = []
time = []
In der zweiten Zeile werden die Anfangswerte gesetzt. Differentialgleichungen, die nicht nur in der Theorie benötigt werden, kommen fast immer mit einem Anfangswerptoblem. Dabei ist das hier so zu lesen:
Die Werte [10, 0] sollen zu Zeitpunkt 0 von dem Gleichungssystem erfüllt werden. Ein Blick zurück auf Gleichung 3 verrät uns, dass der erste Wert im Vektor dem Ort, der zweite der Geschwindigkeit entspricht. Demnach startet unsere Simulation mit einem Ball der aus 10 Metern ohne Anfangsgeschwindigkeit fallen gelassen wird.
dt steht hier für unseren Zeitschritt, end entspricht unserer Endzeit der Simulation. Demnach simulieren wir den Fall des Balles für 10 Sekunden in 0.01 Sekunden Schritten.
Die zwei Listen position und time dienen lediglich dazu, die Zeiten bzw die Positionen des Balles zu speichern.
Weiter gehts:
Code:
while solver.successful() and solver.t < end: solver.integrate(solver.t + dt) position.append(solver.y[0]) time.append(solver.t) if solver.y[0] <= 0: solver.set_initial_value([solver.y[0], -0.9*solver.y[1]], solver.t)
Also läuft die auf Zeile 1 folgende Schleife bis der Solver nicht mehr lösen kann oder unsere Simulationsende erreicht wurde.
Die Zeile darauf wird die Gleichung für den nächsten Zeitpunkt t + dt gelöst (dazu am Schluss noch eine kleine Anmerkung) und in den folgenden 2 Zeilen die Ergebnisse in die vorher genannten Listen verfrachtet.
Die darauf folgende Bedingung
Code:
if solver.y[0] <= 0:
Code:
solver.set_initial_value([solver.y[0], -0.9*solver.y[1]], solver.t)
Um zum Ende zu kommen, wollen wir das Ergebnis nun auch grafisch darstellen:
Code:
plt.plot(time, position)
plt.grid()
plt.xlabel('time (s)', fontsize=17)
plt.ylabel('height (m)', fontsize=17)
plt.title('Bouncing ball', fontsize=23)
plt.show()
Code:
plt.show()
Den kompletten Code gibt es noch einmal hier:
Viel Spaß beim rumspielen damit und falls jemand Fehler findet: immer her damit.
