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

Die bisher diskutierten Lösungsansätze beruhen allesamt darauf, aus dem zuletzt berechneten Zustand $n$ eines Mehrkörpersystems zu einem Zeitpunkt $t_{n}$ in einem einzigen Schritt $\Delta t$ auf den zu einem neuen Zeitpunkt $t_{n+1}$ geltenden Zustand $n+1$ 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 $n$, $n-1$, $n-2$ bis $n-k+1$ in Betracht, kann man daraus ebenfalls den neuen Zustand $n+1$ 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 $n-k+1$ bis $n$ dienen dazu, mittels einer Extrapolationsformel näherungsweise den neuen Zustand $n+1$ 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

Um die Entwicklung eines Systems sukzessive durch Extrapolation beschreiben zu können, müssen zumindest die letzten beiden Zustände $n$ und $n-1$ 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:

${\vec {r}}_{n+1}={\vec {r}}_{n}+{\vec {v}}_{n}\Delta t+{\frac {1}{2}}{\vec {a}}_{n}(\Delta t)^{2}$

Die Beschleunigung kann gemäß der Definition ${\vec {a}}=\Delta {\vec {v}}/\Delta t$ auf eine Geschwindigkeitsdifferenz zurückgeführt werden, wie untenstehende Abbildung zeigt:

${\vec {a}}_{n}={\frac {{\vec {v}}_{n}-{\vec {v}}_{n-1}}{\Delta t}}$

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 $n$, sondern auch den vorletzten $n-1$:

${\vec {r}}_{n+1}={\vec {r}}_{n}+{\frac {1}{2}}(3{\vec {v}}_{n}-{\vec {v}}_{n-1})\Delta t$

Prinzip der Extrapolation der Bewegung eines Massenpunktes nach Adams und Bashforth

Völlig analog dazu kann man aus den zurückliegenden Beschleunigungen ${\vec {a}}_{n}$ und ${\vec {a}}_{n-1}$ eine neue Geschwindigkeit ${\vec {v}}_{n+1}$ gewinnen.

${\vec {v}}_{n+1}={\vec {v}}_{n}+{\frac {1}{2}}(3{\vec {a}}_{n}-{\vec {a}}_{n-1})\Delta t$

Der neue Ort ${\vec {r}}_{n+1}$ liefert mittels des Gravitationsgesetzes auch eine neue Beschleunigung ${\vec {a}}_{n+1}$. Schließlich betrachtet man die Zustände $n$ und $n+1$ als vergangen und gelangt so zum Zustand $n+2$ 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 $n$ insgesamt $k$ sukzessive Größen Geschwindigkeit, Beschleunigung, Ruck usw. betrachten, sofern man anhand von $k$ zurückliegenden Zuständen $n-k+1$ bis $n$ auf den Zustand $n+1$ schließen will. Die entsprechenden Formeln für 3-5 zurückliegende Zustände sind nachfolgend zusammengestellt:

$k=3:{\vec {r}}_{n+1}={\vec {r}}_{n}+{\frac {1}{12}}(23{\vec {v}}_{n}-16{\vec {v}}_{n-1}+5{\vec {v}}_{n-2})\Delta t$

$k=4:{\vec {r}}_{n+1}={\vec {r}}_{n}+{\frac {1}{24}}(55{\vec {v}}_{n}-59{\vec {v}}_{n-1}+37{\vec {v}}_{n-2}-9{\vec {v}}_{n-3})\Delta t$

$k=5:{\vec {r}}_{n+1}={\vec {r}}_{n}+{\frac {1}{720}}(1901{\vec {v}}_{n}-2774{\vec {v}}_{n-1}+2616{\vec {v}}_{n-2}-1274{\vec {v}}_{n-3}+251{\vec {v}}_{n-4})\Delta t$

Die entsprechenden Beziehungen für neue Geschwindigkeiten ${\vec {v}}_{n+1}$ 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

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 ${\vec {v}}_{n+1}$ lässt sich aber dazu verwenden, die ebenso extrapolierte Position ${\vec {r}}_{n+1}$ durch eine Interpolationsformel zu korrigieren. Entsprechend erlaubt die neue Beschleunigung ${\vec {a}}_{n+1}$, die Geschwindigkeit ${\vec {v}}_{n+1}$ 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 $n-1$ bis $n+1$ muss als höchster Term der Ruck berücksichtigt werden:

${\vec {r}}_{n+1}={\vec {r}}_{n}+{\vec {v}}_{n}\Delta t+{\frac {1}{2}}{\vec {a}}_{n}(\Delta t)^{2}+{\frac {1}{6}}{\vec {j}}_{n}(\Delta t)^{3}$

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

${\vec {a}}_{n}={\frac {{\vec {v}}_{n+1}-{\vec {v}}_{n-1}}{2\Delta t}}$

Der Ruck ${\vec {j}}_{n}$ wird gemäß der Definition ${\vec {j}}=\Delta {\vec {a}}/\Delta t$ zunächst auf die Zwischenbeschleunigungen ${\vec {a}}_{n+1/2}$ und ${\vec {a}}_{-1/2}$ zurückgeführt:

${\vec {j}}_{n}={\frac {{\vec {a}}_{n+1/2}-{\vec {a}}_{n-1/2}}{\Delta t}}$

Diese beiden Beschleunigungen lassen sich wiederum durch bereits bekannten Geschwindigkeiten darstellen. Es gilt nämlich ${\vec {a}}_{n+1/2}=({\vec {v}}_{n+1}-{\vec {v}}_{n})/\Delta t$ und entsprechend ${\vec {a}}_{n-1/2}=({\vec {v}}_{n}-{\vec {v}}_{n-1})/\Delta t$. Somit kann der Ruck letztendlich wie die Beschleunigung vollständig durch Geschwindigkeiten wiedergegeben werden:

${\vec {j}}_{n}={\frac {{\vec {v}}_{n+1}-2{\vec {v}}_{n}+{\vec {v}}_{n-1}}{(\Delta t)^{2}}}$

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

${\vec {r}}_{n+1,korr}={\vec {r}}_{n}+{\frac {1}{12}}(5{\vec {v}}_{n+1}+8{\vec {v}}_{n}-{\vec {v}}_{n-1})\Delta t$

Entsprechend gilt für die verbesserte Geschwindigkeit:

${\vec {v}}_{n+1,korr}={\vec {v}}_{n}+{\frac {1}{12}}(5{\vec {a}}_{n+1}+8{\vec {a}}_{n}-{\vec {a}}_{n-1})\Delta t$

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 $k$ zurückliegenden Zuständen muss man für den Zustand $n$ jetzt $k+1$ 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:

$k=3:{\vec {r}}_{n+1}={\vec {r}}_{n}+{\frac {1}{24}}(9{\vec {v}}_{n+1}+19{\vec {v}}_{n}-5{\vec {v}}_{n-1}+{\vec {v}}_{n-2})\Delta t$

$k=4:{\vec {r}}_{n+1}={\vec {r}}_{n}+{\frac {1}{720}}(251{\vec {v}}_{n+1}+646{\vec {v}}_{n}-264{\vec {v}}_{n-1}+106{\vec {v}}_{n-2}-19{\vec {v}}_{n-3})\Delta t$

$k=5:{\vec {r}}_{n+1}={\vec {r}}_{n}+{\frac {1}{1440}}(475{\vec {v}}_{n+1}+1427{\vec {v}}_{n}-798{\vec {v}}_{n-1}+482{\vec {v}}_{n-2}-173{\vec {v}}_{n-3}+27{\vec {v}}_{n-4})\Delta t$

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 $t_{n}$ und 4 den ältesten Wert zum Zeitpunkt $t_{n-4}$ bedeutet. Für die Position müssen nur die Werte zum Zeitpunkt $t_{n}$ vorgehalten werden, so dass eindimensionale Zeigerarrays genügen. Die für einen Zeitpunkt $t_{n+1}$ 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.