Der Barnes-Hut-Algorithmus
Eine Einführung in die Simulation von N-KörperproblemenDas N-Körperproblem
Angenommen in einem Raum befinden sich N Körper. Jeder Körper hat ein Potential, das bedeuted er besitzt ein Kraftfeld. Es könnte sich zum Beispiel um eine elektrische Kraft (Coloumbgesetz) oder um eine Gravitationskraft handeln. Die Kräfte der einzelnen Körper überlagern sich zu einem gemeinsamen Kraftfeld. Dieses Feld wirkt wiederrum auf die Körper, die dadurch eine Beschleunigung erfahren (1. Newtonsches Grundgesetz). Wir haben also ein Problem, in dem das Kraftfeld von den Positionen der Körper abhängt und die Position bzw. die Bewegung der Körper wiederum vom Kraftfeld. In aller Regel hat man eine bekannte Anfangsverteilung der Körper mit bekannten Positionen und Geschwindigkeiten und möchte ausgehend davon berechnen, wie sich das System mit der Zeit entwickelt. Typische Anwendungsfälle für derartige Probleme sind zum Beispiel in der Astronomie die Berechnung von Galaxienentwicklung oder in der Biologie das Verhalten von Molekülen.
Die Bewegungsgleichungen sind in aller Regel recht einfach aufgestellt. Wenn man Beispielsweise eine Galaxie simulieren möchte, dann lautet die Gleichung um die Kraft auf ein Teilchen der Gesamtmenge zu berechnen:
Angenommen Fij ist die Kraft zwischen dem Teilchen i und j. Da ein Körper keine Kraft auf sich selbst ausüben kann, ist die Gesamtzahl aller möglichen Kräftepaare im System N*(N-1). Berücksichtigt man weiterhin, das Nach dem Wechselwirkungsgesetz jede Kraft eine gleich große Gegenkraft hat, also Fij = -Fji ist, so reduziert sich die Gesamtanzahl der Kraftberechnungen letztendlich auf:
Der Aufwand steigt also mit der Ordnung O(N2) an. Wenn sich die Teilchenanzahl verdoppelt, verfierfacht sich der Rechenaufwand, verzehnfacht man die Teichenanzahl, dann verhundertfacht sich der Rechenaufwand. Es ist absehbar, daß diese Art der Lösung des N-Körperproblemes für eine große Teilchenanzahl viel zu aufwendig ist. Möchte man mit sehr großen Teilchenanzahlen rechen, dann benötigt man einen besseren Algorithmus. Daran ändert sich im übrigen auch nichts, wenn man berücksichtigt, daß diese Art von Problem seit einiger Zeit sehr gut auf Spezialhardware (z.B. GPU / CUDA) durch massive parallelisierung lösen kann. Das Problem verschiebt sich dadurch lediglich um einige Zahntausend Teilchen nach "hinten" und auch der Barnes-Hut Algorithmus kann auf der GPU implementiert werden.
Der Barnes-Hut Algorithmus
Der Barnes-Hut-Algorithmus ist ein 1986 von Josh Barnes und Piet Hut [1] veröffentlichtes Näherungsverfahren, das eine effektive Berechnung der Kräfte in einem N-Körper-Problem ermöglicht. Im Gegensatz zur direkten Aufsummierung der Kräfte, deren Rechenaufwand mit O(N2) ansteigt reduziert sich der Aufwand beim Barnes-Hut-Algorithmus auf O(N log N).
Funktionsweise
Der Barnes-Hut-Algorithmus verringert die Anzahl zu berechnenden Kräfte durch geeignetes Zusammenfassen von Teilchengruppen zu Pseudoteilchen. Die Grundidee ist dabei, dass die von einer Partikelgruppe auf ein Einzelpartikel ausgeübte Kraft in sehr guter Näherung durch die Wirkung einer einzelnen Masse im Massenschwerpunkt der Partikelgruppe approximiert werden kann. So kann man beispielsweise die Kraft, die von der Andromeda-Galaxie auf die Sonne ausgeübt wird, sehr gut durch eine Punktmasse nähern, die sich im Zentrum der Andromeda-Galaxie befindet und die deren Masse hat. Diese Näherung ist allerdings nur zulässig, wenn der Abstand der Gruppe vom Einzelteilchen groß und der Gruppendurchmesser im Verhältnis zum Abstand klein ist. Das Verhältnis von Gruppenentfernung (d) zu Gruppendurchmesser (r) wird als Multipol-Akzeptanz Kriterium (engl.: MAC; Multipole-Acceptance-Criterion) bezeichnet:
Unterschreitet θ eine bestimmten Schwellwert, so sollte die Näherung nicht mehr angewendet werden um größere Fehler zu vermeiden. Barnes und Hut verwendeten in ihrer Simulation einen Wert von θ=1. Prinzipiell gilt: Je kleiner θ, um so genauer werden die Simulationsergebnisse. Der Barnes-Hut-Algorithmus wendet dieses Prinzip rekursiv an, so dass es auf verschiedenen Längenskalen wirkt (Bild 2). Wie genau das funktioniert wird im nächsten Abschnitt erläutert.
| Bild 1: Die von einer Partikelverteilung auf ein Einzelpartikel ausgeübte Kraft kann durch die von einer Punktmasse ausgehenden Kraft angenähert werden, wenn die Entfernung der Partikelgruppe zum Einzelteilchen groß im Verhältnis zum Gruppendurchmesser ist. | Bild 2: Die gravitative Wirkung von Sternhaufen und Stern B auf Stern A kann infolge der großen Entfernung als Punktmasse approximiert werden. Doch auch innerhalb des Sternhaufens kann die Gravitationskraft des Sternhaufens auf Stern B durch eine Punktmasse angenähert werden, da Stern B weit genug vom Sternhaufen entfernt ist. |
Die Implementierung
In den folgenden Kapiteln wird die Implementierung des Algorithmus erläutert. Alle Erklärungen und Beispiele beziehen sich auf den zweidimensionalen Fall. Eine verallgemeinerung auf die dritte Dimension ist unkompliziert soll hier aber nicht näher Erläutert werden.Die Baumstruktur
Ausgangspunkt ist eine gegebenen Verteilung von Teilchen im zweidimensionalen Raum. Der erste Schritt des Barnes-Hut Algorithmus ist das einsortieren der Teilchen in eine hierarchische Struktur, den Quadtree (bzw. Octtree in 3D). Dafür wird der Raum in vier Quadranten unterteilt. Diese werden im folgenden nach den Himmelsrichtungen NW, NO, SW, SO benannt. Um ein Teilchen in den Baum einzusortieren wird ermittelt, in welchem Quadranten es sich befindet. Wenn sich bereits ein anderes Teilchen in diesem Quadranten befindet, wird der Quadrant erneut unterteilt und beide Teilchen werden in ihre jeweiligen Unterquadranten einsortiert. Dieser Vorgang wird solange wiederholt, bis alle Teilchen einsortiert sind und sich jedes Teilchen in einem eigenen Quadranten befindet. In Bild 3 ist das Ergebnis dieses Vorganges für eine Beispielverteilung dargestellt.
| Bild 3: rechts: Räumliche Verteilung der Partikel mit ihren jeweiligen Quadranten; links: Zur Verteilung gehörende Baumstruktur. |
Die Erstellung des Quadtree ist ein rekursiver Vorgang, der durch folgenden Pseudocode beschrieben werden kann:
Function MainApp::BuildTree
{
ResetTree
for all particles
rootNode->InsertToNode(particle)
end for
}
Function TreeNode::InsertToNode(newParticle)
{
if number of particles in this node > 1
{
quad = GetQuadrant(newParticle);
if subnode(quad) does not exist
create subnode(quad)
subnode(quad)->InsertToNode(newParticle)
}
else if number of particles in this node == 1
{
quad = GetQuadrant(existingParticle);
if subnode(quad) does not exist
create subnode(quad)
subnode(quad)->InsertToNode(existingParticle)
quad = GetQuadrant(newParticle);
if subnode(quad) does not exist
create subnode(quad)
subnode(quad)->InsertToNode(newParticle);
}
else if node is empty
{
store newParticle in node as existingParticle
}
Increase number of particles
}
Später in der Simulation muß der Baum muß für jeden Zeitschritt neu berechnet werden. Die Komplexität der
Baumerstellung hängt von der Anordnung der Partikel ab (siehe auch [2]).
Der Aufwand um ein Teilchen in den Baum einzufügen ist direkt proportional zur Entfernung des Zielknotens
vom Baumursprung. Gibt es viele dicht beieinander liegende Körper dauert das einsortieren länger, da alle
Teilchen in einem kleinen Gebiet konzentriert und der Baum dementsprechend oft unterteilt werden muß.
Berechnen der Masseverteilung
Bis jetzt enthält der Baum nur informationen über die Räumliche Verteilung der Teilchen. Der nächste Schritt in der Simulation ist die Berechnung der Masseverteilung im Quadtree (bzw. Octtree). Es dabei speziell um die Berechnung der Gesamtmasse aller in den Knoten enthaltenen Teilchen, sowie die Berechnung ihrer Massenschwerpunkte. Diese Daten müßen für alle Knoten des Baumes die Teilchen enthalten berechnet werden die Berechnung erfolgt wieder rekursiv durch abscannen der Baumknoten. Folgender Pseudocode beschreibt den prinzipiellen Ablauf:
Function TreeNode::ComputeMassDistribution()
{
if number of particles in this node equals 1
{
CenterOfMass = Particle.Position
Mass = Particle.Mass
}
else
{
// Compute the center of mass based on the masses of all child quadrants
// and the center of mass as the center of mass of the child quadrants
// weightes with their mass
for all child quadrants that have particles in them
{
Quadrant.ComputeMassDistribution
Mass += Quadrant.Mass
CenterOfMass += Quadrant.Mass * Quadrant.CenterOfMass
}
CenterOfMass /= Mass
}
}
Der Aufwand für diese Prozedur ist O(N log N) [2].
Die Kraftberechnung
Der letzte aber entscheidene Schritt ist die Berechnung der Kraft auf ein beliebiges Teilchen. Dieser Schritt ist es, der den Barnes-Hut-Algorithmus schneller macht als die direkte Aufsummierung der Kräfte. Das wäre jetzt ein passender Moment sich nochmal die Bilder 1 und 2 ins Gedächtnis zu rufen, denn hier wird nun die Implementierung erläutert. Um die Kraft auf ein Teilchen zu berechnen wird der Quadtree erneut rekursiv durchiteriert. Für jeden Knoten wird das Verhältnis θ berechnent:
mit:
- d - Seitenlänge Größe der Box
- r - Entfernung des Teilchens vom Massenschwerpunkt des Knotens
|
| Animation 1: Gezeigt wird eine Verteilung von 5000 Körpern. Es werden nur die Knoten angezeigt, die in die Kraftberechnung auf das Teilchen in der Bildmitte eingehen (Koordinatenursprung). Umso größer θ ist, je weniger Knoten werden für die Berechnung benötigt, umso größer ist allerdings auch der Fehler des berechneten Kraftvektors. (wie so oft: Schelligkeit kontra Genauigkeit!) |
Function MainApp::CalcForce
for all particles
force = RootNode.CalculateForceFromTree(particle)
end for
end
Function force = TreeNode::CalculateForce(targetParticle)
force = 0
if number of particle equals 1
force = Gravitational force between targetParticle and particle
else
r = distance from nodes center of mass to targetParticle
d = height of the node
if (d/r < θ)
force = Gravitational force between targetParticle and node using the nodes center of mass and the total mass of the node
else
for all child nodes n
force += n.CalculateForce(particle)
end for
end if
end
end
Galaxienkollision als Beispiel
Kombiniert man nun den Barnes-Hut-Algorithmus mit einem geeigneten Integrationsverfahren (z.Bsp. Adams-Bashforth), so kann man mit relativ einfachen mitteln in erträglicher Zeit eine 2D Galaxienkollision errechnen:
Quellenverzeichnis
Der Inhalt dieses Artikels basiert auf folgenden Veröffentlichungen:- J. Barnes und P. Hut: ''A hierarchical O(N log N) force-calculation algorithm'' in Nature,324(4) Dezember 1986
- Jim Demmel '' Fast Hierarchical Methods for the N-body Problem'', Part 1, Applications of Parallel Computers (CS267): Lecture 24, April 11, 1996
- Tom Ventimiglia and Kevin Wayne '' Barnes-Hut Galaxy Simulator'', Programming Assignment (Computer Science; COS126), 2003
Autor: Ingo Berg; Dieser Text und die dazugehörigen Mediendateien stehen, sofern nicht anderweitig angegeben, unter der Creative Commons Namensnennung 3.0 Deutschland Lizenz.