Mathesis Wiki

01.06.17:

1. Am Anfang stand die Frage im Raum ob eine bereits geschriebene Physik-Engine genutzt wird oder ob wir selbst eine schreiben. Wir haben uns dazu entschieden selbst eine Engine zu schreiben damit wir ein Verständnis vom Ablauf der Berechnungen haben und uns nicht nachträglich in einer Blackbox einarbeiten müssen.
   
2. Wir haben damit begonnen eine Normalverteilung von Massenpunkten in einer zwei dimensionalen Ebene zu erzeugen. Aus der Verteilung wird die Kraft von jedem Teilchen auf jedes andere berechnet und in eine Matrix eingetragen.

# Normalverteilung in der x-y-Ebene
x = np.random.normal()
y = np.random.normal()

08.06.17:

1. Wir beschäftigen uns mit dem Aufsummieren der Kräfte bezüglich eines Objektes. Da wir bis jetzt alle Kräftebeziehungen in einem Tensor gespeichert haben, gibt es Probleme, die richtigen Daten zu kombinieren. Nachdem wir die richtige Bezeichnung gefunden hatten konnten wir die richtigen Werte aud´s dem Tensor extrahieren.
2. Wir befassen uns nun mit dem physikalischen Prozess der Beschleunigung. Mit Hilfe der Kraft, der Masse und dem Anfangspunkt, müssen wir den Aufenthaltsort nach einem Zeitschritt x feststellen.
3. Lösung: Mit Hilfe des Euler-Verfahrens schreiben wir eine Funktion die uns schlussendlich realistische Werte liefert.
4. Wir kümmern uns nun um eine erste grafische Umsetzung. Das Problem besteht darin die Teilchen nicht bei jedem Schritt neu zu plotten, sondern eine „fließende“ Darstellung zu schaffen.
5. Lösung: Die Position der Teilchen im Graphen wird geupdatet nachdem mit dem Euler-Verfahren die neuen Position der Teilchen berechnet werden.

15.06.2017

1. Am Anfang gaben wir den Teilchen die neuen Attribute: Masse, Volumen und Dichte. Hierbei verglichen wir Daten von reellen Asteroiden um ungefähre Bereiche für die Zufallsverteilung zu erhalten. Nachdem wir nun unterschiedlich große, schwere und dichte Teilchen hatten, ging es an die Feinarbeit.
2. Zu nächst erkannten wir in der Simulation das Problem, dass Teilchen, wenn sie sich zu nah kamen, unkontrolliert aus dem Bild geschleudert wurden. Dies rührt daher, dass nach dem Gravitationsgesetz, die Kraft für einen unendlich kleinen Abstand unendlich wird. Als ersten Lösungsansatz begrenzten wir die Kraftberechnung auf einen minimalen Abstand. Falls der Abstand kleiner wurde, wurde mit dem minimalen Abstand weitergerechnet. Dies sah um einiges besser aus.
3. Nun machten wir uns an das Problem des Verschmelzens verschiedener Körper. Hierzu ignorierten wir vorerst die Möglichkeit eines Abprallens oder eines Zersplitterns bei einer Kollision. Wir schrieben eine Funktion die die Massen und Volumina der Teilchen addiert und dem ersten Teilchen zuschreibt. Nun gab es das Problem das zweite Teilchen zu löschen. Durch Stefans Hilfe fanden wir heraus, dass ein Listenobjekt zu löschen Zeitaufwendig ist. Statt das Teilchen zu löschen, gaben wir den Teilchen ein Attribut „aktiv/inaktiv“. So wurde das Teilchen nicht gelöscht aber zählt in der weiteren Berechnung nicht mit.

22.6.2017

1. Aufteilung in zwei seperate Aufgaben. Zum einen die Einarbeitung in PyGame um die grafische Darstellung zu realisieren und von Graphen als Darstellung wegzukommen. Zum anderen die Betrachtung des bisher vernachlässigten Drehimpulses. Bei der jetzigen zwei-dimensionalen Darstellung und Berechnung der Kräfte mit Massenpunkte ist der Drehimpuls nicht essentiel für genaue Berechnung und Darstellung.

29.07.2017

1. Weitere Einbindung (Anpassung der Daten) von PyGame für die Darstellung.
2. Behandlung des Drehimpuls und einer möglich Einbindung in die Berechnung der Kräft, Positionen, Geschwindigkeiten und Absorptionen und in einander.
3. Die Move- und Velocity-Funktionen in die Klasse Teilchen verschoben. Macht den Code übersichtlicher und das aufrufen der Funktionen einfacher im Quellcode.
4. Einarbeitung in das Leapfrog-Verfahren zur Berechnung der Kräfte, Position und Geschwindigkeit. Es ist eine Ordnung besser als das Euler-Verfahren in der Genauigkeit.

06.07.2017

1. Überprüfung der Physik-Berechnung am „Spezialfall“ des Kreisens eines Planeten um einen Stern (für diesen Zweck wurde die Erde und die Sonne als Referenz genommen)
2. Normalverteilung von Punkten auf einer Kreisbahn um einen Punkt

13.07.2017

1. Kreisbahnverteilung, Impulsberechnung berichtigt, Verschmelzung umgeschrieben
2. Strukturierung in Scripte für Variablen, Methoden, Klassen und Darstellung

20.07.2017

1. Dynamische Farbskala in Abhängigkeit von Masse (leichtestes Teilchen bis theoretisch schwerst möglichen Teilchen)
2. Loggen der Berechnungszeiten aller Funktionen und Darstellung. Resultat: get_force, get_force_matrix und total_force_per_part sollten für eine Beschleunigung in Cython geschrieben werden.

23.08 - 25.08.2017

1. Funktionen get_force, get_force_matrix und total_force_per_part wurden in Cython geschrieben.
2. Kontrolle der Teilchenbewegung um ein Zentralgestirn. Es wurde festgestellt, dass die Teilchen zu schnell abdriften (selbst mit einer Berechnungsungenauigkeit). Darauf hin wurde die Kraftberechnung dahingehend angepasst, dass eine Wechselwirkung zwischen Stern und Teilchen ohne Eintrag in die Liste der Teilchen möglich ist.
3. Definieren von festen Regeln für den Zusammenstoß von Teilchen und wie sich verhalten (elastische Stoß, plastischer Stoß/Verschmelzen und Elemenierung), basierend auf klaren physikalischen Gesichtspunkten.
4. Weitere Ausarbeitung der Farbzuweisung, mit einem feineren Verlauf und der Möglichkeit schnell zwischen dynamischen und statischen Grenzen zu wechseln.

Für das Projekt hat sich die die Gruppe gegen die Verwendung einer bereits fertigen Physik-Engine entschieden. Dagegen sprach vor allem die Einarbeitung in die verwendete Software, da sie sonst eine Blackbox ist.
Die selbst geschriebene Engine besitzt mehrere Klassen von Objekten und Funktionen auf die sie aufbaut und die Berechnungen ausführt.

Von Beginn an war die Klasse Teilchen definiert. Sie erzeugt Teilchen die in der xy-Ebene normalverteilt sind und weist ihnen weitere Attribute zu. Alle Teilchen werden in eine Liste eingetragen die auf die einzelnen Teilchen und ihre Attribute verweist. Das Eintragen der Liste ist aber nicht Bestandteil der Klasse.

# Die Klasse erstellt ein Teilchen als Massenpunkt im Raum mit den Attributen
# x-Position, y-Position, Masse, Geschwindigkeit (Vektor) und Kraft (Vektor).
class Teilchen(object):
	def __init__(self,m = 10**12):
		self.x = 100*np.random.normal()
		self.y = 100*np.random.normal()
		self.mass = m
		self.velocity = np.array([0.,0.])
		self.force = np.array([0.,0.])

Mit dem Fortschreiten des Projekts wurde die Klasse um die Attribute Dichte, Volumen, Radius und Impuls ergänzt. Die Masse ist aus einer Normalverteilung zufällig gewählt und beeinflusst zusammen mit der Masse das Volumen und den Radius des Teilchens.

# Der Code ist in die Klasse ergänzt worden.
		self.density = np.random.normal(950,50)
		self.volume = self.mass/self.density
		self.radius = np.cbrt(3*self.volume/(4*np.pi))
		self.momentum = [0.,0.]

Die Berechnung der auf die Teilchen wirkenden Kräfte erfolgt aus dem Gravitationsgesetz. Ein für die Anwendung geeignetes Berechnungsmethode ist das Euler-Verfahren. Für die Berechnung der Kraft muss der Abstand der Teilchen zueinander bestimmt werden. Mit dem bekannten Abstand zueinander und den Massen der Teilchen lässt sich die Kraft nach dem Gravitationsgesetz berechnen. Für den effizienten Umgang mit den Daten werden Kräfte in einer oberen rechten Dreiecksmatrix eingetragen und an der Diagonalen gespiegelt. Aus einer solchen Matrix lässt sich die Gesamtkraft auf ein Teilchen einfach bestimmen (nicht im Ausschnitt dargestellt).

# Aus der x- und y-Position zweier Teilchen wird der Abstand berechnet.
def get_abstand(part_1,part_2):
	return np.sqrt((part_1.x-part_2.x)**2 + (part_1.y-part_2.y)**2)
 
# Die Kraft wird aus den Massen und den Abstand der Teilchen
# komponentenweise berechnet und als Vektor ausgegeben.
def get_force(part_1,part_2):
	l = get_abstand(part_1,part_2)
	if l != 0:
		f_x = const.G*(((part_1.mass * part_2.mass)/l**3)
		*(part_2.x - part_1.x))
		f_y = const.G*(((part_1.mass * part_2.mass)/l**3)
		*(part_2.y - part_1.y))
		f_ges = [f_x,f_y]
		return np.array(f_ges)
	else: # Kann keine Kraft aus sich selbst berechnen
		return np.array([0.,0.])
 
# Die Kraft wird für jedes Teilchen auf jedes Teilchen berechnen.
# Wenn die Kraft schon (in umgekehrter Richtung) berechnet worden 
# ist wird auf die Berechnung verzichtet um Leistung einzusparen.
# Anschließend werden die nicht berechnet Werte aus berechneten
# Werten gespiegelt.
def get_force_matrix(liste_part):
	temp = init_matrix(len(liste_part),object)
	for i,part_1 in enumerate(liste_part):
		for j,part_2 in enumerate(liste_part):
			if i < j:
				temp[i,j] = get_force(part_1,part_2)
				temp[j,i] = -temp[i,j]
			elif i == j:
				temp[i][j] = np.array([0.,0.])
	return temp

Die Berechnung der Kraft ist bis auf geringe Abstände relativ genau. Für kleine Abstände wachsen die Kräfte gegen Unendlich (+∞), da die Engine mit den Teilchen als Massenpunkte rechnet und somit Radien und Volumina vernachlässigt. Zum umgehen dieses Problem muss eine Ausnahme für geringe Abstände geschrieben werden.

Für eine kleine Anzahl von Teilchen ist das Euler-Verfahren genau und schnell genug. Ab ca. 75 Teilchen wird die Berechnung langsamer. Desweiteren kommt hinzu das mit zunehmender Simulationsdauer die Ungenauigkeiten in der Berechnung zunhemen.
Soll die Simulation für längere Zeiten genau sein muss auf die Verwendung des Leapfrog-Verfahren gegriffen werden. Dies führt auch zu einigen Änderungen im Code, da die Berechnung anders abläuft als bei Euler. Als Lösung dafür wurde das Euler-Verfahren durch das Leapfrog-Verfahren ersetzt. Dazu muss die move-Funktion in die Klassen der Objekte geschrieben werden.

Eine Darstellung der Teilchen und ihrer Bewegungen ist anfänglich über die Ausgabe der Kraftmatrix in die Kommandozeile gemacht worden. Da dies wenig intuitiv ist wurde auf die Matplot-Bibliothek zugegriffen. Die Darstellung erfolgte nun über die Übergabe der Positionen in der xy-Ebene und der Aktualisierung des Graphen. Im Hintergrund wird weiterhin mit dem Euler-Verfahren die Änderung von Position, Kraft und Geschwindigkeit in diskreten Zeitschritten berechnet.
Eine Darstellung der Teilchen ist auf vielen wegen möglich. Die wohl einfachste Lösung (weil unkompliziert), die auch anfänglich verwendet wurde, ist die Ausgabe eine Kraftmatrix in die Kommandozeile bzw. in das Terminal. Es eignet sich gut um die Werte ungefähr auf Plausibilität zu prüfen und eine grobe Vorstellung von der Bewegen zu bekommen, ist aber davon abgesehen wenig intuitiv. Eine schönere und relativ unkomplizierte Darstellung ist mit Graphen möglich. Dafür wird die Matrix mit den Teilchenpositionen ausgelesen und in ein großen Koordinatensystem geschrieben. Eine Bewegung der Teilchen im Graphen wird über die Darstellung eines neuen Graphen mit anderen Positionswerten der Teilchen bewerkstelligt. Der folgende Ausschnitt zeigt die Aktualisierung eines Matplot-Graphen. Das Beispiel ist nicht vollständig in dem benötigten Code.

# Für 1000 Zeitschritte wird die Simulation über einen Matplot-Graphen dargestellt.
for _ in xrange(1000):
	force_matrix = get_force_matrix(liste_teilchen)
	for a,i in enumerate(liste_teilchen):
		i.force = total_force_per_part(force_matrix,a)
	get_velocity(liste_teilchen)
	move(liste_teilchen)
	for i,p in zip(liste_teilchen,plots):
		p[0].set_data(i.x,i.y)
	fig.canvas.draw()
	plt.pause(0.002)

Die im Ausschnitt erwähnte move-Funktion berechnet die Teilchenbewegung für einen bestimmten Zeitschritt nach dem Euler-Verfahren folgender Maßen.

# Aus einer übergebenen Liste von Teilchen wird die neue Position aus
# der alten Position plus Zeitschritt mal Geschwindigkeit berechnet.
def move(liste_part):
	for i in liste_part:
		i.x += zeitschritt*i.velocity[0]
		i.y += zeitschritt*i.velocity[1]

Für eine noch ansprechende Darstellung lassen sich andere Bibliotheken und Programme zur graphischen Darstellung nutzen. Eine dieser Bibliotheken ist PyGame, was auch in diesem Projekt für die Darstellung gegen Ende genommen wurde. Die grafischen Darstellung mit PyGame erfordert eine Anpassung der Werte, da die Darstellung im Fenster nur auf ein Pixel genau sein kann. Dafür müssen die Positionswerte mit einem Faktor auf die Größe des Bildschirm angepasst werden und anschließend auf ganzzahlige Werte (integers) gerundet werden. Um die Genauigkeit der Berechnung nicht zu beeinträchtigen übernimmt das Python-Script für die Darstellung die Anpassung der Werte. Daraus resultiert, dass alle Teilchen in die Berechnung mit einbezogen werden, aber einige Teilchen außerhalb des Darstellungsfensters und dessen Ausschnitts liegen und nicht angezeigt werden.