1.3 Zeitlicher Verlauf der Bahn des Planeten

Da der zeitliche Verlauf der Bahn nicht exakt berechnet werden kann, schlagen wir eine Berechnung vor, die nur eine numerische Nullstellenberechnung (Abschnitt 1.3.1) und die numerische Berechnung eines bestimmten Integrals benötigt (Abschnitt 1.3.2).

1.3.1 Berechnung von r(t)

Aus (1.4) folgt:

\dot{φ} (t) = \frac{L}{m r^{2} (t)}

Einsetzten in (1.3') ergibt eine gewöhnliche Differentialgleichung (DGL) ersten Grades für r(t):

E = \frac{m}{2} ({\dot{r}}^{2} (t) + {(r (t) \frac{L}{m r^{2} (t)})}^{2}) - G M m \frac{1}{r (t)}

\Rightarrow \frac{m}{2} {\dot{r}}^{2} (t) = E + G M m \frac{1}{r (t)} - \frac{L^{2}}{2 m r^{2} (t)}

Das ergibt folgende DGL:

\dot{r} (t) = \pm \sqrt{\frac{2 E}{m} + \frac{2 G M}{r (t)} - \frac{L^{2}}{m^{2} r^{2} (t)}}

= \pm \frac{1}{r (t)} \sqrt{\frac{2 E}{m} r^{2} (t) + 2 G M r (t) - \frac{L^{2}}{m^{2}}}

Im Folgenden betrachten wir nur die Bahn des Planeten vom Perihel zum Aphel. In diesem Bereich ist $\dot{r} (t)$ stets positiv.

Mit

$a = \frac{2 E}{m}, b = 2 G M, c = \frac{L^{2}}{m^{2}}$

ergibt sich daraus folgende DGL:

$(1.8) \dot{r} (t) = \frac{1}{r (t)} \sqrt{a r^{2} (t) + b r (t) + c}$

Da die Bahnkurve des Planeten eine Ellipse ist, hat $\dot{r} (t)$ genau dann eine Nullstelle r_p oder r_a, wenn sich der Planet im Perihel oder im Aphel befindet.

Für die Nullstellen gilt:

$(1.8') r_{p} = r (t = 0) = \frac{- b + \sqrt{b^{2} - 4 a c}}{2 a} \Leftrightarrow \frac{2 a r_{p} + b}{\sqrt{b^{2} - 4 a c}} = 1$

$(1.8'') r_{a} = r (t = \frac{T}{2}) = \frac{- b - \sqrt{b^{2} - 4 a c}}{2 a} \Leftrightarrow \frac{2 a r_{a} + b}{\sqrt{b^{2} - 4 a c}} = - 1$

Hierbei ist T die Umlaufzeit des Planeten (genau ein Jahr für die Erde).

Gesucht ist eine Lösung der DGL (1.8) für die Anfangsbedingung

r (t = 0) = r_{p}

Die DGL (1.8) ist eine DGL mit getrennten Variablen und kann, zumindest implizit, elementar gelöst werden. Allgemein gilt:

Sei f eine Funktion in Abhängigkeit von einer Variablen t und g eine Funktion in Abhängigkeit von einer Variablen r. Dann heißt

\dot{r} = f (t) g (r)

eine DGL mit getrennten Variablen.

Wir definieren Funktionen F und G durch:

F (t) = \int_{0}^{t} f (x) d x, G (r) = \int_{r_{p}}^{r} \frac{1}{g (x)} d x

Dann gilt für eine Lösung r der obigen DGL mit $r (0) = r_{p}$ :

$G (r (t)) = F (t)$

Existenz und Eindeutigkeit der Lösung wird hier nicht weiter betrachtet.

Für DGL (1.8) gilt somit:

F (t) = \int_{0}^{t} 1 d x = t, G (r) = \int_{r_{p}}^{r} \frac{x}{\sqrt{a x^{2} + b x + c}} d x

Daraus ergibt sich:

$(1.9) \int_{r_{p}}^{r (t)} \frac{x}{\sqrt{a x^{2} + b x + c}} d x = t$

Für das Integral in (1.9) existiert eine Stammfunktion, siehe diverse Integraltabellen:

Sei $X = a x^{2} + b x + c$ und $Δ = 4 a c - b^{2}$ , dann gilt für $a < 0$ und $Δ < 0$ :

$\int \frac{x}{\sqrt{X}} d x = \frac{\sqrt{X}}{a} - \frac{b}{2 a} \int \frac{1}{\sqrt{X}} d x = \frac{\sqrt{X}}{a} - \frac{b}{2 a} (- \frac{1}{\sqrt{- a}}) \arcsin (\frac{2 a x + b}{\sqrt{- Δ}})$

Damit folgt aus (1.9) und dem Hauptsatz der Differential- und Integralrechnung:

$(1.10) \frac{\sqrt{a r^{2} (t) + b r (t) + c}}{a} + \frac{b}{2 a \sqrt{- a}} \arcsin (\frac{2 a r (t) + b}{\sqrt{b^{2} - 4 a c}}) -$

\frac{\sqrt{a r_{p}^{2} + b r_{p} + c}}{a} - \frac{b}{2 a \sqrt{- a}} \arcsin (\frac{2 a r_{p} + b}{\sqrt{b^{2} - 4 a c}}) = t

Mit (1.8') ergibt sich:

$(1.10') \frac{\sqrt{a r^{2} (t) + b r (t) + c}}{a} + \frac{b}{2 a \sqrt{- a}} \arcsin (\frac{2 a r (t) + b}{\sqrt{b^{2} - 4 a c}}) - \frac{b}{2 a \sqrt{- a}} \frac{π}{2} = t$

In dieser Gleichung kann r(t) nicht explizit bestimmt werden, jedoch kann r(t) für jedes t>0 durch eine Nullstellensuche numerisch berechnet werden.

Dazu definieren wir die Funktion z(x) für ein festes t > 0 und t < T/2 wie folgt:

z (x) = \frac{\sqrt{a x^{2} + b x + c}}{a} + \frac{b}{2 a \sqrt{- a}} \arcsin (\frac{2 a x + b}{\sqrt{b^{2} - 4 a c}}) - \frac{b}{2 a \sqrt{- a}} \frac{π}{2} - t

Dann folgt aus (1.9) und der Tatsache, dass der Planet vom Perihel zum Aphel die halbe Umlaufzeit benötigt:

$z (r_{p}) = - t < 0$ und $z (r_{a}) = \frac{T}{2} - t > 0$

Daher existiert eine Nullstelle x₀ von z, d.h. z(x₀) = 0, im offenen Intervall (r_p, r_a) mit der Eigenschaft:

r(t) = x₀

Die Nullstellensuche kann z.B. mit der Matlab Funktion „fzero“ durchgeführt werden. Als Anfangswert für die Suche kann r_p benutzt werden. Alternativ können natürlich auch klassische Nullstellensuchverfahren wie das Bisektionsverfahren benutzt werden.

1.3.2 Berechnung von φ(t)

Wie oben folgt aus (1.4) :

\dot{φ} (t) = \frac{L}{m r^{2} (t)}

Aus dem Hauptsatz der Differential- und Integralrechnung folgt daraus:

φ (t) - φ (0) = \int_{0}^{t} \dot{φ} (x) d x = \int_{0}^{t} \frac{L}{m r^{2} (x)} d x

Setzt man φ(0) = 0 so ergibt sich:

$(1.11) φ (t) = \int_{0}^{t} \frac{L}{m r^{2} (x)} d x$

Das obige Integral kann z.B. mit der Matlab Funktion „integral“ numerisch berechnet werden, hierzu werden die in Abschnitt 1.3.1 berechneten Werte für r(t) verwendet.

1.3.3 Berechnung von x(t) und y(t)

Aus der Darstellung der Bahn des Planeten in Polarkoordinaten (Abschnitt 1.3.1 und 1.3.2) ergibt sich für die Darstellung in kartesischen Koordinaten:

$x (t) = r (t) \cos (φ (t))$

$y (t) = r (t) \sin (φ (t))$

top

back