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 zweiter Ordnung, d.h. genauer. Die Ordnung ist gleichzusetzen mit der höchsten Ableitung der DGL. Bei Euler wird nur der Ort betrachtet, innerhalb der DGL die Ausgangsfunktion, und die Geschwindigkeit, die Ableitung des Ortes. Während bei Leapfrog ebenfalls noch die Beschleunigung, Ableitung der Geschwindigkeit betrachtet wird.
Grundlegend sind numerische Integrationsverfahren Annäherungen an Funktionen mit kontinuierlichen, veränderliche Werten, da für Simulationen diskrete, feste Werte z.B. die Koordinaten, benötigt werden. Die jeweiligen Methoden zur Annäherung an den realen verlauf sind aber nicht perfekt, so kommt es zu Abweichungen. Erklärbar ist dieser Fehler durch den Zeitschritt. Anschaulich erklärbar ist dies durch Vergleich mit dem Integrieren, würde eine feste Unterteilung vorgenommen werden würde die Berechnete Fläche der angenäherten Stufen sichtbar von der tatsächlichen abweichen. Die Verkleinerung der Zeitschritte würde dabei den Fehler nicht einmal unbedingt verringern, da sich die kleineren Abweichungen aufsummieren würden. Ebenso bedeuten kleinere Zeitschritte auch mehr Berechnungen und somit höhere Laufzeiten für die Simulation. Somit bleibt für genauere Werte nur eine Genauere Annäherung bei beliebigen Zeitschritten. Die Genauigkeit dieser Verfahren steigt mit ihrer Ordnung, wobei dies gut durch das Taylor Polymon vorstellbar ist. Das Taylor Polymon erster Ordnung besteht aus einer Funktion am Entwicklungspunkt addiert mit der ersten Ableitung dieser, beim T.P. zweiter Ordnung wird die zweite Ableitung noch dazu addiert usw.. Beide liefern eine Annäherung des Funktionsverlaufes am Entwicklungspunkt, wobei das T.P. zweiten Grades genauer ist als das ersten Grades.
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. Dadurch wird der Ort genauer bestimmt, da der Fehler, welcher bei dem gewähltem Zeitschritt auftritt, durch die Aktualisierung verringert wird. $$\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 nahen Code zu schreiben. Ein kurzes Beispiel anhand der Fibonacci-Folge soll das verdeutlichen:

# Fibonacci in Python
def fib(n):
    a, b = 0, 1
    for i in range(n):
        a, b = a + b, a
    return a

// Fibonacci in 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;
}

# Fibonacci in Cython
def cyfib(int n):
    cdef double a = 0
    cdef double b = 1
    cdef int i
 
    for i in range(n):
        a, b = a + b, a
    return a

Wie man sieht, sind sich der Python- und Cython-Code relativ ähnlich, der C-Code ist hingegen umständlicher und benötigt eine zusätzliche Variable in der Schleife.
Der einzige Unterschied zwischen Python und Cython besteht, in diesem Beispiel, darin, dass in Cython der Typ der Variablen explizit deklariert wird, außerdem wird Cython-Code in .pyx nicht in .py Dateien gespeichert. Während der Python-Code nun von einem Python-Skript importiert oder direkt ausgeführt werden kann, muss für die Benutzung unseres Cython-Codes das Cython-Modul „pyximport“ importiert und installiert werden.

import pyximport; pyximport.install()
from time import time
from cyfib import cyfib
 
def fib(n):
    a, b = 0, 1
    for i in range(n):
        a, b = a + b, a
    return a
 
start = time()
x = fib(1000)
end = time()
print "erfolg python %f" % (end - start)
 
start = time()
y = cyfib(1000)
end = time()
print "erfolg cython %f" % (end - start)

Dieses Skript zeigt, wie man unsere einfache Cython-Datei in Python nutzt und direkt mit dem Python-Pendant vergleicht. Zuerst importieren wir das Modul pyximport und benutzen die Funktion install um Cython einzurichten. Danach können wir bequem unsere Datei mit dem Namen cyfib.pyx importieren und ausführen. Der time-Import erlaubt es uns die Zeit, die die jeweilige Funktion benötigt, um das 1000te Glied der Fibonacci-Folge zu berechnen. Wenn man das Skript ausführt, ist die Cython-Funktion mehr als 10-mal schneller als die reine Python Funktion.

erfolg python 0.000105
erfolg cython 0.000005

Der Grund, warum Cython so viel schneller ist, liegt darin, dass beim ersten Skriptaufruf der Cython-Code von Cython in direkten C-Code übersetzt und dann vom C-Compiler in Maschinencode kompiliert wird, welcher dann zur Verfügung steht. Aus dem Cython-Code wird also indirekt der C-Code aus dem Beispiel.
Es ist sehr beeindruckend, dass durch so eine kleine Änderung bereits ein großer Zeitgewinn stattfindet. Insbesondere lohnenswert dabei, auch wenn in diesem Beispiel nicht vorkommend, ist zu versuchen, alle möglichen Schleifen in Cython umzuschreiben.

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