Das Mehrkörperproblem in der Astronomie/ Allgemeine Lösungsmethoden/ Euler-Verfahren

Konstruktion

Bereits der Mathematiker Euler entwickelte das Konzept, ein zu komplexes System in diskreten Zeitschritten $\Delta t$ zu beobachten, anstatt nach einer kontinuierlichen, für beliebige Zeiten $t$ geltenden Lösung zu suchen. Sind diese Schritte klein genug, darf man annehmen, dass die Beschleunigungen und Geschwindigkeiten der Massen unterdessen sich kaum ändern und so während des Intervalls $\Delta t$ näherungsweise als konstant betrachten. Man geht also davon aus, dass die zu einem bestimmten Zeitpunkt $t$ gegebenen Momentanbeschleunigungen und -geschwindigkeiten den mittleren Werten für den gesamten Zeitrahmen $[t,t+\Delta t]$ entsprechen. Die Entwicklung des Systems von einem Zeitpunkt zum nächsten lässt sich dann mit nachfolgendem Schema beschreiben. Mit der neuen Position kennt man auch die neue Beschleunigung und kann so den nächsten Schritt vollziehen.

Schrittweise Berechnung der Positionen und Geschwindigkeiten eines Massenpunktes in einem Kraftfeld nach dem Euler-Verfahren

Obiger Ansatz bedeutet, dass man die wahre Geschwindigkeitskurve durch eine Stufenfunktion annähert und dementsprechend die darunterliegende Fläche (welche wie im letzten Kapitel besprochen die zurückgelegte Strecke darstellt) durch Rechtecke. Die Breite eines solchen Rechtecks ist durch den Zeitschritt $\Delta t$ gegeben, seine Höhe durch die Geschwindigkeit jeweils zu Beginn eines solchen Schritts. Untenstehende Abbildung lässt bereits erahnen, dass insbesondere bei plötzlichen Änderungen der Geschwindigkeit auf diese Weise die tatsächliche Fläche unter der Kurve nur schlecht ausgefüllt wird, d.h. die entsprechende Vorhersage für die Bahn eines Massenpunktes sehr ungenau ist.

Berechnung der zurückgelegten Strecke in Abh. von der Geschwindigkeit durch das Euler-Verfahren. Der wahre Geschwindigkeitsverlauf wird durch eine Stufenfunktion angenähert

Auf analoge Weise werden unter der Beschleunigungskurve Rechtecke platziert, um die Änderung der Geschwindigkeit während eines Zeitschritts zu ermitteln. Es zeigt sich die gleiche Schwäche wie bei der Betrachtung der zurückgelegten Entfernung.

Berechnung der Geschwindigkeitsänderung in Abh. von der Beschleunigung durch das Euler-Verfahren. Der wahre Beschleunigungsverlauf wird durch eine Stufenfunktion angenähert

Genauigkeit

Die unzureichende Genauigkeit des Euler-Verfahrens lässt sich anhand der Kreisbahn leicht demonstrieren. Abermals sei eine sehr kleine Masse $m$ betrachtet, welche sich auf einem Kreis mit Radius $R$ um eine viel größere Masse $M$ bewegt.

Behandlung einer Kreisbahn durch das Euler-Verfahren

Gemäß der Zeichnung soll sich $m$ zu Beginn genau rechts von $M$ aufhalten. Die Geschwindigkeit der kleinen Masse weist dann exakt nach oben (oder unten). Es ist einleuchtend, dass $m$ sich schon nach dem ersten Schritt außerhalb der tatsächlichen Kreisbahn befindet, in einer zu großen Entfernung $R+\Delta R$. Der Fehler $\Delta R$ ist beträchtlich, wie man z.B. für das System Erde-Sonne zeigen kann. Wählt man als Schrittweite 1 Tag, so liegt der Fehler etwa bei 0.015 %. des tatsächlichen Bahnradius. Da die Fehler der einzelnen Schritte dazu tendieren, sich aufzusummieren (wie im letzten Unterkapitel präsentierte praktische Beispiele zeigen), liegt die gesamte relative Abweichung nach 1 Jahr schon etwa bei 5 %. Es liegt auf der Hand, dass ein solches Ergebnis völlig inakzeptabel ist.

Mit kleineren Zeitschritten lässt sich die Genauigkeit des Verfahrens steigern, doch erwächst aus einem solchen Versuch ein neues Dilemma. Eine detaillierte Betrachtung zeigt nämlich für den Fehler pro Schritt folgenden Zusammenhang:

${\frac {\Delta R}{R}}\sim (\Delta t)^{2}$

Der Fehler des Verfahrens hängt nur quadratisch von $\Delta t$ ab, mit z.B. halb so großen Zeitintervallen lässt sich nur die vierfache Genauigkeit pro Schritt erzielen. Verringert man $\Delta t$ allgemein um einen Faktor $n$, so weist jeder Einzelschritt eine um $n^{2}$ verbesserte Genauigkeit auf. Jedoch werden dann auch $n$ Mal so viele Schritte benötigt, um die Entwicklung des Systems über eine bestimmte Gesamtzeit zu verfolgen. Die Genauigkeit der Gesamtsimulation wächst so lediglich linear um das $n$-fache. Um etwa die Erdbahn mit einem jährlichen relativen Fehler von 0.1 % statt 5 % wiederzugeben, wäre eine auf 1/50 reduzierte Schrittweite nötig, d.h. von 1/50 Tag ≈ 1/2 Stunde. Ein Fehler von nur ein Millionstel des wirklichen Radius pro Jahr erforderte gar noch 1000 Mal kleinere Schritte von ungefähr 2 Sekunden. Für eine Simulation mit akzeptablen Rechenzeit wächst die Genauigkeit des Euler-Verfahrens viel zu langsam mit kleinerer Schrittweite.

Der prozentuale Fehler hängt nicht nur von $\Delta t$ selbst, sondern auch dem Abstand der beiden Massen ab. Hierbei gilt:

${\frac {\Delta R}{R}}\sim {\frac {1}{R^{3}}}$

Er verhält sich dem Kubus von $R$ umgekehrt proportional, nimmt also dramatisch zu, wenn zwei Massenpunkte sich sehr nahe kommen. Beträgt ihr Abstand die Hälfte des ursprünglichen Wertes, liegt der Fehler des Verfahrens 8 Mal höher als zu Beginn. Bei einem Drittel der anfänglichen Distanz ist er schon auf das 27-fache des Ausgangswertes angewachsen. Die außerordentlich schnell wachsende Unsicherheit der Lösung im Falle einer engen Begegnung zweier Körper stellt eine grundsätzliche Schwierigkeit dar, welche auch bei genaueren Näherungsverfahren auftritt. Die daraus resultierende Instabilität eines simulierten Sternsystems kommt später noch ausführlich zur Sprache.

Mit Hilfe des 3. Keplerschen Gesetzes lässt sich eine weitere wichtige Aussage über den Fehler des Euler-Verfahrens bei der Kreisbahn gewinnen. Aus der Proportionalität $T^{2}\sim R^{3}$ folgt für den relativen Fehler:

${\frac {\Delta R}{R}}\sim \left({\frac {\Delta t}{T}}\right)^{2}$

Für die Genauigkeit der Simulation ist das Verhältnis zwischen zeitlicher Schrittweite und Umlaufszeit maßgeblich. Dies gilt auch für präzisere Methoden.

Einschub für Fortgeschrittene: Wiedergabe einer Kreisbahn durch das Euler-Verfahren

Die Startposition der kleinen Masse $m$ sei:

${\begin{pmatrix}x\\y\end{pmatrix}}={\begin{pmatrix}R\\0\end{pmatrix}}$

Für ihre Geschwindigkeit an dieser Stelle gilt:

${\begin{pmatrix}v_{x}\\v_{y}\end{pmatrix}}={\begin{pmatrix}0\\\omega R\end{pmatrix}}$

$\omega$ bezeichnet die konstante Winkelgeschwindigkeit, welche man für $m$ auf einer Kreisbahn erwartet. Aus der Forderung, dass die Anziehungskraft $GMm/R^{2}$ gleich der Zentripetalkraft $m\omega ^{2}R$ sein muss, folgt:

$\omega ={\sqrt {\frac {GM}{R^{3}}}}$

Nach dem Euler-Verfahren bewegt sich die kleine Masse während der Zeit $\Delta t$ jedoch nicht auf einem Kreis, sondern einfach geradeaus in y-Richtung. Damit gelangt diese an folgende Position:

${\begin{pmatrix}x\\y\end{pmatrix}}={\begin{pmatrix}R\\R\omega \Delta t\end{pmatrix}}$

Der Abstand $d$ dieser geschätzten Position zur großen Masse $M$ beträgt:

$d=R{\sqrt {1+(\omega \Delta t)^{2}}}$

Für kleine Zeiten $\Delta t$ wird $\omega \Delta t<<1$, so dass man ${\sqrt {1+x}}\approx (1+x/2)$ verwenden darf. Damit vereinfacht sich der Abstand zu $d\approx R(1+(\omega \Delta t)^{2}/2)$. Der Fehler $\Delta R$ beläuft sich also auf $R(\omega \Delta t)^{2}/2$ und damit der relative Fehler $\Delta R/R$ auf:

${\frac {\Delta R}{R}}={\frac {(\omega \Delta t)^{2}}{2}}={\frac {GM(\Delta t)^{2}}{2R^{3}}}$

Mit der Sonnenmasse $M$ von 1.989 10³⁰ kg, dem Radius der Erdbahn $R$ von 1.496 10¹¹ m sowie der Schrittweite $\Delta t$ von 1 Tag bzw. 86400 s erhält man die schon skizzierten Zahlenwerte für die Bahn der Erde um die Sonne.

Der relative Fehler nimmt die gerade erwähnte besonders prägnante Form an, wenn man das 3. Keplersche Gesetz $T^{2}=4\pi ^{2}R^{3}/GM$ heranzieht. Einsetzen liefert unmittelbar:

${\frac {\Delta R}{R}}=2\pi ^{2}\left({\frac {\Delta t}{T}}\right)^{2}$

Dass der relative Fehler mit jedem Zeitschritt tatsächlich immer weiter zunimmt, wird klar, wenn man auch die Geschwindigkeit betrachtet. Gemäß der Näherung wird die kleine Masse während des gesamten Zeitraums $\Delta t$ nur in x-Richtung beschleunigt (nämlich mit $GM/R^{2}$), jedoch nicht in y-Richtung. An der neuen Position stellt sich so folgende Geschwindigkeit ein:

${\begin{pmatrix}v_{x}\\v_{y}\end{pmatrix}}={\begin{pmatrix}-{\frac {GM}{R^{2}}}\Delta t\\\omega R\end{pmatrix}}={\begin{pmatrix}-\omega R\omega \Delta t\\\omega R\end{pmatrix}}$

Der Betrag der Geschwindigkeit $v$ liegt damit bei:

$v=\omega R{\sqrt {1+(\omega \Delta t)^{2}}}$

Er weist also exakt den gleichen relativen Fehler auf wie der Abstand der beiden Massen. Man sieht zudem, dass die genäherten Orts- und Geschwindigkeitsvektoren ebenso senkrecht aufeinander stehen wie die tatsächlichen. Jede neue Position kann als Element einer Kreisbahn betrachtet werden, die sich mit jedem Zeitschritt weitet. Die kleine Masse bewegt sich im Rahmen der Eulerschen Näherung somit auf einer Spiralbahn nach außen.

C-Code: Mehrkörpersimulation mit dem Euler-Verfahren

Die Implementierung einer Mehrkörpersimulation auf Grundlage des Euler-Verfahrens mit festen Zeitschritten ist vergleichsweise einfach. Um von einem Zeitpunkt $t$ zum nächsten $t+\Delta t$ zu gelangen, berechnet man zuerst anhand der Positionen zu Beginn des Zeitschritts die entsprechenden Beschleunigungen. Dann werden die neuen Positionen mittels der anfänglichen Geschwindigkeiten bestimmt und zuletzt die neuen Geschwindigkeiten auf Grundlage der zu Beginn gültigen Beschleunigungen.

Die bislang definierten Variablen und Prozeduren bleiben unverändert gültig. Hinzukommen als Double die Gesamtzeit T, über welche sich die Simulation erstrecken soll, der aktuelle Zeitpunkt t zu Anfang jeweils eines Schritts und die Schrittweite dt.

/* Globale Variablen */

unsigned int N;
double *m;
double **r;
double **v;
double **a;

unsigned int i,k;

double t,dt;
double T;

for (t = 0;t < T;t += dt)
{

/* Beschleunigungen zum Zeitpunkt t */

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

/* Positionen und Geschwindigkeiten zum Zeitpunkt t + dt */

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