Astronomische Berechnungen für Amateure/ Himmelsmechanik/ Lösungsverfahren

Wie wir festgehalten haben ist die Keplergleichung eine transzendente Gleichung. Dies bedeutet, dass es abgesehen von trivialen Fällen keine analytische Lösung in der Form E = f(e,M) gibt, sondern nur numerische Lösungen für konkrete Werte. Die trivialen Ausnahmen betreffen die Fälle, wo M ein ganzzahliges Vielfaches von π ist. Ist nämlich M = k∙π (k = …, –2, –1, 0, 1, 2, …), dann ist E = M, denn sin(k∙π) = 0, und daraus folgt v = E. Anschaulich bedeutet dies, dass sowohl der wahre als auch der auf den Hilfskreis projizierte und der mittlere Planet gleichzeitig durch Perihel und Aphel gehen.


In allen anderen Fällen muss die Keplergleichung durch ein numerisches Verfahren genähert gelöst werden. Dabei stehen im Prinzip zwei grundsätzlich verschiedene Verfahren mit Untervarianten zur Auswahl:

  • Iterationsverfahren
  • Reihenentwicklungen

Wir beschränken uns im folgenden auf die erste Variante. Als Iterationsverfahren bezeichnet man ein Verfahren, bei dem schrittweise aus berechneten Werten neue Werte solange berechnet werden, bis sich zwei aufeinander folgende Werte um weniger als eine vorgegebene Schranke unterscheiden. Die Untervarianten unterscheiden sich gelegentlich auch durch die Wahl des Startwertes.

Das erste Verfahren beginnt mit E0 = M als Startwert. Man berechnet dann schrittweise:


Die Iteration wird abgebrochen, wenn die Differenz Δn+1 kleiner als eine vorgegebene Schranke ist. Will man z.B. E auf 1" = 1/3600° = 4.848 ∙ 10–6 genau bestimmen, dann wird die Iteration abgebrochen, wenn die Differenz 4 ∙ 10–6 oder sogar 1 ∙ 10–6 [1] unterschreitet. Wird die Schranke kleiner gewählt, dann wird das Resultat genauer, benötigt aber mehr Schritte, bis es erreicht wird. Solange nur ein einzelner Wert berechnet werden soll, ist das bei der Rechenleistung der heutigen Computer kein Problem. Sollen aber z.B. Ephemeriden für 50 oder 100 Asteroiden für ein ganzes Jahr im Tagesabstand berechnet werden, dann kommen ordentlich viele Rechnungen zusammen und können einen schwächeren PC schon ganz schön fordern.

Für kleine numerische Exzentrizitäten e konvergiert[2] das vorgestellte Verfahren recht gut. Ab Werten von e, die grösser als 0.4 oder 0.5 sind, ist das Konvergenzverhalten sehr langsam. Es kann dann schneller sein, einen komplizierteren Ausdruck für die nächste Näherung zu berechnen, dafür deutlich weniger Iterationsschleifen zu durchlaufen. Ein solches Verfahren addiert zu einem Näherungswert En einen Korrekturterm, der mit M, e und En errechnet wird. Das Verfahren wird wiederum so lange wiederholt, bis die Differenz zweier aufeinander folgender Werte – also eigentlich der Korrekturterm – eine vorgegebene Schranke unterschreitet. Man berechnet darum schrittweise[3]



Beispiele:

Wir lösen zunächst die Keplergleichung nach der ersten Methode für M = 15° = 0.261 7994 und e = 0.0934[4] mit der Schranke 0.000 0001. Wir erhalten dann die Reihe:

E0 = 0.261 7994
E1 = 0.285 9731     Δ1 = 0.024 1737
E2 = 0.288 1467     Δ2 = 0.002 1736
E3 = 0.288 3414     Δ3 = 0.000 1947
E4 = 0.288 3589     Δ4 = 0.000 0174
E5 = 0.288 3604     Δ5 = 0.000 0016
E6 = 0.288 3606     Δ6 = 0.000 0001
E7 = 0.288 3606     Δ7 = 0.000 0000

Die Tabellenkalkulation OOo Calc zeigt an, dass mit der siebten Iteration die Genauigkeitsschranke unterschritten wird. Damit haben wir das Resultat gefunden:

Wenn M = 15°, dann ist E = 16.521 844° und v = 18.118 566° (e = 0.0934).

Bei gleicher Exzentrizität sind es immer etwa 7 Iterationsschritte, bis die Schranke unterschritten wird, nahezu unabhängig davon, wie gross M ist. Wird die Schranke um den Faktor 10 verkleinert, sind es 8 Iterationsschritte, wird sie um den Faktor 10 vergrössert, sind es 6 Iterationsschritte.

Betrachten wir nun aber M = 15° und e = 0.967[5]. Es sind dann 21 Iterationsschritte notwendig, bis die Schranke 0.000 0001 unterschritten wird:

E0 = 0.261 7994
E1 = 0.512 0774     Δ1 = 0.250 2780
E2 = 0.735 6190     Δ2 = 0.223 5416
E3 = 0.910 7011     Δ3 = 0.175 0821
E4 = 1.025 6654     Δ4 = 0.114 9643
E5 = 1.088 6419     Δ5 = 0.062 9764
E6 =   …

Nach dem 21. Iterationsschritt findet man als Resultat: E = 65.360 217° und daraus v = 157.169 691°.

Bei M = 175° und e = 0.967 wird aber die Schranke erst nach der Iteration mit der Nummer 396 unterschritten. Das Ergebnis ist E = 177.457 645° und daraus v = 179.670 647°.


Setzen wir dafür aber das zweite Iterationsverfahren ein, so wird die Schranke bereits nach dem 3. Iterationsschritt erreicht:

E0 = 3.054 3262
E1 = 3.097 2533      Δ1 = 042 9271
E2 = 3.097 2202      Δ2 = 0.000 0330
E3 = 3.097 2202      Δ3 = 0.000 0000

Das zweite Verfahren konvergiert also sehr schnell zur Lösung E = 177.457 649° und daraus v = 179.670 648°. Es hat allerdings zwei Nachteile und sollte darum nicht einfach bedenkenlos eingesetzt werden: zum einen ist es mit zwei zu berechnenden Winkelfunktionen und einem Bruch für den Computer deutlich aufwendiger als das erste Verfahren. Es kann damit seine Vorteile erst dann ausspielen, wenn das erste Verfahren zu einer sehr grossen Anzahl Iterationsschritte führt. Dies ist bei Exzentrizitäten nahe 1 der Fall. Zum anderen kann ein neues Problem auftreten: statt stetig zum Lösungswert zu konvergieren, können die Werte anfänglich oszillieren. Dies ist speziell bei Exzentrizitäten ganz nahe bei 1 der Fall.


Beispiel:

Nehmen wir nochmals e = 0.967, aber M = 5°. Dann wird das Ergebnis E = 42.258 779° mit dem 8. Iterationsschritt erreicht. Der erste Iterationsschritt liefert allerdings einen völlig unsinnigen Wert von E1 = 2.384…rad = 136.649 441° (!):

E0 = 0.087 2665
E1 = 2.384 9827     Δ1 = 2.297 7162
E2 = 1.425 6490     Δ2 = 0.959 3337
E3 = 0.982 0547     Δ3 = 0.443 5943
E4 = 0.786 3955     Δ4 = 0.195 6591
E5 = 0.740 0885     Δ5 = 0.046 3071
E6 = 0.737 5622     Δ6 = 0.002 5263
E7 = 0.737 5548     Δ7 = 0.000 0073
E8 = 0.737 5548     Δ8 = 0.000 0000

Dieses Verhalten wird besonders deutlich bei grossen Werten für e:


Beispiel:

Wählen wir e = 0.999 und M = 7°, dann finden wir das Resultat nach Iterationsschritt 29 – aber sehen Sie sich die ersten Iterationswerte an. Ausnahmsweise sind auch die Zwischenwerte in Grad umgerechnet:

E0 = 0.122 7.000° ... ...
E1 = 14.536 832.869° Δ1 = 14.414
E2 = 4.816 275.955° Δ2 = 9.720
E3 = -1.529 -87.611° Δ3 = 6.345
E4 = -0.848 -48.562° Δ4 = 0.682
E5 = -0.196 -11.225° Δ5 = 0.652
E6 = 5.951 340.963° Δ6 = 6.147
E7 = -104.664 -5'996.812° Δ7 = 110.615
E8 = -36.381 -2'084.498° Δ8 = 68.283
E9 = 13.586 778.411° Δ9 = 49.967
E10 = -12.872 -737.536° Δ10 = 26.458
E11 = 254.789 14'598.350° Δ11 = 267.662
E12 = 123.909 7'099.442° Δ12 = 130.881
E13 = 18.444 1'056.786° Δ13 = 105.464
E14 = -210.125 -12'039.297° Δ14 = 228.570
E15 = -101.624 -5'822.647° Δ15 = 108.501
E16 = 84.813 4'859.434° Δ16 = 186.437
E17 = 42.450 2'432.227° Δ17 = 42.363
E18 = -2.626 -150.451° Δ18 = 45.076
E19 = -1.419 -81.313° Δ19 = 1.207
E20 = -0.767 -43.943° Δ20 = 0.652
E21 = -0.069 -3.961° Δ21 = 0.698
E22 = 36.049 2'065.432° Δ22 = 36.118
E23 = 1.847 105.842° Δ23 = 34.201
E24 = 1.247 71.445° Δ24 = 0.600
E25 = 0.986 56.518° Δ25 = 0.261
E26 = 0.917 52.558° Δ26 = 0.069
E27 = 0.912 52.272° Δ27 = 0.005
E28 = 0.912 52.270° Δ28 = 0.000
E29 = 0.912 52.270° Δ29 = 0.000

Das Verhalten hängt damit zusammen, dass der Nenner nahezu Null werden kann. Mit einer einfachen Massnahme kann man das Konvergenzverhalten stark beeinflussen: man wählt als Startwert fix E0 = π. Dann ist bereits nach dem 7. Iterationsschritt das Ergebnis erreicht, und es gibt keine oszillierenden Werte.

Zum Schluss sei angemerkt, dass die Ergebnisse stark von der eingesetzten Soft- und Hardware abhängen, besonders davon, mit wievielen Stellen das Programm intern rechnet. Es ist ohne weiteres möglich, dass Sie mit Ihrer Software auf Ihrem Rechner andere Zwischenergebnisse erhalten – das Schlussresultat sollte davon allerdings unbeeinflusst bleiben. Ein guter Rat: testen Sie die Implementation des Algorithmus sorgfältig, bevor Sie Ihre Software produktiv einsetzen!



Übungen

  • Im voran stehenden Kapitel haben Sie die Merkurpositionen für eine heliozentrische Kreisbahn berechnet. In diesem Kapitel sollen Sie nun die Merkurpositionen unter zwei anderen Bedingungen rechnen: 1. Merkur läuft gleichmässig auf einer exzentrischen Kreisbahn, wobei Kreis- und Ellipsenbahn Perihel- und Aphelposition gemeinsam haben. Der Kreis entspricht dem Hilfskreis, den wir zur Herleitung der Keplergleichung gezeichnet haben. Der zu berechnende Winkel entspricht M. 2. Merkur läuft auf einer elliptischen Bahn mit a = 0.387 AE und e = 0.2056. Für diese Aufgabe müssen Sie die Keplergleichung lösen und die wahre Anomalie v berechnen. Berechnen Sie wiederum die x-y-Koordinaten für ein System mit dem Ursprung in der Sonne. Vergleichen Sie die Lösungen der drei Aufgaben.
  • Um die Bedeutung des 2. Keplergesetzes etwas besser zu verstehen, berechnen Sie die wahre Anomalie v und den Abstand r von der Sonne für den Halley'schen Kometen, wenn M = 0° / 45° / 90° / 120° / 150° / 170° / 175° / 179° misst. Die Bahn des Halley'schen Kometen hat eine grosse Halbachse von a = 17.834 AE und eine numerische Exzentrizität von e = 0.967.



Nachweise:

  1. Was faktisch eine Genauigkeit von besser als 0.2" bedeuten würde.
  2. Konvergenz bedeutet in diesem Zusammenhang: ab einem bestimmten Wert m werden die Differenzen Δm immer kleiner, bis sie schliesslich die vorgegebene Schranke unterschreiten.
  3. Für Mathematiker sei angemerkt, dass es sich hierbei um das Newtonsche Verfahren zur Bestimmung von Nullstellen transzendenter, differenzierbarer Funktionen handelt.
  4. e = 0.0934 ist die numerische Exzentrizität der Marsbahn.
  5. e = 0.967 ist die numerische Exzentrizität für die Bahn des Halley'schen Kometen.