Himmelsmechanik

Abbildung: System Erde-Sonne

Im Folgenden betrachten wir die Berechnung der Planetenbewegung, zunächst einmal der Erde, um die Sonne. Entscheidend hierbei ist die Gravitationskraft zwischen der Masse der Sonne ($m_{S}$) und der Masse der Erde ($m_{E}$). Die Formel für diese Kraft lautet :

$${\mathbf f}=-\frac{Gm_{S}m_{E}}{r^{2}}*\frac{{\mathbf r}}{r}$$ (5.66)

Der Faktor $G$ ist die Gravitationskonstante und der Vektor ${\mathbf r} = (x, y)$ ist der Abstandsvektor zwischen den Massenschwerpunkten von Sonne und Erde. Zur Vereinfachung betrachten wir die Bewegung nur im Zweidimensionalen. Für die Berechnung benutzen wir die Formel in der Komponentenschreibweise für die $x$-Richtung $f_{x}$ und in die $y$-Richtung $f_{y}$ mit

$$f_x=-\frac{Gm_Sm_Ex}{r^3} f_y=-\frac{Gm_Sm_Ey}{r^3} r=\sqrt{x^2+y^2}$$ (5.67)

Um angenehmere Zahlen zu bekommen, verwenden wir hier nicht das MKS-System (Meter-Kilo-Sekunde), sondern die Astronomische Einheit (1 AE $=1.5*10^{11}$m$\approx$ Abstand Erde - Sonne) für Längen und das Jahr ($1$a) für die Zeit. Die Bahngeschwindigkeit bei der Kreisbewegung ist somit $v_0=2\pi *\frac{1 {\rm AE}}{1 {\rm a}}$. Da die Gravitationskraft betragsmäßig gleich der Zentrifugalkraft ist, haben wir

$$f=\frac{m_Ev_0^2}{r}=\frac{Gm_Sm_E}{r^2} \Rightarrow Gm_S=v_0^2r=4\pi ^2\frac{{\rm AE}^3}{{\rm a}^2}$$ (5.68)

Für die numerische Berechnung verwenden wir wieder die Newtonsche Bewegungsgleichung

$$\ddot{x}=\frac{f_x}{m_E}=-\frac{4\pi ^2x}{r^3} \ddot{y}=\frac{f_y}{m_E}=-\frac{4\pi ^2y}{r^3}$$ (5.69)

Wir benutzen in diesem Fall wieder die Verlet-Methode für unser Programm und lösen für die numerische Rechnung die beiden Gleichungen

$$\begin{eqnarray}x(t+\Delta t) =-x(t-\Delta t)+2x(t)-4\pi ^2x(t)r_3(t)\Delta t^2 \\ y(t+\Delta t) =-y(t-\Delta t)+2y(t)-4\pi ^2y(t)r_3(t)\Delta t^2 \end{eqnarray}$$

Dabei ist die Abkürzung $r_3(t) = \left( \sqrt{x(t)^2+y(t)^2}\right) ^{-3}$ verwendet worden.

int main()
{
   double p0x, p1x, p2x, dp0x, dp1x, dp2x, t;
   double p0y, p1y, p2y, dp0y, dp1y, dp2y, r3;
   double alpha, t_max, dt;
 
// read parameters
   cout << "Anf.x  ";  cin >> p1x;
   cout << "Anf.y  ";  cin >> p1y;
   cout << "Anf.vx ";  cin >> dp1x;
   cout << "Anf.vy ";  cin >> dp1y;
   cout << "Alpha  ";  cin >> alpha;
   cout << "Gesamtzeit     ";  cin >> t_max;
   cout << "Zeitintervall  ";  cin >> dt;
 
   const double beta=(alpha)/2.;
   const double dt2=dt*dt;
   const double pi24=4.*M_PI*M_PI;
 
// setup initial conditions
   p0x=p1x-dt*dp1x;
   p0y=p1y-dt*dp1y;
   dp0x=dp1x;
   dp0y=dp1y;
 
// loop from t=0 to t=t_max
   t=0.0;
   ofstream outfile("orbit.dat");
   for( int i=0; i<t_max/dt; i++ ) {
        r3=1./pow(p1x*p1x+p1y*p1y,beta);
        p2x=2.*p1x-p0x-pi24*dt2*r3*p1x;  // Verlet algorithm
        p2y=2.*p1y-p0y-pi24*dt2*r3*p1y;
        outfile << t << ' ' << p1x << ' ' << p1y << '\n';
        t=t+dt;
        p0x=p1x;                         // change new->old
        p0y=p1y;
        p1x=p2x;
        p1y=p2y;
   }
   outfile.close();
}

Für ideale Anfangsbedingungen, das wäre $x_0=1$, $y_0=0$, $vx_0=0$ und $vy_0=2\pi$, erhält man für die Erdbewegung eine Kreisbahn, siehe Abbildung (links). Verändert man die Anfangsgeschwindigkeit auf $vy_0=5$, so ergibt sich eine elliptische Bahn (rechts).

Eine weitere Variationsmöglichkeit besteht in der Änderung des Gravitationsgesetzes. Anstatt die Gravitationskraft proportional zu $1/r^2$ abfallen zu lassen, kann man sie auch proportional zu $1/r^{\alpha -1}$ machen. Allgemein wird im obigen Programm r=1/sqrt(pow(r, alpha)) verwendet, wobei pow(r, alpha) der Potenzierung $r^\alpha$ entspricht. Abbildung zeigt vier verschiedene Beispiele mit Werten für $\beta$ zwischen $3.00$ und $2.01$, wobei $\beta=\alpha-1$. Die Rotation der Ellipse bei den Beispielen für $\beta \approx 2$ nennt man Präzession.

Für $\beta =3.00$ fällt auf, daß die Erde, nachdem sie in die Nähe der Sonne ``geflogen'' ist, sich mit hoher Geschwindigkeit aus dem System entfernt. Dies ist allerdings ein numerisches Problem. Für diesen Effekt ist die Wahl eines zu großen $\Delta t$ verantwortlich, daß bei großen Geschwindigkeitsänderungen kleiner sein sollte. Allerdings würde bei einem kleineren $\Delta t$ die Simulation zu lange dauern. Man benötigt also ein automatisches Verfahren, das die Größe von $\Delta t$ regelt. Eine Möglichkeit ist es, die Änderung der Energie für verschiedene $\Delta t$ zu vergleichen und den Zeitschritt dementsprechend zu verändern (adaptives Verfahren). Dazu bietet sich z.B. das Runge-Kutta Verfahren an, da dort sowieso Werte zur Zeit $t + \Delta t/2$ berechnet werden. Die Energie der Erde ist

$$E_{Erde}=-\frac{Gm_Sm_E}{r}+\frac{m_E}{2}v^2$$ (5.70)

Die Energie ergibt sich also aus den Koordinaten $x$ und $y$ für $r=\sqrt{x^2+y^2}$ direkt und für ${\mathbf v}=\Delta {\mathbf r}/\Delta t$ indirekt. Wir berechnen zuerst die Energie $E_1$ aus

$$\begin{eqnarray}x_1=x(t+\Delta t)=f_{\Delta t}(x(t),y(t)) \\ y_1=y(t+\Delta t)=g_{\Delta t}(x(t),y(t)) \end{eqnarray}$$

und die Energie $E_2$ mit halb so großem $\Delta t$ durch

$$\begin{eqnarray}x_2=f_{\frac{\Delta t}{2}}(x(t+\frac{\Delta t}{2}),y(t+\frac{\De... ...ac{\Delta t}{2}}(x(t+\frac{\Delta t}{2}),y(t+\frac{\Delta t}{2})) \end{eqnarray}$$

Die Funktionen $f$ und $g$ berechnen die neuen $x$- bzw. $y$-Werte für $\Delta t$ oder $\Delta t/2$. Das entscheidende Kriterium ist die Differenz

$$\Delta E=\vert E_2-E_1\vert$$ (5.71)

Ist $\Delta E$ groß, so ist $\Delta t$ zu klein gewählt, ist $\Delta E$ sehr klein, so kann $\Delta t$ größer gewählt werden. Im Programm könnte z.B. bei $\Delta E>10^{-5}\rightarrow \Delta t'=\Delta t/2$ und bei $\Delta E<10^{-7}\rightarrow \Delta t'=\Delta t*2$ gewählt werden.