Ziel dieses Projekts war es, eine funktionierende n-Teilchen Gravitationssimulation in Python zu programmieren, in der jedes Objekt mit jedem wechselwirkt. Dazu mussten wir uns einerseits in Python mit zugehörigen passenden Datenstrukturen (Klassen) und Funktionen, wie auch dem numerischen Lösen von Differentialgleichungen 2. Ordnung, auseinandersetzen.

Fabian Heisinger
Erik Kriesch
Thomas Oppermann
Ludwig Skuras

A Gentle Introduction To Numerical Python
Computational Many Particle Physics
Computational Physics Python
Computing Particle Properties
Gravitational N-Body Simulation
Numerische Simulation 001
Gilbert Strang: Wissenschaftliches Rechnen, ISBN 978-3-540-78494-4
Kurt W. Smith: Cython

Außerdem natürlich auch die Numpy- und Pygame-Dokumentationen

Für das Projekt haben wir uns als 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.

Die Grundlage der Berechnungen bildet das Gravitationsgesetz: $$\vec{F} = G\cdot\frac{m_1\cdot m_2}{r^3}\cdot\vec{r}$$ damit folgt für die Gesamtkraft auf de i-te Teilchen bei $n$ Teilchen: $$\vec{F}_i=G\cdot\sum_{j=1;j\neq i}^{n}\frac{m_i\cdot m_j}{r_{ij}^3}\cdot\vec{r}_{ij}$$ Die Gravitationskonstante ist mit $G = 6,673\cdot 10^{-11} m^3\cdot kg^{-1}\cdot s^{-2}$ bekannt und die anderen Größen sind als Eigenschaften der Objekte gespeichert. Daraus wird der Abstand zwischen zwei Körpern über die Distanz zwischen ihren Postionen im Raum bestimmt. Aus diesen Werten lässt sich die Kraft von einem Objekt zu allen anderen Objekten berechnen und in eine Matrix eingetragen. Dies wird für alle Objekte wiederholt, bis die Matrix voll ist. Um die Gesamtkraft auf ein Objekt zu bestimmen müssen alle Kräfte die in der passenden Zeile der Matrix stehen aufsummiert werden. Daraus resultieren die auf die Objekte wirkenden Kräfte für einen bestimmten Zeitschritt. Je nach verwendetem Verfahren ist die Berechnung über diesen Zeitschritt einem numerischen Fehler unterlegen. Im Folgendem werden das Euler-Verfahren und die Leapfrog-Methode vorgestellt,da sie beide zur Benutzung kamen bzw. kommen.
Integration, Euler, Leapfrog, Numerik

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. Eine 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 (+$\infty$), 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ßes Koordinatensystem geschrieben. Eine Bewegung der Teilchen im Graphen wird über die Darstellung eines neuen Graphen mit anderen Positionswerten der Teilchen bewerkstelligt (Es wird ein Graph über den anderen gelegt). 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 vorherigen 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]

Im weiterem verlauf wurde diese Verfahren durch die Leapfrog Methode ausgetauscht, dieses ist ähnlich einfach wie die Eulersche, aber von dritter Ordnung, d.h. genauer.
Bei Leapfrog findet die Berechnung von Ort und Beschleunigung, sowie der Geschwindigkeit einen halben Zeitschritt voneinander entfernt statt. Zuerst wird der Ort während des halben Zeitschrittes berechnet, danach die Geschwindigkeit durch die Beschleunigung aktualisiert und zuletzt der Ort durch die neue Geschwindigkeit im nächst ganzen Zeitschritt bestimmt. $$\vec{x} \left( t + \frac{1}{2} \right) = \vec{x} \left( t \right) + \frac{1}{2} \vec{v} \left( t \right) dt $$ $$\vec{v} \left( t + \frac{1}{2} \right) = \vec{x} \left( t \right) + \frac{F}{m} dt $$ $$\vec{x} \left( t + 1 \right) = \vec{x} \left( t + \frac{1}{2} \right) + \frac{1}{2} \vec{v} \left( t + \frac{1}{2} \right) dt$$ Im Code wird dies noch auf die Achsen aufgeteilt.

	def move(self):
		self.x += 0.5 * self.velocity[0] * dt
		self.y +=  0.5 * self.velocity[1] * dt
		self.velocity[0] += dt * self.force[0]/self.mass
		self.velocity[1] += dt * self.force[1]/self.mass
		self.x += 0.5 * self.velocity[0] * dt
		self.y +=  0.5 * self.velocity[1] * dt

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.

Einer der wesentlichen Aspekte von Python ist die Tatsache, dass in Python eine dynamische Typisierung zu Grunde liegt. Das bedeutet, dass man Variablen, im Gegensatz zu Sprachen wie C oder Fortran, nicht explizit einen Typ zuweist, die Typisierung nimmt der Python-Interpreter selber vor. Einerseits ist das sehr nützlich, da der Quellcode dadurch leichter verständlich ist und Variablen auch einfach umgeändert werden können. Der Nachteil besteht allerdings darin, dass die Erschliesung des Typs zur Laufzeit des Programms geschieht, dadurch ist Python, vorallem bei vielen numerischen Berechnungen, wie wir sie als Kraftberechnungen haben, sehr viel langsamer als Programmiersprachen, die statische Typdeklerationen benutzen und kompiliert werden.

Als reines Python-Programm, läuft die Simulation bis für ca. 100 Teilchen recht flüssig. Da wir in der Funktion get_force_matrix $n^2$ Einträge haben, müssten wir theoretisch auch $n^2$ Einträge berechnen.
Da wir aber wissen, dass alle Hauptdiagonaleinträge der Matrix $0$ sind, verringert sich die Anzahl der Berechnungen auf $n^2-n=n(n-1)$. Zusätzlich wissen wir, dass die Kraft, die von Teilchen $i$ auf Teilchen $j$ wirkt, gleich der negativen Kraft von $j$ auf $i$ entspricht. Die Matrix ist also schiefsymetrisch und daher reicht es die obere oder untere Dreiecksmatrix zu berechen, die restlichen Werte können durch eine einfach Spiegelung bestimmt werden, der Rechenaufwand verkürzt sich damit weiter auf $\frac{1}{2}n(n-1)\approx \mathcal O(n^2)$
Trotz unserer Vereinfachungen wächst der Rechenaufwand also annähernd quadratisch, so müssen bei 10 mal so vielen Teilchen ungefähr 100 mal so viele Rechnungen durchgeführt werden. Diese Kombination aus exponentiell-wachsendem Rechenaufwand und der dynamischen Typisierung, welche natürlich auch immer laufzeithungriger, da immer mehr Variablen erschlossen werden müssen, wird, lassen das Programm also schnell an seine Grenzen stoßen.

Es muss also eine Möglichkeit gefunden werden, damit die Simulation auch für große $n$ flüssig läuft.
Die Antwort auf dieses Problem lautet Cython. Das Modul Cython erlaubt es, gewisse Programmteile als C-Code zu schreiben und vor der Laufzeit zu kompilieren, um das gesamte Programm, durch eben kompilierte Funktionen und damit verbundender statischer Typisierung, um ein vielfaches zu beschleunigen.

Als erstes ist es wichtig, herauszufinden, welche Funktionen besonders laufzeitintensiv sind und eine Implementierung in Cython sich somit lohnen würde. Das macht man zum Beispiel mit cProfile:

python -m cProfile -o analyse.stat run_project.py

Unser Python-Programm wird mit dem Parameter cProfile aufgerufen und während der Laufzeit analysiert. Die Analysedaten werden in der Datei analyse.stat gespeichert und können mit folgendermaßen Python-Skript ausgegeben werden:

import pstats
 
p = pstats.Stats("analyse.stat")
p.strip_dirs().sort_stats("time").print_stats(20)

Dieser Aufruf erzeugt eine Liste von den 20, am häufigsten genutzten Funktionen unseres Programms. Der Output kann z.B so aussehen:

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
  5781600   23.011    0.000   52.819    0.000 functions.py:25(get_force)
  5781600   22.357    0.000   22.357    0.000 functions.py:68(get_abstand)
     1168   16.701    0.014   69.762    0.060 functions.py:41(get_force_matrix)
   116800    8.430    0.000    8.430    0.000 {method 'reduce' of 'numpy.ufunc' objects}
  5898674    7.619    0.000    7.619    0.000 {numpy.core.multiarray.array}
        1    2.197    2.197   84.441   84.441 run_project.py:1(<module>)
     1168    1.132    0.001    1.132    0.001 {pygame.display.update}
   116800    0.801    0.000    0.801    0.000 classes.py:20(move)
     1168    0.527    0.000    0.527    0.000 {method 'fill' of 'pygame.Surface' objects}
   116800    0.245    0.000    8.809    0.000 fromnumeric.py:1730(sum)
   116800    0.233    0.000    0.576    0.000 graphic_functions.py:4(draw)
   117968    0.171    0.000    0.171    0.000 {pygame.draw.circle}
   233600    0.145    0.000    0.145    0.000 graphic_functions.py:15(norm_pos)
   116800    0.133    0.000    8.942    0.000 functions.py:59(total_force_per_part)
   123061    0.078    0.000    0.078    0.000 {isinstance}
     1168    0.068    0.000    0.068    0.000 {numpy.core.multiarray.zeros}
   116800    0.059    0.000    8.489    0.000 _methods.py:31(_sum)
   116800    0.051    0.000    0.051    0.000 graphic_functions.py:18(norm_rad)
        2    0.043    0.021    0.095    0.048 __init__.py:45(<module>)
        1    0.040    0.040    0.043    0.043 case.py:1(<module>)

Die Zahl ncalls ist die Anzahl der Aufrufe, tottime die insgesamte Laufzeit der Funktion oder des Moduls und percall gibt die Zeit pro Aufruf an.
Wie zu erwarten, sind vor allem die Funktionen, die die Grundlage unserer Simulation bilden, für über 70% der Laufzeit verantwortlich (also get_way, get_force und get_force_matrix). Diese Funktionen in Cython umzuschreiben ist daher naheliegend.

Ein großer Vorteil von Cython ist, dass es möglich ist in sehr Python nahem Code zu schreiben. Ein kurzes Beispiel anhand der Fibonacci-Folge soll das verdeutlichen:

<code c>

double cfib(int n){

  double a = 0;
  double b = 1;
  double tmp;
  int i;
  
  for(i = 0; i < n; i++){
      tmp = a;
      a = a + b;
      b = tmp;
  }
  return a;

} <\code>

Hier geht es zu den Protokollen, die Einblick in den Projektverlauf/Stand und die einzelnen Labor- und Blocktermine geben.

Hier geht es zu verschiedenen Testprogrammen, mit denen wir vereinfacht geprüft haben, ob alles richtig gerechnet wird oder welche Methode zur Berechnung besser geeignet ist.

Hier geht es zu unserem Github-Repository, dort findet ihr unseren gesamten Quellcode, einmal in puren Python und einmal mit Cython-Verbesserungen