Das Mehrkörperproblem in der Astronomie/ Hierarchische Algorithmen/ Räumliche Hierarchie
Idee
Um die in einem Mehrkörpersystem herrschenden Kräfte exakt zu bestimmen, muss man wie schon mehrfach erwähnt sämtliche Massenpaare individuell betrachten, wodurch der Rechenaufwand quadratisch mit der Anzahl der Mitglieder wächst. Dieser rapide Anstieg wird jedoch vermieden, wenn man stattdessen nur die in der näheren Umgebung eines Körpers sich befindlichen Nachbarn exakt behandelt, weiter entfernte Massenpunkte hingegen zu Gruppen zusammenfasst und anstelle der Einzelobjekte nur noch die Schwerpunkte dieser Cluster zur Kräfteberechnung heranzieht. Je weiter eine solche Gruppe entfernt ist, umso größer darf deren räumliche Ausdehnung sein, ohne die von ihr ausgehende Gravitation zu ungenau abzuschätzen. Auf Grundlage dieser Idee lässt sich ein Algorithmus entwickeln, dessen Rechenaufwand in der gewünschten Größenordnung liegt.
Um Gruppen von Massenpunkten zu definieren, die je nach Entfernung unterschiedlich groß gehalten werden können, wird gemäß dem Verfahren von Barnes und Hut (1986) [1] das Ensemble in Würfel (Oktanten) verschiedener Kantenlänge eingeteilt, wobei nun unterschiedliche Raumebenen darstellt. Im Gegensatz zur Konstruktion der Zeithierarchie startet man jetzt nicht auf der niedrigsten, sondern der höchsten Ebene (welche nun aber umgekehrt durch gekennzeichnet ist und nicht durch den höchsten vorkommenden Wert von ). Um das System auf der obersten Ebene mit 8 Würfeln der Kantenlänge sicher zu überdecken, muss sich nach dem Abstand des fernsten Körpers vom Gesamtschwerpunkt des Ensembles richten.
Für jeden dieser Würfel wird die Anzahl der darin sich befindlichen Objekte ermittelt. Ist mehr als ein solches in einem Oktanten vorhanden, wird dieser in 8 kleinere Würfel der Größe unterteilt, es tritt die Raumebene hinzu. Im nächsten Schritt werden diese kleineren Kuben überprüft. Findet sich in einem solchen trotz der kleineren Kantenlänge erneut mehr als ein Massenpunkt, erfolgt für diesen eine neuerliche Aufteilung in noch kleinere Würfel mit einer Ausdehnung von nunmehr entsprechend der Ebene . Dieses Verfahren wird solange fortgesetzt, bis in jedem Oktanten nur noch ein Körper anzutreffen ist.

Rücken im Verlauf einer Simulation Massenpunkte eng zusammen, so werden dort u.U. sehr kleine Würfel benötigt, um auf ein Objekt pro solchem zu kommen, d.h. eventuell müssen gleich mehrere Raumebenen angelegt werden. Entfernen sie sich wieder voneinander, sind solch kleine Oktanten und entsprechend tiefe Ebenen nicht mehr erforderlich.
Zwischen Raum- und Zeithierarchie besteht ein fundamentaler Unterschied. Bei letzterer genügt es, die Zeitebenen und die dazugehörigen dynamischen Zeitschritte der einzelnen Körper zu kennen. Hingegen ist für die räumliche Hierarchie die Kenntnis allein der Ebenen und der Dazugehörigkeiten der Massenpunkte zu den entsprechenden Würfeln unzureichend. Wie nachfolgend diskutiert wird, müssen für jeden belegten Kubus auch die Positionen und sofern die Simulation auf dem Hermite-Polynome-Verfahren beruht, auch die Geschwindigkeiten ihrer Schwerpunkte bekannt sein, welche sich mit jedem Simulationsschritt ändern.
Nach jeder Ausführung des Prädiktors muss daher die Raumhierarchie aktualisiert werden. Zunächst überprüft man, ob durch den Prädiktor zumindest ein Massenpunkt bezogen auf die niedrigste Raumebene in einen anderen Würfel überwechselt. Ist das nicht der Fall, können die momentan definierten Oktanten beibehalten werden, nur die Positionen und gegebenenfalls die Geschwindigkeiten ihrer Schwerpunkte müssen neu berechnet werden. Ansonsten muss das hierarchische Gitter selbst neu aufgebaut werden.
Prozeduren im Würfelgitter
Einordnung der Massenpunkte
Die Einstufung der Massenpunkte in die Raumhierarchie ist viel anspruchsvoller als deren Zuordnung zur Zeithierarchie. Bei letzterer genügt es, für jeden Körper nur die Zeitebene vorzuhalten, auf welcher er sich tatsächlich befindet. Im Falle der Raumhierarchie jedoch muss für jedes Objekt nicht nur die Lage bezüglich der tatsächlichen, sondern auch aller darüber liegenden Ebenen bekannt sein. Befindet sich nämlich ein Körper nahe an einem Massenpunkt, dessen Beschleunigung bestimmt werden soll, muss jener als Einzelobjekt oder bestenfalls als zu einer Gruppe geringer Ausdehnung (die schon mit einem kleinen Würfel überdeckt werden kann) zugehörig betrachtet werden. Ist der gleiche Körper weit von einem anderen Massenpunkt entfernt, darf er von diesem aus gesehen als Mitglied eines weiträumigen Clusters (der nur von einem großen Oktanten abgedeckt wird) aufgefasst werden. Betrachtet man in obiger Abbildung z.B. den Massenpunkt oben rechts, so ist dessen nächster Nachbar als Einzelobjekt zu behandeln. Aus der Sicht des Massenpunktes unten links dürfen die beiden Körper oben rechts hingegen zu einer Gruppe zusammengefasst werden. Je nach Perspektive ist ein Objekt somit unterschiedlichen Raumebenen zugehörig.
Die Lage eines Massenpunktes in der Raumhierarchie kann durch ganzzahlige kartesische Koordinaten angegeben werden. Auf der höchsten Ebene sind nur 2*2*2 Würfel vorhanden, die möglichen Werte pro Koordinate also 0 oder 1. Auf der nächsten Ebene existieren 4*4*4 Kuben, so dass sich der mögliche Wertebereich von 0 bis 3 erstreckt. Allgemein weist eine Ebene Oktanten auf, so dass jede Koordinate Werte von 0 bis annehmen kann. In obigem zweidimensionalen Beispiel hat der Massenpunkt rechts oben auf der Ebene 0 die Position (1,1), auf der Ebene 1 die Lage (3,2) und auf seiner niedrigsten Ebene 2 die Koordinaten (7,5).
Die Einordnung der Massenpunkte in das Würfelgitter beruht auf folgenden Algorithmus, für welchen zunächst nur die Struktur erläutert wird. Der dazugehörige C-Code wird im letzten Abschnitt dieses Unterkapitels angegeben.
Start auf höchster Ebene e = 0
while (zumindest zwei Massenpunkte im gleichen Würfel)
{
Berechnung der Positionen aller Körper i im Würfelgitter bezüglich der aktuell betrachteten Ebene e
for (i = 0;i < N;i ++)
{
if (Ebene des Körpers i = aktuell betrachtete Ebene e)
{
Neuer Würfel für Körper i
for (j = i + 1;j < N;j ++)
{
if (Körper j im gleichen Würfel wie Körper i)
{
Aufnahme des Körpers j in den Würfel von Körper i
Verschieben beider Körper i und j auf die nächsttiefere Ebene
}
}
}
}
if (zumindest ein Paar i,j mit gemeinsamem Würfel auf aktuell betrachteter Ebene gefunden)
Analyse auf nächsttieferer Ebene
e++
}
Zu Beginn auf der obersten Ebene, wird für den ersten Massenpunkt sogleich ein Würfel der Kantenlänge angelegt. Die Positionen aller anderen Objekte auf der aktuell untersuchten Ebene werden mit derjenigen des Körpers verglichen. Bei gleicher Position (d.h. Zugehörigkeit zum gleichen Oktanten wie Massenpunkt ) wird zum Kubus von hinzugefügt, gleichzeitig beide Körper und der nächsttieferen Ebene zugeordnet. Solche Objekte werden von der laufenden -Schleife nicht mehr bearbeitet, da diese nur Massenpunkte berücksichtigt, welche der aktuell analysierten Ebene angehören. Dies bedeutet aber, dass nicht alle möglichen Paare geprüft werden müssen, um die Mitglieder des Ensembles in die Raumhierarchie einzuordnen.
Bleibt ein Körper in seinem Würfel allein, wird er nicht auf die nächsttiefere Ebene geschoben. Geht der Algorithmus zur nächsten Ebene über, so wird ein solches Einzelobjekt nicht mehr von der -Schleife beachtet, da es im Vergleich zur aktuellen Ebene um eine Stufe zurückgeblieben ist. Hingegen erfasst der Algorithmus nun solche Massenpunkte, welche auf der Ebene zuvor von der -Schleife um eine Stufe nach unten versetzt wurden.
Auf jeder Ebene werden für alle belegten Würfel laufende Nummer, Masse, die Position und eventuell auch die Geschwindigkeit des Schwerpunkts festgehalten, zudem auch die laufende Nummer des jeweils übergeordneten Würfels auf der nächsthöheren Ebene. Um die Schwerpunkte der Oktanten schnell neu berechnen zu können, wird für jeden auch eine Liste mit den Nummern der enthaltenen Massenpunkte angelegt.
Bestimmung von Beschleunigung und Ruck
Um die auf einen Körper ausgeübte Beschleunigung und für das Hermite-Polynome-Verfahren auch den Ruck zu ermitteln, werden jetzt nur die Massenpunkte in dessen unmittelbarer Umgebung als Einzelobjekte betrachtet. Von einer gewissen Entfernung an werden alle in einem Würfel sich befindlichen Objekte kollektiv behandelt, indem die Gesamtmasse des Oktanten sowie die Position und gegebenenfalls auch die Geschwindigkeit dessen Schwerpunkts zur Berechnung herangezogen werden anstatt die entsprechenden Größen der individuellen Mitglieder. Diese Näherung ist zulässig, sofern die Kantenlänge eines solchen Würfels nicht zu groß ist im Vergleich zu seiner Entfernung von dem soeben untersuchten Körper, d.h. das Verhältnis zwischen den beiden Längenskalen einen gewissen Wert nicht überschreitet:
Je kleiner dieser Schwellwert gesetzt wird, umso genauer werden die Kräfte zwischen den Mitgliedern eines Mehrkörpersystems bestimmt, umso größer ist jedoch auch der Rechenaufwand. Der hierarchische Algorithmus greift dann erst auf recht tiefer Ebene, da ein im Vergleich zur Kantenlänge eines Kubus erheblicher Abstand desselben verlangt wird. Als Folge davon werden weiterhin viele in der Umgebung eines Körpers sich befindliche Massenpunkte individuell betrachtet, und selbst zu größeren Entfernungen hin zumeist nur kleine Würfel, die lediglich wenige Objekte zu einer Gruppe zusammenfassen. Mit = 0 wird auf eine Gruppenbildung komplett verzichtet.
Umgekehrt bedeutet ein großer Schwellwert, dass nur noch wenige Massenpunkte in der unmittelbaren Umgebung eines Körpers einzeln behandelt werden. Mit zunehmendem Abstand geht das Verfahren zügig zu Würfeln großer Kantenlänge über, die viele Objekte enthalten können. Die nun sehr effiziente Gruppenbildung und damit kürzere Rechenzeit ist jedoch mit einem größeren Fehler bei der Berechnung der Kräfte im Ensemble erkauft.
Gemäß Barnes und Hut (1986) [1] sollte maximal auf 1 gesetzt werden. Den Autoren zufolge werden dann die auf die einzelnen Mitglieder eines Mehrkörpersystems einwirkenden Kräfte mit einer relativen Genauigkeit von 1% bestimmt.
Das Vorgehen bei der Berechnung der auf einen Körper einwirkenden Beschleunigung und des Rucks wird wie der Algorithmus zur Einordnung der Massenpunkte in die räumliche Hierarchie zunächst rein qualitativ vorgestellt.
for (e = 0;e < Anzahl der Ebenen;e ++)
{
for (o = 0;o < Anzahl der belegten Oktanten auf Ebene e;o ++)
{
if (übergeordneter Oktant bereits abgearbeitet)
dann auch aktueller Oktant abgearbeitet
else
{
Betrachte Abstand r zwischen aktuellem Oktant und Körper i
if (Kantenlänge (a / 2^e) / r < alpha oder Oktant enthält nur 1 Massenpunkt)
{
Bestimme Beschleunigung und Ruck, welche der aktuelle Oktant auf Körper i ausübt
aktueller Oktant abgearbeitet
}
}
}
}
Aus der Sicht eines Körpers werden beginnend mit der höchsten Ebene die gesamte Hierarchie und innerhalb jeder Ebene alle belegten Würfel abgearbeitet. Zuerst wird für den gerade betrachteten Oktanten geprüft, ob der darüber liegende bereits in die Berechnung eingegangen ist. In diesem Fall kann der aktuelle Kubus einfach übergangen und sogleich als abgearbeitet gekennzeichnet werden. Ansonsten wird der Abstand zwischen dessen Schwerpunkt und dem Körper ermittelt. Ist im Vergleich zur Kantenlänge des Würfels groß genug, so werden alle darin vorhandenen Massenpunkte wie soeben beschrieben als Kollektiv behandelt und dieser wiederum als abgearbeitet markiert. Kann der Oktant wegen einer im Vergleich zu zu großen Kantenlänge nicht abgearbeitet werden, erfasst der Algorithmus auf der nächsttieferen Ebene dessen Teilwürfel. Enthält ein Kubus nur 1 Massenpunkt, kann er natürlich unabhängig von seiner Größe und Entfernung berücksichtigt werden. Auf diese Weise werden auch diejenigen Objekte einbezogen, welche sich nicht zu Gruppen zusammenfassen lassen, weil sie sich entweder in der unmittelbaren Umgebung des Körpers oder in Bereichen geringer Dichte aufhalten.
Die Aussagen von Barnes und Hut (1986) [1] über die Genauigkeit ihres Verfahrens gelten nur im Mittel. Im Einzelfall hängt der durch die Gruppenbildung bewirkte Fehler stark von der Verteilung der Massenpunkte in einem Würfel ab. Als Beispiele seien die in folgender Skizze dargestellten Fälle diskutiert, in denen jeweils drei Massenpunkte der Masse durch ihren Schwerpunkt ersetzt werden.

Im ersten Fall stehen diese hintereinander und üben die Beschleunigungen , und auf. Mit lauten diese , und . Ersetzt man die Körper durch ihren Schwerpunkt, so ergibt sich die genäherte Beschleunigung . Vergleicht man den exakten Wert mit der Näherung , so gilt folgendes Verhältnis:
Im zweiten Fall sind die Massenpunkte nebeneinander postiert. Die Beschleunigungen sind nun , und wiederum bzw. , und nochmals . Der Vergleich mit der Näherung liefert nun:
Stehen die Massenpunkte hintereinander, so wird durch eine Gruppenbildung die von diesen ausgeübte Beschleunigung mit einer im Vergleich zur Entfernung zunehmenden Größe des Würfels weit unterschätzt. Mit = 0.5 liegt der relative Fehler bei knapp 20%, mit = 1 aber schon bei etwa 80%. Die Ursache für diese erhebliche Abweichung liegt in der -Abhängigkeit der Gravitation. Liegt der Schwerpunkt eines Oktanten weiter entfernt als ein individueller Massenpunkt, so wird mit der Gruppenbildung die vom einzelnen Körper ausgehende Kraft unterschätzt. Liegt der Schwerpunkt näher als ein Einzelobjekt, so verhält es sich umgekehrt. Die -Abhängigkeit sorgt jedoch dafür, dass der Fehler im ersten Fall deutlich schwerer wiegt als im letzteren.
Liegen die Massenpunkte nebeneinander, so verhält sich das Verfahren wesentlich stabiler. Mit = 0.5 und 1 finden sich nun relative Fehler von einigen % und etwas mehr als 10%. Mit einer Gruppenbildung wird die ausgeübte Beschleunigung leicht überschätzt, da die Entfernung zum Schwerpunkt etwas geringer ist als zu den Mittelpunkten der seitlichen Flächen des Würfels.

Gemäß den hier erörterten Beispielen scheint der Fehler des Barnes-Hut-Algorithmus deutlich höher zu liegen als von den Autoren ursprünglich selbst genannt. In einer weiteren Arbeit aus dem Jahr 1989 [2] geben diese in der Tat an, dass bei einer geringen Anzahl von Massenpunkten pro Würfel (weniger als 10) mit einem relativen Fehler der Kraft von bis zu 10% gerechnet werden muss. Für die Genauigkeit des Verfahrens ist nicht allein entscheidend, sondern auch die Zahl der Objekte pro Würfel. Je größer diese ist, umso geringer ist der durch Inhomogenitäten bewirkte Fehler bei der Kräfteberechnung. Vor allem im dichten Zentralbereich eines Mehrkörpersystems sind somit zuverlässige Ergebnisse zu erwarten. Man kann zudem davon ausgehen, dass die dort sich befindlichen Massenpunkte von allen Seiten von Würfeln mit relativ vielen Einzelkörpern umgeben sind und sich die pro Oktant eingeführten Fehler so in erheblichem Maße gegenseitig aufheben. Finden sich mehr als 1000 Massenpunkte pro Würfel, ist selbst mit einem von 1 ein relativer Fehler der Kraft von nur einigen 0.1% gegeben.
Wird die Anziehungskraft gemäß der Methode von Aarseth geglättet, wirkt sich dies zusätzlich stabilisierend auf den Barnes-Hut-Algorithmus aus. Dessen relativer Fehler hängt jetzt auch davon ab, wie weit ein Oktant im Vergleich zum Plummerradius entfernt ist, also vom Verhältnis . Betrachtet man obiges Beispiel abermals für den ungünstigsten Fall hintereinander stehender Massenpunkte, so gilt:
Beträgt der Abstand des Oktanten das doppelte des Plummerradius, so tritt bei einem von 0.5 ein relativer Fehler von nur noch circa 3% Prozent in Erscheinung. Wird ein gleich 1 zugelassen, steigt bei unveränderter Glättung der Gravitation die Abweichung auf immerhin noch etwa 15%. Ist der Distanz des Würfels gleich, so liegen die entsprechenden relativen Fehler nur noch bei circa 2% und 10%. Zudem kehrt der durch die Gruppenbildung begangene Fehler bei der Kräfteberechnung jetzt seine Richtung um: Diese liefert jetzt größere Beschleunigungen als die exakte Betrachtung, da für Abstände kleiner als der Plummerradius Aarseth's Glättungsformel eine um so geringere Kraft liefert, je geringer die Entfernung zwischen zwei Massenpunkten ist. Gerade dann, wenn eine stark inhomogene Verteilung von Massenpunkten gleichzeitig von sehr kleinen Abständen zwischen einzelnen Körpern (d.h. in der Größenordnung des Plummerradius) begleitet wird, erweist sich für ein geglättetes Gravitationsgesetz das Barnes-Hut-Verfahren also als recht stabil.

Insgesamt wird die Genauigkeit einer Mehrkörpersimulation mittels hierarchischer Kräfteberechnung zumeist von deren Fehlern dominiert und nicht von der Wahl der Lösungsformel für die sukzessive Berechnung von zurückgelegten Strecken und Geschwindigkeiten. Auch dies rechtfertigt, das nominell ungenauere Leapfrog-Verfahren den Hermite-Polynomen vorzuziehen.
Abschließend sei als praktisches Beispiel ein weiteres Mal das schon wiederholt besprochene Drei-Körper-System herangezogen. Der Parameter wird auf 0.1 gesetzt, d.h. aus der Sicht eines Massenpunktes die beiden übrigen Objekte erst dann zu einer Gruppe zusammengefasst, wenn der überdeckende Würfel mindestens 10 Mal kleiner ist als dessen Entfernung. Obwohl damit bei der Berechnung der von einem solchen Paar ausgehenden Gesamtkraft lediglich ein geringer Fehler in der Größenordnung von zumeist nur 0.1% begangen wird, sind die Auswirkungen der Näherung bereits drastisch. Die originalen und die genäherten Bahnen weichen im Laufe der Zeit immer mehr voneinander ab, bis schließlich jede Ähnlichkeit verschwunden ist. Für die Bearbeitung von Problemen, wo es auf eine möglichst exakte Vorhersage der einzelnen Orbits ankommt (z.B. der Frage nach der Langzeitstabilität der Planetenbahnen), ist der Barnes-Hut-Algorithmus also völlig ungeeignet. Erfolgreich kann er dagegen bei aus sehr vielen Einzelobjekten bestehenden Systemen eingesetzt werden, wo die genaue individuelle Bahn eines Körpers in der Regel nicht von Bedeutung ist. Bei Kugelsternhaufen oder Galaxien beispielsweise ist es meist ausreichend, wenn großräumige Eigenschaften wie z.B. die Verteilung der Dichte oder die Häufigkeit von Geschwindigkeiten richtig modelliert werden. Welche Auswirkungen hierbei unterschiedliche Einstellungen von haben, wird im Praxiskapitel gezeigt.

Aktualisierung der Raumhierarchie
Nach jeder Anwendung des Prädiktors muss die Raumhierarchie aktualisiert werden. Da die Massenpunkte durch eine einmalige Anwendung der Prognoseformel nur kurze Strecken zurücklegen, bleiben diese selbst auf der niedrigsten Raumebene oft allesamt in ihren aktuellen Würfeln. Für viele Simulationsschritte genügt es somit, die Schwerpunkte der momentan definierten Oktanten neu zu bestimmen. Dazu dient das folgende Verfahren, dessen Code wiederum im letzten Abschnitt dieses Unterkapitels gegeben wird.
for (e = Anzahl der Ebenen - 1;e > 0;e --)
{
for (o = 0;o < Anzahl der belegten Oktanten auf Ebene e;o ++)
{
Summiere für Schwerpunktberechnung über alle im aktuellen Oktanten vorhandenen Massenpunkte i
sofern diese nicht schon abgearbeitet sind
Markiere jeden vorhandenen Massenpunkt i als abgearbeitet
}
if (e > 0)
{
for (o = 0;o < Anzahl der belegten Oktenten auf Ebene e;o ++)
{
Betrachte die soeben aktualisierten Oktanten auf Ebene e selbst als Massenpunkte
Addiere sie zu ihren übergeordneten Oktanten auf Ebene e - 1 hinzu
}
}
}
Die Neuberechnung der Würfelschwerpunkte startet im Gegensatz zu den bisher vorgestellten Algorithmen auf der untersten Raumebene. Zuerst werden die kleinsten Kuben aktualisiert, indem über die dort vorhandenen Massenpunkte (auf der niedrigsten Ebene zwangsläufig immer nur einer pro Würfel) aufsummiert wird, wobei jeder Körper als abgearbeitet markiert wird. Die Oktanten können nun selbst als Massenpunkte angesehen und zu den ihnen übergeordneten Würfeln addiert werden. Dieses Prinzip wird auf allen Ebenen wiederholt. Massenpunkte , welche tiefer eingeordnet sind als die momentan abgearbeitete Ebene, müssen somit nicht nochmals individuell berücksichtigt werden. Auf jeder Ebene wird der Inhalt eines Oktanten jeweils an seinen übergeordneten Würfel weitergereicht.