Das Mehrkörperproblem in der Astronomie/ Allgemeine Lösungsmethoden/ Zwischenschritt-Verfahren: Leapfrog und Runge-Kutta

Ein entscheidender Schwachpunkt des Euler-Verfahrens besteht darin, dass man die Momentanbeschleunigungen und -geschwindigkeiten zu einem Zeitpunkt $t$ benutzt, um den Zustand des Systems zu einer Zeit $t+\mathrm{\Delta }t$ zu bestimmen. Es ist naheliegend, dass man den tatsächlichen Verlauf der Geschwindigkeiten und Beschleunigungen während der Zeitspanne $\mathrm{\Delta }t$ genauer wiedergeben kann, wenn man zusätzlich dabei auftretende Zwischenwerte berücksichtigt.

Im Folgenden werden zwei Methoden zur Konstruktion solcher Zwischenwerte vorgestellt. Das Leapfrog-Verfahren kommt mit einem einzigen Zwischenschritt aus, vermag die Genauigkeit einer Mehrkörpersimulation im Vergleich zur Euler-Methode dennoch um etliche Größenordnungen zu steigern. Das Runge-Kutta-Verfahren erfordert mehrere Zwischenstufen, liefert dafür aber noch präzisere Ergebnisse.

Leapfrog-Verfahren

Herangehensweise

Der einfachste Weg hin zu Zwischenwerten besteht darin, die Momentangeschwindigkeiten zur Zeit $t$ nicht für einen ganzen, sondern nur einen halben Zeitschritt gelten zu lassen. Damit gewinnt man Positionen zum Zeitpunkt $t+\mathrm{\Delta }t/2$. Aus solchen Zwischenpositionen folgen unmittelbar Zwischenwerte für die Momentanbeschleunigungen. Diese lassen genauer auf die Momentangeschwindigkeiten zum Zeitpunkt $t+\mathrm{\Delta }t$ schließen als die im Rahmen des Euler-Verfahrens herangezogenen Werte zu Beginn eines Zeitschritts. Die neuen Momentangeschwindigkeiten gestatten es schließlich, aus den Zwischenpositionen zum Zeitpunkt $t+\mathrm{\Delta }t/2$ die für $t+\mathrm{\Delta }t$ gültigen Endpositionen zu berechnen. Das Schema sieht also jetzt folgendermaßen aus:

Schrittweise Berechnung der Positionen und Geschwindigkeiten eines Massenpunktes in einem Kraftfeld nach dem Leapfrog-Verfahren. Im Gegensatz zum Euler-Verfahren erfolgt diese in zwei Etappen

Um von einer alten zu einer neuen Position zu gelangen, arbeitet man nun mit dem Mittel der alten und neuen Momentangeschwindigkeiten. Dieses stellt eine viel bessere Näherung der für die Zeitspanne zwischen $t$ und $t+\mathrm{\Delta }t$ geltenden Durchschnittsgeschwindigkeit dar als die von dem Euler-Verfahren ausschließlich benutzte alte Momentangeschwindigkeit zu Beginn des Zeitintervalls.

Die Verwendung des Mittels der Momentangeschwindigkeiten zu Anfang und Ende eines Zeitschritts bedeutet, dass man die wahre Fläche unter der Geschwindigkeitskurve nicht mehr durch Rechtecke, sondern Trapeze ersetzt. Es leuchtet sofort ein, dass dieses Vorgehen wesentlich zuverlässiger ist.

Berechnung der zurückgelegten Strecke in Abh. von der Geschwindigkeit durch das Leapfrog-Verfahren. Der wahre Geschwindigkeitsverlauf wird im Gegensatz zum Euler-Verfahren nicht durch eine Stufen-, sondern eine stückweise lineare Funktion angenähert

Als kritischer Aspekt bleibt bestehen, dass man die Momentangeschwindigkeit zum Zeitpunkt $t+\mathrm{\Delta }t$ bereits kennen muss, bevor man die entsprechende Position berechnet hat. Dies ist nur nach dem Schema des Euler-Verfahrens möglich, d.h. erneut wird die Fläche unter der Beschleunigungskurve durch Rechtecke angenähert. Die Höhe eines solchen wird nun aber jeweils durch die Beschleunigung in der Mitte eines Zeitintervalls festgelegt und nicht durch den Wert zu Beginn. Somit wird auch für die Beschleunigung die Fläche unter der Kurve genauer wiedergegeben.

Berechnung der Geschwindigkeitsänderung in Abh. von der Beschleunigung durch das Leapfrog-Verfahren. Der wahre Beschleunigungsverlauf wird durch eine Stufenfunktion angenähert. Im Gegensatz zum Euler-Verfahren werden die Beschleunigungen jedoch nicht zu Beginn der Zeitintervalle genommen, sondern in deren Mitte.

Genauigkeit

Wieder sei die Güte der Lösung am Beispiel der Kreisbahn demonstriert. Nach dem ersten Halbschritt scheint keine signifikante Verbesserung eingetreten zu sein, da dieser genau dem Schema des Euler-Verfahrens folgt. Die Beschleunigung an der Zwischenposition wirkt jedoch derart auf die kleine Masse $m$, dass mit der für den zweiten Halbschritt geltenden neuen Geschwindigkeit der Fehler des ersten Halbschritts fast vollständig kompensiert wird.

Behandlung einer Kreisbahn durch das Leapfrog-Verfahren

Eine genaue Analyse der Erdbahn liefert für einen Zeitschritt von 1 Tag (zwei Halbschritten von je 1/2 Tag) einen erstaunlich kleinen Fehler von nur etwa 5 Milliardstel des wahren Bahnradius. Nimmt man wieder den ungünstigsten Fall an, dass die Fehler der einzelnen Schritte sich aufaddieren (die später skizzierten Tests zeigen, dass dies für die Leapfrog-Methode nicht zutrifft), erhält man auf 1 Jahr hochgerechnet einen relativen Fehler von etwa 2 Millionstel! Die auf dem Leapfrog-Verfahren beruhende Lösung ist um mehrere Größenordnungen genauer als diejenige nach Euler. Ein Fehler von 2 Millionstel des Radius pro Jahr ist allerdings immer noch viel zu hoch, um beispielsweise die Stabilität der Erdbahn (welche ja nicht isoliert ist, sondern von den übrigen Planeten des Sonnensystems gestört wird) auf geologischer Zeitskala zu prüfen.

Das im Vergleich zur Euler-Methode gute Resultat ist kein Zufall, denn für den prozentualen Fehler pro Schritt gilt:

Dieser ist jetzt der 4. Potenz des Zeitschritts $\mathrm{\Delta }t$ proportional, d.h. eine Halbierung der Zeitintervalle erbringt pro Schritt das 16-fache an Genauigkeit. Verringert man allgemein $\mathrm{\Delta }t$ um einen Faktor $n$, steigert man die Präzision eines Einzelschritts um $n^4$ bzw. der gesamten Simulation um $n^3$. Kleinere Zeitschritte erhöhen nun sehr effizient die Genauigkeit der Lösung. Umgekehrt muss man aber sehr darauf achten, zwecks kürzerer Rechenzeit die Schrittweite nicht allzu sehr zu vergröbern, da dann die Qualität des Verfahrens sehr rasch abnimmt.

Die Effizienz des Leapfrog-Verfahrens darf keineswegs zu dem Schluss verleiten, man könne mit genügend kleiner Schrittweite eine im Prinzip beliebige Genauigkeit erzielen. Die bei Computersimulationen verwendeten Gleitkommazahlen weisen selbstverständlich eine endlich Anzahl von Nachkommastellen auf, so dass jeder Zeitschritt zwangsläufig mit Rundungsfehlern behaftet ist. Kleinere Schritte bedeuten mehr Einzelberechnungen und damit Rundungsfehler, wodurch die höhere nominelle Genauigkeit der Näherung zunehmend unterlaufen wird. Auch hier handelt es sich um eine prinzipielle, jedem nominell noch so guten Verfahren anhaftende Einschränkung.

Dem beachtlichen Abschneiden der Leapfrog-Methode bezüglich der zeitlichen Schrittweite steht leider eine enorme Instabilität hinsichtlich des Abstandes der beiden Massen gegenüber. Es schält sich nämlich folgende Beziehung heraus:

Der prozentuale Fehler steigt sage und schreibe mit der 6. Potenz des Kehrwertes von $R$ an – bei der Hälfte des ursprünglichen Abstands schon auf das 64-fache, bei einem Drittel desselben gar bereits auf das 729-fache! Man könnte daher vermuten, dass sehr enge Vorübergänge zweier Punktmassen mit dem Euler-Verfahren zuverlässiger behandelt würden. In der Tat gibt es einen kritischen Abstand, bei dem letzteres einen kleineren prozentualen Fehler pro Zeitschritt liefert. Allerdings ist der Fehler selbst eines einzigen Schrittes an dieser Stelle schon so hoch, dass keine der beiden Näherungen mehr auch nur ansatzweise tauglich ist. Beide Verfahren brechen zusammen, lange bevor sich zwei Massenpunkte bis auf den relevanten Abstand genähert haben. Damit aber ist dieser in der Praxis völlig belanglos.

Um ein Mehrkörpersystem auch unter den Bedingungen von Beinahezusammenstößen korrekt behandeln zu können, sind spezielle Techniken notwendig, denen das ganze nächste Kapitel gewidmet ist. Dazu gehört unter anderem eine dynamische zeitliche Schrittweite. Dort wo die Körper eng zusammenstehen, werden die Berechnungen mit kleineren $\mathrm{\Delta }t$ durchgeführt als in Gebieten geringer Dichte.

Analog zum Euler-Verfahren kann man auch für die Leapfrog-Methode das 3. eplersche Gesetz auf den relativen Fehler anwenden. Abermals erweist sich das Verhältnis zeitliche Schrittweite / Umlaufdauer als entscheidend, denn man erhält:

Vorlage:Kasten

Vorlage:Kasten

Aufgrund des vergleichsweise geringen Rechenaufwandes wird das Leapfrog-Verfahren häufig in der Praxis eingesetzt. Mit der hier geschilderten Darstellung als Ziwschenschritt-Methode ist es aber nicht möglich, dynamische zeitliche Schrittweiten anzuwenden. Jedoch lassen sich die Beziehungen zwischen Positionen, Geschwindigeiten und Beschleunigungen zu diesem Zweck relativ leicht umformulieren, wie das nächste Unterkapitel zeigt.

Klassisches Runge-Kutta-Verfahren

Heun-Verfahren als Vorstufe

Von allen Methoden, mehrere Zwischenschritte von alten nach neuen Positionen zu bilden, ist das Runge-Kutta-Verfahren bei weitem am gebräuchlichsten. Es leitet sich jedoch nicht von dem Leapfrog-Algorithmus ab, sondern dem folgenden, auf den Mathematiker Heun zurückgehenden Schema. Auf den ersten Blick ist dieses dem Leapfrog-Verfahren unterlegen, weil für die Berechnung der Geschwindigkeit am Ende des Zeitintervalls $\mathrm{\Delta }t$ die Beschleunigung zu dessen Beginn verwendet wird und nicht der Zwischenwert in dessen Mitte.

Schrittweise Berechnung der Positionen und Geschwindigkeiten eines Massenpunktes in einem Kraftfeld nach dem Heun-Verfahren. Im Gegensatz zum Leapfrog-Verfahren wird die (neue) Geschwindigkeit am Ende eines Schritts aus der (alten) Beschleunigung zu Beginn eines Zeitintervalls abgeleitet und nicht aus der nach einem halben Schritt herrschenden Zwischenbeschleunigung

Als Folge dessen zeigt der relative Fehler eines Zeitschritts lediglich eine Proportionalität $\mathrm{\Delta }R/R\sim (\mathrm{\Delta }t{)}^3$. Die Genauigkeit eines solchen Schritts steigt so nur mit $n^3$ und diejenige einer kompletten Simulation nur mit $n^2$, wenn man $\mathrm{\Delta }t$ um den Faktor $n$ verringert. Das Verfahren lässt sich jedoch, wie die Mathematiker Runge und Kutta erkannten, stringenter zu mehreren Zwischenschritten ausbauen als die Leapfrog-Methode.

Erweiterung zum klassischen Verfahren

Um eine solche Verbesserung zu erreichen, startet man zunächst wie bei dem Leapfrog-Verfahren mit der Geschwindigkeit ${\overrightarrow{v}}_1$ und Beschleunigung ${\overrightarrow{a}}_1$ an der Position ${\overrightarrow{r}}_1$ zu Beginn eines Zeitintervalls und geht damit einen halben Schritt zur Position ${\overrightarrow{r}}_2$, wo neue Werte für Geschwindigkeit ${\overrightarrow{v}}_2$ und Beschleunigung ${\overrightarrow{a}}_2$ gelten. Mit diesen startet man nun abermals von der Ausgangsposition ${\overrightarrow{r}}_1$ aus und gelangt so wiederum mit einem Halbschritt an die Position ${\overrightarrow{r}}_3$, wo die Geschwindigkeit ${\overrightarrow{v}}_3$ und Beschleunigung ${\overrightarrow{a}}_3$ herrschen. Man durchläuft das Zeitintervall $[t,t+\mathrm{\Delta }t/2]$ einmal unter den Bedingungen zu seinem Beginn (Werte 1) und im Gegensatz zur Leapfrog-Methode ein zweites Mal unter den ungefähren Verhältnissen zu dessen Ende (Werte 2). Dementsprechend stellt das Resultat (Werte 3) im Vergleich zum Leapfrog-Verfahren einen genaueren Zwischenwert für den Zeitpunkt $t+\mathrm{\Delta }t/2$ bereit.

Hiermit gelangt man wiederum von der Startposition ausgehend nach einen ganzen Zeitschritt zur Position ${\overrightarrow{r}}_4$, welche eine Geschwindigkeit ${\overrightarrow{v}}_4$ und Beschleunigung ${\overrightarrow{a}}_4$ liefert. Diese Werte 4 stellen jedoch noch nicht das Endergebnis der Prozedur dar. Vielmehr bildet man aus den Anfangswerten 1, den Zwischenwerten 2 und 3 sowie den ungefähren Endwerten 4 Mittelwerte für Geschwindigkeit und Beschleunigung während des ganzen Zeitintervalls, wobei die Zwischenwerte doppelt gewichtet werden. Diese Mittelwerte führen nun tatsächlich von der Ausgangs- zur Endposition.

Schrittweise Berechnung der Positionen und Geschwindigkeiten eines Massenpunktes in einem Kraftfeld nach dem klassischen Runge-Kutta-Verfahren

Das Runge-Kutta-Verfahren ist zu aufwendig, um dessen Güte anhand einer expliziten algebraischen Rechnung für die Kreisbahn darzulegen. Man kann jedoch allgemein zeigen, dass es die Leapfrog-Methode noch um eine Potenz übertrifft:

Der prozentuale Fehler eines Zeitschritts ist der 5. Potenz, derjenige einer gesamten Simulation der 4. Potenz von $\mathrm{\Delta }t$ proportional. Es gilt jedoch zu beachten, dass dies mit wesentlich mehr Einzelberechnungen erkauft ist und damit die Problematik der Rundungsfehler umso stärker ins Gewicht fällt. Wieder gibt eine optimale Schrittweite, unterhalb derer aufgrund der endlichen Genauigkeit der Gleitkommazahlen die Qualität des Verfahrens sich nicht weiter verbessert, im Gegenteil bei zugleich längerer Rechenzeit sogar abnimmt.

Vorlage:Kasten

Vorlage:Kasten