Das Mehrkörperproblem in der Astronomie/ Hierarchische Algorithmen/ Zeitliche Hierarchie

Konstruktion

Wie soeben angedeutet, lassen sich die dynamischen Zeitskalen der Massenpunkte eines Mehrkörpersystems in eine hierarchische Struktur einordnen, indem man einen kleinsten elementaren Zeitschritt $\Delta t_{E}$ definiert und diesen sukzessive verdoppelt.

Prinzip der Simulation eines N-Körper-Systems mit hierarchischen Zeitschritten pro Massenpunkt. Jedem Körper ist je nach seiner individuellen Schrittweite eine hierarchische Schrittweite

\Delta t=2^{n}\Delta t_{E}

auf einer Zeitebene

n

zugeordnet, wobei

\Delta t_{E}

die kleinste im Ensemble vorkommende Schrittweite bedeutet

Um $\Delta t_{E}$ zu Beginn der Simulation festzulegen, muss man sich an der kleinsten anfänglich vorliegenden dynamischen Zeitskala $\Delta t_{min}$ orientieren. $\Delta t_{E}$ repräsentiert die Ebene 0, die niedrigste der Hierarchie. Durch Verdoppelung $2\Delta t_{E}$ legt man die nächsthöhere Ebene 1 an, durch weitere Verdoppelung $4\Delta t_{E}$ die Ebene 2 usw.

Fällt im Laufe der Modellierung die dynamische Zeitskala eines Körpers unter den Wert, welcher seine Zeitebene definiert, so muss er um eine Ebene nach unten geschoben werden. Tritt z.B. für ein Objekt auf Ebene 1 ein $\Delta t<2\Delta t_{E}$ auf, so wechselt dieses zur Ebene 0. Wird sogar ein $\Delta t<\Delta t_{E}$ gefunden, so ist die bisherige elementare Schrittweite nicht länger gültig. An ihrer Stelle muss ein neuer Elementarschritt $\Delta t_{E}/2$ definiert werden, welcher fortan die Ebene 0 vertritt. Dementsprechend werden alle anderen Ebenen um 1 Niveau höher eingestuft. Aus der bisherigen Ebene 0 wird die Ebene 1, die bisherige Ebene 1 ist nun die Ebene 2 usw. Um zu vermeiden, dass diese Situation schon kurz nach dem Start der Simulation eintritt, empfiehlt es sich, zu Anfang $\Delta t_{E}$ so festzulegen, dass $\Delta t_{min}$ mit einem Betrag von $3/2\Delta t_{E}$ genau zwischen $\Delta t_{E}$ und $2\Delta t_{E}$ liegt.

Nimmt umgekehrt die dynamische Zeitskala so stark zu, dass sie den Wert der nächsthöheren Ebene überschreitet, gelangt der betroffene Massenpunkt um eine Stufe nach oben. Hat etwa ein Körper auf Ebene 1 ein $\Delta t>4\Delta t_{E}$ , so steigt er zur Ebene 2 auf. Ein solches Höhersteigen ist nur zu solchen Zeitpunkten möglich, welche einem Vielfachen des auf der Zielebene gültigen Zeitschritts entsprechen, ein Wechsel z.B. von Ebene 1 zu Ebene 2 nur zu Zeitpunkten $0$ , $4\Delta t_{E}$ , $8\Delta t_{E}$ usw., nicht aber zu Zeitpunkten $2\Delta t_{E}$ , $6\Delta t_{E}$ usw. Ein aufsteigender Massenpunkt soll ja gleichzeitig mit den auf der Zielebene schon vorhandenen Objekten bewegt werden.

Ist auf der Ebene 0 kein Körper mehr vorhanden, so wird an Stelle von $\Delta t_{E}$ die auf der Ebene 1 gültige Schrittweite $2\Delta t_{E}$ als neue elementare Schrittweite angesetzt. Sämtliche Ebenen werden dadurch um eine Stufe niedriger eingeordnet. Aus der bisherigen Ebene 1 wird eine Ebene 0, aus der Ebene 2 eine Ebene 1 usw.

Die höchste erforderliche Ebene richtet sich nach der größten im Ensemble auftretenden dynamischen Zeitskala $\Delta t_{max}$ . $\Delta t_{E}$ muss solange verdoppelt werden, bis die so erzeugte Zeitskala größer als die Hälfte von $\Delta t_{max}$ ist.

Wie die niedrigste ist auch die höchste Ebene dynamisch. Überschreitet für einen Massenpunkt $\Delta t$ das Doppelte der aktuell größten Zeitskala, so muss für diesen eine neue oberste Ebene angelegt werden. Wird die höchste Ebene nicht mehr benötigt, weil die dynamischen Zeitskalen aller Objekte unter der entsprechenden Schrittweite liegen, kann diese einfach entfernt werden.

Lösungsschema

Die Steigerung der Effizienz des Prädiktor-Korrektor-Algorithmus beruht darauf, nicht mehr die individuellen Zeitschritte $\Delta t_{i}$ der einzelnen Massenpunkte zu verwenden, sondern diejenigen $\Delta t_{E}$ , $2\Delta t_{E}$ , $4\Delta t_{E}$ usw. der Zeitebenen, auf welchen sich diese befinden. Da diese Schrittweiten nur einige wenige unterschiedliche Werte annehmen können, welche untereinander durch Verdoppelung hervorgehen, ist quasi automatisch sichergestellt, dass zumeist mehrere Körper den gleichen niedrigsten Wert $t+\Delta t$ aufweisen. Die sowieso immer erforderliche Extrapolation aller Positionen (und eventuell auch der Geschwindigkeiten) hin zu diesem Zeitpunkt kann so fast immer genutzt werden, um gleichzeitig mehrere Objekte zu bewegen. Das Lösungsschema sieht nun folgendermaßen aus:

Prinzip der Simulation eines N-Körper-Systems mit hierarchischen Zeitschritten pro Massenpunkt. Das Schema ähnelt demjenigen für individuelle Zeitschritte pro Körper, doch können nun, da die Schrittweite nur noch wenige diskrete Werte

\Delta t=2^{n}\Delta t_{E}

annehmen kann, fast immer mehrere Massenpunkte mit einem gemeinsamen kleinsten

t+\Delta t

gleichzeitig bearbeitet werden

1) Jeder Körper hat zum gegenwärtigen Zeitpunkt $t$ einen Zeitschritt $2^{n}\Delta t_{E}$ . Die Ebene $n$ muss so gewählt werden, dass der individuelle Zeitschritt $\Delta t_{i}$ zwischen $2^{n}\Delta t_{E}$ und $2^{n+1}\Delta t_{E}$ liegt. Man greife diejenigen Massenpunkte heraus, welche gemeinsam für $t+2^{n}\Delta t_{E}$ den niedrigsten Wert aufweisen (mittlere rote Kreise). Findet sich ausnahmsweise nur ein Objekt mit diesem Minimalwert, so wird wie bisher nur dieses als Einziges betrachtet.

Die aktuell zu bearbeitenden Massenpunkte lassen sich sehr einfach identifizieren, indem man die Simulation mit der zeitlichen Schrittweite $\Delta t_{E}$ ablaufen lässt. Objekte auf der Zeitebene 0 werden bei jedem Schritt bewegt, solche auf der Ebene 1 jeden 2. Schritt, solche auf der Ebene 2 jeden 4. Schritt usw.

2) Man wende unverändert den Prädiktor auf die Positionen ${\vec {x}}_{n}$ aller Massenpunkte an, um die zum Zeitpunkt $t+2^{n}\Delta t_{E}$ gültigen Positionen ${\vec {x}}_{n+1}$ zu ermitteln (türkise Kreise). Mit den Hermite-Polynomen als grundlegendes Verfahren bestimme man entsprechend auch die Geschwindigkeiten ${\vec {v}}_{n+1}$ .

3) Aus den Positionen ${\vec {x}}_{n+1}$ gewinne man wie üblich die dazugehörigen Beschleunigungen ${\vec {a}}_{n+1}$ für die ausgewählten Körper, im Falle der Hermite-Polynome als fundamentale Lösungsformel aus den ${\vec {x}}_{n+1}$ und ${\vec {v}}_{n+1}$ auch die Rucks ${\vec {j}}_{n+1}$ .

4) Wie gehabt verbessere man die Positionen ${\vec {x}}_{n+1}$ der gerade untersuchten Massenpunkte mit Hilfe des Korrektors und anschließend deren Geschwindigkeiten ${\vec {v}}_{n+1}$ . Zuletzt aktualisiere man die Zeitebenen, auf welchen diese sich befinden. Falls erforderlich, führe man neue Ebenen ein oder entferne nicht mehr benötigte.

5) Abermals werden nur die in Schritt 1 ausgewählten Körper versetzt, alle anderen verbleiben auf ihren Positionen ${\vec {x}}_{n}$ zum Zeitpunkt $t$ . Wieder kehre man zu Schritt 1 zurück.

Wie individuelle Zeitskalen pro Massenpunkt können auch hierarchische Schrittweiten auf jedes Prädiktor-Korrektor-Verfahren angewandt werden. Beispiele hierzu finden sich z.B. bei Makino (1991) ^[1] sowie Dehnen und Read (2011) ^[2]. Im Folgenden werden nur die auf die Hermite-Polynome sowie die Leapfrog-Methode beruhenden Lösungsschemata diskutiert. Das wiederum extrem unhandliche Mehrschritt-Verfahren ist zwar erneut etwas genauer, liefert aber keine grundlegend neuen Erkenntnisse über den Algorithmus.

Prädiktor-Korrektor-Verfahren mit hierarchischen Zeitschritten

Beispiel

Eine Simulation mit hierarchischen Zeitschritten ist von gleicher Güte wie eine solche mit individuellen, wobei die Hermite-Polynome erwartungsgemäß erneut besser abschneiden als die Leapfrog-Methode. Betrachtet man abermals das schon mehrfach vorgeführte Drei-Körper-System, so befinden sich die jeweils nahestehenden Objekte zumeist beide auf der Zeitebene 0 und das dritte Mitglied auf einer höheren Ebene. Kommen sich zwei Massenpunkte sehr nahe, so wird der elementare Zeitschritt immer kleiner, was den dritten Körper (dessen dynamische Zeitskala absolut betrachtet nahezu unverändert bleibt) auf eine immer höhere Zeitebene gelangen lässt. Entfernen sich die Stoßpartner wieder voneinander, wird $\Delta t_{E}$ erneut größer, wodurch der dritte Massenpunkt auf niedrigere Ebenen zurückkehrt. Da $\Delta t_{E}$ falls erforderlich stets um den gleichbleibenden Faktor 2 geändert wird, erscheinen diese Anpassungen in logarithmischer Darstellung als Schritte konstanter Größe.

Simulation eines Drei-Körper-Systems mit dem Hermite-Polynome-Verfahren über einen Zeitraum von 100 Jahren mit hierarchischen Zeitschritten pro Massenpunkt. Die Einstufung in die jeweilige Zeitebene beruht auf dem Kriterium, dass die entsprechende Schrittweite 0.5% des Verhältnisses Beschleunigung / Ruck nicht übersteigen soll.. Für jeden Massenpunkt ist die Variation der Zeitebene im Verlauf der Simulation dargestellt. Das Auftreten hoher Ebenen korreliert mit dem Erscheinen sehr kleiner minimaler Schrittweiten infolge sehr enger Begegnungen von jeweils zwei Körpern

Simulation eines Drei-Körper-Systems mit dem Hermite-Polynome-Verfahren über einen Zeitraum von 100 Jahren mit hierarchischen Zeitschritten pro Massenpunkt. Die Einstufung in die jeweilige Zeitebene beruht auf dem Kriterium, dass die entsprechende Schrittweite 0.5% des Verhältnisses Beschleunigung / Ruck nicht übersteigen soll. Dargestellt ist die Variation des kleinsten im Ensemble auftretenden Zeitschritts

\Delta t_{E}

im Verlauf der Simulation, welcher die niedrigste Zeitebene definiert

C-Code: Hierarchische Zeitschritte mit Hermite-Polynome-Verfahren

Fast alle Variablen können dem Code für die Simulation mit individuellen Zeitschritten entnommen werden. Neu ist das eindimensionale Array ebene_zeit, welches die den dynamischen Zeiten tdyn entsprechenden Zeitebenen für jeden Massenpunkt angibt. Da andauernd Potenzen von 2 benötigt werden, werden diese mittels des Arrays potenz für die Werte von 2⁰ bis 2¹⁵ tabelliert (d.h. es werden maximal 15 Zeitebenen zugelassen). Das ebenfalls eindimensionale Array status_zeit schließlich zeigt für jeden Körper an, wann er zum letzten Mal bewegt wurde, damit die aktuell erforderliche zeitliche Schrittweite dt korrekt bestimmt werden kann.

tdynmin bezeichnet die kleinste im Ensemble vorkommende dynamische Zeitskala, dtmin den davon abgeleiteten elementaren Zeitschritt. Z zählt die im Verlauf der Simulation abgearbeiteten Zeitschritte ab. Ein Massenpunkt auf der Zeitebene $n$ muss dann bewegt werden, wenn Z durch $2^{n}$ teilbar ist. Ein Aufstieg zur nächsthöheren Ebene $n+1$ erfordert zwingend, dass Z durch $2^{n+1}$ teilbar ist.

Ändert sich der elementare Zeitschritt dtmin, muss Z entsprechend angepasst werden. Eine Halbierung von dtmin bedeutet eine Verdoppelung von Z und umgekehrt, wobei letzteres nur möglich ist, wenn Z gerade ist.

/* Globale Variablen */
 
unsigned int i,k;
unsigned int N;
 
double *m;
double ***r;
double ***v;
double ***a;
double ***j;
double s[3];
double c[3];
double *tdyn;
int *ebene_zeit;
unsigned int potenz[16];
double *status_zeit;
 
double t,dt,dt2,dt3,dt4,dt5;
double dtmin;
unsigned long Z;
double tdynmin;
double T;
 
double gamma;
double epsilon;

/* Initialisierung aller Massenpunkte aus Anfangspositionen und -geschwindigkeiten  */
/* Beschleunigungen und Rucks                                                       */
/* Dynamische Zeitskalen                                                            */

tdynmin = T;
for (i = 0;i < N;i ++)
{
  ruck (i,N,epsilon,m,r[0],v[0],a[0][i],j[0][i]);
  tdyn[i] = gamma * betrag (a[0][i]) / betrag (j[0][i]);
  if (tdyn[i] < tdynmin)
  tdynmin = tdyn[i];
}
dtmin = tdynmin / 1.5;

/* Zeitebenen                      */
/* Zeitpunkte der letzten Bewegung */

for (i = 0;i < N;i ++)
{
  k = 0;
  while (pow (2,k+1) * dtmin < tdyn[i])
  k ++;
  if (k <= 14)
  ebene_zeit[i] = k;
  else
  ebene_zeit[i] = 14;
  status_zeit[i] = 0;
}

/* Tabelle der Zweier-Potenzen */

potenz[0] = 1;
for (i = 1;i <= 15;i ++)
potenz[i] = 2 * potenz[i-1];

t = 0;
Z = 1;

while (t < T)
{

  /* Prädiktor-Schritt                                                                            */
  /* Wie bisher für alle Körper                                                                   */
  /* Zeitschritt = aktueller Zeitpunkt - Zeitpunkt der letzten Bewegung + elementarer Zeitschritt */

  for (i = 0;i < N;i ++)
  {
    dt = t - status_zeit[i] + dtmin;
    dt2 = pow (dt,2);
    dt3 = pow (dt,3);
    for (k = 0;k < 3;k ++)
    {
      r[1][i][k] = r[0][i][k] + v[0][i][k] * dt + a[0][i][k] * dt2 / 2 + j[0][i][k] * dt3 / 6;
      v[1][i][k] = v[0][i][k] + a[0][i][k] * dt + j[0][i][k] * dt2 / 2;
    }
  }

  /* Korrektor-Schritt                                                                             */
  /* Für alle Massenpunkte, für deren Zeitebene n gilt                                             */
  /* Anzahl bisher erfolgter Simulationsschritte durch 2^n teilbar                                 */
  /* Beschleunigungen und Rucks an vorläufig neuen Positionen                                      */
  /* Verbesserung der neuen Positionen und Geschwindigkeiten durch Einsetzen in die Lösungsformeln */
  /* Zeitschritt = aktueller Zeitpunkt - Zeitpunkt der letzten Bewegung + elementarer Zeitschritt  */

  for (i = 0;i < N;i ++)
  {
    if (Z % potenz[ebene_zeit[i]] == 0)
    ruck (i,N,epsilon,m,r[1],v[1],a[1][i],j[1][i]);
  }

  for (i = 0;i < N;i ++)
  {
    if (Z % potenz[ebene_zeit[i]] == 0)
    {
      dt = t - status_zeit[i] + dtmin;
      dt2 = pow (dt,2);
      dt3 = pow (dt,3);
      dt4 = pow (dt,4);
      dt5 = pow (dt,5);
      for (k = 0;k < 3;k ++)
      {
        s[k] = (- 6 * (a[0][i][k] - a[1][i][k]) - (4 * j[0][i][k] + 2 * j[1][i][k]) * dt) / dt2;
        c[k] = (12 * (a[0][i][k] - a[1][i][k]) + 6 * (j[0][i][k] + j[1][i][k]) * dt) / dt3;
        r[1][i][k] += s[k] * dt4 / 24 + c[k] * dt5 / 120;
        v[1][i][k] += s[k] * dt3 / 6 + c[k] * dt4 / 24;
      }
    }
  }

  /* Aktualisierung der soeben simulierten Objekte */

  tdynmin = T;
  for (i = 0;i < N;i ++)
  {
    if (Z % potenz[ebene_zeit[i]] == 0)
    {

  /* Verbesserte Positionen und Geschwindigkeiten */
  /* Beschleunigungen und Rucks                   */

      for (k = 0;k < 3;k ++)
      {
        r[0][i][k] = r[1][i][k];
        v[0][i][k] = v[1][i][k];
        a[0][i][k] = a[1][i][k];
        j[0][i][k] = j[1][i][k];
      }

  /* Dynamische Zeitskalen           */
  /* Zeitpunkte der letzten Bewegung */

      tdyn[i] = gamma * betrag (a[0][i]) / betrag (j[0][i]);
      if (tdyn[i] < tdynmin)
      tdynmin = tdyn[i];
      status_zeit[i] = t + dtmin;

  /* Abstieg um 1 Ebene falls neue dynamische Zeitskala < 2^n elementare Zeitschritte */
  /* Zu jedem Zeitpunkt möglich                                                       */
  /* Körper auf bisheriger Ebene 0 vorläufig nun auf Ebene -1                         */

      if (tdyn[i] < potenz[ebene_zeit[i]] * dtmin)
      ebene_zeit[i] --;

  /* Aufstieg um 1 Ebene falls neue dynamische Zeitskala > 2^(n+1) elementare Zeitschritte */
  /* Nur möglich, wenn Anzahl bisher erfolgter Simulationsschritte durch 2^(n+1) teilbar   */

      if (tdyn[i] > potenz[ebene_zeit[i]+1] * dtmin && Z % potenz[ebene_zeit[i]+1] == 0 && ebene_zeit[i] < 14)
      ebene_zeit[i] ++;
    }
  }

  /* Falls kleinste neue dynamische Zeitskala < elementarer Zeitschritt */
  /* Halbierung des elementaren Zeitschritts                            */
  /* Verdoppelung der Anzahl bisher erfolgter Simulationsschritte       */
  /* Aufstieg aller Massenpunkte um 1 Zeitebene                         */
  /* Damit Körper auf vorläufiger Ebene -1 wieder auf Ebene 0           */

  if (tdynmin < dtmin)
  {
    dtmin /= 2;
    Z *= 2;
    for (i = 0;i < N;i ++)
    {
      if (ebene_zeit[i] < 14)
      ebene_zeit[i] ++;
    }
  }

  /* Falls kleinste neue dynamische Zeitskala > 2 * elementarer Zeitschritt                    */
  /* Verdoppelung des elementaren Zeitschritts                                                 */
  /* Halbierung der Anzahl bisher erfolgter Simulationsschritte (möglich wenn durch 2 teilbar) */
  /* Abstieg aller Massenpunkte um 1 Zeitebene                                                 */

  if (tdynmin > 2 * dtmin && Z % 2 == 0)
  {
    dtmin *= 2;
    Z /= 2;
    for (i = 0;i < N;i ++)
    ebene_zeit[i] --;
  }

  t += dtmin;
  Z ++;
}

C-Code: Hierarchische Zeitschritte mit Leapfrog-Verfahren

Der Code für die Leapfrog-Methode bleibt mit demjenigen für die Hermite-Polynome weitgehend identisch, da abermals nur andere Berechnungsvorschriften für Positionen und Geschwindigkeiten gebraucht werden. Die Zuordnung einzelner Massenpunkte zu Zeitebenen geschieht ebenso auf gleiche Weise, jedoch dient als Kriterium erneut das Verhältnis $v/a$ anstelle von $a/j$ . Die Korrektor-Formeln werden weiterhin als Iteration drei Mal hintereinander ausgewertet.

/* Globale Variablen */
 
unsigned int i,k;
unsigned int N;
 
double *m;
double ***r;
double ***v;
double ***a;
double *tdyn;
int *ebene_zeit;
unsigned int potenz[16];
double *status_zeit;
 
double t,dt,dt2;
double dtmin;
unsigned long Z;
double tdynmin;
double T;
 
double gamma;
double epsilon;

/* Initialisierung aller Massenpunkte aus Anfangspositionen und -geschwindigkeiten  */
/* Beschleunigungen                                                                 */
/* Dynamische Zeitskalen                                                            */

tdynmin = T;
for (i = 0;i < N;i ++)
{
  beschleunigung (i,N,epsilon,m,r[0],a[0][i]);
  tdyn[i] = gamma * betrag (v[0][i]) / betrag (a[0][i]);
  if (tdyn[i] < tdynmin)
  tdynmin = tdyn[i];
}
dtmin = tdynmin / 1.5;

/* Zeitebenen                      */
/* Zeitpunkte der letzten Bewegung */

for (i = 0;i < N;i ++)
{
  k = 0;
  while (pow (2,k+1) * dtmin < tdyn[i])
  k ++;
  if (k <= 14)
  ebene_zeit[i] = k;
  else
  ebene_zeit[i] = 14;
  status_zeit[i] = 0;
}

/* Tabelle der Zweier-Potenzen */

potenz[0] = 1;
for (i = 1;i <= 15;i ++)
potenz[i] = 2 * potenz[i-1];

t = 0;
Z = 1;

while (t < T)
{

  /* Prädiktor-Schritt                                                                            */
  /* Für Positionen wie bisher für alle Körper                                                    */
  /* Für Geschwindigkeiten nur für die Körper mit dem gemeinsamen kleinsten neuen Zeitpunkt       */
  /* Zeitschritt = aktueller Zeitpunkt - Zeitpunkt der letzten Bewegung + elementarer Zeitschritt */

  for (i = 0;i < N;i ++)
  {
    dt = t - status_zeit[i] + dtmin;
    dt2 = pow (dt,2);
    for (k = 0;k < 3;k ++)
    {
      r[1][i][k] = r[0][i][k] + v[0][i][k] * dt + a[0][i][k] * dt2 / 2;
      if (Z % potenz[ebene_zeit[i]] == 0)
      v[1][i][k] = v[0][i][k] + a[0][i][k] * dt + j[0][i][k] * dt2 / 2;
    }
  }

  /* Korrektor-Schritt                                                                             */
  /* Für alle Massenpunkte, für deren Zeitebene n gilt                                             */
  /* Anzahl bisher erfolgter Simulationsschritte durch 2^n teilbar                                 */
  /* Beschleunigungen an vorläufig neuen Positionen                                                */
  /* Verbesserung der neuen Positionen und Geschwindigkeiten durch Einsetzen in die Lösungsformeln */
  /* Zeitschritt = aktueller Zeitpunkt - Zeitpunkt der letzten Bewegung + elementarer Zeitschritt  */

  for (z < 0;z < 3;z ++)
  {
    for (i = 0;i < N;i ++)
    {
      if (Z % potenz[ebene_zeit[i]] == 0)
      beschleunigung (i,N,epsilon,m,r[1],a[1][i]);
    }

    for (i = 0;i < N;i ++)
    {
      if (Z % potenz[ebene_zeit[i]] == 0)
      {
        dt = t - status_zeit[i] + dtmin;
        for (k = 0;k < 3;k ++)
        {
          r[1][i][k] = r[0][i][k] + (v[0][i][k] + v[1][i][k]) * dt / 2;
          v[1][i][k] = v[0][i][k] + (a[0][i][k] + a[1][i][k]) * dt / 2;
        }
      }
    }
  }

  /* Aktualisierung der soeben simulierten Objekte */

  tdynmin = T;
  for (i = 0;i < N;i ++)
  {
    if (Z % potenz[ebene_zeit[i]] == 0)
    {

  /* Verbesserte Positionen und Geschwindigkeiten */
  /* Beschleunigungen                             */

      for (k = 0;k < 3;k ++)
      {
        r[0][i][k] = r[1][i][k];
        v[0][i][k] = v[1][i][k];
        a[0][i][k] = a[1][i][k];
      }

  /* Dynamische Zeitskalen           */
  /* Zeitpunkte der letzten Bewegung */

      tdyn[i] = gamma * betrag (v[0][i]) / betrag (a[0][i]);
      if (tdyn[i] < tdynmin)
      tdynmin = tdyn[i];
      status_zeit[i] = t + dtmin;

  /* Abstieg um 1 Ebene falls neue dynamische Zeitskala < 2^n elementare Zeitschritte */
  /* Zu jedem Zeitpunkt möglich                                                       */
  /* Körper auf bisheriger Ebene 0 vorläufig nun auf Ebene -1                         */

      if (tdyn[i] < potenz[ebene_zeit[i]] * dtmin)
      ebene_zeit[i] --;

  /* Aufstieg um 1 Ebene falls neue dynamische Zeitskala > 2^(n+1) elementare Zeitschritte */
  /* Nur möglich, wenn Anzahl bisher erfolgter Simulationsschritte durch 2^(n+1) teilbar   */

      if (tdyn[i] > potenz[ebene_zeit[i]+1] * dtmin && Z % potenz[ebene_zeit[i]+1] == 0 && ebene_zeit[i] < 14)
      ebene_zeit[i] ++;
    }
  }

  /* Falls kleinste neue dynamische Zeitskala < elementarer Zeitschritt */
  /* Halbierung des elementaren Zeitschritts                            */
  /* Verdoppelung der Anzahl bisher erfolgter Simulationsschritte       */
  /* Aufstieg aller Massenpunkte um 1 Zeitebene                         */
  /* Damit Körper auf vorläufiger Ebene -1 wieder auf Ebene 0           */

  if (tdynmin < dtmin)
  {
    dtmin /= 2;
    Z *= 2;
    for (i = 0;i < N;i ++)
    {
      if (ebene_zeit[i] < 14)
      ebene_zeit[i] ++;
    }
  }

  /* Falls kleinste neue dynamische Zeitskala > 2 * elementarer Zeitschritt                    */
  /* Verdoppelung des elementaren Zeitschritts                                                 */
  /* Halbierung der Anzahl bisher erfolgter Simulationsschritte (möglich wenn durch 2 teilbar) */
  /* Abstieg aller Massenpunkte um 1 Zeitebene                                                 */

  if (tdynmin > 2 * dtmin && Z % 2 == 0)
  {
    dtmin *= 2;
    Z /= 2;
    for (i = 0;i < N;i ++)
    ebene_zeit[i] --;
  }

  t += dtmin;
  Z ++;
}

Einzelnachweise

↑ Makino J., A modified Aarseth code for GRAPE and vector processors, in: Publications of the Astronomical Society of Japan Band 43, S.859 ff, 1991
↑ Dehnen W., Read J.I., N-body simulations of gravitational dynamics, in: the European Physical Journal Plus, Band 126, 2011

[1] Makino J., A modified Aarseth code for GRAPE and vector processors, in: Publications of the Astronomical Society of Japan Band 43, S.859 ff, 1991

[2] Dehnen W., Read J.I., N-body simulations of gravitational dynamics, in: the European Physical Journal Plus, Band 126, 2011

[1]

[2]