Das Mehrkörperproblem in der Astronomie/ Grundlagen/ Das Newtonsche Gravitationsgesetz

Kernaussagen

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

${\vec {F}}_{ji}={\frac {Gm_{i}m_{j}}{r_{ji}^{2}}}{\vec {e}}_{ji}$

$G$ bezeichnet dabei die Gravitationskonstante mit dem Wert 6.67384 10^-11 m³ kg^-1 s^-2. ${\vec {e}}_{ji}$ ist der Einheitsvektor auf der Verbindungslinie der beiden Massen, er zeigt von $i$ nach $j$.

In der Praxis interessiert nicht die auf einen Massenpunkt einwirkende Kraft, sondern nur dessen Beschleunigung. Aus dem Newtonschen Kraftgesetz folgt, dass die Masse $m_{i}$ folgende Beschleunigung ${\vec {a}}_{ji}$ durch die Masse $m_{j}$ erfährt:

${\vec {a}}_{ji}={\frac {Gm_{j}}{r_{ji}^{2}}}{\vec {e}}_{ji}$

Entsprechend erleidet die Masse $m_{j}$ eine Beschleunigung ${\vec {a}}_{ij}$ durch die Masse $m_{i}$ in entgegengesetzter Richtung.

${\vec {a}}_{ij}=-{\frac {Gm_{i}}{r_{ji}^{2}}}{\vec {e}}_{ji}$

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 zunehmendem 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 $m_{i}$ sei so durch einen Vektor ${{\vec {r}}_{i}}=(x_{i},y_{i},z_{i})$ und diejenige der Masse $m_{j}$ durch einen solchen ${{\vec {r}}_{j}}=(x_{j},y_{j},z_{j})$ 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 ${{\vec {r}}_{ji}}$ an. In kartesischen Koordinaten lautet dieser:

${{\vec {r}}_{ji}}={\begin{pmatrix}x_{j}-x_{i}\\y_{j}-y_{i}\\z_{j}-z_{i}\end{pmatrix}}$

Den Einheitsvektor ${{\vec {e}}_{ji}}$ gewinnt man, indem man den Abstandsvektor mit seinem eigenen Betrag dividiert:

${{\vec {e}}_{ji}}={\frac {1}{|{{\vec {r}}_{ji}}|}}{\begin{pmatrix}x_{j}-x_{i}\\y_{j}-y_{i}\\z_{j}-z_{i}\end{pmatrix}}$

Einsetzen in die oben gegebenen Beschleunigungen liefert:

${{\vec {a}}_{ji}}={\frac {Gm_{j}}{|{{\vec {r}}_{ji}}|^{3}}}{\begin{pmatrix}x_{j}-x_{i}\\y_{j}-y_{i}\\z_{j}-z_{i}\end{pmatrix}}$

${{\vec {a}}_{ij}}=-{\frac {Gm_{i}}{|{{\vec {r}}_{ji}}|^{3}}}{\begin{pmatrix}x_{j}-x_{i}\\y_{j}-y_{i}\\z_{j}-z_{i}\end{pmatrix}}$

Der Betrag des Abstandsvektors folgt aus dem Satz des Pythagoras:

$|{{\vec {r}}_{ji}}|={\sqrt {(x_{j}-x_{i})^{2}+(y_{j}-y_{i})^{2}+(z_{j}-z_{i})^{2}}}$

Elementare Schlussfolgerungen

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 $m$ betrachtet werden, welche um eine große Masse $M$ auf einer Kreisbahn mit Radius $r$ umläuft.

Dabei sei zuerst die Bahngeschwindigkeit $v_{B}$ betrachtet. Man erhält sie durch Gleichsetzen der Gravitationsbeschleunigung $GM/r^{2}$ mit der Zentripetalbeschleunigung $v^{2}/r$:

$v_{B}={\sqrt {\frac {GM}{r}}}$

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 $v_{B}$ 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 $2\pi r$ der Bahn folgt die Umlaufdauer $T$. Es gilt $T=2\pi r/v_{B}$, woraus durch Einsetzen und Quadrieren folgt:

$T^{2}={\frac {4\pi ^{2}}{GM}}r^{3}$

Das Quadrat der Umlaufdauer 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 Umlaufdauer um den Faktor ${\sqrt {8}}$ entsprechend etwa 2.8. Bei dreifacher Entfernung nimmt $T$ um das ${\sqrt {27}}$-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 Umlaufdauern verglichen werden (eine Gegenüberstellung der Absolutentfernungen zur Sonne und der absoluten Werte der Umlaufdauern 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 Umlaufdauer, Bahngeschwindigkeit und Gravitationsbeschleunigung $a_{B}$ nochmals genau an, so stellt man fest, dass sie folgenden Zusammenhang erfüllen:

$T=2\pi {\frac {v_{B}}{a_{B}}}$

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 Anweisung. 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.