Das Mehrkörperproblem in der Astronomie/ Druckversion

Aus Wikibooks
Zur Navigation springen Zur Suche springen

Mehrkörpersysteme in der Astronomie[Bearbeiten]

Inhaltsverzeichnis

Die kosmische Hierarchie[Bearbeiten]

Sirius in künstlerischer Darstellung

Wie bereits seit langem bekannt, sind die Himmelskörper nicht regellos im Universum verteilt, sondern innerhalb von Systemen höchst unterschiedlicher Größe und Struktur angeordnet. Monde umkreisen Planeten, diese wiederum bewegen sich um Zentralgestirne. Die Sterne selbst sind ebenfalls keine Einzelgänger. Sehr häufig umkreisen sich zwei Sterne gegenseitig, bilden also ein sogenanntes Doppelsternsystem. Ein Beispiel hierfür ist Sirius, der hellste Stern am Nachthimmel, welcher einen sehr lichtschwachen Begleiter in Gestalt eines weißen Zwerges besitzt. Doch auch Mehrfachsysteme kommen vor. So besteht Castor, das hellste Objekt im Sternbild Zwillinge, in Wahrheit aus 6 Sternen: drei äußerst enge Paare, die um einen gemeinsamen Schwerpunkt laufen, wie Struve (1952) [1] erkannte.

Die Plejaden, ein offener Sternhaufen

Selbst solche Objekte stehen in der Hierarchie der Sternsysteme freilich noch fast ganz unten. Mehr als 1000 sogenannte offene Sternhaufen sind bekannt, ziemlich dichte Ansammlungen von typischerweise mehreren 100 Sternen. Ein typischer Vertreter sind die Plejaden, das Siebengestirn, für welches Jones (1970) [2] 340 zugehörige Objekte bis zur 15.Größenklasse (entsprechend Hauptreihensternen des Spektraltyps K) auflistet. Aktuelle Beobachtungen von Bouy und anderen (2015) [3], welche auch sehr lichtschwache Rote und Braune Zwerge berücksichtigen, lassen sogar 2100 Mitglieder erkennen. Besonders spektakulär erscheinen dem Beobachter die Kugelhaufen, wo hunderttausende von Sternen auf so engem Raum angeordnet sind, dass diese selbst mit leistungsstarken Instrumenten häufig nur am Randbereich des Systems noch als Individuen erkannt werden können. 47 Tucanae im Sternbild Tukan ist eines der schönsten Beispiele.

47 Tucanae, ein Kugelsternhaufen

Einzel-, Doppel- und Mehrfachsterne sowie Kugel- und offene Sternhaufen finden sich wiederum zu Galaxien unterschiedlicher und oft sehr komplexer Gestalt zusammen. Spiralgalaxien wie die Andromedagalaxie bestehen aus drei Komponenten: einer im Vergleich zu ihrem Durchmesser sehr dünnen Scheibe (welche die Spiralarme und offenen Sternhaufen enthält), einer linsenförmigen Verdickung im Zentrum (auch in der deutschen Literatur oft mit dem englischen Begriff Bulge bezeichnet) sowie einem annähernd kugelsymmetrischen Halo (wo sich die Kugelhaufen befinden). Balkenspiralgalaxien ähneln den einfachen Spiralgalaxien, doch unterscheiden sie sich von diesen dadurch, dass die Spiralarme an einer balkenförmigen Struktur ansetzen. Unsere eigene Milchstraße, lange Zeit als gewöhnliche Spiralgalaxie betrachtet, könnte laut Benjamin (2008) [4] ebenfalls ein Balkensystem darstellen. Sowohl gewöhnliche als auch Balkenspiralen können mehrere hundert Milliarden Sterne enthalten. Darüber hinaus enthalten sie große Mengen an interstellarer Materie in Form von Gas und Staub. Den Löwenanteil an ihrer Masse aber scheint dunkle Materie zu stellen, die sich allein durch ihrer Gravitationswirkung verrät.

Der Andromedanebel, eine Spiralgalaxie

Zum Ensemble der großen Sternsysteme gehören auch elliptische Galaxien. Wie der Name bereits andeutet, weisen sie die Gestalt eines Ellipsoids auf, ohne dabei wie Spiralgalaxien auffällige Strukturen zu zeigen. Elliptische Systeme sind oft noch massereicher als Spiralgalaxien, wobei wiederum dunkle Materie zu dominieren scheint. Der Anteil interstellarer Materie ist oft sehr gering.

Mit der Verfügbarkeit immer größerer und damit lichtempfindlicherer Teleskope sind in den letzten Jahrzehnten neben den großen verstärkt die Zwerggalaxien ins Blickfeld der Forschung gerückt, wobei als bekannteste Vertreter die Magellanschen Wolken zu nennen sind. Hinsichtlich ihrer Größe und Sternanzahl stehen sie zwischen den Kugelhaufen und den großen Galaxien. Meist sind sie von elliptischer oder unregelmäßiger Gestalt, wohingegen Zwergspiralen selten sind. Zwergsysteme weisen im Vergleich zu großen Galaxien meist einen geringeren Anteil interstellarer Materie auf, werden dafür aber um so mehr von dunkler Materie beherrscht.

Trotz ihrer gewaltigen Dimensionen stellen selbst die Galaxien noch nicht das Ende der Stufenleiter dar, sie ordnen sich vielmehr in eine noch höhere Hierarchie ein. Große Galaxien wie die Milchstraße halten sich einen Hofstaat von mehreren Zwerggalaxien (so sind die eben erwähnten Magellanschen Wolken gravitativ an die Milchstraße gebunden). Mehrere solche von kleinen Satelliten umgebene sowie einzelne Galaxien schließen sich weiter zu Galaxienhaufen zusammen, zum Beispiel bilden die von der Milchstraße und der Andromedagalaxie definierten Satellitenysteme zusammen mit einigen weiteren Galaxien die Lokale Gruppe. Galaxienhaufen können bis zu einige 1000 Galaxien enthalten.

NGC 1300, eine Balkenspiralgalaxie
ESO 325-G004, eine elliptische Galaxie, eingebettet in den Galaxienhaufen Abell S740

Die Galaxienhaufen schließlich bilden die größten bekannten Systeme des Universums, die sogenannten Superhaufen. Der uns am nächsten stehende ist der Virgo-Superhaufen, zu welchem auch die Lokale Gruppe gerechnet wird. Selbst auf dieser höchsten Ebene liegt keine irreguläre Verteilung der Materie vor. Haufen und Superhaufen scheinen vielmehr sich entlang von Filamenten anzuordnen, zwischen denen sich riesige, kaum Galaxien enthaltende Leerräume befinden.

Die Sagittarius-Zwerggalaxie

Um die Entstehung und Dynamik all dieser Strukturen zu verstehen, hat sich die Simulation von Mehrkörpersystemen als sehr nützliches Werkzeug erwiesen. Im Prinzip stellt man hierbei das zu untersuchende System als ein Ensemble von Massepunkten dar und berechnet mit bestmöglicher Genauigkeit, wie sich diese infolge ihrer wechselseitigen Anziehungskraft bewegen. Welches reale Objekt dabei einem Massepunkt entspricht, hängt von der betrachteten Hierarchieebene ab. Interessiert man sich für die Vorgänge innerhalb eines Sternhaufens, so repräsentiert jeder Massepunkt einen Stern. Gegebenenfalls um einen Stern laufende Planeten sind dagegen aufgrund ihrer im Vergleich zum Zentralgestirn geringen Masse in der Regel irrelevant. Zudem ist die Ausdehnung ihrer Bahnen viel kleiner als der typische Abstand selbst zum nächsten Stern. Von der Ebene eines Sternhaufens betrachtet, darf ein Stern mitsamt seines Planetensystems somit als ein einziger Körper behandelt werden. Will man großräumige Strukturen des Kosmos untersuchen, so entspricht einem Massepunkt eine Galaxie, wenn nicht gar ein ganzer Galaxienhaufen. Die exakte Verteilung der Sterne innerhalb der Galaxien ist für die relative Bewegung der Galaxien und erst recht der Galaxienhaufen untereinander zumeist bedeutungslos.

R136, ein extrem dichter Sternhaufen

Oft reicht es zur Erklärung bestimmter Strukturen jedoch nicht aus, lediglich mit einem Ensemble von Massepunkten als Modell zu arbeiten, die unter wechselseitigem Einfluß der Schwerkraft stehen. Auch anderer Phänomene, wie etwa hydrodynamische Vorgänge oder Strahlung, spielen häufig eine wichtige Rolle. Man nimmt beispielsweise an, dass ein offener Sternhaufen aus einem durch Masseanziehung bedingten Kollaps einer großen, kalten Wolke aus molekularem Wasserstoff hervorgeht. Im Verlauf eines solchen Zusammenbruchs zerfällt die Wolke zu kleineren Klumpen, aus denen schließlich einzelne, z.T. sehr massereiche Sterne entstehen. Solche geben jedoch eine intensive ultraviolette Strahlung ab und ionisieren dadurch ihre Umgebung. Zusätzlich verdrängen sie, durch Strahlungsdruck und Sternwind, das sie umgebende Gas. Nach wenigen Millionen Jahren zerplatzen die ersten sehr massereichen Sterne durch Supernova-Explosionen. Durch die dabei entstehenden Schockwellen wird viel Gas aus dem noch jungen Sternhaufen herausgeschleudert. Die anfänglich sehr aktive Sternentstehung kommt durch diese Verluste an Rohmaterial mehr und mehr zum Erliegen.

Die Simulation von Mehrkörpersystemen wird auch, aber nicht nur, durch den Einfluss solch zusätzlicher Effekte häufig sehr komplex. In bestimmten Fällen ist die Ersetzung realer Himmelskörper durch Massepunkte auch schlicht nicht zulässig. Ein frappantes Beispiel ist der Sternhaufen R136 im Tarantelnebel, welcher die Zentralregion der Großen Magellanschen Wolke darstellt. Dieser Haufen enthält Sterne so hoher Leuchtkraft, dass man von Einzelmassen von bis zu 300 Sonnenmassen ausgehen muss, was nach der aktuellen Theorie der Sternentstehung eigentlich unmöglich ist. Jedoch ist die Sterndichte in diesem Haufen so hoch, dass gemäß Fuji und Portegies Zwart (2013) [5] tatsächlich einzelne Mitglieder miteinander kollidieren und so zu Objekten weit höherer Masse verschmelzen können.


Die Simulation unterschiedlicher Typen von Mehrkörpersystemen im Kosmos soll in diesem Buch detailliert behandelt werden. Dabei wird ausschließlich das Modell unter der Schwerkraft sich bewegender Massepunkte zum Einsatz kommen, da es trotz aller hier skizzierten Vereinfachungen viele Beobachtungen korrekt wiedergeben kann. Zunächst werden die dafür unentbehrlichen physikalischen Grundlagen erläutert, anschließend allgemeine Verfahren für das Lösen von Bewegungsgleichungen vorgestellt. Ein großer Nachteil des auf Massepunkten beruhenden Modells besteht darin, dass sich solche Punkte beliebig nahe kommen könnten und dann im Prinzip unbegrenzt große Anziehungskräfte aufeinander ausüben würden. Wie dieses Problem umgangen werden kann, wird in einem eigenen Kapitel dargelegt. Eine weitere Schwierigkeit für astronomische Mehrkörpersimulationen erwächst aus der Tatsache, dass wegen der großen Reichweite der Gravitation jedes Mitglied eines solchen Systems auf alle übrigen Angehörigen einwirkt. Bei einer auf einem Standard-Verfahren beruhenden Lösung der Bewegungsgleichungen wächst der Rechenaufwand daher mit dem Quadrat der Anzahl der Massepunkte. In einem weiteren Kapitel werden deshalb Algorithmen vorgestellt, welche eine wesentlich effizientere Behandlung von Mehrkörpersystemen gestatten. Im letzten Kapitel wird erläutert, wie sich diese optimierten Lösungsverfahren auf tatsächlich beobachtete Ensembles vom Planetensystem bis herauf zum Universum als Ganzes anwenden lassen.


Grundlagen[Bearbeiten]

Im Folgenden soll ungeachtet aller hier angedeuteten Komplikationen die Schwerkraft als einzige in einem Mehrkörpersystem wirkende Kraft berücksichtigt werden. Zudem sollen - abgesehen von einigen Diskussionen in speziellen Abschnitten - alle Körper als ideale Massenpunkte betrachtet werden. Damit wird indirekt angenommen, dass der Abstand zweier Sterne stets sehr viel größer ist als ihr Durchmesser. Bevor mögliche Lösungsansätze im Detail besprochen werden, sollen zuvor einige für deren Verständnis unentbehrliche Grundbegriffe skizziert werden.

Vektoren und Skalare[Bearbeiten]

Bei vielen der in einem Sternsystem interessierenden Größen, z.B. den Geschwindigkeiten der einzelnen Mitglieder, kommt es nicht nur auf den Betrag an, sondern auch auf die Richtung, in welche die Massenpunkte sich bewegen - wie bei einer Autofahrt, wo man sich nicht nur Gedanken darüber macht, wie viele Kilometer man pro Stunde zurücklegen will, sondern auch darüber, welchen Weg man einschlagen muss, um an das Ziel zu gelangen. Solche Größen werden als Vektoren bezeichnet. Dagegen kann man einer Größe wie der Masse keine Richtung zuordnen. Eine solche Größe nennt man Skalar.

Vektoren werden in der Regel durch einen Pfeil über dem Größensymbol gekennzeichnet. Diskutiert man z.B. die Geschwindigkeit einschließlich ihrer Richtung, so schreibt man .

Komponentendarstellung und Betrag eines Vektors[Bearbeiten]

Die Tatsache, dass Vektoren eine Richtungsangabe beinhalten, führt zwangsläufig dazu, dass solche aus mehreren Komponenten bestehen. Will man z.B. von Nürnberg nach Frankfurt gelangen, muss man sich in nordwestliche Richtung bewegen. Die Angabe "Nordwesten" aber schließt bereits zwei Komponenten ein, nämlich "Westen" und "Norden". Man könnte diese durch kartesische Koordinaten wiedergeben, in dem man z.B. die West-Ost-Richtung als x-Achse und die Süd-Nord-Richtung als y-Achse darstellt. Der von Nürnberg nach Frankfurt zielende Vektor hätte dann eine in negative x-Richtung zeigende Komponente und eine etwa gleich große in positive y-Richtung weisende.

Im Alltag kommt man meist mit zwei Richtungskomponenten aus. In einem Sternsystem dagegen mit voller räumlicher Bewegung der einzelnen Massenpunkte sind drei Komponennten erforderlich, zu der x- und y-Achse tritt die z-Achse als weitere kartesische Komponente hinzu. Will man die einzelnen Komponenten eines Vektors auflisten, so verwendet man oft folgende Schreibweise:

, und sind dabei die Geschwindigkeiten in Richtung der drei kartesischen Achsen.

Nicht nur die Geschwindigkeit, auch der Ort eines Körpers ist ein Vektor, denn für eine Ortsangabe werden gleichfalls mehrere Komponennten benötigt. Auf der Erde genügen dazu der Längen- und Breitengrad. In einem Sternsystem müssen abermals drei Komponenten angegeben werden, wobei man als Bezugspunkt zumeist den Schwerpunkt wählt (dieser ist, wie im letzten Abschnitt dieses Unterkapitels dargelegt wird, der auf die Verteilung der einzelnen Massen bezogene Mittelpunkt des Systems). Man schreibt dann z.B.

Der Ortsvektor bezeichnet demgemäß die Richtung, in welcher man vom Ursprung aus zu dem gewünschten Objekt schauen muss.

Den Betrag eines Vektors kann man aus seinen Komponenten mittels des Satzes des Pythagoras berechnen. Um kenntlich zu machen, dass man den Betrag meint ohne Richtungsangaben, lässt man entweder den Pfeil über dem Größensymbol weg (schreibt also z.B. einfach ) oder umgibt die als Vektor markierte Größe mit Betragsstrichen (z.B. ). Auf den Geschwindigkeitsvektor angewandt liefert der Satz des Pythagoras den Absolutwert der Gesamtgeschwindigkeit:

Im Falle des Ortsvektors gewinnt man den Abstand eines Massenpunktes vom Schwerpunkt des Systems.

Vektoren stellt man sich anschaulich oft als Pfeile vor. Die Komponenten geben die Richtung an, nach welcher der Pfeil orientiert ist, der Betrag dessen Länge.

Im Verlauf einer Mehrkörpersimulation müssen Betragsberechnungen sehr häufig durchgeführt werden, so dass es angebracht ist, dafür eine eigene Prozedur zu definieren. Übergeben wird dieser der vektor, dessen Betrag zu bestimmen ist. vektor ist ein Array, das drei Double stellvertretend für die drei kartesischen Komponenten enthält. Der Rückgabewert ist ebenfalls vom Typ Double.


double betrag (double *vektor)
{
  return (sqrt (pow (vektor[0],2) + pow (vektor[1],2) + pow (vektor[2],2)));
}


Die wichtigste Anwendung dieser Prozedur besteht in der Berechnung des Abstands zweier Körper. Zunächst bestimmt man dafür den Differenzvektor der Ortsvektoren der beiden Objekte und anschließend dessen Betrag.

Nicht nur für obige, sondern auch für alle nachfolgend in diesem Buch besprochenen Prozeduren werden für reelle Zahlen (Gleitkommazahlen) Variablen des Datentyps Double verwendet. Zwar werden mit diesem Datentyp für jede Gleitkommazahl 8 Byte Speicherplatz benötigt, doch werden mit einer Genauigkeit von 15-16 Dezimalstellen Rundungsfehler minimiert, was wegen der bei einer Mehrkörpersimulation auftretenden großen Zahl von Rechenoperationen von großer Wichtigkeit ist. Mit dem Datentyp Float kommt man mit 4 Byte pro reelle Zahl aus, doch kann mit nur 7-8 vorgehaltenen Dezimalstellen keine ausreichende Rechengenauigkeit gewährleistet werden.

Addition und Subtraktion von Vektoren[Bearbeiten]

Auf Vektoren können die Grundrechenarten Addition und Subtraktion gleichermaßen angewandt werden wie auf Skalare. Dies geschieht, indem man die einzelnen Komponenten addiert bzw. voneinander subtrahiert. Betrachtet man zwei Vektoren und , so gilt für die Vektorsumme :

Amalog gilt für die Vektordifferenz :

Graphisch kann man die Vektoraddition durchführen, in dem man an die Spitze des Vektors den Fuß von setzt. Der Summenvektor erstreckt sich dann von dem Fuß von bis zur Spitze von .


Addition zweier Vektoren +


Im Falle der Subtraktion wird ebenfalls der Fuß von an die Spitze von gesetzt, dabei aber zugleich in seine Gegenrichtung umgeklappt. Der Differenzvektor weist wiederum vom Fuß von zur Spitze des umgeklappten Vektors .


Subtraktion zweier Vektoren -


Bei der Simulation eines Sternsystems wird die Vektoraddition vor allem zur Berechnung der auf ein Mitglied wirkenden Gesamtkraft benötigt. Jeder andere Massenpunkt wirkt auf einen beliebig herausgegriffenen Körper eine gewisse Anziehungskraft aus. Um die gesamte Gravitation zu kennen, müssen alle diese einzelnen Kraftvektoren aufsummiert werden.

Um die Anziehung zwischen zwei Massenpunkten zu bestimmen, muss deren Abstand bekannt sein. Dazu wiederum müssen wie schon angedeutet deren Ortsvektoren voneinander subtrahiert werden, woraus der sogenannte Abstandsvektor resultiert.

Multiplikation eines Vektors mit einem Skalar und Einheitsvektor[Bearbeiten]

Wie Skalare können auch Vektoren mit weiteren Skalaren multipliziert werden. Um einen Vektor mit einem Skalar zu multiplizieren, bildet man für jede seiner Komponenten das Produkt mit :

Multipliziert man einen Vektor mit einem Skalar > 1, so wird dieser einfach in der ursprüngliche Richtung gestreckt. Liegt der Skalar zwischen 0 und 1, resultiert eine Stauchung wiederum in ursprünglicher Richtung. Bei Multiplikation mit einem negativen Skalar wird der Vektor nicht nur gestreckt (falls < -1) oder gestaucht (für zwischen -1 und 0), sondern zugleich in seine Gegenrichtung umgeklappt.

Die Division mit einem Skalar ist gleichbedeutend mit der Multiplikation mit . Wird ein Vektor durch Multiplikation mit einem Skalar > 1 gestreckt, so wird er durch Division mit demselben gestaucht (die Division wirkt dann wie eine Multiplikation mit einem Skalar < 1).


Multiplikation eines Vektors mit einem Skalar


Eine besondere, aber in der Praxis oft benötige Rechenoperation ist die Division eines Vektors mit seinem eigenen Betrag. Der dadurch gewonnene Vektor hat exakt die Länge 1 und wird als Einheitsvektor bezeichnet; er zeigt in die gleiche Richtung wie der Ausgangsvektor:

Der Einheitsvektor stellt eine reine Richtungsangabe da. So bezeichnet beispielsweise der Einheitsvektor der Geschwindigkeit die Richtung, in welche sich ein Objekt gerade bewegt. Der Einheitsvektor des Abstandsvektors gibt die Verbindungslinie zweier Massenpunkte an.

Der einheitsvektor eines vektor lässt sich mit folgender Prozedur gewinnen. Zunächst liefert die Betragsprozedur in Form der Double laenge den Betrag von vektor. Danach wird jede Komponente von vektor mit laenge dividiert. einheitsvektor ist wie vektor ein Array mit 3 Double.

Der Nullvektor, dessen Komponenten und somit auch dessen Betrag gleich 0 sind, ist nicht normierbar. Um diesen Sonderfall programmiertechnisch abfangen zu können, gibt die Prozedur für den Nullvektor wiederum den Nullvektor zurück.

void einheit (double *vektor, double *einheitsvektor)
{
  double laenge;
  unsigned int k;

  laenge = betrag (vektor);
  for (k = 0;k < 3;k ++)
  {
    if (laenge > 0)
    einheitsvektor[k] = vektor[k] / laenge;
    else
    einheitsvektor[k] = 0;
  }
}


Skalarprodukt und Kreuzprodukt[Bearbeiten]

Das "Produkt" zweier Vektoren hat im Gegensatz zur Multiplikation eines Vektors mit einem Skalar mit einer gewöhnlichen Multiplikation nur wenig zu tun. Es existieren zwei verschiedene Formen, das Skalarprodukt und das Kreuzprodukt.

Das Skalarprodukt zweier Vektoren und wird so genannt, weil es nicht wie die bisher skizzierten Rechenoperationen wieder einen Vektor, sondern einen Skalar liefert. Es wird in der Regel durch den gewöhnlichen Multiplikationspunkt zwischen den beiden betrachteten Vektoren kenntlich gemacht, man schreibt also . Gegeben ist es durch die Beträge von und sowie den Cosinus des von den beiden Vektoren eingeschlossenen Winkels :

Graphisch lässt sich das Skalarprodukt folgendermaßen darstellen. Man projiziert den Vektor senkrecht auf und erhält so den Vektor . Das Skalarprodukt folgt dann durch Multiplikation der Beträge von und . Ist der Winkel spitz, zeigt in die gleiche Richtung wie und das Skalarprodukt ist positiv. Stehen die beiden zu "multiplizierenden" Vektoren senkrecht aufeinander, so verschwindet die Projektion von , d.h. das Skalarprodukt ist gleich 0. Im Falle eines stumpfen Zwischenwinkels zeigt der Projektionsvektor in die Gegenrichtung von , wodurch das Skalarprodukt negativ wird.


Projektion eines Vektors auf einen Vektor zur Berechnung des Skalarprodukts


Schließlich sei auch die kartesische Darstellung des Skalarprodukts gegeben. Man multipliziert jeweils die Komponenten der beiden Vektoren miteinander und summiert die einzelnen Produkte auf.

Das skalarprodukt zweier Vektoren vektor1 und vektor2 lässt sich durch folgende Prozedur berechnen. Beide sind Arrays mit je 3 Double, der Rückgabewert ebenfalls vom Typ Double.


double skalarprodukt (double *vektor1, double *vektor2)
{
  return (vektor1[0] * vektor2[0] + vektor1[1] * vektor2[1] + vektor1[2] * vektor2[2]);
}


Das Skalarprodukt taucht in mehreren in diesem Buch verwendeten Formeln auf, wobei jedoch keine herausragende Anwendung genannt werden kann. In der Geometrie wird das Skalarprodukt z.B. verwendet, um nachzuprüfen, ob zwei Geraden aufeinander senkrecht stehen. Dies ist dann der Fall, wenn das Skalarprodukt ihrer Richtungsvektoren gleich 0 ist.

Das Kreuzprodukt erinnert noch weniger an die gewöhnliche Multiplikation als das Skalarprodukt. Es liefert wieder einen Vektor und wird durch ein "x" zwischen den zu kombinierenden Vektoren gekennzeichnet (z.B. ). Der resultierende Vektor steht senkrecht auf den beiden Ausgangsvektoren, sein Betrag entspricht der Fläche des von diesen aufgespamnnten Parallelogramms:


Kreuzprodukt zweier Vektoren x


Die kartesische Darstellung des Kreuzprodukts ist relativ kompliziert:

Im Gegensazu zu Vektoraddition und Skalarprodukt folgt das Kreuzprodukt dem Kommutativgesetz nicht, es gilt:

Auch für das Kreuzprodukt sei eine Prozedur gegeben. Übergeben werden dieser die zu untersuchenden Vektoren vektor1 und vektor2, zurückgegeben wird das Ergebnis vektor3. Sämtliche Vektoren sind Arrays mit je 3 Double.


void kreuzprodukt (double *vektor1, double *vektor2, double *vektor3)
{
  vektor3[0] = vektor1[1] * vektor2[2] - vektor1[2] * vektor2[1];
  vektor3[1] = vektor1[2] * vektor2[0] - vektor1[0] * vektor2[2];
  vektor3[2] = vektor1[0] * vektor2[1] - vektor1[1] * vektor2[0];
}


Eine klassische Anwendung des Kreuzprodukts in der Astronomie ist die Bestimmung des Bahndrehimpulses der Planetenbewegung, welcher durch das Kreuzprodukt des Orts- und Geschwindigkeitsvektors des Planeten, multipliziert mit dessen Masse, gegeben ist. Der Bahndrehimpuls verändert sich während des Umlaufs um die Sonne nicht, woraus das Zweite Keplersche Gesetz folgt, wonach der Ortsvektor (als Verbindungslinie zwischen Sonne und Planet) in gleichen Zeiten gleiche Flächen überstreicht. Dieses Gesetz wird oft auch als Flächensatz bezeichnet.


Zweites Keplersche Gesetz


Schwerpunkt eines Systems von Massenpunkten[Bearbeiten]

Für die Simulation eines Mehrkörpersystems wählt man wie bereits erwähnt als Koordinatenursprung zumeist dessen Schwerpunkt. Man erhält dessen Ortsvektor , indem man für die Ortsvektoren der einzelnen Mitglieder deren gewichtetes Mittel bildet, wobei als Gewichtsfaktoren die individuellen Massen dienen:

Jede Komponente von stellt eine Art Mittelwert der entsprechenden Komponenten der Einzelobjekte dar. Im Gegensatz zum gewöhnlichen (arithmetischen) Mittel tragen die individuellen Körper jedoch nicht in gleichem Maße dazu bei, sondern mehr oder weniger je nach ihrer Masse. Ist z.B. doppelt so groß wie , so wird das Mitglied 1 bei der Mittelwertbildung gegenüber dem Mitglied 2 quasi doppelt gezählt. Der Schwerpunkt stellt somit den Mittelpunkt des Systems dar, wobei aber berücksichtigt wird, wie die einzelnen Massen verteilt sind.

Ein alltägliches Beispiel für den Schwerpunkt als gewichtetes Mittel unterschiedlicher Massen ist der Hebel. Um mit einer kleinen Masse z.B. eine doppelt so große Masse aufwiegen zu können, muss der Hebelarm für doppelt so lang sein wie für . Der Schwerpunkt des Sonnensystems liegt sehr nahe am Sonnenmittelpunkt (nahe deren Oberfläche), weil die Masse der Sonne sehr viel größer ist als diejenige selbst des massereichsten Planeten, des Jupiter.

Um die Ortsvektoren der einzelnen Massenpunkte auf den Schwerpunkt zu beziehen, zieht man einfach dessen Ortsvektor von diesen ab:

Lässt man eine Mehrkörpersimulation ablaufen, so ist es wünschenswert, dass der Schwerpunkt selbst ruht, also im Koordinatenursprung verharrt. Dies gelingt, indem man auch die Geschwindigkeiten in das Schwerpunktsystem umrechnet. Man geht dabei vor wie für die Ortsvektoren, d.h. bestimmt erst die Geschwindigkeit des Schwerpunkts und zieht diese dann von den ursprünglichen Geschwindigkeiten der Einzelobjekte ab.

Die Geschwindigkeit des Schwerpunkts folgt in Analogie zu seinem Ortsvektor aus dem gewichteten Mittel der Geschwindigkeiten der individuellen Körper:

Entsprechend gewinnt man die korrigierten Einzelgeschwindigkeiten aus den originalen durch Subtraktion der Schwerpunktgeschwindigkeit:

Die Frage, wie die anfänglichen Positionen und Geschwindigkeiten der Massenpunkte festzulegen sind, wird im Praxiskapitel ausführlich behandelt. Es wird dort aufgezeigt, dass die dafür erforderliche Herangehensweise sehr stark von dem Typ des zu untersuchenden Systems abhängt.

Die Umrechnung von Orten und Geschwindigkeiten ins Schwerpunktsystem kann mit folgender Prozedur bewerkstelligt werden. Übergeben werden dieser die Anzahl der Körper N, die Gesamtmasse M des Systems, die Einzelmassen m sowie die Vektoren r im ursprünglichen System. N ist vom Typ unsigned Integer, M eine Double. An für sich könnte die Gesamtmasse auch innerhalb der Prozedur bestimmt werden. Jedoch wird diese auch für andere Anwendungen benötigt, so dass M für die später diskutierten Praxisbeispiele vorab berechnet wird.

m ist ein eindimensionales Array, welches aus N Double besteht entsprechend der Masse jedes einzelnen Objekts. r stellt ein zweidimensionales Array dar, welches N Vektoren enthält, die wiederum je 3 Double auflisten. Diese Vektoren können sowohl Orte als auch Geschwindigkeiten repräsentieren, je nachdem was man umrechnen will.

Die Prozedur bestimmt zuerst den Zähler obiger Ausdrücke. Anschließend wird dieser mit M dividiert. Der so gewonnene Schwerpunkt wird durch das Array s repräsentiert. Im letzten Schritt erfolgt für r der Übergang ins Schwerpunktsystem, d.h. dieses Array nimmt am Ende auch die Rückgabewerte auf. Die Ausgangswerte werden dabei überschrieben, d.h. nicht weiter verwendet.


void schwerpunkt (unsigned int N, double M, double *m, double **r)
{
  unsigned int i,k;
  double s[3];

  for (k = 0;k < 3;k ++)
  s[k] = 0;

  for (i = 0;i < N;i ++)
  {
    for (k = 0;k < 3;k ++)
    s[k] += m[i] * r[i][k];
  }
  for (k = 0;k < 3;k ++)  
  s[k] /= M;

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


Die Koordinatentransformation in das Schwerpunktsystem muss nur ein einziges Mal zu Beginn der Simulation durchgeführt werden. Begründet ist dies durch die Erhaltung der Energie (welche im Unterkapitel Energieerhaltung detailliert besprochen wird) sowie die Unabhängigkeit der physikalischen Gesetze vom Beobachtungsort und der Blickrichtung. Die fundamentale Annahme, dass überall im Kosmos die gleichen physikalischen Gesetze gelten, bezeichnet man als Homogenität des Raumes. Die nicht weniger bedeutsame Annahme, dass (auf genügend großer Längenskala) das Universum in jeder Blickrichtung die gleichen Eigenschaften aufweist, als Isotropie des Raumes.

Mittlere und momentane Geschwindigkeit[Bearbeiten]

Definition[Bearbeiten]

Um die Bahnen der einzelnen Mitglieder eines Sternsystems vorhersagen zu können, braucht man ein Verfahren, das deren momentane Geschwindigkeiten zu jedem Zeitpunkt liefert (es soll zunächst nur vom Betrag der Geschwindigkeit die Rede sein). Tatsächlich messen kann man jedoch nur die mittlere Geschwindigkeit , welche sich auf die in einem gewissen Zeitraum zurückgelegte Strecke bezieht.

Kommt man z.B. mit dem Auto in 2 Stunden 200 km weit, so beträgt die mittlere Geschwindigkeit 200 km / 2 h = 100 km/h. Dem Ideal der momentanen Geschwindigkeit kann man aber zumindest nahe kommen, indem man den durch die Messung überspannten Zeitraum sehr klein macht. Dies geschieht z.B. durch den Tachometer, welcher die Drehzahl (und damit mit Hilfe des Abrollumfangs der Räder die gefahrene Strecke) auf einer Zeitskala von Zehntelsekunden erfasst (gemessen wird dabei allein der Betrag der Geschwindigkeit, nicht jedoch deren Richtung).


Zurückgelegte Strecke als Fläche unter der Kurve momentane Geschwindigkeit gegen die Zeit


Bei einer für jeden beliebigen Zeitpunkt bekannten momentanen Geschwindigkeit lässt sich die während eines bestimmten Zeitraums zurückgelegte Strecke bestimmen, indem man gegen aufträgt und die Fläche unter der entsprechenden Kurve ermittelt. Bei konstanten Geschwindigkeit ist diese Fläche ein Rechteck entsprechend der bekannten Regel Strecke = Geschwindigkeit ∙Zeit.

Für eine Mehrkörpersimulation müssen Positionen und Geschwindigkeiten wie zu Beginn dieses Kapitels diskutiert natürlich als Vektoren betrachtet werden. Man erhält den mittleren Geschwindigkeitsvektor, indem man für jede Komponente des Ortsvektors deren Änderung mit der Zeit betrachtet:

Bei konstantem Geschwindigkeitsvektor gilt wiederum für die Änderung des Ortsvektors während einer bestimmten Zeit:


Einschub für Fortgeschrittene: Definition der Geschwindigkeit

Die mittlere Geschwindigkeit folgt aus der während eines Zeitraums zurückgelegten Strecke gemäß:

Je kleiner man wählt, umso besser stimmt die mittlere Geschwindigkeit mit der momentanen Geschwindigkeit überein. Im Grenzfall gehen die beiden Geschwindigkeiten ineinander über.

Aus dieser Definition von als Grenzwert von folgt, dass die momentane Geschwindigkeit durch die erste Zeitableitung der Strecke gegeben ist.

Umgekehrt gewinnt man die während eines endlichen Zeitraums bewältige Distanz durch Integration über .

Durch Division mit erhält man den exakten Zusammenhang zwischen momentaner und mittlerer Geschwindigkeit.

Mittlere und momentane Beschleunigung[Bearbeiten]

Definition[Bearbeiten]

Die Kenntnis der momentanen Geschwindigkeiten - wiederum sollen zunächst nur Beträge statt Vektoren diskutiert werden - setzt wiederum diejenige der momentanen Beschleunigungen voraus. Diese stellt wie die momentane Geschwindigkeit eine Idealisierung dar, denn auch eine Änderung der Geschwindigkeit lässt sich nur für einen bestimmten von Null verschiedenen Zeitraums erfassen. Daraus resultiert eine mittlere Beschleunigung

Diese Mittelung kommt auch in Redewendungen wie „in 10 s von 0 auf 100“ zum Ausdruck. Um die entsprechende mittlere Beschleunigung zu bestimmen, müssen für Zeit und Geschwindigkeit natürlich analoge Einheiten benutzt werden. 100 km/h entsprechen 100 / 3.6 ≈ 28 m/s, woraus eine mittlere Beschleunigung von 28 m/s / 10 s = 2.8 m/s2 folgt – pro Sekunde nimmt die Geschwindigkeit im Mittel um 2.8 m/s zu.

Wie für die Geschwindigkeit kann man auch für die Beschleunigung zumindest näherungsweise Momentanwerte gewinnen, indem man nur über sehr kurze Zeiträume misst. Fährt man mit hoher Geschwindigkeit in eine Kurve, so erfährt man plötzlich eine senkrecht zur Fahrtrichtung nach „außen“ wirkende Beschleunigung.

In der technischen Praxis wird die Beschleunigung nicht mit Hilfe obiger Definition, sondern auf Grundlage des Newtonschen Kraftgesetzes bestimmmt, in dem die auf eine bekannte Masse ausgeübte Kraft gemessen wird. Die Definition der Beschleunigung als Änderung der Geschwindigkeit pro Zeitintervall ist jedoch fundamental für die im nächsten Kapitel vorgestellten Lösungsmethoden des Mehrkörperproblems.

Auf das Newtonsche Kraftgesetz wird im folgenden Unterkapitel ausführlich eingegangen. Anwendungen wie Trägheitsnavigation oder Seismologie belegen, dass man anhand dieses Prinzips Beschleunigungen auch für sehr kurze Zeiträume messen und damit dem Ideal der Momentanbeschleunigung sehr nahe kommen kann.


Änderung der Geschwindigkeit als Fläche unter der Kurve momentane Beschleunigung gegen die Zeit


Die während eines gewissen Zeitraums eintretende Änderung der Geschwindigkeit lässt sich auf ähnliche Weise berechnen wie die zurückgelegte Strecke. Man trägt nun die Momentanbeschleunigung gegen die Zeit auf und betrachtet wiederum die Fläche unter der dazugehörigen Kurve. Bei gleichbleibender Beschleunigung erhält man abermals ein Rechteck in Übereinstimmung mit der Vorschrift Änderung der Geschwindigkeit = Beschleunigung ∙ Zeit.

Ebenso wie Ort und Geschwindigkeit stellt auch die Beschleunigung einen Vektor dar. Um den mittleren Beschleunigungsvektor zu bestimmen, muss man für den Geschwindigkeitsvektor die zeitliche Änderung einer jeden Komponente untersuchen:

Im Falle eines konstanten Beschleunigungsvektors lautet dementsprechend die Änderung des Geschwindigkeitsvektors im Verlauf eines gewissen Zeitintervalls:


Einschub für Fortgeschrittene: Definition der Beschleunigung

Die folgende Diskussion ist derjenigen über die Geschwindigkeit vollkommen analog. Die mittlere Beschleunigung ist gegeben durch die während eines Zeitraums sich einstellenden Änderung der Geschwindigkeit .

Erneut ist der momentane Wert als Grenzwert des Mittelwerts definiert.

Dies bedeutet, dass die momentane Beschleunigung mit der ersten Zeitableitung der Geschwindigkeit identisch ist.

Die während eines endlichen Zeitraums beobachtbare Änderung der Geschwindigkeit folgt aus dem Integral über .

Wieder liefert die Division mit die Beziehung zwischen momentanem und mittlerem Wert.

Beschleunigung und Strecke[Bearbeiten]

Oft ist nicht nur der Zusammenhang zwischen Beschleunigung und Geschwindigkeit, sondern auch direkt zwischen Beschleunigung und zurückgelegter Strecke von Interesse (wobei abermals zunächst nur Beträge diskutiert werden). Für den Spezialfall einer konstanten Beschleunigung lässt er sich wiederum graphisch recht einfach herleiten.


Zurückgelegte Strecke als Dreiecksfläche unter der Kurve momentane Geschwindigkeit gegen die Zeit bei Bewegung mit konstanter Beschleunigung


Wird ein Körper aus der Ruhelage heraus konstant beschleunigt, so nimmt seine Geschwindigkeit direkt proportional zur verflossenen Zeit zu. Die dabei zurückgelegte Strecke ist dann als Dreiecksfläche unter der Geschwindigkeitskurve gegeben. Die Grundlinie des Dreiecks ist durch die verstrichene Zeit , dessen Höhe durch gegeben. Aus der wohlbekannten Beziehung Dreiecksfläche = Grundlinie ∙ Höhe / 2 folgt so:

Diese Beziehung gibt z.B. an, welche Strecke man nach einer gewissen Zeit beim freien Fall unter Vernachlässigung des Luftwiderstands zurückgelegt hat. Für ist dann der Wert von 9.81 m/s2 für die Fallbeschleunigung einzusetzen. Nach 1 s ist man 4.9 m, nach 2 s aber schon das vierfache, nämlich 19.6 m gefallen.

Setzt man bzw. in obige Formel ein, gewinnt man folgende von der Zeit unabhängige Beziehung zwischen Strecke und Geschwindigkeit:

Damit kann man z.B. berechnen, mit welcher Geschwindigkeit man bei einer bestimmten Fallhöhe auf dem Boden ankommt. Bei einem Sprung vom 5-Meter-Brett beträgt die Aufprallgeschwindigkeit 9.9 m/s entsprechend 36 km/h. Ein vom 10-Meter-Brett fallender Turmspringer bringt es auf 14.0 m/s bzw. 50 km/h.

Weist ein konstant beschleunigter Körper zu Beginn bereits eine gewisse Geschwindigkeit auf, so muss zur nach einer Zeit bewältigten Strecke der Anteil hinzuaddiert werden, so dass gilt.

Betrachtet man Ort, Geschwindigkeit und Beschleunigung abermals als Vektoren, so ist obige Formel für jede Komponente einzeln anzuwenden:


Einschub für Fortgeschrittene: Beschleunigung und zurückgelegte Strecke

Aus den Definitionen von momentaner Beschleunigung und Geschwindigkeit folgt, dass erstere durch die zweite Zeitableitung der Strecke gegeben ist.

Dadurch ist die Berechnung der zurückgelegten Distanz aus der momentanen Beschleunigung allgemein etwas aufwendiger. Die momentane Geschwindigkeit nach einem Zeitraum lässt sich schreiben als:

Die nochmalige Integration über obigen Ausdruck liefert für :

Für einen aus dem Ruhezustand heraus konstant beschleunigten Körper erhält man daraus wiederum die bereits oben graphisch abgeleitete klassische Beziehung.

Das Newtonsche Kraftgesetz[Bearbeiten]

Als nächstes stellt sich die Frage, wie die Beschleunigungen zu bestimmen sind, welche auf die Mitglieder eines Mehrkörpersystems einwirken. Den Schlüssel hierzu liefert das Newtonsche Kraftgesetz, welches anschaulich oft als Kraft = Masse • Beschleunigung formuliert wird.

Um ein Auto der Masse von 1000 kg in 10 Sekunden von 0 auf 100 km/h zu beschleunigen, ist demnach eine Kraft von 1000 kg • 2.8 m/s2 = 2800 Newton erforderlich (wobei hier Luftwiderstand und Haftreibung vernachlässigt sind). Kennt man die auf einen Massenpunkt einwirkende Kraft, so folgt aus dem Kraftgesetz die Beziehung Beschleunigung = Kraft / Masse.

Die Beschleunigung, welche ein Körper erfährt, ist also der auf ihn ausgeübten Kraft direkt und seiner Masse umgekehrt proportional. Um eine doppelt so große Masse der gleichen Beschleunigung zu unterwerfen, ist die doppelte Kraft erforderlich.

Das Kraftgesetz betrachtet nicht nur die Beträge von Kraft und Beschleunigung, sondern auch deren Richtungen. Es gilt die einleuchtende Regel, dass die Beschleunigung in dieselbe Richtung erfolgt, in welcher die Kraft gerichtet ist. Dies kommt auch schon in obiger Vektorschreibweise zum Ausdruck, wonach der Beschleunigungsvektor lediglich mit einem Skalar > 0 (nämlich der Masse) multipliziert werden muss, um den Kraftvektor zu erhalten. Bei einer solchen Rechenoperation aber bleibt - wie schon ganz zu Anfang dieses Kapitels erläutert - die ursprüngliche Richtung des Vektors erhalten.

Wie schon erwähnt, ermöglicht es das Newtonsche Kraftgesetz, Beschleunigungen zu messen, indem man die auf eine bekannte Testmasse einwirkenden Kräfte betrachtet. Man stelle sich z.B. einen Aufzug vor, in welchem an einer Federwaage eine Masse von 1 kg hängt. Ruht der Aufzug, wirkt allein die Schwerebeschleunigung (etwa 9.8 m/s2) auf die Masse ein. Die Federwaage zeigt dann eine Kraft von 9.8 Newton an. Setzt sich der Aufzug nach oben in Bewegung, so wird die Testmasse einer zusätzlichen Beschleunigung unterworfen. Die Feder wird dadurch stärker ausgelenkt, sie zeigt jetzt eine größere Kraft an, woraus die Beschleunigung des Aufzugs gemäß bestimmt werden kann. Zeigt die Federwage z.B. 12 Newton hat, so entspricht das einer Beschleunigung von (12 / 1 - 9.8) m/s2 = 2.2 m/s2.

Himmelskörper kann man natürlich nicht als Testmassen benutzen, um die auf sie ausgeübten Beschleunigungen zu messen. Wie im nächsten Unterkapitel gezeigt, liefert für solche das Newtonsche Gravitationsgesetz belastbare Kriterien für die wechselseitig herrschenden Kräfte und damit Beschleunigungen.

Das Newtonsche Gravitationsgesetz[Bearbeiten]

Kernaussagen[Bearbeiten]

Die Berechnung der in einem Mehrkörpersystem herrschenden Kräfte selbst beruht auf dem Newtonschen Gravitationsgesetz. Dieses besagt, dass die zwischen zwei Massen und herrschende Anziehungskraft diesen direkt und dem Quadrat ihres Abstandes umgekehrt proportional ist, und die Kraft in Richtung ihrer Verbindungslinie wirkt. soll dabei die Kraft bezeichnen, welche der Körper auf den Körper ausübt.

bezeichnet dabei die Gravitationskonstante mit dem Wert 6.67384 10-11 m3 kg-1 s-2. ist der Einheistvektor auf der Verbindungslinie der beiden Massen, er zeigt von nach .

In der Praxis interessiert nicht die auf einen Massenpunkt einwirkende Kraft, sondern nur dessen Beschleunigung. Aus dem Newtonschen Kraftgesetz folgt, dass die Masse folgende Beschleunigung durch die Masse erfährt:

Entsprechend erleidet die Masse eine Beschleunigung durch die Masse in entgegengesetzter Richtung.

Aus der umgekehrt quadratischen Abhängigkeit der Anziehungskraft vom Abstand der Massen folgt, dass bei doppeltem Abstand diese auf ein Viertel des Vergleichswertes abgefallen ist, bei dreifachem Abstand auf ein Neuntel usw. Die Schwerkraft fällt also nur langsam mit zunehmenden Abstand ab. Im Prinzip besitzt sie eine unendliche Reichweite, was bei der Modellierung von Mehrkörperensembles erhebliche Schwierigkeiten nach sich ziehen kann. Andererseits nimmt die wechselseitige Anziehungskraft enorm zu, wenn zwei Mitglieder sich sehr nahe kommen. Da solche zumeist als Massenpunkte behandelt werden, können im Verlauf einer Simulation grundsätzlich beliebig kleine Abstände und damit unbegrenzt große Kräfte auftreten. Auch dieser Effekt erschwert eine zuverlässige Behandlung astronomischer Systeme enorm.

Um die gesamte Beschleunigung eines Massenpunkts zu kennen, müssen sämtliche Einzelbeschleunigungen aufaddiert werden, welche alle übrigen Mitglieder des Ensembles auf diesen ausüben. Eine solche Gesamtbetrachtung ist für alle Körper notwendig, d.h. es müssen alle möglichen Abstandspaare betrachtet werden, die man für die individuellen Massenpunkte bilden kann. Daraus folgt, dass der Aufwand zur Berechnung aller Anziehungskräfte quadratisch mit der Anzahl der Mitglieder eines Mehrkörpersystems ansteigt. Eine exakte Erfassung der Kräfte ist mit vertretbarem Rechenaufwand somit nur für eine Simulation von bis zu einigen 1000 Massenpunkten möglich. Für komplexere Systeme ist man auf Näherungsverfahren angewiesen, welche später im Detail erläutert werden.

Wie bereits erörtert, werden für praktische Simulationen oft kartesische Koordinaten mit dem Schwerpunkt des zu untersuchenden Ensembles als Ursprung verwendet. Die momentane Position der Masse sei so durch einen Vektor und diejenige der Masse durch einen solchen gegeben.


Berechnung der wechselseitigen Beschleunigung zweier Massen durch deren Anziehungskraft in kartesischen Koordinaten


Für die Berechnung der wechselseitigen Beschleunigungen der beiden Massen kommt es allein auf den Abstandsvektor an. In kartesischen Koordinaten lautet dieser:

Den Einheitsvektor gewinnt man, indem man den Abstandsvektor mit seinem eigenen Betrag dividiert:

Einsetzen in die oben gegebenen Beschleunigungen liefert:

Der Betrag des Abstandsvektors folgt aus dem Satz des Pythagoras:

Elementare Schlussfolgerungen[Bearbeiten]

An dieser Stelle sollen bereits einige Schlussfolgerungen gezogen werden, welche für die Simulation von Mehrkörpersystemen von Bedeutung sind. Dazu soll das sehr einfache Beispiel einer kleinen Masse betrachtet werden, welche um eine große Masse auf einer Kreisbahn mit Radius umläuft.

Dabei sei zuerst die Bahngeschwindigkeit betrachtet. Man erhält sie durch Gleichsetzen der Gravitationsbeschleunigung mit der Zentripetalbeschleunigung :

Je weiter ein Körper von seiner Zentralmasse entfernt ist, um so geringer ist seine Bahngeschwindigkeit, was auch für Ellipsenbahnen zutrifft. Im Falle einer Kreisbahn nimmt umgekehrt proportional zur Quadratwurzel der Entfernung ab, d.h. bei vierfacher Entfernung beträgt die Bahngeschwindigkeit nur noch die Hälfte, bei neunfachem Abstand noch ein Drittel des Bezugswertes usw.

Aus der Bahngeschwindigkeit und dem Umfang der Bahn folgt die Umlaufsdauer . Es gilt , woraus durch Einsetzen und Quadrieren folgt:

Das Quadrat der Umlaufsdauer ist dem Kubus des Bahnradius proportional, wobei die Proportionalitätskonstante durch die Zentralmasse bestimmt ist. Diese Beziehung gilt unverändert auch für Ellipsen, wobei an die Stelle des Bahnradius die große Halbachse tritt. Es handelt sich hierbei um das 3.Keplersche Gesetz.

Verdoppelt man die Entfernung, so steigt die Umlaufsdauer um den Faktor entsprechend etwa 2.8. Bei dreifacher Entfernung nimmt um das -fache zu, entsprechend ungefähr dem Faktor 5.2.

Das 3.Keplersche Gesetz bietet eine sehr anschauliche Prüfmöglichkeit für die Gültigkeit des Gravitationsgesetzes und damit der daraus abgeleiteten Kräfte und Beschleunigungen. Die Entfernung der Erde von der Sonne kann aus Parallaxenmessungen bestimmt werden, die Abstände zu den übrigen Planeten des inneren Sonnensystems auch aus Radarmessungen. Die Verhältnisse der Abstände zur Sonne sind also mit hoher Genauigkeit bekannt und können mit den ebenfalls genau bekannten Verhältnissen der Umlaufsdauern verglichen werden (eine Gegenüberstellung der Absolutentfernungen zur Sonne und der absoluten Werten der Umlaufsdauern erforderte auch noch eine unabhängige Kenntnis der Sonnenmasse).

Man muss bei dieser Argumentation allerdings beachten, dass das 3.Keplersche Gesetz streng genommen nur für ein aus lediglich zwei Körpern bestehendes System gültig ist. Allerdings wird die auf einen Planeten einwirkende Schwerkraft fast völlig von der Sonne dominiert, so dass dessen Bahn als eine von den übrigen Planeten nur geringfügig gestörte Ellipse betrachtet werden darf. Auf diese Problematik wird im Praxiskapitel ausführlich eingegangen.

Sieht man sich Umlaufsdauer, Bahngeschwindigkeit und Gravitationsbeschleunigung nochmals genau an, so stellt man fest, dass sie folgenden Zusammenhang erfüllen:

Wie später gezeigt wird, ist dieses Verhältnis von Geschwindigkeit und Beschleunigung eines Massenpunktes ein fundamentales Kriterium für die zeitliche Schrittweite, mit welcher dessen Bewegung in einem Mehrkörpersystem modelliert werden soll. Es wird sehr häufig dynamische Zeit genannt.

C-Code: Beschleunigungen im Mehrkörpersystem

Nun stehen alle erforderlichen Kenntnisse für die praktische Berechnung der auf einen Massenpunkt einwirkenden Beschleunigung bereit. Angesichts der immer gleichartigen Berechnungen für jede Vektorkomponente bietet es sich an, die Positionen (x,y,z) zu einem Array r und die Beschleunigungen (ax,ay,az) zu einem Array a für den Variablentyp Double zusammenzufassen. Da r und a für jeden Körper gebraucht werden, müssen diese Größen tatsächlich als zweidimensionale Arrays deklariert werden. Die Masse m muss ebenfalls für jedes Mitglied des Ensembles gegeben sein. Als Skalar erfordert diese nur ein eindimensionales Array wiederum für den Variablentyp Double.

Die Unsigned Integer i und k dienen dazu, die einzelnen Massenpunkte und deren Vektorkomponenten abzuzählen. Die Unsigned Integer N gibt die Gesamtzahl der Körper an. Für den Abstandsvektor dr zweier Objekte genügt ein einfaches Array für den Variablentyp Double mit 3 Elementen. Die Double d3 schließlich bezeichnet die 3. Potenz dessen Betrags.

Programmiertechnisch werden zweidimensionale Vektorarrays wie r und a im Folgenden als eindimensionale Arrays von Zeigern behandelt. Für jeden Massenpunkt gibt es einen Zeiger, der auf jede der 3 Komponenten des entsprechenden Vektors deuten kann. Eindimensionale Skalararrays wie m werden jeweils durch einen einzigen Zeiger realisiert. Ein solcher Zeiger kann auf jeden der N Körper deuten.

Die Berechnung der Beschleunigung ist ein während einer Mehrkörpersimulation sich ständig wiederholender Vorgang. Somit ist es naheliegend, diese als Prozedur zu definieren, welcher die Nummer objekt des gerade zu bearbeitenden Massenpunkts, N, m und r übergeben werden. Zurückgegeben wird die Beschleunigung a[objekt]. Mit dieser Konstruktion werden i, k, dr und d3 zu lokalen Variablen. Aufgerufen wird die Prozedur durch die Anweisung beschleunigung (objekt,N,m,r,a[objekt]).


/* Globale Variablen */

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

void beschleunigung (unsigned int objekt, unsigned int N, double *m, double **r, double *a)
{

/* Lokale Variablen */

  unsigned int i,k;
  double dr[3];
  double d3;

/* Initialisierung der Beschleunigung */

  for (k = 0;k < 3;k ++)
  a[k] = 0;

/* Berechnung der Beschleunigung */

  for (i = 0;i < N;i ++)
  {
    if (i != objekt)
    {

/* Abstandsvektor zwischen einem beliebigen Körper i und dem untersuchten Massenpunkt objekt */

      for (k = 0;k < 3;k ++)
      dr[k] = r[i][k] - r[objekt][k];

/* 3.Potenz dessen Betrags */

      d3 = pow (betrag (dr),3);

/* Beitrag der Beschleunigung durch den Körper i */

      for (k = 0;k < 3;k ++)
      a[k] += G * m[i] * dr[k] / d3;
    }
  }
}


Vor allem für mehrdimensionale Arrays ist eine sorgfältige Speicheradressierung erforderlich, die hier zumindest kurz angedeutet werden soll. Benötigt werden hierzu zumindest die Anweisungen malloc (um Speicher zu adressieren) und free (um solchen wieder freizugeben).

Für eindimensionale Skalararrays wie die Masse ist die Speicheradressierung sehr einfach. Mit der Anweisung malloc wird für jeden Körper Speicher für eine Double bereitgestellt.

m = malloc (N * sizeof (double));

Bei zweidimensionalen Arrays wie hier etwa der Position erfolgt die Speicheradressierung in zwei Schritten. Zunächst wird für jedes Objekt Speicher für einen Zeiger (Vektorobjekt) reserviert (deshalb sizeof (double *) und nicht sizeof (double)). Anschließend wird für jeden dieser Zeiger Platz für jeweils 3 Double (Vektorkomponenten) geschaffen.

r = malloc (N * sizeof (double *));
for (i = 0;i < N;i ++)
r[i] = malloc (3 * sizeof (double));

Um reservierten Speicher wieder freizugeben, genügt für Skalararrays abermals eine einzige Anweisungen. Mit free wird der gesamte durch die N Double beanspruchte Platz wieder frei.

free (m);

Zweidimensionale Arrays verlangen auch bei der Speicherfreigabe ein zweistufiges Vorgehen, wobei aber im Vergleich zur Speicheradressierung die Reihenfolge umgekehrt wird. Zuerst muss für jeden Körper der durch die Vektorkomponenten belegte Speicher freigegeben werden. Danach kann der durch die Vektorobjekte selbst belegte Platz geräumt werden.

for (i = 0;i < N;i ++)
free (r[i]);
free (r);

Für eine detaillierte Behandlung der Speicherverwaltung mittels C oder C++ sei auf die einschlägige Literatur (auch unter Wikibooks selbst) verwiesen.

Energieerhaltung[Bearbeiten]

Jedes noch so komplexe System weist unabhängig von seinem genauen Aufbau einige allgemeingültige Eigenschaften auf. So ist in jedem isolierten System die Gesamtenergie konstant. Allerdings sind die meisten natürlichen Systeme nicht isoliert, z.B. tauscht jedes Lebewesen Energie und Stoffe mit seiner Umgebung aus. Auch Sternsysteme sind es streng genommen nie, da die von einer Masse ausgehende Schwerkraft wie bereits erwähnt zwar mit zunehmender Entfernung immer schwächer wird, aber niemals ganz verschwindet.

Ist ein Mehrkörperensemble im Vergleich zum Abstand zu Körpern ähnlicher Masse klein genug, dürfen die von diesen fernen Nachbarn ausgeübten Kräfte meist vernachlässigt werden. Eine solche Entscheidung bedarf allerdings eines gewissen Fingerspitzengefühls. Selbst die nächsten Sterne sind zu weit entfernt, um die Bahnen der Planeten nachweisbar zu beeinflussen. In diesem Sinne ist das Sonnensystem faktisch isoliert. Die in der sogenannten Oortschen Wolke sich aufhaltenden Kometen jedoch sind so weit von der Sonne entfernt, dass nahe Sterne sehr wohl deren Orbits massiv stören können. Auf diese Weise gelangen immer wieder einige von ihnen in das innere Sonnensystem, das sich so doch in Wechselwirkung mit den Sternen seiner Umgebung befindet.

Das abstrahierten Massenpunkte einer Computersimulation dürfen selbstverständlich immer als isoliert betrachtet werden. Die große Schwierigkeit besteht darin, zu entscheiden, ob eine solche Sichtweise den tatsächlichen Gegebenheiten angemessen ist.

Kinetische Energie[Bearbeiten]

Ein allein der Schwerkraft unterworfenes System weist zwei Formen von Energie auf, kinetische Energie und potentielle Energie. Unter der kinetischen Energie (auch Bewegungsenergie genannt) versteht man die Energie, welche man aufwenden muss, um eine Masse vom Ruhezustand aus auf eine bestimmte Geschwindigkeit zu beschleunigen. Sie ist der beschleunigten Masse und dem Quadrat der Endgeschwindigkeit direkt proportional. Um eine doppelt so hohe Geschwindigkeit zu erreichen, muss das vierfache an Energie zur Verfügung stehen. Um die gesamte kinetische Energie eines Mehrkörpersystems zu kennen, müssen die Energien aller einzelnen Massenpunkte aufaddiert werden.

Der Zusammenhang zwischen Bewegungsenergie, Masse und Geschwindigkeit lässt sich leicht aus dem Prinzip Arbeit = Kraft • Weg herleiten, wobei hier wieder nur Beträge betrachtet werden sollen. Man betrachte einen Körper der Masse , welcher aus dem Ruhezustand auf einer Strecke einer konstanten Beschleunigung unterworfen wird. Nach dem Newtonschen Kraftgesetz ist dazu eine Kraft

erforderlich. Das Arbeitsprinzip besagt, dass dazu wiederum eine Energie

aufgewandt werden muss. Setzt man die schon aus dem Abschnitt "mittlere und momentane Beschleunigung" bekannte Beziehung zwischen Geschwindigkeit, Beschleunigung und Wegstrecke ein, so erhält man:


Einschub für Fortgeschrittene: Kinetische Energie bei beliebiger Beschleunigung

Obige Formel lässt sich auch für beliebige Beschleunigungen beweisen. sei nun eine nur für eine beliebig kleine Strecke wirksame Momentanbeschleunigung. Um den Weg zurückzulegen, wird die Arbeit

benötigt. Mit Hilfe der Definitionen und resultiert daraus

und durch Integration aus dem Ruhezustand


Potentielle Energie[Bearbeiten]

Die potentielle Energie gibt die Energie an, die man benötigt, um eine Masse vollständig aus dem Anziehungsbereich einer anderen zu entfernen. Bei dieser Beschreibung wird man sofort wieder mit dem Problem der in Wahrheit nie vorhandenen Isolation von Himmelskörpern konfrontiert. Da die Schwerkraft in ihrer Reichweite nicht beschränkt ist, kann ein Körper strenggenommen niemals der Anziehung eines anderen entkommen. Jedoch ist es möglich, wie die Raumfahrt anschaulich gemacht hat, eine (kleine) Masse auf eine Bahn zu bringen, auf welcher sie sich immer weiter von der (großen) Masse entfernt, von wo aus sie gestartet ist. Die Energie, die man für eine derartige Bahn aufwenden muss, ist den beiden Massen direkt und ihrem ursprünglichen Abstand umgekehrt proportional. Je näher sich die Massen beim Start zueinander befinden, um so mehr Energie ist erforderlich, um die beiden effektiv voneinander zu trennen, bei z.B. halbem Startabstand die doppelte Energie. Dies ist freilich nicht überraschend, da bei kleineren Anfangsabstand eine entsprechend größere Anziehungskraft überwunden werden muss.


Potentielle Energie im Schwerefeld als Fläche unter der Kurve Anziehungskraft gegen Entfernung


Den exakten Wert der potentiellen Energie erhält man, indem man die Anziehungskraft vom Startabstand an gegen die Entfernung aufträgt und die Fläche unter der entsprechenden Kurve berechnet. Man folgt also abermals der Regel Arbeit = Kraft • Weg, wobei aber nun sich die Kraft entlang der zurückgelegten Strecke ständig ändert. Obwohl nie gänzlich verschwindet, hat diese Fläche einen endlichen Inhalt. Man erhält als Ergebnis, was leider aber nicht elementar beweisbar ist:

In einem aus vielen Massenpunkten bestehendem System existiert für jedes Mitglied gegenüber allen anderen potentielle Energie. Um ein solches aus dem Gesamtsystem zu entfernen, müssen die Anziehungskräfte jeweils aller anderen Körper überwunden werden. Um die gesamte potentielle Energie zu bestimmen, muss über alle möglichen Abstandspaare addiert werden. Sie wird als negative Größe angesetzt, da sie ja aufgewandt werden muss, um das Ensemble aufzulösen.


Einschub für Fortgeschrittene: Potentielle Energie im Schwerefeld

Der Beweis für obige Beziehung lässt sich folgendermaßen skizzieren. Nach dem Gravitationsgesetz wirkt auf eine Kraft

Um von um eine kleine Wegstrecke zu entfernen, muss die Arbeit

verrichtet werden. Um gänzlich entkommen zu lassen, wird das Integral über die Arbeit vom ursprünglichen Abstand bis ins Unendliche fortgeführt. Daraus ergibt sich:


Kinetische und potentielle Energie zusammen bilden die konstante Gesamtenergie des Systems, doch besteht zwischen beiden ein ständiges Wechselspiel. Entfernen sich zwei Massenpunkte voneinander, werden sie durch ihre gegenseitige Anziehung abgebremst. Sie verlieren kinetische Energie, gewinnen aber an potentieller Energie. Nähern sie sich, verhält es sich umgekehrt. Die Massen werden auf höhere Geschwindigkeiten beschleunigt, aus potentieller Energie wird kinetische Energie.

Ein alltägliches Beispiel für die Umwandlung der einen in die andere Energieform ist das Wasserkraftwerk. Das ins Tal fließende Wasser nimmt an Geschwindigkeit, d.h. an kinetischer Energie zu, welche zum Antrieb einer Turbine genutzt werden kann. Gleichzeitig verliert es potentielle Energie. Um etwa den See eines Pumpspeicherkraftwerks wieder aufzufüllen, muss ja eine andere Energiequelle, z.B. ein Windkraftwerk, angezapft werden.

Elementare Schlussfolgerungen[Bearbeiten]

Betrachtet man die beiden Energieformen für ein System aus 2 Körpern, erwächst aus dem beständigen Hin und Her von kinetischer und potentieller Energie folgende wichtige Regel. Um die beiden Massen für immer voneinander zu entfernen, muss die kinetische Energie des Systems dem Betrage nach mindestens gleich dessen potentieller Energie sein. Überwiegt der Betrag der potentiellen Energie, was einer negativen Gesamtenergie entspricht, bleiben sie einander auf endlichen Bahnen gebunden.

Dies lässt sich anhand einer Kreisbahn mit Radius leicht demonstrieren. Die Bahngeschwindigkeit eines Körpers der Masse , welcher um ein Zentralobjekt der Masse läuft, beträgt wie im letzten Unterkapitel gezeigt und damit seine kinetische Energie . Für die potentielle Energie gilt . Sofern sehr viel größer ist als , kann man davon ausgehen, dass die Zentralmasse sich fast unbeweglich am Schwerpunkt beider Massen aufhält, so dass deren kinetische Energie vernachlässigbar ist. Damit beläuft sich die Gesamtenergie auf:

Die Gesamtenergie eines auf einer Kreisbahn laufenden Körpers ist somit dem Bahnradius umgekehrt proportional. Im Falle einer Ellipsenbahn tritt an die Stelle des Kreisradius die große Halbachse. Bei einer unendlich weiten Bahn ist die Gesamtenergie in der Tat gleich Null

Aus diesem Grenzfall, dass kinetische und potentielle Energie betragsmäßig genau gleich sind, folgt die sogenannte Fluchtgeschwindigkeit . Diese gibt die Mindestgeschwindigkeit an, auf welche die Masse beschleunigt werden muss, um vom Abstand aus die Masse zu verlassen. Das Gleichsetzen der Ausdrücke für die kinetische und potentielle Energie liefert unmittelbar:

Die Fluchtgeschwindigkeit unterschiedet sich von der Bahngeschwindigkeit allein um den Faktor und zeigt somit die gleiche Abhängigkeit von .

Bei Simulationen von Sternsystemen finden sich nach einer gewissen Zeit stets einzelne Massenpunkte, die im Vergleich zu der Mehrzahl der Mitglieder sehr weit vom Schwerpunkt des Ensembles entfernt sind. Für solche Körper können die übrigen Objekte weitgehend als eine im Schwerpunkt ruhende Zentralmasse aufgefasst werden. Zeigen diese Ausreißer eine mit zunehmender Entfernung abfallende Geschwindigkeit kleiner als , so sind sie weiterhin an das Gesamtsystem gebunden. Anderenfalls handelt es sich um Massenpunkte, welche bereits im Begriff sind, das Ensemble zu verlassen.

Der Hinauswurf einzelner Mitglieder aus einem Mehrkörpersystem widerspricht nicht der Energieerhaltung. Die für 2 Massenpunkte geltende Regel, dass ein System mit einem Überschuss an potentieller Energie stabil ist, lässt sich nämlich nicht auf komplexere Ensembles übertragen. Für ein aus 3 oder mehr Körpern bestehendes System garantiert eine negative Gesamtenergie nur eine relative Stabilität. Dort können Massen durch enge Vorübergänge ohne weiteres so stark beschleunigt werden, dass sie dem Gesamtsystem entkommen. Die dazu erforderliche kinetische Energie wird durch eine Verringerung der potentiellen Energie gewonnen. Die zurückgebliebenen Körper rücken im Mittel etwas näher zusammen.

Diese Eigenschaft eines Mehrkörpersystems wird z.B. für die interplanetare Raumfahrt genutzt. Durch nahe Vorübergänge an großen Planeten können Raumsonden erheblich beschleunigt und so die Flugzeiten ins äußere Sonnensystem enorm verkürzt werden.


C-Code: Gesamtenergie im Mehrkörpersystem

Die für die Berechnung der Beschleunigung eines Massenpunktes eingeführten Variablen m, r, dr, i, k und N können auch für die Bestimmung der Energie verwendet werden. Für die kinetische Energie muss die Geschwindigkeit v analog zur Position r als Vektorarray für den Variablentyp Double deklariert werden, die Energie selbst als einfache Double Ekin. Für die potentielle Energie ist das Vektorarray r schon vorhanden, zusätzlich benötigt werden als einfache Double Epot für die Energie selbst und d für den Betrag des Abstandsvektors, ferner als weiterer Zähler die Unsigned Integer j.

Für die Berechnung der Bahnen der Massenpunkte wird die Energie nicht benötigt, wohl aber zur Kontrolle der Stabilität einer Mehrkörpersimulation. Die Bestimmung der Energie erfolgt daher ebenfalls im Rahmen einer regelmäßig aufgerufenen Prozedur. Die Übergabeparameter sind N, m, r und v, die Rückgabevariablen Ekin und Epot. Der Aufruf geschieht durch energie (N,m,r,v,&Ekin,&Epot).

/* Globale Variablen */

unsigned int N;
double *m;
double **r;
double **v;
double Ekin,Epot;

void energie (unsigned int N, double *m, double **r, double **v, double *Ekin, double *Epot)
{

/* Lokale Variablen */
 
  unsigned int i,j,k;
  double dr[3];
  double d;

/* Berechnung der kinetischen Energie */

  *Ekin = 0;
  for (i = 0;i < N;i ++)
  *Ekin += m[i] * pow (betrag (v[i]),2) / 2;

/* Berechnung der potentiellen Energie */

 *Epot = 0;
  for (i = 0;i < N;i ++)
  {

/* Beitrag pro Abstandspaar zur potentiellen Energie */

    for (j = i + 1;j < N;j ++)
    {

/* Abstandsvektor zwischen zwei Körpern j und i und dessen Betrag */

      for (k = 0;k < 3;k ++)
      dr[k] = r[j][k] - r[i][k];
      d = betrag (dr);

/* Potentielle Energie pro Abstandspaar */

      Epot -= G * m[i] * m[j] / d;
    }
  }
}

Ruck[Bearbeiten]

Definition[Bearbeiten]

Wie in den ersten Abschnitten dieses Kapitels gezeigt wurde, ist die mittlere Geschwindigkeit definiert als das Verhältnis zurückgelegte Strecke / Zeit , die mittlere Beschleunigung wiederum als das Verhältnis Geschwindigkeitsänderung / Zeit . Die entsprechenden momentanen Größen gewinnt man, indem man das Zeitintervall so klein wie möglich hält. Die sich anschließende Diskussion von Kraft und Energie macht deutlich, dass zum physikalischen Verständnis eines Systems, in welchem allein die Schwerkraft wirkt, diese Größen ausreichen.

Für bestimmte Lösungsverfahren des Mehrkörperproblems - der Gegenstand des nächsten Kapitels - ist es jedoch erforderlich, über Strecke, Geschwindigkeit und Beschleunigung hinaus weitere Größen durch fortgesetzte Betrachtung der zeitlichen Änderung einzuführen. Als dabei Wichtigste ist der mittlere Ruck zu nennen, dessen Symbol j auf seine englische Bezeichnung "jerk" zurückgeht. Er gibt die pro Zeitintervall auftretende Änderung der Beschleunigung an.

Die mathematische Definition entspricht durchaus der Alltagserfahrung. Das Gefühl, bei einer Notbremsung durchgeschüttelt zu werden, ist nicht allein auf die hohen Absolutwerte der dabei auftretenden Beschleunigungen zurückzuführen, sondern auch auf deren abrupte Änderungen. Tritt man heftig auf die Bremse, so erfährt man sofort eine starke der Fahrtrichtung entgegengesetzte Beschleunigung. Kommt das Auto zum Stillstand, verschwindet diese ebenso rasch wieder. Tatsächlich ist die Fahrdynamik eine der wichtigsten Anwendungen des Rucks.

Wie für die Geschwindigkeit und Beschleunigung kann man durch extrem kleine Momentanwerte auch für den Ruck ableiten und mit den Methoden der Differential- und Integralrechnung allgemeine Zusammenhänge zwischen den momentanen Größen , und aufstellen. Für das Verständnis von auf dem Ruck basierenden Lösungsmethoden genügt aber die Betrachtung des mittleren bzw. eines für kleine Zeitintervalle konstanten Rucks , so dass auf eine weitere Diskussion der momentanen Größen verzichtet wird. Wird ein Körper einem konstanten Ruck unterworfen, so ändert sich seine Beschleunigung gemäß obiger Definition linear mit der Zeit.

Ruck und Geschwindigkeit[Bearbeiten]

Für eine auf dem Ruck beruhende Lösung eines Mehrkörperproblems wird selbstverständlich ein Zusammenhang zwischen diesem und der Geschwindigkeit benötigt. Für den Fall eines konstanten Rucks lässt sich ein solcher elementar ableiten, das Vorgehen entspricht völlig der Beziehung zwischen konstanter Beschleunigung und zurückgelegter Strecke.


Änderung der Geschwindigkeit als Dreiecksfläche unter der Kurve momentane Beschleunigung gegen die Zeit bei Bewegung mit konstantem Ruck


Wirkt der Ruck aus der Ruhelage heraus, so liegt nach einer Zeit eine zu dieser direkt proportionale Beschleunigung vor. Trägt man die Beschleunigung über die Zeit auf, so erhält man als Fläche unter der Kurve (analog zum Problem der bei konstanter Beschleunigung bewältigten Distanz) ein Dreieck. Somit besteht folgender Zusammenhang (man vergleiche mit für konstante Beschleunigung).

Weist ein Körper schon zu Anfang eine Beschleunigung auf, so gilt wegen des Zusatzbeitrags

Ruck und Strecke[Bearbeiten]

Schließlich muss auch die Beziehung zwischen Ruck und zurückgelegter Strecke betrachtet werden. Selbst für den einfachen Fall eines konstanten Rucks lässt sich diese leider ohne höhere Mathematik nicht angeben. Da die Geschwindigkeit wie soeben gezeigt nun quadratisch von der Zeit abhängt, liegt die gesuchte Strecke im Diagramm gegen jetzt als Fläche unter einer Parabel vor, nicht als Dreieck wie im Fall einer konstanten Beschleunigung. Für einen zu Beginn ruhenden Körper liefert die Integralrechnung folgendes Ergebnis.

Existieren schon zu Anfang eine Beschleunigung und eine Geschwindigkeit , so müssen deren Beiträge und berücksichtigt werden. Damit gilt allgemein


Einschub für Fortgeschrittene: Geschwindigkeit, Beschleunigung und Ruck als Taylor-Entwicklung der Strecke

Gemäß des Taylorschen Satzes lässt sich eine reelle Funktion in der Umgebung einer Stelle als Potenzreihe darstellen, sofern diese dort stetig und beliebig oft differenzierbar ist.

Setzt man sowie und betrachtet die ersten 3 Glieder der Summe, so erhält man

Man erkennt leicht, dass dies dem Zusammenhang zwischen konstantem Ruck und der während einer gewissen Zeit bewältigten Distanz entspricht. Der Vergleich liefert unmittelbar, dass Geschwindigkeit, Beschleunigung und Ruck der 1., 2. und 3.Ableitung der Strecke nach der Zeit entsprechen.


Ruck im Schwerefeld[Bearbeiten]

Im Zusammenhang mit einem astronomischen Mehrkörpersystem interessieren natürlich vor allem die Rucks, den zwei Massen und aufgrund der Schwerkraft aufeinander ausüben. Man findet, was sich wiederum einer elementaren Betrachtung entzieht, dass diese gegenseitigen Rucks im Gegensatz zu den Beschleunigungen nicht nur vom Abstandsvektor , sondern auch der vektoriellen Geschwindigkeitsdifferenz zwischen und abhängen. Im ersten Term erscheint diese allein, im zweiten als Bestandteil des Skalarprodukts mit dem Abstandsvektor.


Einschub für Fortgeschrittene: Zeitableitung der Gravitationsbeschleunigung

Obige Formeln erhält man, in dem man die im Abschnitt "Das Newtonsche Gravitationsgesetz" hergeleiteten wechselseitigen Schwerebeschleunigungen zweier Massen nach der Zeit ableitet. Betrachtet man z.B. die x-Komponente der durch auf ausgeübten Beschleunigung, so liefert die Quotientenregel für die entsprechende Komponente des Rucks

Die Ableitung des Abstandsvektors nach der Zeit liefert unmittelbar die Geschwindigkeitsdifferenz. Weniger offensichtlich ist die Ableitung dessen Betrags. Schreibt man diesen aber ausführlich gemäß des Satzes von Pythagoras hin (), so erkennt man, dass das Skalarprodukt durch zweimaliges Nachdifferenzieren innerhalb der Wurzel entsteht.

Elementare Schlussfolgerungen[Bearbeiten]

Auf die schon wiederholt erörterte Kreisbahn soll nun auch der Ruck angewandt werden. Bei unbeweglicher Zentralmasse sind vektorieller Abstand und Geschwindigkeitsdifferenz zwischen und direkt durch die entsprechenden Größen der umlaufenden Masse gegeben. Auf einem kreisförmigen Orbit stehen Orts- und Geschwindigkeitsvektor zudem stets senkrecht aufeinander, so dass das Skalarprodukt der beiden verschwindet. Damit ist der Betrag des auf ausgeübten Rucks einfach durch gegeben. Setzt man noch die Bahngeschwindigkeit ein, so lautet das Endergebnis:

Der Ruck nimmt sehr rasch mit zunehmendem Bahnradius ab - bei doppeltem Radius um einen Faktor (etwa 11.3), bei dreifachem um einen Faktor (ungefähr 46.8) usw.

Die starke Abhängigkeit des Rucks von der Ausdehnung der Bahn lässt sich auch im Alltag nachvollziehen. Bei einer kurvenreichen Fahrt kann nicht nur das häufige Bremsen und Gas geben Schwindelgefühl hervorrufen, sondern auch die heftigen Rucks beim Passieren enger Kurven.

Betrachtet man das Verhältnis von Beschleunigung und Ruck für eine Kreisbahn, so stellt man fest, dass es mit dem Verhältnis von Umlaufgeschwindigkeit und Beschleunigung exakt identisch ist und somit ebenfalls als Maß für die dynamische Zeit dienen kann:


C-Code: Rucks im Mehrkörpersystem

Die Berechnung des in einem Mehrkörpersystem auf einen Massenpunkt ausgeübten Rucks folgt der gleichen Struktur wie diejenige der Beschleunigung. Benötigt werden jetzt aber nicht nur die Positionen r, sondern auch die Geschwindigkeiten v der einzelnen Mitglieder, für welche wie für die Rucks j Zeigerarrays für den Variablentyp Double definiert werden müssen. Zusätzlich zur 3. wird die 5.Potenz d5 des Abstandvektors dr gebraucht, hinzu tritt die vektorielle Geschwindigkeitsdifferenz dv zwischen zwei Körpern. skalar schließlich bezeichnet das aus Abstand und Geschwindigkeitsdifferenz gebildete Skalarprodukt.

Wie die Beschleunigung wird auch der Ruck im Rahmen einer ständig wiederkehrenden Prozedur bestimmt. Zu den bisherigen Übergabeparametern objekt, N, m und r tritt das Vektorarray v hinzu. Da der Ruck nie allein, sondern immer zusammen mit der Beschleunigung betrachtet wird, werden für den aktuell untersuchten Körper beide Größen a[objekt] und j[objekt] zurückgegeben. Der Aufruf der Prozedur geschieht durch ruck(objekt,N,m,r,v,a[objekt],j[objekt]).

/* Globale Variablen */

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

void ruck (unsigned int objekt, unsigned int N, double *m, double **r, double **v, double *a, double *j)
{

/* Lokale Variablen */

  unsigned int i,k;
  double dr[3];
  double d3,d5;
  double dv[3];
  double skalar;

/* Initialisierung von Beschleunigung und Ruck */

  for (k = 0;k < 3;k ++)
  {
    a[k] = 0;
    j[k] = 0;
  }

/* Berechnung von Beschleunigung und Ruck */

  for (i = 0;i < N;i ++)
  {
    if (i != objekt)
    {

/* Abstandsvektor und Geschwindigkeitsdifferenz zwischen einem beliebigen Körper i */
/* und dem untersuchten Massenpunkt objekt                                         */

      for (k = 0;k < 3;k ++)
      {
        dr[k] = r[i][k] - r[objekt][k];
        dv[k] = v[i][k] - v[objekt][k];
      }

/* 3.und 5.Potenz des Betrags des Abstandsvektors */

      d3 = pow (betrag (dr),3);
      d5 = pow (betrag (dr),5);

/* Skalarprodukt von Abstandsvektor und Geschwindigkeitsdifferenz */

      skalar = skalarprodukt (dr,dv);

/* Beitrag der Beschleunigung und des Rucks durch den Körper i */

      for (k = 0;k < 3;k ++)
      {
        a[k] += G * m[i] * dr[k] / d3;
        j[k] += G * m[i] * (dv[k] / d3 - 3 * skalar * dr[k] / d5);
      }
    }
  }
}

Allgemeine Lösungsmethoden[Bearbeiten]

Wie im vorausgegangenen Kapitel diskutiert, hängt die Beschleunigung eines jeden Massenpunkts von dessen eigener und den Positionen aller anderen ab, so dass sämtliche Beschleunigungen untereinander gekoppelt sind. Eine exakte Lösung dieses Problems, welche die Positionen und Geschwindigkeiten aller Massenpunkte für jeden beliebigen Zeitpunkt und somit deren genaue Bahnen angibt, ist nur für ein aus zwei Körpern bestehendes System möglich. Schon ab drei Massenpunkten ist man auf Näherungen angewiesen, die im Kern auf folgenden, bereits von Euler ausgesprochenen Gedanken zurückgehen.

Euler-Verfahren[Bearbeiten]

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 nachfolgenden 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 kleineres 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äziseren 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;
    }
  }
}

Zwischenschritt-Verfahren: Leapfrog und Runge-Kutta[Bearbeiten]

Ein entscheidender Schwachpunkt des Euler-Verfahrens besteht darin, dass man die Momentanbeschleunigungen und -geschwindigkeiten zu einem Zeitpunkt benutzt, um den Zustand des Systems zu einer Zeit zu bestimmen. Es ist naheliegend, dass man den tatsächlichen Verlauf der Geschwindigkeiten und Beschleunigungen während der Zeitspanne 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[Bearbeiten]

Herangehensweise[Bearbeiten]

Der einfachste Weg hin zu Zwischenwerten besteht darin, die Momentangeschwindigkeiten zur Zeit nicht für einen ganzen, sondern nur einen halben Zeitschritt gelten zu lassen. Damit gewinnt man Positionen zum Zeitpunkt . Aus solchen Zwischenpositionen folgen unmittelbar Zwischenwerte für die Momentanbeschleunigungen. Diese lassen genauer auf die Momentangeschwindigkeiten zum Zeitpunkt 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 die für 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 und geltenden Durchschnittsgeschwindigkeit dar als die von dem Euler-Verfahren auschließ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 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[Bearbeiten]

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 , 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 Fehlen 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 proportional, d.h. eine Halbierung der Zeitintervalle erbringt pro Schritt das 16-fache an Genauigkeit. Verringert man allgemein um einen Faktor , steigert man die Präzision eines Einzelschritts um bzw. der gesamten Simulation um . 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 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 durchgeführt als in Gebieten geringer Dichte.

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


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

Die kleine Masse startet erneut von der Position aus. Die Anfangsgeschwindigkeit soll jetzt aber nur für einen halben Zeitschritt gelten. Für die Zwischenposition und ihren Abstand von der großen Masse gilt somit:

Daraus folgt mittels des Gravitationsgesetzes die entsprechende Zwischenbeschleunigung , wobei der immer wieder im Nenner auftauchende etwas unhandliche Ausdruck durch abgekürzt wird.

Um auf die Geschwindigkeit nach dem ganzen Schritt zu schließen, wird obige Beschleunigung für das volle Zeitintervall festgehalten. Das Ergebnis lautet:

Mit der ebenfalls nur für eine Zeit angewandten Endgeschwindigkeit gelangt man von der Zwischen- zur Schlussposition. Für diese gilt:

Nun stehen alle für die Bestimmung der relativen Fehler erforderlichen Ergebnisse bereit. Da wieder angenommen wird, muss man bei der Berechnung des Endabstandes und des Betrags der Endgeschwindigkeit nur die niedrigsten Potenzen von berücksichtigen. Dabei ist zu beachten, dass auch der Term solche Potenzen enthält, denn wegen gilt die Näherung:

Unter Verwendung aller Terme lautet zunächst:

Bei dem quadratischen Glied fällt nach dem Einsetzen des genäherten Wertes von die 1 weg, so dass sich dieses effektiv als ein Term der 4.Potenz mit dem Beitrag erweist. Bei dem nächsten Glied genügt es, mit zu rechnen, da alle weiteren Terme nur höhere Potenzen von liefern. Es lautet so . Das letzte Glied kann einfach gestrichen werden. All diese Überlegungen führen zu:

Die Wurzel kann auf die gleiche Art und Weise genähert werden wie bereits unter dem Euler-Verfahren besprochen, womit sich folgender relativer Fehler ergibt:

Die abermalige Anwendung der Sonnenmasse und des Bahnradius der Erde bestätigt die bereits diskutierte Genauigkeit des Leapfrog-Verfahrens. Setzt man wiederum das 3.Keplersche Gesetz ein, so gewinnt man den einprägsamen Ausdruck:

Setzt man die Fehler von Leapfrog- und Euler-Methode gleich, so erhält man einen kritischen Abstand , unterhalb dessen das Euler-Verfahren anscheinend besser abschneidet.

Setzt man jedoch wiederum in den relativen Fehler ein, so gilt einfach:

Für ist der relative Fehler größer als der Abstand zwischen den beiden Massen selbst! Damit ist der Zusammenbruch beider Verfahren schon lange vor Erreichen des kritischen Abstands gezeigt.

Zuletzt sei auch der relative Fehler der Geschwindigkeit erörtert. Der Betrag der Endgeschwindigkeit lautet unter Mitnahme aller Terme:

Beim quadratischen Glied fällt analog zum Abstand die 1 weg, so dass sich dieses erneut als Term der 4.Potenz entpuppt. Damit ist für den Vorfaktor erlaubt, alle anderen Terme würden nur höhere Potenzen von darstellen. Es ergibt sich so der Beitrag . Für das Glied 4.Potenz gilt wiederum , so dass es den Wert beisteuert. Somit lautet das Endergebnis:

Wie bei dem Euler-Verfahrens weist der Geschwindigkeitsbetrag absolut betrachtet den gleichen relativen Fehler auf wie der Abstand. Doch nun wird die Geschwindigkeit leicht unter- anstatt überschätzt, was angesichts des leicht überschätzten Abstandes eine zusätzliche Stabilisierung der Simulation bedeutet. Die Euler-Methode überschätzt sowohl als auch , was eine Unterschätzung der Gravitation, aber eine Überschätzung der Zentripetalkraft nach sich zieht. Dementsprechend driftet die kleine Masse mit jedem Zeitschritt rasch nach „außen“. Nun aber werden beide Kräfte unterschätzt, so dass sie weiterhin in guter Näherung im Gleichgewicht bleiben. Die Fehler der Einzelschritte weisen im Vergleich zum Euler-Verfahren so eine viel schwächere Tendenz auf, sich aufzusummieren.


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

Für eine effiziente Programmierung ist es hilfreich, wenn man gemäß dem zu Beginn dieses Abschnitts gezeigten Schemas einige Schritte hintereinander aufschreibt. Man erkennt dann, dass das Leapfrog-Verfahren folgendem Algorithmus folgt:


Wechselseitiges Überspringen von Positionen und Geschwindigkeiten bei der Implementierung des Leapfrog-Verfahrens


Zu Beginn der Simulation zum Zeitpunkt muss man sich mit der Geschwindigkeit um einen halben Zeitschritt zur Position bewegen. Mit der dort herrschenden Beschleunigung gewinnt man die Geschwindigkeit , d.h. die Geschwindigkeit zum Zeitpunkt kann man überspringen. Mit gelangt man nicht nur bis zur Positionen , sondern bis nach , kann also überspringen. Die Beschleunigung liefert wiederum die Geschwindigkeit und jene die Position usw. Dieses permanente wechselseitige Überspringen von Positionen und Geschwindigkeiten gibt dem Verfahren seinen Namen, welches im Englischen Bockspringen bedeutet. Für den nachfolgenden Code werden im Vergleich zum Euler-Verfahren keine neuen Variablen und Prozeduren benötigt.

Ein Vergleich der Codes für die beiden Methoden zeigt, dass für das Leapfrog-Verfahren kein höherer Rechenaufwand erforderlich ist. Seine im Vergleich zum Euler-Verfahren viel bessere Genauigkeit verdankt es allein der geschickten Wahl der Zeitpunkte, für welche Positionen, Geschwindigkeiten und Beschleunigungen berechnet werden.

/* 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)
{

/* Initialer Halbschritt zur Position r(dt/2) */

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

/* Beschleunigungen a(dt/2), a(3dt/2) usw. */

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

/* Geschwindigkeiten v(t+dt), v(t+2dt) usw. */
/* Positionen r(t+3dt/2), r(t+5dt/2) usw.   */
 
  for (i = 0;i < N;i ++)
  {
    for (k = 0;k < 3;k ++)
    {
      v[i][k] += a[i][k] * dt;
      r[i][k] += v[i][k] * dt;
    }
  }
}

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[Bearbeiten]

Heun-Verfahren als Vorstufe[Bearbeiten]

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 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 . Die Genauigkeit eines solchen Schritts steigt so nur mit und diejenige einer kompletten Simulation nur mit , wenn man um den Faktor 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[Bearbeiten]

Um eine solche Verbesserung zu erreichen, startet man zunächst wie bei dem Leapfrog-Verfahren mit der Geschwindigkeit und Beschleunigung an der Position zu Beginn eines Zeitintervalls und geht damit einen halben Schritt zur Position , wo neue Werte für Geschwindigkeit und Beschleunigung gelten. Mit diesen startet man nun abermals von der Ausgangsposition aus und gelangt so wiederum mit einem Halbschritt an die Position , wo die Geschwindigkeit und Beschleunigung herrschen. Man durchläuft das Zeitintervall 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 bereit.

Hiermit gelangt man wiederum von der Startposition ausgehend nach einen ganzen Zeitschritt zur Position , welche eine Geschwindigkeit und Beschleunigung 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 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.


Einschub für Fortgeschrittene: Runge-Kutta-Verfahren und Simpsonregel

Die nochmalige (nominelle) Steigerung der Genauigkeit kann man sich folgendermaßen klarmachen. Das Leapfrog-Verfahren verwendet pro Zeitschritt nur zwei Geschwindigkeitswerte (nämlich zu dessen Beginn und Ende), so dass die Fläche unter der Geschwindigkeitskurve durch Trapeze ersetzt wird. Jetzt stehen drei Werte zur Verfügung, nämlich zu Beginn, am Ende und zusätzlich die Zwischenwerte und , welche durch Mittelung zusammengefasst werden (da beide aus dem ersten Halbschritt hervorgehen).


Deutung des klassischen Runge-Kutta-Verfahrens durch die Simpson-Regel


Die Fläche unter der Geschwindigkeitskurve wird jetzt durch Parabel-Stücke angenähert, welche jeweils durch die drei pro Simulationsschritt gegebenen Werte gelegt werden. Der Flächeninhalt und damit die während eines Schritts zurückgelegte Strecke kann mit Hilfe der Simpsonregel. berechnet werden. Diese Formel erklärt die doppelte Gewichtung der Werte 2 und 3 im letzten Schritt des Runge-Kutta-Verfahrens, denn es gilt:

Auch für die Beschleunigung kann man mit Parabel-Flächen arbeiten, wohingegen man sich bei der Leapfrog-Methode mit Rechtecken als Ersatz für die tatsächliche Fläche unter der Beschleunigungskurve begnügen muss. Für die Änderung der Geschwindigkeit gilt:

Die Simpson-Regel begründet auch die Abhängigkeit der Genauigkeit des Runge-Kutta-Verfahrens von . Die Abweichung zwischen der Fläche unter einer Ersatzparabel und derjenigen unter der tatsächlichen Kurve ist direkt proportional zur 5.Potenz des Abstands zwischen zwei Stützstellen.


C-Code: Mehrkörpersimulation mit dem Runge-Kutta-Verfahren

Bedingt durch die zahlreichen Zwischenschritte ist der Code komplexer als für die bisher diskutierten Methoden. Für Position, Geschwindigkeit und Beschleunigung müssen nun jeweils 4 Datensätze vorgehalten werden - ein Datensatz pro Zwischenschritt. Diese Größen werden daher jetzt als dreidimensionale Arrays bzw. als zweidimensionale Arrays von Zeigern deklariert. Nicht nur für jeden Massenpunkt, sondern auch für jeden Zwischenschritt existiert damit ein Zeiger, der wie bisher auf jede der drei Vektorkoordinaten deuten kann. Der erste Index eines jeden Zeigerarrays bezieht sich auf den gerade abgearbeiteten Zwischenschriit (läuft aber von 0 bis 3 anstatt von 1 bis 4, da die Indizes für ein Array gemäß der Syntax von C mit 0 und nicht mit 1 beginnen), der zweite jeweils auf einen bestimmten Körper. Beim Aufruf der Beschleunigungs-Prozedur muss für r und a nun jeweils der Index des aktuellen Zwischenschrittes angegeben werden.

/* 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 a1 an den Positionen r1 */

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

/* Positionen r2 = r1 + v1 * dt/2                             */
/* Geschwindigkeiten v2 = v1 + a1 * dt/2 an den Positionen r2 */
/* Beschleunigungen a2 an den Positionen r2                   */

  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 / 2;
      v[1][i][k] = v[0][i][k] + a[0][i][k] * dt / 2;
    }
  }
  for (i = 0;i < N;i ++)
  beschleunigung (i,N,m,r[1],a[1][i]);

/* Positionen r3 = r1 + v2 * dt/2                             */
/* Geschwindigkeiten v3 = v1 + a2 * dt/2 an den Positionen r3 */
/* Beschleunigungen a3 an den Positionen r3                   */

  for (i = 0;i < N;i ++)
  {
    for (k = 0;k < 3;k ++)
    {
      r[2][i][k] = r[0][i][k] + v[1][i][k] * dt / 2;
      v[2][i][k] = v[0][i][k] + a[1][i][k] * dt / 2;
    }
  }
  for (i = 0;i < N;i ++)
  beschleunigung (i,N,m,r[2],a[2][i]);

/* Positionen r4 = r1 + v3 * dt                             */
/* Geschwindigkeiten v4 = v1 + a3 * dt an den Positionen r4 */
/* Beschleunigungen a4 an den Positionen r4                 */

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

/* Endpositionen = Startpositionen + mittlere Geschwimndigkeiten * dt             */
/* Endgeschwindigkeiten = Startgeschwindigkeiten + mittlere Beschleunigungen * dt */

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

Im Prinzip könnte man die pro Zwischenschritt immer gleiche Abfolge von neuer Position, Geschwindigkeit und Beschleunigung durch eine weitere Prozedur zusammenfassen, doch wird zugunsten einer klareren Darstellung der Zusammenhänge darauf verzichtet.

Da in den nachfolgenden Buchkapiteln immer wieder dreidimensionale Arrays (wie hier für die Position) benutzt werden, soll auch für solche die Speicherverwaltung kurz skizziert werden. Die Speicheradressierung geschieht nun in drei Stufen. Im ersten Schritt wird für jedes erforderliche Zeigerarray (Vektorarray-Objekt) Speicher zur Verfügung gestellt, daher die Anweisung sizeof (double **). Um alle Zwischenschritte des Runge-Kutta-Verfahrens abzudecken, werden insgesamt 4 Zeigerrrays benötigt. Im zweiten Schritt wird innerhalb jeden Zeigerarrays für jedes Vektorobjekt Platz geschaffen (was dem ersten Schritt der Speicheradressierung für ein zweidimensionales Array entspricht). Zuletzt wird für jedes Vektorobjekt Speicher für dessen Komponenten reserviert (wie im Falle zweidimensionaler Arrays).

r = malloc (4 * sizeof (double **));
for (j = 0;j < 4;j ++)
{
  r[j] = malloc (N * sizeof (double *));
  for (i = 0;i < N;i ++)
  r[j][i] = malloc (3 * sizeof (double));
}

Die Speicherfreigabe durchläuft jetzt ebenfalls drei Stufen, wobei wie für zweidimensionale Arrays im Vergleich zur Speicheradressierung die Abfolge sich umkehrt. Abermals wird zuerst der Speicher für die Vektorkomponenten geleert, dann für die Vektorobjekte, zuletzt für die Vektorarray-Objekte.

for (j = 0;j < 4;j ++)
{
  for (i = 0;i < N;i ++)
  free (r[j][i]);
  free (r[j]);
}
free (r);

Trotz seiner großen allgemeinen Bedeutung und Genauigkeit wird das Runge-Kutta-Verfahren nur relativ selten auf astronomische Mehrkörpersysteme angewandt, da es sich aufgrund seines komplexen Vorgehens bei der Bildung von Zwischenschritten kaum mit dynamischen Zeitintervallen vereinbaren lässt. Es kann bei Aufgabenstellungen eingesetzt werden, bei denen konstant gehalten werden darf. Als Beispiel sei eine Diplomarbeit von Geisler (2006) [6] genannt, welcher untersuchte, wie die ausgedehnten Gashüllen sehr leuchtkräftiger Sterne durch die von einem Begleitstern ausgehende Schwerkraft gestört werden. Die Gashülle wird durch Massenpunkte jeweils gleicher Masse dargestellt. Da im Vergleich zu den Sternmassen sehr gering ist, üben diese Modellkörper untereinander nur sehr kleine Beschleunigungen aufeinander aus. Die Sterne werden nicht als punktförmig behandelt, sondern deren endlicher Radius wird berücksichtigt, so dass ein verschwindender Abstand zwischen einem Mitglied der Gashülle und einem Stern ausgeschlossen ist. Die Simulation bleibt somit auch ohne dynamische Zeitschritte ausreichend stabil, was man anhand der Energieerhaltung überprüfen kann.


Prädiktor-Korrektor-Verfahren: Leapfrog und Hermite-Polynome[Bearbeiten]

Idee[Bearbeiten]

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

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:

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 und derjenige einer ganzen Simulation proportional zu . Mit Ruck ergeben sich Proportionalitäten zu bzw. .

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

Leapfrog-Verfahrens[Bearbeiten]

Alternative Formulierung[Bearbeiten]

Wie Hut und andere (1995) [7] 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:

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:

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 viemehr 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 . 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[Bearbeiten]

Klassischer Korrektor[Bearbeiten]

Um die auf während eines Zwitintervalls 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 und tragen. Unter dem mittleren Knall versteht man die Änderung des Rucks , unter dem mittleren Knistern die Änderung des Knalls im Verlauf von .

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

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:

Die 4 Formeln für die gesuchten 4 neuen Werte für , , und bilden ein vollständiges lineares Gleichungssystem, wodurch man und eliminieren kann. Das Ergebnis wirkt erstaunlich elegant:

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 jewels einen auf der Beschleunigung bzw. dem Ruck beruhende Zusatztern. Dieser sorgt dafür, dass der relative Fehler eines Simulationsschritts wie für das Runge-Kutta-Verfahren zu proportional ist, und nicht wie für die Leapfrog-Methode nur proportional zu .


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:


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[Bearbeiten]

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

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:

Für den Knall und das Knistern geben Makino und Aarseth (1992) [8] folgende Ausdrücke an:


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];
    }
  }
}