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 E₀ = M als Startwert. Man berechnet dann schrittweise:

$E_{0}\quad =\quad M$

$E_{1}\quad =\quad M\;+\;e\cdot \sin(E_{0})\;{\text{;}}\quad \left|E_{1}-E_{0}\right|\;=\;\Delta _{1}$

$E_{2}\quad =\quad M\;+\;e\cdot \sin(E_{1})\;{\text{;}}\quad \left|E_{2}-E_{1}\right|\;=\;\Delta _{2}$

$E_{3}\quad =\quad M\;+\;e\cdot \sin(E_{2})\;{\text{;}}\quad \left|E_{3}-E_{2}\right|\;=\;\Delta _{3}$

$\ldots \qquad \quad \ldots$

$E_{n+1}\quad =\quad M\;+\;e\cdot \sin(E_{n})\;{\text{;}}\quad \left|E_{n+1}-E_{n}\right|\;=\;\Delta _{n+1}$

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 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 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 E_n einen Korrekturterm, der mit M, e und E_n 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

$E_{0}\quad =\quad M$

$E_{1}\quad =\quad E_{0}\;+\quad {\frac {M+e\cdot \sin(E_{0})-E_{0}}{1-e\cdot \cos(E_{0})}}\;{\text{;}}\quad \left|E_{1}-E_{0}\right|\;=\;\Delta _{1}$

$E_{2}\quad =\quad E_{1}\;+\quad {\frac {M+e\cdot \sin(E_{1})-E_{1}}{1-e\cdot \cos(E_{1})}}\;{\text{;}}\quad \left|E_{2}-E_{1}\right|\;=\;\Delta _{2}$

$\ldots \qquad \quad \ldots$

$E_{n+1}\quad =\quad E_{n}\;+\quad {\frac {M+e\cdot \sin(E_{n})-E_{n}}{1-e\cdot \cos(E_{n})}}\;{\text{;}}\quad \left|E_{n+1}-E_{n}\right|\;=\;\Delta _{n+1}$

Beispiele:

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

E₀ = 0.261 7994

E₁ = 0.285 9731     Δ₁ = 0.024 1737

E₂ = 0.288 1467     Δ₂ = 0.002 1736

E₃ = 0.288 3414     Δ₃ = 0.000 1947

E₄ = 0.288 3589     Δ₄ = 0.000 0174

E₅ = 0.288 3604     Δ₅ = 0.000 0016

E₆ = 0.288 3606     Δ₆ = 0.000 0001

E₇ = 0.288 3606     Δ₇ = 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. Es sind dann 21 Iterationsschritte notwendig, bis die Schranke 0.000 0001 unterschritten wird:

E₀ = 0.261 7994

E₁ = 0.512 0774     Δ₁ = 0.250 2780

E₂ = 0.735 6190     Δ₂ = 0.223 5416

E₃ = 0.910 7011     Δ₃ = 0.175 0821

E₄ = 1.025 6654     Δ₄ = 0.114 9643

E₅ = 1.088 6419     Δ₅ = 0.062 9764

E₆ =   …

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:

E₀ = 3.054 3262

E₁ = 3.097 2533      Δ₁ = 042 9271

E₂ = 3.097 2202      Δ₂ = 0.000 0330

E₃ = 3.097 2202      Δ₃ = 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 E₁ = 2.384…rad = 136.649 441° (!):

E₀ = 0.087 2665

E₁ = 2.384 9827     Δ₁ = 2.297 7162

E₂ = 1.425 6490     Δ₂ = 0.959 3337

E₃ = 0.982 0547     Δ₃ = 0.443 5943

E₄ = 0.786 3955     Δ₄ = 0.195 6591

E₅ = 0.740 0885     Δ₅ = 0.046 3071

E₆ = 0.737 5622     Δ₆ = 0.002 5263

E₇ = 0.737 5548     Δ₇ = 0.000 0073

E₈ = 0.737 5548     Δ₈ = 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 E₀ = π. 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: