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

Aus Wikibooks
Zur Navigation springen Zur Suche springen

Konstruktion[Bearbeiten]

Wie soeben angedeutet, lassen sich die dynamischen Zeitskalen der Massenpunkte eines Mehrkörpersystems in eine hierarchische Struktur einordnen, indem man einen kleinsten elemenaren Zeitschritt 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 auf einer Zeitebene zugeordnet, wobei die kleinste im Ensemble vorkommende Schrittweite bedeutet


Um zu Beginn der Simulation festzulegen, muss man sich an der kleinsten anfänglich vorliegenden dynamsichen Zeitskala ortientieren. repräsentiert die Ebene 0, die niedrigste der Hierarchie. Durch Verdoppelnung legt man die nächsthöhere Ebene 1 an, durch weitere Verdopplelung die Ebene 2 usw.

Fällt im Laufe der Modellierung die dynamische Zeitskala eines Körpers unter den Wert, welcher seine Zeitebende definiert, so muss er um eine Ebene nach unten geschoben werden. Tritt z.B. für ein Objekt auf Ebene 1 ein auf, so wechselt dieses zur Ebene 0. Wird sogar ein gefunden, so ist die bisherige elementare Schrittweite nicht länger gültig. An ihrer Stelle muss ein neuer Elementarschritt 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 so festzulegen, dass mit einem Betrag von genau zwischen und 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 , 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 , , usw., nicht aber zu Zeitpunkten , 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 die auf der Ebene 1 gültige Schrittweite 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 . muss solange verdoppelt werden, bis die so erzeugte Zeitskala größer als die Hälfte von ist.

Wie die niedrigste ist auch die höchste Ebene dynamisch. Überschreitet für einen Massenpunkt 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[Bearbeiten]

Die Steigerung der Effizienz des Prädiktor-Korrektor-Algorithmus beruht darauf, nicht mehr die individuellen Zeitschritte der einzelnen Massenpunkte zu verwenden, sondern diejenigen , , 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 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 siehr 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 annehmen kann, fast immer mehrere Massenpunkte mit einem gemeinsamen kleinsten gleichzeitig bearbeitet werden


1) Jeder Körper hat zum gegenwärtigen Zeitpunkt einen Zeitschritt . Die Ebene muss so gewählt werden, dass der individuelle Zeitschritt zwischen und liegt. Man greife diejenigen Massenpunkte heraus, welche gemeinsam für 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, in dem man die Simulation mit der zeitlichen Schrittweite 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 aller Massenpunkte an, um die zum Zeitpunkt gültigen Positionen zu ermitteln (türkise Kreise). Mit den Hermite-Polynomen als grundlegendes Verfahren bestimme man entsprechend auch die Geschwindigkeiten .

3) Aus den Positionen gewinne man wie üblich die dazugehörigen Beschleunigungen für die ausgewählten Körper, im Falle der Hermite-Polynome als fundamentale Lösungsformel aus den und auch die Rucks .

4) Wie gehabt verbessere man die Positionen der gerade untersuchten Massenpunkte mit Hilfe des Korrektors und anschließend deren Geschwindigkeiten . 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 zum Zeitpunkt . 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 neue Erkenntnisse über den Algorithmus.

Prädiktor-Korrektor-Verfahren mit hierarchischen Zeitschritten[Bearbeiten]

Beispiel[Bearbeiten]

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 erneut größer, wodurch der dritte Massenpunkt auf niedrigere Ebenen zurückkehrt. Da 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 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 20 bis 215 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 muss dann bewegt werden, wenn Z durch teilbar ist. Ein Aufstieg zur nächsthöheren Ebene erfordert zwingend, dass Z durch 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 anstelle von . 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

  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