Das Mehrkörperproblem in der Astronomie/ Allgemeine Lösungsmethoden/ Mehrschrittverfahren

Die bisher diskutierten Lösungsansätze beruhen allesamt darauf, aus dem zuletzt berechneten Zustand eines Mehrkörpersystems zu einem Zeitpunkt in einem einzigen Schritt auf den zu einem neuen Zeitpunkt geltenden Zustand zu schließen (auch wenn dabei je nach Verfahren mehr oder weniger Zwischenschritte benutzt werden). Man nennt diese daher Einschrittverfahren. Zieht man jedoch mehrere vergangene Zustände , , bis in Betracht, kann man daraus ebenfalls den neuen Zustand prognostizieren. Ein solches Vorgehen wird als Mehrschrittverfahren bezeichnet.

Um ein solches zu konstruieren, zieht man wiederum das Prädiktor-Korrektor-Prinzip heran. Die aktuell vorgehaltenen Zustände von bis dienen dazu, mittels einer Extrapolationsformel näherungsweise den neuen Zustand zu bestimmen. Eine Interpolationsformel liefert aus den bisherigen Zuständen und dem genäherten neuen Zustand für letzteren einen verbesserten Wert.

Extrapolation nach Adams-Bashforth

Bearbeiten

Um die Entwicklung eines Systems sukzessive durch Extrapolation beschreiben zu können, müssen zumindest die letzten beiden Zustände   und   berücksichtigt werden. Um aus diesen für die Positionen der Massenpunkte eine Extrapolationsformel zu entwickeln, geht man von der schon mehrfach verwendeten Beziehung zwischen konstanter Beschleunigung und zurückgelegter Strecke aus:

 

Die Beschleunigung kann gemäß der Definition   auf eine Geschwindigkeitsdifferenz zurückgeführt werden, wie untenstehende Abbildung zeigt:

 

Setzt man dieses oben ein, so erhält man einen Ausdruck, welche als höchste Terme nur Geschwindigkeiten enthält, aber nicht nur für den letzten Zustand  , sondern auch den vorletzten  :

 


 
Prinzip der Extrapolation der Bewegung eines Massenpunktes nach Adams und Bashforth


Völlig analog dazu kann man aus den zurückliegenden Beschleunigungen   und   eine neue Geschwindigkeit   gewinnen.

 

Der neue Ort   liefert mittels des Gravitationsgesetzes auch eine neue Beschleunigung  . Schließlich betrachtet man die Zustände   und   als vergangen und gelangt so zum Zustand   usw.

Um mittels einer Extrapolationsformel die gleiche Stabilität zu erreichen wie für eine auf dem Runge-Kutta-Verfahren oder den Hermite-Polynomen beruhende Lösung, müssen mehr als nur die beiden letzten Zustände herangezogen werden. Dies gelingt, indem man auch höhere zur zurückgelegten Strecke beitragende Terme wie Ruck, Knall und Knistern berücksichtigt und durch Differenzenbildung erneut auf Geschwindigkeiten zurückführt. Allgemein muss man für den Zustand   insgesamt   sukzessive Größen Geschwindigkeit, Beschleunigung, Ruck usw. betrachten, sofern man anhand von   zurückliegenden Zuständen   bis   auf den Zustand   schließen will. Die entsprechenden Formeln für 3-5 zurückliegende Zustände sind nachfolgend zusammengestellt:

 
 
 

Die entsprechenden Beziehungen für neue Geschwindigkeiten   sind wiederum völlig analog. Erneut braucht man anstelle der zurückliegenden Geschwindigkeiten nur die Beschleunigungen einzusetzen.

Die hier skizzierte Methode wird nach ihren Entdeckern Adams-Bashforth-Verfahren genannt. Um dieses für eine praktische Simulation in Gang zu setzen, müssen die ersten Zustände jedoch mittels eines der bisher behandelten Einschrittverfahren erzeugt werden. Die Genauigkeit dieser Startwerte sollte in der gleichen Größenordnung liegen wie diejenige der Extrapolationsformel.

Interpolation nach Adams-Moulton

Bearbeiten

Die Extrapolationsformeln nach Adams-Bashforth sind wie jedes Extrapolationsverfahren mit einer gewissen Instabilität behaftet. Die von diesen aus vergangenen Zuständen eines Mehrkörpersystems vorhergesagte neue Geschwindigkeit   lässt sich aber dazu verwenden, die ebenso extrapolierte Position   durch eine Interpolationsformel zu korrigieren. Entsprechend erlaubt die neue Beschleunigung  , die Geschwindigkeit   zu verbessern.

Um anhand alter und neuer Geschwindigkeiten eine Interpolationsvorschrift zu konstruieren, betrachtet man abermals die Strecke als Funktion von Geschwindkeit, Beschleunigung, Ruck usw. Bei drei gegebenen Zuständen von   bis   muss als höchster Term der Ruck berücksichtigt werden:

 

Die Beschleunigung   kann jetzt, da sowohl links als auch rechts von   jeweils ein Zustand zur Verfügung steht, auf "symmetrische" Weise durch Geschwindigkeiten ausgedrückt werden (siehe nachfolgendes Diagramm):

 

Der Ruck   wird gemäß der Definition   zunächst auf die Zwischenbeschleunigungen   und   zurückgeführt:

 

Diese beiden Beschleunigungen lassen sich wiederum durch bereits bekannten Geschwindigkeiten darstellen. Es gilt nämlich   und entsprechend  . Somit kann der Ruck letztendlich wie die Beschleunigung vollständig durch Geschwindigkeiten wiedergegeben werden:

 

Das Einsetzen der auf die Geschwindigkeiten zurückgeführten Beschleunigung und des Rucks liefert für die korrigierte Position:

 

Entsprechend gilt für die verbesserte Geschwindigkeit:

 

Zuletzt liefert die korrigierte Position zudem eine korrigierte Beschleunigung.


 
Prinzip der Interpolation der Bewegung eines Massenpunktes nach Adams und Moulton


Um auch bei der Interpolation mehr vergangene Zustände zu berücksichtigen, drückt man wie bei der Konstruktion der Extrapolationsformeln höhere Terme der Strecke durch Differenzen aus. Bei   zurückliegenden Zuständen muss man für den Zustand   jetzt   Größen Geschwindigkeit, Beschleunigung usw. einbeziehen. Wieder seien für 3-5 zurückliegende Zustände die Resultate vorgestellt, wobei für die korrigierten Geschwindigkeiten abermals einfach die zurückliegenden Beschleunigungen zu verwenden sind:

 
 
 

Die Korrektur mittels obiger Interpolationsformeln trägt den Namen Adams-Moulton-Verfahren. Die Anzahl der berücksichtigten alten Zustände muss selbstverständlich bei Extrapolation und Interpolation gleich sein. Um eine ausreichende Stabilität der Simulation zu erreichen, genügt eine einmalige Anwendung der Interpolation.

C-Code: Mehrkörpersimulation mit dem Mehrschrittverfahren

Der Algorithmus für das Mehrschrittverfahren zerfällt in zwei Teile. Zunächst werden mittels der Runge-Kutta-Methode 5 zurückliegende Geschwindigkeits- und Beschleunigungszustände für das Ensemble berechnet, was durch eine eigene Prozedur geschieht. Anschließend wird auf die initialen Zustände sukzessive das Mehrschrittverfahren angewandt, um immer wieder neue Werte zu generieren.

Die für das Runge-Kutta-Verfahren erforderlichen Prozeduren und Variablen wurden bereits im entsprechenden Abschnitt erläutert. Für das Mehrschrittverfahren werden Geschwindigkeit und Beschleunigung ebenfalls als zweidimensionale Zeigerarrays deklariert, wobei der erste Index nun angibt, auf welchen Zeitpunkt der Wert sich bezieht. Er läuft von 0 bis 4, wobei 0 den aktuellen Wert zum Zeitpunkt   und 4 den ältesten Wert zum Zeitpunkt   bedeutet. Für die Position müssen nur die Werte zum Zeitpunkt   vorgehalten werden, so dass eindimensionale Zeigerarrays genügen. Die für einen Zeitpunkt   gesuchten neuen Werte werden als rneu, vneu und aneu bezeichnet und ebenfalls als einfache Zeigerarrays deklariert. Schließlich ist die Unsigned Integer ordnung zu erwähnen, welche die Zustände für Geschwindigkeit und Beschleunigung abzählt.

/* Globale Variablen */

unsigned int N;
double *m;
double **r;
double **rneu;
double ***v;
double **vneu;
double ***a;
double **aneu;

unsigned int i,k;
unsigned int ordnung;

double t,dt;
double T;

/* Berechnung der ersten Zustände des Systems mittels des Runge-Kutta-Verfahrens */

rungekutta ();

while (t < T)
{

/* Prädiktor-Schritt                                                  */
/* Generieren neuer Zustände durch Extrapolation nach Adams-Bashforth */

  for (i = 0;i < N;i ++)
  {
    for (k = 0;k < 3;k ++)
    {
      rneu[i][k] = r[i][k] + (1901 * v[0][i][k] - 2774 * v[1][i][k] + 2616 * v[2][i][k] -
                              1274 * v[3][i][k] + 251 * v[4][i][k]) / 720 * dt;
      vneu[i][k] = v[0][i][k] + (1901 * a[0][i][k] - 2774 * a[1][i][k] + 2616 * a[2][i][k] -
                                 1274 * a[3][i][k] + 251 * a[4][i][k]) / 720 * dt;
    }
  }

/* Neue Beschleunigungen an neuen Positionen */

  for (i = 0;i < N;i ++)
  beschleunigung (i,N,m,rneu,aneu[i]);

/* Korrektor-Schritt - nur einmal durchlaufen                          */
/* Korrektur der neuen Zustände durch Interpolation nach Adams-Moulton */

  for (i = 0;i < N;i ++)
  {
    for (k = 0;k < 3;k ++)
    {
      rneu[i][k] = r[i][k] + (475 * vneu[i][k] + 1427 * v[0][i][k] - 798 * v[1][i][k] +
                              482 * v[2][i][k] - 173 * v[3][i][k] + 27 * v[4][i][k]) / 1440 * dt;
      vneu[i][k] = v[0][i][k] + (475 * aneu[i][k] + 1427 * a[0][i][k] - 798 * a[1][i][k] +
                                 482 * a[2][i][k] - 173 * a[3][i][k] + 27 * a[4][i][k]) / 1440 * dt;
    }
  }

/* Korrigierte Beschleunigungen an korrigierten Positionen */

  for (i = 0;i < N;i ++)
  beschleunigung (i,N,m,rneu,aneu[i]);

/* Werte(n-4) bekommen die Werte(n-3) zugewiesen, Werte(n-3) die Werte(n-2) usw. */

  for (i = 0;i < N;i ++)
  {
    for (ordnung = 4;ordnung > 0;ordnung --)
    {
      for (k = 0;k < 3;k ++)
      {
        v[ordnung][i][k] = v[ordnung-1][i][k];
        a[ordnung][i][k] = a[ordnung-1][i][k];
      }
    }
  }

/* Werte(n) bekommen die Werte(n+1) zugewiesen */

  for (i = 0;i < N;i ++)
  {
    for (k = 0;k < 3;k ++)
    {
      r[i][k] = rneu[i][k];
      v[0][i][k] = vneu[i][k];
      a[0][i][k] = aneu[i][k];
    }
  }

  t += dt;
}

Die für die Zwischenwerte des Runge-Kutta-Verfahrens erforderlichen Variablen werden nun lokal deklariert und ihre Verwendung nur für diesen Zweck durch die Bezeichnung "rk" kenntlich gemacht.

void rungekutta (void)
{
  double ***rrk;
  double ***vrk;
  double ***ark;

/* Vorhalten der Positionen x(n) und Geschwindigkeiten v(n-4) */

  for (i = 0;i < N;i ++)
  {
    for (k = 0;k < 3;k ++)
    {
      rrk[0][i][k] = r[i][k];
      vrk[0][i][k] = v[4][i][k];
    }
  }

  for (ordnung = 4;ordnung > 0;ordnung --)
  {
    for (i = 0;i < N;i ++)
    beschleunigung (i,N,m,rrk[0],ark[0][i]);

/* Vorhalten der Beschleunigungen von a(n-4) bis a(n-1) */

    for (i = 0;i < N;i ++)
    {
      for (k = 0;k < 3;k ++)
      a[ordnung][i][k] = ark[0][i][k];
    }

/* Rechenschritte des Runge-Kutta-Verfahrens */

    for (i = 0;i < N;i ++)
    {
      for (k = 0;k < 3;k ++)
      {
        rrk[1][i][k] = rrk[0][i][k] + vrk[0][i][k] * dt / 2;
        vrk[1][i][k] = vrk[0][i][k] + ark[0][i][k] * dt / 2;
      }
    }
    for (i = 0;i < N;i ++)
    beschleunigung (i,N,m,rrk[1],ark[1][i]);

    for (i = 0;i < N;i ++)
    {
      for (k = 0;k < 3;k ++)
      {
        rrk[2][i][k] = rrk[0][i][k] + vrk[1][i][k] * dt / 2;
        vrk[2][i][k] = vrk[0][i][k] + ark[1][i][k] * dt / 2;
      }
    }
    for (i = 0;i < N;i ++)
    beschleunigung (i,N,m,rrk[2],ark[2][i]);

    for (i = 0;i < N;i ++)
    {
      for (k = 0;k < 3;k ++)
      {
        rrk[3][i][k] = rrk[0][i][k] + vrk[2][i][k] * dt;
        vrk[3][i][k] = vrk[0][i][k] + ark[2][i][k] * dt;
      }
    }
    for (i = 0;i < N;i ++)
    beschleunigung (i,N,m,rrk[3],ark[3][i]);

    for (i = 0;i < N;i ++)
    {
      for (k = 0;k < 3;k ++)
      {
        rrk[0][i][k] += (vrk[0][i][k] + 2 * vrk[1][i][k] + 2 * vrk[2][i][k] + vrk[3][i][k]) * dt / 6;
        vrk[0][i][k] += (ark[0][i][k] + 2 * ark[1][i][k] + 2 * ark[2][i][k] + ark[3][i][k]) * dt / 6;
      }
    }

/* Vorhalten der Geschwindigkeiten von v(n-3) bis v(n) */

    for (i = 0;i < N;i ++)
    {
      for (k = 0;k < 3;k ++)
      v[ordnung-1][i][k] = vrk[0][i][k];
    }
  }

  for (i = 0;i < N;i ++)
  beschleunigung (i,N,m,rrk[0],ark[0][i]);

/* Vorhalten der Positionen x(n) und der Beschleunigungen a(n) */

  for (i = 0;i < N;i ++)
  {
    for (k = 0;k < 3;k ++)
    {
      r[i][k] = rrk[0][i][k];
      a[ordnung][i][k] = ark[0][i][k];
    }
  }
}

Die Mehrschrittmethode spielt in der Praxis astronomischer Mehrkörpersimulationen nur eine geringe Rolle. Als Prädiktor-Korrektor-Verfahren gestattet sie zwar die Anwendung dynamischer Zeitschritte, doch im Vergleich zur Leapfrog- oder Hermite-Polynome-Methode mit weit höherem mathematischen Aufwand.