Unterschiede

Hier werden die Unterschiede zwischen zwei Versionen gezeigt.

--- ss17:viele_dinge_fliegen_im_weltall_durcheinander [2017/08/31 21:22]
mtoppermann [Die Physik-Engine]
+++ ss17:viele_dinge_fliegen_im_weltall_durcheinander [2017/09/07 15:16] (aktuell)
studierenderfh98 [Die Teilchenbewegung]
@@ Zeile 25: / Zeile 25: @@
 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: $$F = G\cdot\frac{m_1\cdot m_2}{r^3}\cdot\overrightarrow{r}
+Die Grundlage der Berechnungen bildet das Gravitationsgesetz:
-$$ 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.\\
+$$\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//
@@ Zeile 114: / Zeile 117: @@
 	plt.pause(0.002)
 </code>
-Die im Ausschnitt erwähnte move-Funktion berechnet die Teilchenbewegung für einen bestimmten Zeitschritt nach dem Euler-Verfahren folgender Maßen.
+\\
+==== Die Teilchenbewegung ===
+Die im vorherigen Ausschnitt erwähnte move-Funktion berechnet die Teilchenbewegung für einen bestimmten Zeitschritt nach dem Euler-Verfahren folgender Maßen.
 <code python>
 # Aus einer übergebenen Liste von Teilchen wird die neue Position aus
@@ Zeile 122: / Zeile 127: @@
 		i.x += zeitschritt*i.velocity[0]
 		i.y += zeitschritt*i.velocity[1]
+</code>
+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.
+<code python>
+	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
 </code>
@@ Zeile 133: / Zeile 154: @@
 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, besitzt unser Programm eine potentielle Laufzeit von $O(n^2)$.\\
+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.\\
-Bei $n=100$ Teilchen müssen als $100^2 = 10000$ Einträge bestimmt werden, bei $n=500$ sind es bereits $500^2=250000$ Einträge.\\
+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)$\\
-Die 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.\\
+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.\\
-====Analyse des Programms====
+====Programmanalyse====
-Als erstes ist es wichtig, herauszufinden, welche Funktionen besonders laufzeitintensiv sind und eine Implementierung in Cython sich somit lohnen würde. Die Programmanalyse macht man zum Beispiel mit dem Modul cProfile:
+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:
 <code shell>
 python -m cProfile -o analyse.stat run_project.py
@@ Zeile 179: / Zeile 201: @@
 </code>
-Wie zu erwarten, sind vor allem die Funktionen, die die Grundlage unserer Simulation bilden, für über 70% der Laufzeit verantwortlich (namentlich get_way, get_force und get_force_matrix). Diese Funktionen in Cython zu schreiben ist daher naheliegend.
+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.
+====Cython schreiben====
+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:
+<code python>
+# Fibonacci in Python
+def fib(n):
+    a, b = 0, 1
+    for i in range(n):
+        a, b = a + b, a
+    return a
+</code>
+<code c>
+// 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;
+}
+</code>
+<code python>
+# 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
+</code>
+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.
+<code python>
+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)
+</code>
+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.
+<code>
+erfolg python 0.000105
+erfolg cython 0.000005
+</code>
+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.
 =====Links=====
 [[ss17:viele_dinge_fliegen_im_weltall_durcheinander:protokolle|Hier]] geht es zu den Protokollen, die Einblick in den Projektverlauf/Stand und die einzelnen Labor- und Blocktermine geben.\\

Mathesis Wiki

Benutzer-Werkzeuge

Webseiten-Werkzeuge

Unterschiede

Seiten-Werkzeuge