Das Mehrkörperproblem in der Astronomie/ Allgemeine Lösungsmethoden/ Prädiktor-Korrektor-Verfahren: Leapfrog und Hermite-Polynome

Idee

Die Schwäche des Euler-Verfahrens, während eines Zeitintervalls $\Delta t$ Geschwindigkeiten und Beschleunigungen konstant zu halten, lässt sich nicht nur durch diskrete Zwischenwerte überwinden. Man kann stattdessen diese Größen während $\Delta t$ sich auch kontinuierlich ändern lassen, ohne Zwischenschritte definieren zu müssen. Wie schon im Grundlagenkapitel dargelegt, gilt unter Annahme einer konstanten Beschleunigung:

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

${\vec {v}}_{neu}={\vec {v}}_{alt}+{\vec {a}}_{alt}\Delta t$

Damit lässt man zumindest eine lineare Änderung der Geschwindigkeit während eines Zeitintervalls zu. Ein konstanter Ruck gestattet für die Geschwindigkeit eine quadratische und für die Beschleunigung eine lineare Abhängigkeit von der Zeit:

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

${\vec {v}}_{neu}={\vec {v}}_{alt}+{\vec {a}}_{alt}\Delta t+{\frac {1}{2}}{\vec {j}}_{alt}(\Delta t)^{2}$

Wie nun aufgezeigt wird, weisen diese Versuche tatsächlich auf den richtigen Weg, sind aber so noch nicht präzise genug. Ohne Ruck bleibt wie für das Euler-Verfahren der relative Fehler eines Zeitschritts proportional zu $(\Delta t)^{2}$ und derjenige einer ganzen Simulation proportional zu $\Delta t$. Mit Ruck ergeben sich Proportionalitäten zu $(\Delta t)^{3}$ bzw. $(\Delta t)^{2}$.

Die von obigen Formeln gelieferten neuen Positionen und Geschwindigkeiten stellen für die vollständigen Lösungen dennoch korrekte erste Prognosen dar. Dementsprechend werden diese Ausdrücke als Prädiktoren bezeichnet. Die vorläufigen Resultate für die ${\vec {x}}_{neu}$ und ${\vec {v}}_{neu}$ werden anschließend in Korrekturformeln eingesetzt, woraus sich verbesserte Werte ergeben. Diese den eigentlichen Kern der Verfahren ausmachenden Beziehungen werden Korrektoren genannt.

Leapfrog-Verfahrens

Alternative Formulierung

Wie Hut und andere (1995) zeigten, lässt sich aus dem Leapfrog-Verfahren recht einfach eine Korrektur für die auf eine pro Zeitschritt konstante Beschleunigung basierende Lösung ableiten. Wie im letzten Unterkapitel erläutert, lautet die klassische Formulierung:

${\vec {r}}_{zwischen}={\vec {r}}_{alt}+{\frac {1}{2}}{\vec {v}}_{alt}\Delta t$

${\vec {v}}_{neu}={\vec {v}}_{alt}+{\vec {a}}_{zwischen}\Delta t$

${\vec {r}}_{neu}={\vec {r}}_{zwischen}+{\frac {1}{2}}{\vec {v}}_{neu}\Delta t$

Setzt man die erste Gleichung in die letzte ein, so lässt sich die Zwischenposition eliminieren, wodurch man sofort den Korrektor für die Position erhält. Die Zwischenbeschleunigung in der zweiten Gleichung kann man nicht durch algebraische Umstellungen loswerden, doch darf man diese näherungsweise durch das Mittel der Beschleunigungen an alter und vorläufiger neuer Position ersetzen. Damit ergibt sich auch ein Korrektor für die Geschwindigkeit, welcher hinsichtlich der Form der Korrektur der Position völlig analog ist:

${\vec {r}}_{neu,korr}={\vec {r}}_{alt}+{\frac {1}{2}}({\vec {v}}_{alt}+{\vec {v}}_{neu})\Delta t$

${\vec {v}}_{neu,korr}={\vec {v}}_{alt}+{\frac {1}{2}}({\vec {a}}_{alt}+{\vec {a}}_{neu})\Delta t$

Im letzten Unterkapitel vorgestellte Testrechnungen lassen auf den ersten Blick vermuten, dass das Leapfrog-Verfahren in Prädiktor-Korrektor-Gestalt wesentlich ungenauer ist als in Zwischenschritt-Form. Man kann die Resultate der Korrekturformeln jedoch wiederum in diese einsetzen, d.h. sie im Rahmen einer Iteration mehrmals durchlaufen und so die Genauigkeit im Vergleich zum Zwischenschritt-Algorithmus vielmehr enorm steigern. Mit drei Durchläufen ist die optimale Balance zwischen nomineller Präzision und Rundungsfehlern erreicht.

C-Code: Mehrkörpersimulation mit den Leapfrog-Verfahren

Durch den Wegfall aufwendiger Konstruktionen von Zwischenschritten ist die Programmierung von Prädiktor-Korrektor-Methoden von Natur aus recht übersichtlich. Da neue Positionen und Geschwindigkeiten durch Iteration ermittelt werden müssen, werden aber für alle Größen von der Position bis zur Beschleunigung generell Datensätze sowohl für alte als auch neue Werte benötigt. Wie für das Runge-Kutta-Verfahren müssen sie daher als zweidimensionale Zeigerarrays deklariert werden, wobei der erste Index jetzt aber lediglich die beiden Werte 0 und 1 für alt und neu annimmt. Die wiederholt erforderlichen Berechnungen von Beschleunigungen werden wie üblich durch die Beschleunigung-Prozedur ausgeführt.

Gegenüber den bisherigen Programmierabschnitten treten als neue Variablen die Double dt2 und die Unsigned integer z auf. dt2 steht für das Quadrat einer Zeitspanne $\Delta t$. z stellt einen weiteren Zähler dar, welcher angibt, wie oft die Korrekturformeln durchlaufen wurden.

/* Globale Variablen */
 
unsigned int N;
double *m;
double ***r;
double ***v;
double ***a;
 
unsigned int i,k,z;
 
double t,dt,dt2;
double T;

dt2 = pow (dt,2);
for (t = 0;t < T;t += dt)
{
  /* Beschleunigungen an alten Positionen */

  for (i = 0;i < N;i ++)
  beschleunigung (i,N,m,r[0],a[0][i]);  

  /* Prädiktor-Schritt                                                                           */
  /* Erste Näherung für neue Positionen nach dem Schema rneu = ralt + valt * dt + aalt * dt2 / 2 */ 
  /* Erste Näherung für neue Geschwindigkeiten nach dem Schema vneu = valt + aalt * dt           */ 

  for (i = 0;i < N;i ++)
  {
    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;
      v[1][i][k] = v[0][i][k] + a[0][i][k] * dt;
    }
  }

  /* Korrektor-Schritt - drei Mal durchlaufene Iteration                                           */
  /* Beschleunigung an vorläufig neuen Positionen                                                  */
  /* Verbesserung der neuen Positionen und Geschwindigkeiten durch Einsetzen in die Lösungsformeln */

  for (z = 0;z < 3;z ++)
  {
    for (i = 0;i < N;i ++)
    beschleunigung (i,N,m,r[1],a[1][i]);

    for (i = 0;i < N;i ++)
    {
      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;
      }
    }
  }

  /* Übernahme der verbesserten neuen Positionen und Geschwindigkeiten */

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

Mit der Formulierung als Prädiktor-Korrektor-Algorithmus ist das Leapfrog-Verfahren als Hilfsmittel astronomischer Mehrkörperbetrachtungen weit verbreitet. Auch mit einer mehrmaligen Iteration bleibt der Rechenaufwand akzeptabel. Vor allem aber lässt sich, wie im folgenden Kapitel dargelegt, die Leapfrog-Methode nun ohne Verlust an Stabilität einer Simulation auf dynamische Zeitschritte anwenden.

Hermite-Polynome

Klassischer Korrektor

Um die auf während eines Zeitintervalls $\Delta t$ konstanten Rucks basierende Erstprognose zu korrigieren, muss man auch höhere Terme berücksichtigen, welche noch über den Ruck hinausgehen. Diese sind der Knall und das Knistern, welche aufgrund ihrer englischen Benennungen "snap" und "crackle" die Symbole $s$ und $c$ tragen. Unter dem mittleren Knall $<{\vec {s}}>$ versteht man die Änderung des Rucks $\Delta {\vec {j}}$, unter dem mittleren Knistern $<{\vec {c}}>$ die Änderung des Knalls $\Delta {\vec {s}}$ im Verlauf von $\Delta t$.

$<{\vec {s}}>={\frac {\Delta {\vec {j}}}{\Delta t}}$

$<{\vec {c}}>={\frac {\Delta {\vec {s}}}{\Delta t}}$

Mit Hilfe dieser Definitionen lässt sich folgenden Ansatz verwenden, welcher jedoch nur mit höherer Mathematik begründet werden kann:

${\vec {r}}_{neu}={\vec {r}}_{alt}+{\vec {v}}_{alt}\Delta t+{\frac {1}{2}}{\vec {a}}_{alt}(\Delta t)^{2}+{\frac {1}{6}}{\vec {j}}_{alt}(\Delta t)^{3}+{\frac {1}{24}}{\vec {s}}_{alt}(\Delta t)^{4}$

${\vec {v}}_{neu}={\vec {v}}_{alt}+{\vec {a}}_{alt}\Delta t+{\frac {1}{2}}{\vec {j}}_{alt}(\Delta t)^{2}+{\frac {1}{6}}{\vec {s}}_{alt}(\Delta t)^{3}+{\frac {1}{24}}{\vec {c}}_{alt}(\Delta t)^{4}$

Schon die Berechnung des von gegenseitig sich anziehenden Massen ausgeübten Rucks ist relativ kompliziert. Zwar lassen sich auch für den aus dem Gravitationsgesetz folgenden Knall und das Knistern Formeln herleiten, doch sind diese so umfangreich, dass eine Auswertung mit vertretbarem Rechenaufwand nicht mehr möglich ist. Doch lassen sich diese äußerst unhandlichen Größen eliminieren, wenn man zusätzlich folgende Beziehungen verwendet:

${\vec {a}}_{neu}={\vec {a}}_{alt}+{\vec {j}}_{alt}\Delta t+{\frac {1}{2}}{\vec {s}}_{alt}(\Delta t)^{2}+{\frac {1}{6}}{\vec {c}}_{alt}(\Delta t)^{3}$

${\vec {j}}_{neu}={\vec {j}}_{alt}+{\vec {s}}_{alt}\Delta t+{\frac {1}{2}}{\vec {c}}_{alt}(\Delta t)^{2}$

Die 4 Formeln für die gesuchten 4 neuen Werte für ${\vec {r}}$, ${\vec {v}}$, ${\vec {a}}$ und ${\vec {j}}$ bilden ein vollständiges lineares Gleichungssystem, wodurch man ${\vec {s}}$ und ${\vec {c}}$ eliminieren kann. Das Ergebnis wirkt erstaunlich elegant:

${\vec {r}}_{neu,korr}={\vec {r}}_{alt}+{\frac {1}{2}}({\vec {v}}_{alt}+{\vec {v}}_{neu})\Delta t+{\frac {1}{12}}({\vec {a}}_{alt}-{\vec {a}}_{neu})(\Delta t)^{2}$

${\vec {v}}_{neu,korr}={\vec {v}}_{alt}+{\frac {1}{2}}({\vec {a}}_{alt}+{\vec {a}}_{neu})\Delta t+{\frac {1}{12}}({\vec {j}}_{alt}-{\vec {j}}_{neu})(\Delta t)^{2}$

Um von alten zu neuen Positionen zu gelangen, braucht man nur Geschwindigkeiten und Beschleunigungen. Für neue Geschwindigkeiten werden lediglich Beschleunigungen und Rucks benötigt, welche noch mit akzeptablem Aufwand ermittelt werden können.

Die neuen Korrekturformeln schließen diejenigen des Leapfrog-Verfahrens mit ein, enthalten aber darüber hinaus jeweils einen auf der Beschleunigung bzw. dem Ruck beruhenden Zusatzterm. Dieser sorgt dafür, dass der relative Fehler eines Simulationsschritts wie für das Runge-Kutta-Verfahren zu $(\Delta t)^{5}$ proportional ist, und nicht wie für die Leapfrog-Methode nur proportional zu $(\Delta t)^{4}$.

Einschub für Fortgeschrittene: Knall und Knistern als Taylor-Entwicklung der Strecke

Nicht nur Geschwindigkeit, Beschleunigung und Ruck, sondern auch Knall und Knistern lassen sich als Zeitableitungen der Strecke deuten. Dazu wendet man abermals den Taylorschen Satz an und berücksichtigt jetzt die ersten 5 Glieder der Reihe. Der Knall stellt sich dann als die 4., das Knistern als die 5. Zeitableitung der Strecke heraus. Die Taylor-Entwicklung erklärt auch die Ansätze zur Ableitung der Lösungsformeln, denn es gilt:

$f(t)={\frac {df}{dt}}t+{\frac {1}{2}}{\frac {d^{2}f}{dt^{2}}}t^{2}+{\frac {1}{6}}{\frac {d^{3}f}{dt^{3}}}t^{3}+{\frac {1}{24}}{\frac {d^{4}f}{dt^{4}}}t^{4}+{\frac {1}{120}}{\frac {d^{5}f}{dt^{5}}}t^{5}+...$

Die hier vorgestellten Lösungsformeln gehören einer speziellen Klasse von Polynomen, den sogenannten Hermite-Polynome an. Sie wurden ursprünglich von dem Mathematiker Hermite zur Lösung der nach ihm benannten hermiteschen Differentialgleichung eingeführt, haben aber seitdem vor allem in der Physik weitere Anwendungen gefunden.

Erweiterung durch Makino und Aarseth

Das soeben vorgestellte, in der Fachliteratur oft beschriebene und somit schon klassisch zu nennende Verfahren lässt sich gemäß Makino und Aarseth (1992) noch verbessern, indem man das Knistern nicht nur für die Bestimmung neuer Geschwindigkeiten, sondern auch neuer Positionen heranzieht:

${\vec {r}}_{neu}={\vec {r}}_{alt}+{\vec {v}}_{alt}\Delta t+{\frac {1}{2}}{\vec {a}}_{alt}(\Delta t)^{2}+{\frac {1}{6}}{\vec {j}}_{alt}(\Delta t)^{3}+{\frac {1}{24}}{\vec {s}}_{alt}(\Delta t)^{4}+{\frac {1}{120}}{\vec {c}}_{alt}(\Delta t)^{4}$

${\vec {v}}_{neu}={\vec {v}}_{alt}+{\vec {a}}_{alt}\Delta t+{\frac {1}{2}}{\vec {j}}_{alt}(\Delta t)^{2}+{\frac {1}{6}}{\vec {s}}_{alt}(\Delta t)^{3}+{\frac {1}{24}}{\vec {c}}_{alt}(\Delta t)^{4}$

Dies hat etwas kompliziertere Korrekturformeln zur Folge, welche dafür aber nur einmal durchlaufen werden müssen, um eine selbst über die Runge-Kutta-Methode noch hinausgehende Genauigkeit zu erzielen. Sie lauten:

${\vec {r}}_{neu,korr}={\vec {r}}_{neu}+{\frac {1}{24}}{\vec {s}}(\Delta t)^{4}+{\frac {1}{120}}{\vec {c}}(\Delta t)^{5}$

${\vec {v}}_{neu,korr}={\vec {v}}_{neu}+{\frac {1}{6}}{\vec {s}}(\Delta t)^{3}+{\frac {1}{24}}{\vec {c}}(\Delta t)^{4}$

Für den Knall ${\vec {s}}$ und das Knistern ${\vec {c}}$ geben Makino und Aarseth (1992) folgende Ausdrücke an:

${\vec {s}}={\frac {-6({\vec {a}}_{alt}-{\vec {a}}_{neu})-(4{\vec {j}}_{alt}+2{\vec {j}}_{neu})\Delta t}{(\Delta t)^{2}}}$

${\vec {c}}={\frac {12({\vec {a}}_{alt}-{\vec {a}}_{neu})+6({\vec {j}}_{alt}+{\vec {j}}_{neu})\Delta t}{(\Delta t)^{3}}}$

C-Code: Mehrkörpersimulation mit den Hermite-Polynomen

Die Programmierung des Hermite-Polynome-Verfahrens ist zur Leapfrog-Methode völlig analog, im Wesentlichen müssen nur erweiterte Formeln für Erstprognose und Korrektur eingesetzt werden. Letztere muss lediglich ein einziges Mal durchgeführt werden, so dass keine Iterationsschleife erforderlich ist. Beschleunigungen und Rucks werden wie üblich berechnet.

Die neuen Größen Knall und Knistern werden jeweils nur für den gerade abgearbeiteten Massenpunkt gebraucht, so dass einfache Vektorrrays s und c für Double ausreichen. Für die immer wiederkehrenden Potenzen der zeitlichen Schrittweite dt werden die Double dt3 bis dt5 verwendet.

/* Globale Variablen */
 
unsigned int N;
double *m;
double ***r;
double ***v;
double ***a;
double ***j;
double s[3];
double c[3];
 
unsigned int i,k;
 
double t,dt,dt2,dt3,dt4,dt5;
double T;

dt2 = pow (dt,2);
dt3 = pow (dt,3);
dt4 = pow (dt,4);
dt5 = pow (dt,5);
for (t = 0;t < T;t += dt)
{
  /* Beschleunigungen und Rucks an alten Positionen */

  for (i = 0;i < N;i ++)
  ruck (i,N,m,r[0],v[0],a[0][i],j[0][i]);  

  /* Prädiktor-Schritt                                                                                            */
  /* Erste Näherung für neue Positionen nach dem Schema rneu = ralt + valt * dt + aalt * dt2 / 2 + jalt * dt3 / 6 */ 
  /* Erste Näherung für neue Geschwindigkeiten nach dem Schema vneu = valt + aalt * dt + jalt * dt2 / 2           */ 

  for (i = 0;i < N;i ++)
  {
    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 - nur einmal durchlaufen                                                    */
  /* Beschleunigung und Ruck an vorläufig neuen Positionen                                         */
  /* Verbesserung der neuen Positionen und Geschwindigkeiten durch Einsetzen in die Lösungsformeln */

  for (i = 0;i < N;i ++)
  ruck (i,N,m,r[1],v[1],a[1][i],j[1][i]);

  for (i = 0;i < N;i ++)
  {
    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;
    }
  }

  /* Übernahme der verbesserten neuen Positionen und Geschwindigkeiten */

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

Die Hermite-Polynome-Methode wird trotz ihrer Genauigkeit seltener verwendet als das Leapfrog-Verfahren. Um zu vermeiden, dass der Aufwand zur Berechnung der Anziehungskräfte quadratisch mit der Anzahl der Massenpunkte wächst, werden in der Praxis Näherungen angewandt, welche der hohen Präzision der Lösungsformeln zuwiderlaufen. Der durch den Ruck bedingte doch noch erhebliche Rechenaufwand ist damit nicht mit einem Mehrwert verbunden.