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

Konstruktion

Bearbeiten

Bereits der Mathematiker Euler entwickelte das Konzept, ein zu komplexes System in diskreten Zeitschritten   zu beobachten, anstatt nach einer kontinuierlichen, für beliebige Zeiten   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   näherungsweise als konstant betrachten. Man geht also davon aus, dass die zu einem bestimmten Zeitpunkt   gegebenen Momentanbeschleunigungen und -geschwindigkeiten den mittleren Werten für den gesamten Zeitrahmen   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   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

Bearbeiten

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


 
Behandlung einer Kreisbahn durch das Euler-Verfahren


Gemäß der Zeichnung soll sich   zu Beginn genau rechts von   aufhalten. Die Geschwindigkeit der kleinen Masse weist dann exakt nach oben (oder unten). Es ist einleuchtend, dass   sich schon nach dem ersten Schritt außerhalb der tatsächlichen Kreisbahn befindet, in einer zu großen Entfernung  . Der Fehler   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:

 

Der Fehler des Verfahrens hängt nur quadratisch von   ab, mit z.B. halb so großen Zeitintervallen lässt sich nur die vierfache Genauigkeit pro Schritt erzielen. Verringert man   allgemein um einen Faktor  , so weist jeder Einzelschritt eine um   verbesserte Genauigkeit auf. Jedoch werden dann auch   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  -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   selbst, sondern auch dem Abstand der beiden Massen ab. Hierbei gilt:

 

Er verhält sich dem Kubus von   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   folgt für den relativen Fehler:

 

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   sei:

 

Für ihre Geschwindigkeit an dieser Stelle gilt:

 

  bezeichnet die konstante Winkelgeschwindigkeit, welche man für   auf einer Kreisbahn erwartet. Aus der Forderung, dass die Anziehungskraft   gleich der Zentripetalkraft   sein muss, folgt:

 

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

 

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

 

Für kleine Zeiten   wird  , so dass man   verwenden darf. Damit vereinfacht sich der Abstand zu  . Der Fehler   beläuft sich also auf   und damit der relative Fehler   auf:

 

Mit der Sonnenmasse   von 1.989 1030 kg, dem Radius der Erdbahn   von 1.496 1011 m sowie der Schrittweite   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   heranzieht. Einsetzen liefert unmittelbar:

 

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   nur in x-Richtung beschleunigt (nämlich mit  ), jedoch nicht in y-Richtung. An der neuen Position stellt sich so folgende Geschwindigkeit ein:

 

Der Betrag der Geschwindigkeit   liegt damit bei:

 

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   zum nächsten   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;
    }
  }
}