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