Homework: The Keppler problem¶
- Simulate the motion of Earth in the solar system as a two body problem (taking into account only the Sun and the Earth). Derive the Newton's equation for the relative coordinate:
Turn the above equation into atronomical units (AU). The astronomical units are: - length is meassured in units of distance between Earth and Sun $R\approx 1.5\;10^{11}$m - time is meassured in years.
Plot the Earth's orbit $(x,y)$ for 5000 years, and verify it is stable (the orbit is a circle with no time dependence). Also plot $x(t)$ and $y(t)$ for the last 5 out of 5000 years to see that period is what is expected.
- Simulate the three body problem, which consists of Sun, Earth and Jupiter. The latter has mass $m_J/m_S\approx 9.55\; 10^{-4}$, and distance $5.2\; AU$. Note that the mass of Earth is $m_E/m_S\approx 3 10^{-6}$ and distance $1\; AU$.
Simulate 5000 Earth years and plot the orbits of Sun, Eart, and Jupiter ($x(t)$ versus $y(t)$ for all three objects). Are the orbits stable?
Check how strong is the influence of the Jupiter on motion of the Earth. Plot $x(t)$ for the last 5 years of 5000 years for the case with and withouth Jupiter.
- It turns out that our solar system has uneven distribution of asteroids, as demonstrated by the picture below. We are plotting the number of asteroids as a function of the distance from the Sun. There are many gaps in the distribution plot, which are now named Kirkwood gaps, after Daniel Kirkwood, who discovered them. He showed that gaps are associated with Jupiter, because the orbits are in resonance with Jupiter's motion. For example, the 2/1 gap is such that an asteroid placed there would complete two orbits every time Jupiter completes one. Similarly there are 3/1, 5/2, and 7/3 resonance, all related to Jupiter.
We would like to simulate Kirkwood gaps. To simulate 2/1 gap, we need to find the distance from the Sun that such asteroid would be placed at, and his initial velocity. Than we need to check the long term stability of such orbit.
We first recall that all orbits in the solar system satisfy $R_i^3/T_i^2=const$. This is because centrifugal force ($m_i\omega_i^2 R_i$) and gravitational force ($G M m_i/R_i^2$) have to be balanced, hence $m_i (2\pi/T_i)^2 R_i = G M m_i/R_i^2$, where $M$ is the Solar's mass. We hence see that $R_i^3/T_i^2=G M/(4\pi^2)$. In AU units this equation is $R_i^3/T_i^2=1$
For Asteroid that completes the orbit in half the Saturn's year, it must satisfy $T_{asteroid}=R_{Saturn}^{3/2}/2$, and hence $R_{asteroid}=R_{Saturn}/2^{2/3}$, which is $R_{asteroid}=5.2/2^{2/3}\approx 3.2758$. It's starting velocity should be $v_{asteroid}=2\pi/\sqrt{R_{asteroid}}\approx 3.4715$.
For 3/1 gap, we should similarly have $R_{asteroid}=5.2/3^{2/3}\approx 2.5$ and $v_{asteroid}=2\pi/\sqrt{R_{asteroid}}\approx 3.974$.
For homework, simulate the three body problem: Sun, Jupiter and Asteroid. You can assume that the mass of Asteroid is vanishingly small. Here it is also safe to ignore Earth and its influence on Asteroid. Simulate 5000 Earth years of motion, and plot the orbits of the three objects $x(t)$ versus $y(t)$. Make sure you subtract the center of motion movement $\vec{R}_{cm}=m_1\vec{r}_1+m_2\vec{r}_2+m_3\vec{r_3}$ when plotting the orbits.
Do you see any change of the orbit of Asteroid over 5000 years?
Sketch of the solution¶
- For Keppler problem we use $\vec{r}_1$ and $\vec{r}_2$ for vector positions of the two bodies, like the Earth and the Sun. The Newton's Eq require:
The center of mass can be fixed at the origin $m_1 \vec{r}_1+m_2 \vec{r}_2=0$.
The relative vector $$\vec{r}=\vec{r}_1-\vec{r}_2$$ then satisfies the equation \begin{equation} \ddot{\vec{r}}=G(m_1+m_2)\frac{\vec{r}}{|\vec{r}|^3} \end{equation} The gravitational constant $G$ can be obtained from information about the Earth's orbit. We know that Earth's orbit is almost circular. For circular orbits we know that centrifugal force and gravitational force have to be balanced, which means $m_E \omega^2 R= G m_S m_E/R^2$, where $R$ is Sun-Earth distance $m_S$ and $m_E$ are Sun and Earth mass, and $\omega$ frequency of the Earth's rotation, which is $\omega=2\pi/T$, with $T$ being one year. We thus have $$G=\frac{R^3}{m_S} \left(\frac{2\pi}{T}\right)^2 .$$
We want to meassure time in years $T$, and distance in Sun-Earth distance $R$. We thus define $$\vec{r}=R \vec{r_d}$$ $$t = T t_d $$ and in this AU units the above equation for relative vector is $$\ddot{\vec{r}_d} = -4\pi^2 (1+\frac{m_E}{m_S})\frac{\vec{r}_d}{r_d^3}$$ The ratio of the mass $m_E/m_S= 3\, 10^{-6}$ is safe to neglect.
The initial conditions will be choosen so that we have a circular orbit. For example $\vec{r}_d=[1,0]$ and $\dot{\vec{r}_d}=[0,2\pi]$. The velocity for a circular motion is $v=\omega R$, which is $v=2\pi R/T$, and hence in AU units is just $\dot{y}_0=2\pi$.
We can choose the following set of variables $[x,y,\dot{x},\dot{y}]$, where $\vec{r}_d=[x,y]$. With these we can solve the Kepler equations by: \begin{eqnarray} \begin{bmatrix}\frac{d x}{dt}\\\frac{d y}{dt}\\ \frac{d \dot{x}}{dt}\\ \frac{d \dot{y}}{dt}\\ \end{bmatrix} =\begin{bmatrix} \dot{x}\\\dot{y}\\ -4\pi^2 \frac{x}{(x^2+y^2)^{3/2}}\\ -4\pi^2 \frac{x}{(x^2+y^2)^{3/2}} \end{bmatrix} \end{eqnarray} with initial conditions $y_0=[1,0,0,2\pi]$
- In the Keppler's three body problem, we have Sun, Earth and Jupiter with $$\frac{m_2}{m_1}=3\; 10^{-6}$$ and $$\frac{m_3}{m_1}=9.55\; 10^{-4}.$$
The Newton's equations are \begin{eqnarray} && \ddot{\vec{r}}_1 = -G m_1\left( \frac{m_2}{m_1}\frac{\vec{r}_1-\vec{r}_2}{|\vec{r}_1-\vec{r}_2|^3} +\frac{m_3}{m_1}\frac{\vec{r}_1-\vec{r}_3}{|\vec{r}_1-\vec{r}_3|^3}\right) \\ && \ddot{\vec{r}}_2 = -G m_1\left(\frac{\vec{r}_2-\vec{r}_1}{|\vec{r}_2-\vec{r}_1|^3} +\frac{m_3}{m_1}\frac{\vec{r}_2-\vec{r}_3}{|\vec{r}_2-\vec{r}_3|^3}\right) \\ && \ddot{\vec{r}}_3 = -G m_1\left(\frac{\vec{r}_3-\vec{r}_1}{|\vec{r}_3-\vec{r}_1|^3} +\frac{m_2}{m_1}\frac{\vec{r}_3-\vec{r}_2}{|\vec{r}_3-\vec{r}_2|^3}\right) \\ \end{eqnarray} and in AU units become \begin{eqnarray} && \ddot{\vec{r}}_1 = -4\pi^2\left( \frac{m_2}{m_1}\frac{\vec{r}_1-\vec{r}_2}{|\vec{r}_1-\vec{r}_2|^3} +\frac{m_3}{m_1}\frac{\vec{r}_1-\vec{r}_3}{|\vec{r}_1-\vec{r}_3|^3}\right) \\ && \ddot{\vec{r}}_2 = -4\pi^2\left(\frac{\vec{r}_2-\vec{r}_1}{|\vec{r}_2-\vec{r}_1|^3} +\frac{m_3}{m_1}\frac{\vec{r}_2-\vec{r}_3}{|\vec{r}_2-\vec{r}_3|^3}\right) \\ && \ddot{\vec{r}}_3 = -4\pi^2\left(\frac{\vec{r}_3-\vec{r}_1}{|\vec{r}_3-\vec{r}_1|^3} +\frac{m_2}{m_1}\frac{\vec{r}_3-\vec{r}_2}{|\vec{r}_3-\vec{r}_2|^3}\right) \\ \end{eqnarray}
The initial condistions for Sun can be $\vec{r}_1=0$ and $\dot{\vec{r}}_1=0$, for Earth as before $\vec{r}_2=[1,0]$ and $\dot{\vec{r}}_2=[0,2\pi]$, and for Jupiter $\vec{r}_3=[R_J/R,0]$ and $\dot{\vec{r}}_3=[0,\frac{2\pi}{\sqrt{R_J/R}}]$, where $R_J/R\approx 5.2$. The velocity folows from circular motion of planets, for which we know that gravitational force and centrigunal force are balanced $\omega^2 R_i = G m_S/R_i^2$ or $v_i^2/R_i=G m_S/R_i^2$ and hence $v_i^2 R_i=G m_S = R^3 (2\pi/T)^2$. In AU units each planet with circular orbit should satisfy $v_i=\frac{2\pi}{\sqrt{R_i}}$
- Because asteroids are very light, we can neglect their mass compared to mass of Sun and Jupiter. In this case we can treat each Asteroid as an independent problem. We need to simulate the three bodies: Sun, Jupiter, and Asteroid. The masses are $m_2/m_1=9.55\; 10^{-4}$ and $m_3/m_1=0$.
The equations are identical to previous equations simulating Sun, Eart and Jupiter.
The initial condistions for circular orbits are also $R_{asteroid\;(2/1)}=R_{Saturn}/2^{2/3}$, $R_{asteroid\;(3/1)}=R_{Saturn}/3^{2/3}$ and $v_{asteroid}=2\pi/\sqrt{R_{asteroid}}$.
# last 5 of 5000 years
plt.plot(sol.t, sol.y[0])
plt.xlim([4995,5000]);
plt.xlabel('t(years)')
plt.ylabel('x(Earth-Sun distance)');
# Eart's orbit for 5000 years
fig, ax = plt.subplots(1,1)
ax.set_title('Earts orbit over 5000 years')
ax.plot(sol.y[0],sol.y[1],lw=0.1)
ax.set_aspect('equal', adjustable='box')
