Simulation von physikalischen Problemen (Python)

Da es bez�glich dessen 2 kleine Anfragen gab, trete ich das kleine Beispiel etwas aus und zeige wie man Schritt f�r Schritt zu einer numerischen L�sung eines physikalischen Problems kommt.

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:
[Only registered and activated users can see links. Click Here To Register...]

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 [Only registered and activated users can see links. Click Here To Register...] erledigen.

Ben�tigte Pakete:
Ben�tigt wird Python mit numpy, scipy und matplotlib. Installiert werden diese Paket am einfachsten mit [Only registered and activated users can see links. Click Here To Register...] oder man ladet sich direkt [Only registered and activated users can see links. Click Here To Register...], 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

ode (ordinary differential equations - (Gew�hnliche Differentialgleichungen)) wird zum tats�chlichen L�sen des Problems ben�tigt.
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]

Das t in der Signatur wird im sp�teren Verlauf von dem Solver ben�tigt, weil Funktionen des Systems eben doch auch direkt von t abh�ngen k�nnen (was in unserem vereinfachten Beispiel nicht der Fall ist).
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 Zeile 1 wird unser Objekt, dass sp�ter zum L�sen der Gleichung ben�tigt wird, erzeugt. Das Argument 'vode' teilt dem Objekt mit, den Runge-Kutta-Algorithmus zur L�sung zu verwenden.
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)

Unser Objekt Solver beinhaltet als Eigenschaft eine Variable t, welche den Aktuellen Zeitpunkt repr�sentiert.
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:

sollte selbsterkl�rend sein. Wenn unser Ball den Nulldurchgang durchquert (sprich der Ball ber�hrt den Boden) passiert folgendes:

Code:

solver.set_initial_value([solver.y[0], -0.9*solver.y[1]], solver.t)

Die Anfangswerte werden wie auch zu Beginn neu gesetzt, und zwar wird der Ort �bernommen (also befindet sich unser Ball am Anfang auf dem Boden) und als Geschwindigkeit wird die Aufprallgeschwindigkeit (solver.y[1]) ged�d�mpft mit -0.9 genommen. Das Minus wurde spendiert, weil der Ball die Geschwindigkeit nun in die genau gegengesetzte Richtung besitzt.

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

Der Code sollte sich wie ein Buch lesen k�nnen. Wir tragen die Zeit �ber der Position auf und verpassen der Grafik noch etwas aussehen, in dem wir Titel und Achsenbeschriftung hinzuf�gen. Zu guter letzt wird via

Code:

plt.show()

die Grafik mithilfe von matploblib angezeigt:
[Only registered and activated users can see links. Click Here To Register...]

Den kompletten Code gibt es noch einmal hier:

Spoiler

Viel Spa� beim rumspielen damit und falls jemand Fehler findet: immer her damit.

Quote:

Ein gutes Tutorial, auch gut um jetzt zB. Grafiken auf einer Webseite anzuzeigen.

Das steht hier nicht im Vordergrund. Die grafische Darstellung bietet nur nochmal eine �berpr�fungsm�glichkeit. Hier gehts nicht ums plotten, sondern um die numerische L�sung von Differentialgleichungssystemen.

Quote:

anstatt jetzt zB. wxMaxima(welches darauf ausgelegt ist) zu verwenden um sowas zu l�sen?

Nein. wxMaxima ist ein Computeralgebrasystem (wie auch Maple oder Mathematica). Das benutzt man normalerweise, um Probleme analytisch zu l�sen. Das ist nicht immer m�glich (bzw. ist das in der Realit�t fast nie m�glich). Kann mir sehr gut vorstellen, dass auch wxMaxima numerische Methoden besitzt, dass diese einfacher zu benutzen sind als Python oder Matlab bezweifel ich aber. Du kannst mir gerne zeigen, dass es einfacher geht.

Quote:

wxMaxima ist leicht installiert

Trifft auch auf Python zu.

Quote:

keine zus�tzlichen Libs werden ben�tigt

Wird alles mitinstalliert, wenn man sich f�r pythonxy entscheidet.

Quote:

es ist cross platform

Ist Python ebenso.

Quote:

man lernt es schnell und man hat das Gef�hl man programmiert

Bei Python hat man nicht nur das Gef�hl, dass man programmiert - man tut es.

Quote:

da man Bl�cke(aka. Methoden) erstellen kann, Schleifen, if-Anweisungen, usw. verwenden kann.

Absolute Basics in Python.

Quote:

Grafiken kann man ebenfalls erstellen sowohl als auch diese lokal zu speichern. (in 2D oder 3D)

Absolut auch mit Python machbar. Matplotlibs 3D Module sind afaik noch experimental - haben aber bei ersten anl�ufen gut funktioniert.
Wem das noch nicht reicht, kann zu Matlab (wunderbare M�glichkeiten in 3D zu plotten und als Student oft f�r ca. 20€ zu bekommen) oder der kostenlosen Alternative Octave greifen.