EphemPedia

Anhänger der rechnenden Astronomie

Benutzer-Werkzeuge

Webseiten-Werkzeuge


loesung_der_keplergleichung

Lösung der Keplergleichung

Zur Berechnung der Position eines Himmelskörpers (Planet, Planetoid, periodischer Komet, usw.) zu einem gegebenen Zeitpunkt auf seiner elliptischen Bahn gibt es verschiedene Methoden:

  1. Numerische Integration: Dies liefert die genauesten Daten, diese Berechnungen bestehen aber aus tausenden zu summierenden Termen, weshalb die Berechnung sehr schnell aufwändig wird.
  2. Eine gekürzte Version von Punkt 1, bei der nur eine begrenzte Anzahl an Termen berücksichtigt wird. Die Berechnugnen werden dadurch schneller, aber die Genauigkeit ist nicht mehr ganz so hoch. Für die meisten Anwendungen in der Praxis ist das völlig ausreichend.
  3. Berechnen der heliozentrischen Positionen aus den mittleren Bahnelementen des Körpers (Keplerbahnen). Dies ist zwar die ungenaueste Methode, liefert aber sehr schnelle Ergebnisse und ist für manche Anwendungen genau genug. Man könnte z.B. eine Grafik mit den aktuellen Planetenpositonen erzeugen, wobei natürlich nicht auf Bogensekunden genau gerechnet werden muss, da dies in der Grafik gar nicht darstellbar ist.

In diesem Kapitel geht es um den dritten Fall. Es soll die wahre Anomalie $\nu$ des Objekts berechnet werden. Dies kann entweder durch Lösung der Keplergleichung erreicht werden oder – wenn die Bahnexzentrizität $\epsilon$ nicht zu groß ist – durch eine Reihenentwicklung (Mittelpunktsgleichung) oder eine einfache Näherungsformel.

Begriffe und Zusammenhänge

In Abb. 1 ist die Hälfte einer elliptischen Bahn dargestellt (Ellipsenbogen $PKA$). Die Sonne steht in einem Brennpunkt der Ellipse $S$. Der andere, leere Brennpunkt der Ellipse heißt $F$. Zur Geometrie der Ellipse siehe hier.

Abb. 1: Abmessungen für die Keplergleichung

$S\dots$ Sonne (Brennpunkt)
$P\dots$ Perihel (Sonnennähe)
$A\dots$ Aphel (Sonnenferne)
$K\dots$ umlaufender Körper (Planet, Komet,…)
$a = \overline{ZP}\dots$ große Halbachse
$e = \overline{ZS}\dots$ lineare Exzentrizität
$b\dots$ kleine Halbachse, $b^2 = a^2 - e^2$
$\epsilon = \frac{\overline{ZS}}{\overline{ZP}} = \frac{e}{a}\dots$ numerische Exzentrizität
$\overline{SP} = q = a\cdot(1 - \epsilon)\dots$ Perihelabstand
$\overline{SK} = r \dots$ Radiusvektor (Abstand Sonne-Körper)
$\overline{SA} = Q = a\cdot(1 + \epsilon)\dots$ Aphelabstand
$\overline{PA} = 2\cdot a = q + Q\dots$ große Achse

Die Strecke $\overline{PA} $ ist die große Achse der Bahn. Das Zentrum $Z$ der Ellipse liegt sowohl exakt in der Mitte zwischen dem Perihel $P$ und dem Aphel $A$ als auch in der Mitte zwischen den Brennpunkten.

$K$ (für Körper) ist die Position des Objekts zu einem gegebenen Zeitpunkt. Die Strecke $\overline{SK}$ ist der Radiusvektor $r$ des Objekts zu diesem Zeitpunkt; diese Entfernung wird in Astronomischen Einheiten $AU$ angegeben. Die wahre Anomalie $\nu$ zu diesem Zeitpunkt ist der Winkel zwischen $\overline{SP}$ und $\overline{SK}$; das ist der von dem Objekt seit dem Durchlaufen des Perihels $ P $ beschriebene Winkel.

Die Halbachse $\overline{ZP}$ wird im allgemeinen mit $a$ bezeichnet und auch in Astronomischen Einheiten angegeben. Die numerische Exzentrizität $\epsilon$ der Bahn ist definiert als das Verhältnis der Strecke $\overline{ZS}$ (= lineare Exzentrizität) zur Strecke $\overline{ZP}$:

$$\epsilon = \frac{\overline{ZS}}{\overline{ZP}} = \frac{e}{a}\tag{1}$$

Diese Größe ist dimensionslos und liegt bei einer Ellipse immer zwischen 0 und 1. Die Exzentrizität $\epsilon = 0$ wäre bei einem Kreis gegeben. Der Perihelabstand und Aphelabstand werden mit $ q $ bzw. $ Q $ bezeichnet. Befindet sich der Körper im Perihel, ist $ \nu = 0^{\circ} $ und $ r = q $, während im Aphel $ \nu = 180^{\circ} $ ist und $ r = Q $.

Nun betrachtet man in Abb. 2 einen fiktiven Planeten oder Kometen $K'$, der eine Kreisbahn um die Sonne beschreibt, und sich mit konstanter Geschwindigkeit bewegt und die gleiche Umlaufzeit hat, wie der reale Planet oder Komet $K$ auf der elliptischen Bahn. Darüber hinaus wird angesetzt, dass sich das fiktive Objekt im Punkt $P'$ befindet, wenn der reale Körper im Perihel $P$ steht. Als Hilfskreis wird der Scheitelkreis der Hauptachse $a$ genommen, auf dem sich der fiktive Planet bewege. Einige Zeit später, wenn sich der wahre Körper in $K$ befindet, ist das fiktive Objekt bis nach $K'$ gelaufen. Wie bereits angemerkt ist der Winkel $\nu = \angle PSK$ die wahre Anomalie $\nu$ des Objekts zum gegebenen Zeitpunkt. Zur gleichen Zeit durchläuft $K'$ den Winkel $\angle P'ZK'$, und dies ist die mittlere Anomalie, die im allgemeinen mit $M$ bezeichnet wird.

Abb. 2: Maßgebliche Winkel für die Keplergleichung

Die mittlere Anomalie $M$ ist also der Winkelabstand zum Perihel, den der Planet haben würde, wenn er mit konstanter Geschwindigkeit die Sonne umliefe.

Definitionsgemäß wächst der Winkel $M$ linear (gleichförmig) mit der Zeit. Der Wert von $M$ ist leicht zu finden, da im Perihel $M = 0^{\circ}$ gilt und der Winkel während eines vollständigen Umlaufes um exakt $360^{\circ}$ zunimmt.

Das eigentliche Problem besteht aber darin, die wahre Anomalie $\nu$ zu finden, wenn die mittlere Anomalie $M$ und die Bahnexzentrizität $\epsilon$ bekannt sind. Falls man nicht auf eine Reihenentwicklung zurückgreifen will, ist man gezwungen, die Keplergleichung zu lösen. Hierzu ist es nötig, einen Hilfswinkel $E$, die exzentrische Anomalie, einzuführen. Deren geometrische Definition ist in Abb. 2 angegeben. Der äußere gestrichelte Hilfskreis hat den Durchmesser $2a = \overline{AP}$ der großen Achse. Man zeichnet $\overline{HK}$ rechtwinklig zu $\overline{AP}$. Den Winkel $\angle PZH$ nennt man die exzentrische Anomalie $E$.

Steht der Planet im Perihel, dann sind die Winkel $\nu = E = M = 0^{\circ}$. In der Nähe des Perihels hat der wahre Planet eine größere Geschwindigkeit als der mittlere, fiktive Planet. Dies ist dem 2. Keplerschen Gesetz geschuldet, wonach ein Körper in gleichen Zeiten mit dem Radiusvektor $r$ dieselben Flächen überfährt, und sich der Körper deshalb in Perihelnähe schneller bewegt als in Aphelnähe. Deshalb gilt zwischen Perihel und Aphel (wenn sich der Planet von der Sonne entfernt)

$$\nu \gt M$$

Damit gilt

$$0^{\circ} \lt M \lt E \lt \nu \lt 180^{\circ}$$

Im Aphel ist $\nu = E = M = 180^{\circ}$. Nach dem Durchschreiten des Aphels auf dem Weg zurück zum Perihel bleibt der wahre Planet hinter dem mittleren Planeten zurück. Für diesen Abschnitt gilt dann

$$\nu \lt M$$

Den Radiusvektor $r$, also die Entfernung des Körpers von der Sonne, erhält man über eine der folgenden Beziehungen:

$$r = a\cdot (1 - \epsilon\cdot \cos E)\tag{2}$$ oder

$$r = \large \frac{{a\cdot (1 - {\epsilon^2})}}{{1 +\epsilon\cdot \cos \nu }}\tag{3}$$ oder

$$r = \large \frac{{q\cdot (1 - \epsilon)}}{{1 +\epsilon\cdot \cos \nu }}\tag{4}$$

Methoden zur Bestimmung der wahren Anomalie

Wenn nun die exzentrische Anomalie $E$ bekannt wäre, dann kann man daraus die wahre Anomalie $\nu$ berechnen mit der Barkerschen Gleichung

\[\tan\left(\frac{\nu }{2}\right) = \sqrt{\frac{{1 + \epsilon}}{{1 - \epsilon}}} \cdot \tan\left(\frac{E}{2}\right)\tag{5}\]

Die Keplergleichung lautet:

\[\Large E = M + \epsilon\cdot \sin E \tag{6}\]

Diese Gleichung muss nach $E$ aufgelöst werden. Es handelt sich hier aber um eine transzendente Gleichung, die nicht direkt gelöst werden kann, weil $E$ einmal alleine auftaucht und auch als Argument der Sinusfunktion. Im folgenden sollen nun 4 Methoden zur Ermittelung von $E$ bzw. $\nu$ beschrieben werden.

Methode 1

In die Formel der Keplergleichung müssen die Winkel $M$ und $E$ in Bogenmaß eingesetzt werden. In den Programmiersprachen werden die Argumente von trigonometrischen Funktionen ohnehin im Bogenmaß erwartet. Soll die Berechnung aber im Gradmodus durchgeführt werden, dann muss $\epsilon$ mit dem Umrechnungsfaktor $\tfrac{180}{\pi}$ vom Bogenmaß in Grad, multipliziert werden. Man nennt die umgerechnete Größe $\epsilon_0$ modifizierte Exzentrizität. Die Keplergleichung lautet dann

\[ E = M + \epsilon_0\cdot \sin E \tag{7}\]

$\epsilon_0 = \epsilon\cdot \tfrac{180}{\pi}\dots$ modifizierte Exzentrizität
$E, M\dots$ Winkel in Grad

Damit kann man nun im Gradmaß rechnen. Zur Lösung der Gleichung setzt man auf der rechten Seite einen genäherten Wert für $E$ ein. Die Formel ergibt dann einen besseren Wert für $E$. Dieser Vorgang ist so oft zu wiederholen, bis die geforderte Genauigkeit erreicht ist; in einem Programm kann das automatisch durch eine Schleife erfolgen. Für den ersten Schritt kann man z.B. $E = M$ als Startwert setzen.

Man hat also die Rechenschritte

\[\begin{align} E_0 &= M\\ E_1 &= M + e_0\cdot \sin E_0\\ E_2 &= M + e_0\cdot \sin E_1\\ E_3 &= M + e_0\cdot \sin E_2\\ &\textsf{... usw.} \end{align}\]

$E_1, E_2, E_3, \dots$ sind aufeinanderfolgende bessere Näherungen der exzentrischen Anomalie $E$. Diese schrittweise Annäherung an den korrekten Wert nennt man Iteration.

Beispiel

Man löse die Keplergleichung für $\epsilon = 0.1$ und $M = 5^{\circ}$ mit einer Geanauigkeit von $0\overset{\circ}{.}000001$ (6. Kommastelle) mittels Iteration.


Man findet $\epsilon_0 = \tfrac{180}{\pi}\cdot \epsilon = 5\overset{\circ}{.}72957795$.

Die Keplergleichung lautet damit

$E = 5 + 5\overset{\circ}{.}72957795\cdot \sin E$

Alle Größen sind hier in Grad gegeben. Beginnend mit $M = E = 5^{\circ}$ erhält man schrittweise die Werte

\(\begin{align} &5.499366 \Rightarrow \textsf{erster Wert}\\ &5.549093\\ &5.554042\\ &5.554535\\ &5.554584\\ &5.554589\\ &5.554589 \Rightarrow \textsf{keine Änderung in der 6. Kommastelle mehr} \end{align}\)

Der gesuchte Wert ist also $E = 5\overset{\circ}{.}554589$. Es waren in diesem Fall 6 Iterationen notwendig.

Diese Methode ist sehr einfach und sie konvergiert immer. Bei kleinem $\epsilon$ wird es nie Probleme geben. Die Zahl der nötigen Iterationen steigt jedoch mit zunehmendem $\epsilon$. So sind z.B. mit den Startwerten $\epsilon = 0.99$ und $M = 2^{\circ}$ nicht weniger als 94 Iterationen notwendig, um das Ergebnis mit einer Genauigkeit von $0\overset{\circ}{.}000001$ zu erhalten:

Die folgende Tabelle 1 zeigt die Zwischenergebnisse der 94 Iterationen für $\epsilon = 0.99$ und $M = 2^{\circ}$:

Tabelle 1: Iterationen nach Methode 1
# $M[^{\circ}]$ # $M[^{\circ}]$ # $M[^{\circ}]$ # $M[^{\circ}]$
01 2.000000 25 30.817592 49 32.338768 73 32.360703
02 3.979598 26 31.059475 50 32.342408 74 32.360753
03 5.936635 27 31.264867 51 32.345452 75 32.360795
04 7.866758 28 31.438865 52 32.347998 76 32.360829
05 9.763644 29 31.585971 53 32.350128 77 32.360859
06 11.619294 30 31.710129 54 32.351909 78 32.360883
07 13.424417 31 31.814766 55 32.353398 79 32.360903
08 15.168909 32 31.902843 56 32.354644 80 32.360920
09 16.842404 33 31.976903 57 32.355686 81 32.360935
10 18.434883 34 32.039122 58 32.356557 82 32.360947
11 19.937269 35 32.091354 59 32.357286 83 32.360957
12 21.341978 36 32.135176 60 32.357895 84 32.360965
13 22.643349 37 32.171921 61 32.358405 85 32.360972
14 23.837929 38 32.202720 62 32.358831 86 32.360978
15 24.924579 39 32.228525 63 32.359187 87 32.360983
16 25.904408 40 32.250138 64 32.359485 88 32.360987
17 26.780556 41 32.268237 65 32.359735 89 32.360990
18 27.557863 42 32.283389 66 32.359943 90 32.360993
19 28.242483 43 32.296071 67 32.360117 91 32.360995
20 28.841471 44 32.306685 68 32.360263 92 32.360997
21 29.362399 45 32.315567 69 32.360385 93 32.360999
22 29.813009 46 32.322999 70 32.360487 94 32.361000
23 30.200940 47 32.329216 71 32.360572 95 32.361002
24 30.533515 48 32.334417 72 32.360643 96 32.361002

Die ersten beiden kursiven Zahlen sind die Startwerte, also die ersten beiden berechneten Werte, mit denen die Iterationsschleife gestartet wird. Die restlichen sind die 94 Iterationen, bis eine Genauigkeit von 6 Kommastellen erreicht wird, d.h. die Differenz in der 6. Kommastelle = 0 wird. Nach 50 Iterationen weicht das Ergebnis (rot: $32\overset{\circ}{.}345452$) immer noch um mehr als $0\overset{\circ}{.}01$ vom korrekten Wert ($32\overset{\circ}{.}361002$) ab!

Es zeigt sich, dass derart viele Iterationen nur bei sehr großen Werten für $\epsilon$ auftreten, die sehr nahe an $1$ liegen. Die höchste Exzentrizität der Planeten hat Merkur* mit etwa $\epsilon = 0.205631^{\circ}$ (für Epoche $J2000$). Die Berechnung für alle Planeten mit dieser Methode sollte daher keine Probleme bereiten. Für Kometen, die eine sehr hohe Exzentrizität haben können, empfiehlt sich diese Methode nicht.

* Pluto hat eine größere Exzentrizität als Merkur mit $\epsilon_{Pluto} = 0.2488^{\circ}$, gilt aber seit 2006 nicht mehr als Planet (Zwergplanet).

Beispiel JavaScript

/* Sinusfunktion in Grad */  
const DEGS = Math.PI/180;

function sin(x) {
  return Math.sin(x * DEGS);
} 

/*
Berechnung von E in der Keplergleichung, Methode 1
Parameter:
M... mittlere Anomalie, dezimaler Winkelwert
e... numerische Exzentrizität, dimensionlslos
Rückgabe:
E... exzentrische Anomalie, dezimaler Winkelwert
*/

function solveKepler1 (M, e) {
  let E0, E1;
  // Umrechnen von e in Grad
  e /= DEGS;
  // Startwerte festlegen
  E0 = M;
  E1 = M + e * sin(E0);
  // Schleife für Iteration:
  // Die absolute Differenz soll z.B. nicht
  // größer als 0.000001 sein = 1e-6
  while (Math.abs(E1 - E0) > 1e-6) {
    E0 = E1;
    E1 = M + e * sin(E0);
  }
  return E1;
}

console.log(solveKepler1(5, 0.1));
// 5.554589200183587  (6 Iterationen)

Methode 2

Wenn die Bahnexzentrizität $\epsilon \gt 0.3$ ist, kann die oben beschriebene Methode 1 so langsam konvergieren, dass es ratsam ist, eine bessere Iterationsformel zu benutzen: Die Newton-Raphson Iteration. Einen besseren Wert $E_1$ für $E$ erhält man dabei mit

$${E_1}^{rad} = {E_0} + \frac{{M + \epsilon \cdot \sin({E_0}) - {E_0}}}{{1 - \epsilon \cdot \cos{E_0}}}\tag{8}$$

wobei $E_0$ der zuletzt für $E$ erhaltene Wert ist. In dieser Formel sind die Winkel $M$, $E_0$ und $E_1$ in Bogenmaß ausgedrückt. Will man im Gradmodus arbeiten, dann muss man die Exzentrizität $\epsilon$ nur im Zähler des Bruches durch die modifizierte Exzentrizität $\epsilon_0 = \epsilon \cdot \tfrac{180}{\pi}$ ersetzen:

$${E_1}^{\circ} = {E_0} + \frac{{M + \epsilon \cdot \tfrac{180}{\pi } \cdot \sin({E_0}) - {E_0}}}{{1 - \epsilon \cdot \cos{E_0}}}\tag{9}$$

Der Prozess muss auch hier so oft wiederholt werden, bis die geforderte Genauigkeit erreicht wird. Man beachte den Unterschied zwischen der Formeln der Methode 1 und jener der Methode 2. Die Methode 1 liefert direkt eine neue Näherung für $E$. Während auch die Formel der Methode 2 eine neue Näherung $E_1$ für die exzentrische Anomalie liefert, ist der Bruch im zweiten Teil rechts bereits eine Korrektur zum vorherigen Wert $E_0$. Die Newton-Raphson Iteration konvergiert schneller als in der Methode 1.

Beispiel JavaScript

Hier ein Beispielcode für Methode 2, falls man in Grad rechen möchte.

/* Sinus- und Cosinusfunktion in Grad */

const DEGS = Math.PI/180;

function sin(x) {
  return Math.sin(x * DEGS);
}

function cos(x) {
  return Math.cos(x * DEGS);
}

/*
Berechnung von E in der Keplergleichung, Methode 2
Diese Methode benutzt den Korrekturwert als
Schleifenbedingung
Parameter:
M... mittlere Anomalie, dezimaler Winkelwert
e... numerische Exzentrizität, dimensionlslos
Rückgabe:
E... exzentrische Anomalie, dezimaler Winkelwert
*/

function solveKepler2 (M, e) {
  let E0, E1, corr;
  // Startwert setzen
  E1 = M;
  // Der Start-Korrekturwert ist beliebig, nur größer als
  // die angestebte Genauigkeit, damit die Schleife
  // betreten wird.
  corr = 1;
  // Schleife für Iteration
  // Die absolute Korrekturwert soll z.B. 
  // nicht größer als 0.000001 sein = 1e-6
  while (Math.abs(corr) > 1e-6) {
    E0 = E1;
    corr = (M + e/DEGS * sin(E0) - E0) / (1 - e * cos(E0));
    E1 = E0 + corr;
  }
  return E1;
}
        
console.log(solveKepler2(5, 0.1));
// 5.554589253872384  (3 Iterationen)

Wenn man wieder das Beispiel aus Methode 1 berechnet, also $\epsilon = 0.1$ und $M = 5^{\circ}$, erhält man hier bereits nach 3 Iterationen eine Genauigkeit von 9 Kommastellen.

Tabelle 2
$E_0$ Korrekturwert $E_1$
$5.000000000$ $0.554616193$ $5.554616193$
$5.554616193$ $-0.000026939$ $5.554589254$
$5.554589254$ $-0.000000000$ $5.554589254$

Es sei darauf hingewiesen, dass eine Genauigkeit von $10^{-9}$ absurd ist, das entspricht etwa 0.000004 Bogensekunden! Der Grenzwert von $10^{-9}$ wurde hier zu Anschauungszwecken gewählt. Er lässt sich aber in der Funktion nach Belieben setzen.

In Tabelle 3 sind die Ergebnisse aus der Keplergleichung für verschiedene Werte von $\epsilon$ und $M$ zu finden. Außerdem ist die Anzahl der Iterationen angegeben, die bei Benutzung der ersten und der zweiten Methode zur Lösung mit dem Startwert $E = M$ nötig sind. Die Iterationen wurden so lange ausgeführt, bis sich der neue Wert für $\epsilon$ um weniger als $0\overset{\circ}{.}000001$ vom vorigen unterscheidet.

Tabelle 3: Anzahl der Iterationen
Ergebnis Iterationen
$\epsilon$ $M[^{\circ}]$ $E[^{\circ}]$ Methode 1 Methode 2
$0.1 $ $5^{\circ} $ $5\overset{\circ}{.}554589 $ $6 $ $2 $
$0.2 $ $5^{\circ} $ $6\overset{\circ}{.}246908 $ $9 $ $2 $
$0.3 $ $5^{\circ} $ $7\overset{\circ}{.}134960 $ $12 $ $2 $
$0.4 $ $5^{\circ} $ $8\overset{\circ}{.}313903 $ $16 $ $2 $
$0.5 $ $5^{\circ} $ $9\overset{\circ}{.}950063 $ $21 $ $2 $
$0.6 $ $5^{\circ} $ $12\overset{\circ}{.}356653$ $28 $ $3 $
$0.7 $ $5^{\circ} $ $16\overset{\circ}{.}167990$ $39 $ $3 $
$0.8 $ $5^{\circ} $ $22\overset{\circ}{.}656579$ $52 $ $4 $
$0.9 $ $5^{\circ} $ $33\overset{\circ}{.}344447$ $58 $ $5 $
$0.99$ $5^{\circ} $ $45\overset{\circ}{.}361023$ $50 $ $11$
$0.99$ $1^{\circ} $ $24\overset{\circ}{.}725822$ $150$ $6 $
$0.99$ $33^{\circ}$ $89\overset{\circ}{.}722155$ $6 $ $5 $

Allgemein kann man sagen, dass ein größerer Wert von $\epsilon$ eine höhere Zahl von Iterationen erfordern wird. Dies gilt sowohl für die erste als auch für die zweite Methode. Mit der zweiten Methode ist die Anzahl der Iterationen aufgrund der starken Konvergenz jedoch viel kleiner.

Für kleine Werte der Exzentrizität, also etwa für $\epsilon \lt 0.3$, ist wohl die erste Methode die bessere: Es kann vorteilhafter sein, 5 oder 10 leichte Iterationen anstatt zwei Iterationen mit einer viel komplexeren Formel durchzuführen. Erst für größere Werte der Exzentrizität ist die zweite Methode der ersten vorzuziehen.

Für große $\epsilon$ hat die erste Methode nachteilige Folgen, weil sie zu langsam konvergiert. Man sehe sich die vorletzte Zeile der Tabelle 3 an: Um $E$ für die Werte $\epsilon = 0.99$ und $M = 1^{\circ}$ zu erhalten, sind nicht weniger als 150 Iterationen nötig (für eine Genauigkeit von 6 Kommastellen).

Instabile Zone

Auch die Methode 2 hat bezüglich ihrer Konvergenz einige Fallstricke. Die folgende Abb. 3 ist eine dreidimensionale Darstellung der Anzahl der Iterationen zum Erreichen von $E$ mit einer Genauigkeit von $10^{-9}$ Grad in Abhängigkeit von der Bahnexzentrizität $\epsilon$ und der mittleren Anomalie $M$ unter Benutzung der Methode 2. Wie oben wird $M$ als Startwert für $E$ benutzt. Die linke Ecke in der Nähe von $\epsilon \approx 1$ und kleinem $M$ ist die „instabile Zone“. Man sieht eine große Zahl von eng benachbarten Spitzen; die Anzahl der zum Erreichen des Ergebnisses nötigen Schritte mit der erwähnten Genauigkeit ändert sich beträchtlich, selbst bei kleinsten Veränderungen von $\epsilon$ oder $M$.

Abb. 3: Anzahl der Iterationen zum Lösen der Keplergleichung (Methode 2)

Die „instabile Zone“ wird in Abb. 4 genauer illustriert. Die Höhe (Anzahl der Iterationen) wurde auf 50 begrenzt, auch wenn es Werte gibt, die mehr als 1000(!) Iterationen erfordern.

Abb. 4: Anzahl der Iterationen im Bereich $\epsilon\ge 0.96$ und $M$ von 0 bis 40 Grad (Methode 2)

Für Abb. 4 wurden alle Werte für $0.960 \le \epsilon \le 0.999$ (Intervall 0.001) und $0^{\circ} \le M \le 40^{\circ}$ (Intervall $0\overset{\circ}{.}1$) berechnet, also 16040 Werte. Es zeigt sich, dass die meisten Iterationsschritte für folgende Werte von $\epsilon$ und $M$ benötigt werden:

Iterationen nach Methode 2

Es werden alle Werte berücksichtigt, die mehr als 1000 Iterationen benötigen.

Tabelle 4: Anzahl der Iterationen in der
instabilen Zone
$\epsilon$ $M^{\circ}$ Iterationen
$0.983$ $13\overset{\circ}{.}8$ $1242$
$0.990$ $24\overset{\circ}{.}5$ $1249$
$0.994$ $3\overset{\circ}{.}0 $ $1018$
$0.997$ $5\overset{\circ}{.}4 $ $1179$
$0.997$ $17\overset{\circ}{.}6$ $1317$
$0.997$ $20\overset{\circ}{.}4$ $1689$
$0.997$ $20\overset{\circ}{.}6$ $1366$
$0.998$ $21\overset{\circ}{.}8$ $1032$
$0.999$ $1\overset{\circ}{.}3 $ $1091$
$0.999$ $20\overset{\circ}{.}8$ $2755$(!)

In der letzten Zeile der Tabelle 2 sieht man, dass bei den Werten $\epsilon = 0.999$ und $M = 20\overset{\circ}{.}8$ insgesamt 2755 Iterationen benötigt werden, um das Ergebnis auf 6 Kommastellen genau zu berechnen: $E = 76\overset{\circ}{.}443861$. Es ist bemerkenswert, dass eine winzige Änderung auf $M = 20\overset{\circ}{.}81$ nur mehr 11 Iterationen durchlaufen. Nimmt man hingegen $\epsilon = 0.999$ und $M = 20\overset{\circ}{.}82$, werden sogar 7358 Iterationen gebraucht. Während der Berechnung springen die Werte auf die Größenordnung von $5\cdot 10^{126}$!


Methode 3 - Reihenentwicklung

Wenn die Exzentrizität $\epsilon$ klein ist, kann man sich die iterative Berechnung der Keplergleichung ersparen, indem man eine Reihenentwicklung für die wahre Anomalie verwendet. Die Berechnung von $E$ wird also nicht mehr benötigt, man erhält sofort die wahre Anomalie $\nu$.

W.M. Smart gibt in seinem Buch Textbook on Spherical Astronomy folgenden Zusammenhang an:

\[\begin{align} C &= \nu - M =\\ &+ \left(2 \cdot \epsilon - \frac{1}{4} \cdot {\epsilon^3}\right) \cdot \sin (M)\\ &+ \left(\frac{5}{4} \cdot {\epsilon^2}\right) \cdot \sin (2\cdot M) \\ &+ \left(\frac{{13}}{{12}} \cdot {\epsilon^3}\right) \cdot \sin (3\cdot M) \end{align}\tag{10}\]

$C$ nennt man „Equation of the center“, also Mittelpunktsgleichung.

Ist $M$ und $\epsilon$ gegeben, kann man daraus direkt die wahre Anomalie $\nu$ bestimmen. Wie oben angegeben ist die Gleichung nur verwendbar, wenn die Exzentrizität sehr klein ist.

J. Meeus gibt in seinen Astronomical Algorithms folgende Beziehung an, die bis zur 5. Potenz in der Exzentrizität $\epsilon$ geht:

\[\begin{align} C &= \nu - M =\\ &+ \left(2 \cdot \epsilon - \frac{1}{4} \cdot {\epsilon^3} + \frac{5}{96} \cdot {\epsilon^5}\right) \cdot \sin (M) \\ &+ \left(\frac{5}{4} \cdot {\epsilon^2} - \frac{{11}}{{24}} \cdot {\epsilon^4}\right) \cdot \sin (2\cdot M) \\ &+ \left(\frac{{13}}{{12}} \cdot {\epsilon^3} - \frac{{43}}{{64}} \cdot {\epsilon^5}\right) \cdot \sin (3\cdot M) \\ &+ \left(\frac{{103}}{{96}} \cdot {\epsilon^4}\right) \cdot \sin (4\cdot M) \\ &+ \left(\frac{{1097}}{{960}} \cdot {\epsilon^5}\right) \cdot \sin (5\cdot M) \end{align}\tag{11}\]

Man erkennt: Lässt man die Terme für $\epsilon^4$ und $\epsilon^5$ hier weg, erhält man die obige Formel von Smart.

Meeus gibt folgende Abweichungen an, wenn man anstatt der iterativen Lösung der Keplergleichung die Mittelpunktsgleichung verwendet:

Abweichungen von der Mittelpunktsgleichung

Tabelle 5
$\epsilon$ Terme bis $\epsilon^5$ Terme bis $\epsilon^3$
$0.03$ $0\overset{''}{.}00032$ $0\overset{''}{.}2371$
$0.05$ $0\overset{''}{.}0071$ $1\overset{''}{.}838$
$0.10$ $0\overset{''}{.}45$ $29\overset{''}{.}72$ ($0\overset{\circ}{.}0083$)
$0.15$ $5\overset{''}{.}2$ $151\overset{''}{.}8$ ($0\overset{\circ}{.}0422$)
$0.20$ $29\overset{''}{.}2$ $482\overset{''}{.}7$ ($0\overset{\circ}{.}1342$)
$0.25$ $111\overset{''}{.}3$ $1182\overset{''}{.}8$ ($0\overset{\circ}{.}3286$)
$0.30$ $330\overset{''}{.}5$ $2455\overset{''}{.}1$ ($0\overset{\circ}{.}6822$)

Für größere Werte als $\epsilon \gt 0.3$ sollte die Mittelpunktsgleichung ohnehin nicht verwendet werden.

In Tabelle 6 werden die Differenzen zwischen der Berechnung mittels Iteration (Keplergleichung) und der Berechnung mittels Mittelpunktsgleichung für die Planeten Merkur bis Neptun angegeben. Die Abbruchbedingung für die Iteration wurde hier bei ${10^{-6}}$ Grad angesetzt. Die Exzentrizitäten der Planeten gelten jeweils für die Epoche $J2000$.

Abweichungen der Mittelpunktsgleichung für die Planeten

Tabelle 6
Planet Exzentrizität $\epsilon$ Größte Differenz
Terme bis $\epsilon^5$ Terme bis $\epsilon^3$
Merkur $0.20563175$ $34\overset{''}{.}54$ bei $M = 69\overset{\circ}{.}7$ $539\overset{''}{.}66$ bei $M = 60\overset{\circ}{.}9$
Venus $0.00677192$ $4''\cdot 10^{-8}$ bei $M = 73\overset{\circ}{.}2$ $0\overset{''}{.}000612$ bei $M = 65\overset{\circ}{.}3$
Erde $0.01670863$ $9\overset{''}{.}8\cdot 10^{-6}$ bei $M = 72\overset{\circ}{.}8$ $0\overset{''}{.}022741$ bei $M = 65\overset{\circ}{.}1$
Mars $0.09340065$ $0\overset{''}{.}3016$ bei $M = 71\overset{\circ}{.}6$ $22\overset{''}{.}5943$ bei $M = 65\overset{\circ}{.}3$
Jupiter $0.04849793$ $0\overset{''}{.}0059$ bei $M = 72\overset{\circ}{.}3$ $1\overset{''}{.}6267$ bei $M = 64\overset{\circ}{.}4$
Saturn $0.05554814$ $0\overset{''}{.}0133$ bei $M = 72\overset{\circ}{.}2$ $2\overset{''}{.}8041$ bei $M = 64\overset{\circ}{.}3$
Uranus $0.04638122$ $0\overset{''}{.}0045$ bei $M = 72\overset{\circ}{.}3$ $1\overset{''}{.}3601$ bei $M = 64\overset{\circ}{.}5$
Neptun $0.00945575$ $3''\cdot 10^{-7} $ bei $M = 73\overset{\circ}{.}0$ $0\overset{''}{.}0023$ bei $M = 65\overset{\circ}{.}3$

Am Beispiel der Erde mit einer Exzentrizität von $\epsilon = 0.016709$ sieht man, dass man hier die Terme in der 4. und 5. Potenz ohne Probleme weglassen kann. Nur bei Merkur, der eine Exzentrizität von $0.2056$ aufweist, erhält man eine größere Abweichung, wenn man mit der Mittelpunktsgleichung rechnet. Auffallend ist, dass die größten Abweichungen alle im Bereich $60^{\circ} \lt M \lt 73^{\circ}$ liegen.


Methode 4 - Einfache Näherungsformel

Eine relativ einfache Formel, welche keine Iteration erfordert, lautet

$$\tan E = \large \frac{{\sin M}}{{\cos M - \epsilon}}\tag{12}$$

Sie gibt einen Näherungswert für $E$ an und gilt wieder nur für kleine Werte der Exzentrizität. Für dieselben Daten wie in Beispiel für Methode 1, nämlich $M = 5^{\circ}$ und $\epsilon = 0.1$ ergibt sich der Wert

$$\tan E = \frac{0.08715574}{0.89619470}\; \Rightarrow \; E = 5\overset{\circ}{.}554599$$

Der genaue Wert ist $5\overset{\circ}{.}554589$, der Fehler beträgt also nur $0\overset{\circ}{.}00001 = 0\overset{''}{.}036$. Der Fehler wächst allerdings relativ schnell mit zunehmender Exzentrizität.

Tabelle 7
$\epsilon$ Größter Fehler im Vergleich zur iterativen Lösung
$0.05$ $0\overset{\circ}{.}0012$
$0.10$ $0\overset{\circ}{.}0096$
$0.15$ $0\overset{\circ}{.}0327$
$0.20$ $0\overset{\circ}{.}0783$
$0.25$ $0\overset{\circ}{.}1552$
$0.30$ $0\overset{\circ}{.}2731$
$0.50$ $1\overset{\circ}{.}42$
$0.75$ $6\overset{\circ}{.}43$
$0.95$ $24\overset{\circ}{.}7$

Für die Erdumlaufbahn mit $\epsilon = 0.016709$ beträgt der Fehler weniger als $0\overset{''}{.}2$. In diesem Fall kann diese Formel sicher verwendet werden, es sei denn, eine hohe Genauigkeit ist erforderlich.

Gegebenenfalls kann diese Näherungsformel einen Startwert für die Methoden 1 bzw. 2 liefern.

Grafiken

In Abb. 5 wurde für einige Werte der Exzentrizität $\epsilon$ die exzentrische Anomalie über der mittleren Anomalie aufgetragen. Man sieht die Übereinstimmung der Winkel bei $0^{\circ}$ (Perihel), $180^{\circ}$ (Aphel) und $360^{\circ}$ (wieder Perihel).

  • In der ersten Hälfte eilt der Winkel $E$ dem Winkel $M$ voraus: $E \gt M$
  • In der zweiten Hälfte läuft der Winkel $E$ dem Winkel $M$ hinterher: $E \lt M$

Abb. 5: Mittlere und exzentrische Anomalie

In Abb. 6 wurde für einen halben Umlauf, also für $0^{\circ} \le M \le 180^{\circ}$, die Abweichung der Mittelpunktsgleichung (Terme bis $\epsilon^5$) von der Lösung der Keplergleichung aufgetragen. Da die Mittelpunktsgleichung die wahre Exzentrizität $\nu$ liefert, wurde auch bei der Iteration der Keplergleichung $\nu$ berechnet. Man sieht das rasche Fehlerwachstum, wenn die Exzentrizität größer wird.

Abb. 6: Abweichung der Mittelpunktsgleichung von der Lösung mittels Keplergleichung

In Abb. 7 wurde ebenfalls für einen halben Umlauf die Abweichung der Näherungslösung aus Methode 4 von der Lösung der Keplergleichung aufgetragen. Die Näherungsformel sollte nur für sehr kleine Exzentrizitäten verwendet werden. Hier sind die Fehler für $\epsilon = 0.01\;\textrm{-}\; 0.06$ dargestellt. Die Näherungsformel liefert eigentlich die exzentrische Anomalie, daher wurde die Ausgabe nach Ermittlung von $E$ mit derselben Methode in die wahre Anomalie umgerechnet, die bei der Keplergleichung verwendet wurde.

Abb. 7: Abweichung der Methode 4 von der Lösung mittels Keplergleichung

Es folgen noch zwei Abbildungen, in denen die Näherungslösung (Methode 4) für die Exzentrizitäten der Planeten Merkur bis Neptun berechnet und mit der Keplergleichung (Methode 2) verglichen werden.

Abb. 8: Abweichung der Methode 4 von der Lösung mittels Keplergleichung für die Planeten

Abb. 9: Hier ergibt sich die größe Abweichung für Jupiter mit ca. $4''$. Die Abweichung für die Erde bleibt immer unter $0\overset{''}{.}2$.

Abb. 9: Abweichung der Methode 4 von der Lösung mittels Keplergleichung für die Planeten

Größere Abweichungen gibt es für die Planeten Merkur, Mars und Saturn. Man sieht, dass vor allem Merkur ausreißt, er hat auch die größte Exzentrizität. Die größte Abweichung liegt bei $306\overset{''}{.}91$ beim Winkel $70\overset{\circ}{.}39$. Zur Erstellung einer Grafik wäre das aber immer noch ausreichend, da $300'' = 5' = 0\overset{\circ}{.}08333$ immer noch ein recht kleiner Gradwert ist und in einer Grafik praktisch nicht sichtbar wäre.

loesung_der_keplergleichung.txt · Zuletzt geändert: 2024/05/30 14:41 von hcgreier