Der Barnes-Hut-Algorithmus

Das N-Körperproblem

Angenommen in einem Raum befinden sich N Körper. Jeder Körper hat ein Potenzial, 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).

Das ist ein Problem, bei 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:

wobei:

• \(m_i, m_j\) - Massen der Körper i und j
• \(\vec{r_i}, \vec{r_j}\) - Ortsvektoren der Körper i und j
• \(G\) - Gravitationskonstante

Angenommen F_ij ist die Kraft zwischen den 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 F_ij = -F_ji ist, so reduziert sich die Gesamtanzahl der Kraftberechnungen letztendlich auf:

Der Rechenaufwand steigt mit der Ordnung O(N²) an. Wenn sich die Teilchenanzahl verdoppelt, vervierfacht sich der Rechenaufwand. Verzehnfacht man die Teichenanzahl, dann verhundertfacht sich der Rechenaufwand. Es ist absehbar, dass 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 auch nichts, wenn man berücksichtigt, dass diese Art von Problem seit einiger Zeit sehr gut auf Spezialhardware (z.B. GPU / CUDA) durch massive Parallelisierung gelöst werden kann. Das Problem verschiebt sich dadurch lediglich um einige Zehntausend 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 Näherungsverfahren, dass eine effektive Berechnung der Kräfte in einem N-Körper-Problem ermöglicht. Er wurde erstmalig 1986 von Josh Barnes und Piet Hut in der Zeitschrift Nature beschrieben [1]. Im Gegensatz zum Verfahren der direkten Aufsummierung der Kräfte, dessen Rechenaufwand mit O(N²) ansteigt, reduziert sich der Aufwand beim Barnes-Hut-Algorithmus auf O(N log N).

Der Algorithmus verringert die Anzahl zu berechnenden Kräfte durch geeignetes Zusammenfassen von Teilchengruppen zu Pseudoteilchen. Die Grundidee ist, 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 angenähert 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 Multipol-Akzeptanz Kriterium genannt (engl.: MAC; Multipole-Acceptance-Criterion):

Unterschreitet θ einen bestimmten Schwellwert, so sollte die Näherung nicht angewendet werden um größere Fehler zu vermeiden. Barnes und Hut verwendeten in ihrer Simulation einen Schwellwert von θ=1. Je kleiner θ gewählt wird, umso genauer werden die Simulationsergebnisse. Der Barnes-Hut-Algorithmus wendet dieses Prinzip rekursiv an, so dass es auf verschiedenen Längenskalen wirkt (Bild 2).

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.

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 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 sogenannten 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.

Räumliche Verteilung der Partikel mit ihren jeweiligen Quadranten und der zur Verteilung gehörender Quadtree.

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 muss der Baum 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, so dauert das Einsortieren länger, da alle Teilchen in einem kleinen Gebiet konzentriert und der Baum dementsprechend oft unterteilt werden muss.

Berechnen der Masseverteilung

Bis jetzt enthält der Baum nur Informationen über die räumliche Verteilung der Körper. Der nächste Schritt in der Simulation ist die Berechnung der Masseverteilung im Quadtree (bzw. Octtree für 3-D Verteilungen). Es geht dabei speziell um die Gesamtmasse aller in den Knoten enthaltenen Körper, sowie die Berechnung ihrer Massenschwerpunkte. Diese Daten müssen für alle Knoten des Baumes, die nicht leer sind, ermittelt werden. Die Berechnung erfolgt 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. Um die Kraft auf ein Teilchen zu berechnen wird der Quadtree erneut rekursiv durchlaufen. Für jeden Knoten wird das Verhältnis θ nach folgender Formel berechnet:

• d - Seitenlänge der Box
• r - Entfernung des Teilchens vom Massenschwerpunkt des Knotens

Anschließend wird getestet, ob θ unter oder über einem bestimmten Schwellwert liegt. (Meistens verwendet man eine Schwelle von θ < 1). Liegt θ darunter, so kann die Kraft, die von dem Quadranten ausgeht genähert werden, indem der Masseschwerpunkt des Knotens und die Knotenmasse in das Gravitationsgesetz eingesetzt werden. Ist θ größer, so muss der Quadrant erneut unterteilt und alle Unterquadranten nach der selben Methode geprüft werden.

Iteriert wird solange es noch Unterquadranten gibt. Im ungünstigsten Fall also bis zu den Endknoten, für die Schwerpunkt und Masse identisch mit dem einzigen darin enthaltenen Teilchen sind. Das bedeuted, dass man für einen Schwellwert von θ=0 die gleiche Lösung bekommt wie bei der direkten Aufsummierung der Kräfte (nur viel langsamer, da der Baum ja trotzdem erzeugt werden muss). Die Wahl von θ beeinflußt die notwendige Iterationstiefe und damit die Gesamtanzahl der notwendigen Kraftberechnungen. Die folgende Animation verdeutlicht diesen Zusammenhang exemplarisch.

Barnes Hut Tree für eine Verteilung von 5000 Körpern.

In der Animation werden nur Knoten gezeigt, die in die Kraftberechnung auf das Teilchen in der Bildmitte eingehen (Koordinatenursprung). Je grösser θ ist, desto weniger Knoten werden für die Berechnung benötigt, umso grösser ist allerdings auch der Fehler des berechneten Kraftvektors. Wie so oft gilt auch hier, dass man Schelligkeit mit Genauigkeit "erkaufen" muss! Der Pseudocode zur Kraftberechnung sieht wie folgt aus:

 
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 
    else
      for all child nodes n
        force += n.CalculateForce(particle)
      end for
    end if
  end
end

Galaxienkollision als Beispiel

Kombiniert man den Barnes-Hut-Algorithmus mit einem geeigneten Integrationsverfahren (z.Bsp. Adams-Bashforth), so kann man mit relativ einfachen Mitteln und in erträglicher Zeit eine 2D Galaxienkollision errechnen:

Download

Das Beispielprojekt ist in C++ geschrieben. Zum kompilieren des Projektes werden die OpenGL, SDL und SDL_TTF Bibliotheken und Includefiles benötigt.

Barnes-Hut-Simulator herunterladen

Quellenangaben

1. J. Barnes, P. Hut: "A hierarchical O(N log N) force-calculation algorithm" Nature; 324(4); Seiten: 446-449; Dezember 1986
2. Jim Demmel: "Fast Hierarchical Methods for the N-body Problem", Part 1, Applications of Parallel Computers (CS267): Lecture 24, April 11, 1996
3. Tom Ventimiglia, Kevin Wayne: "Barnes-Hut Galaxy Simulator" Programming Assignment (Computer Science; COS126), 2003