Celestial Mechanics

Kunal Sarkhel, Sam Jugus, and Harnam Arneja

Background on Celestial Mechanics

Celestial Mechanics refer to the motion of objects in space. This project involves mathematical modeling of these mechanics considering the gravational attraction between two objects with mass. The gravational force law betwen two objects is given by the equation below: $$ F_g= \frac{G m_1 m_2}{r^2} $$ In this equation G represents the unviersal gravational constant, m represents the mass of each object and r represents the distance between two objects. Since $F=ma$ the previous equation can be divided by either $m_1$ or $m_2$ to find $a_1$ or $a_2$ respectively. Notice that when doing this $a_1$ depends on $m_2$ and $a_2$ depends on $m_1$. The equations are shown below: $$ a_1= \frac{Gm_2}{r^2}\\ a_2= \frac{Gm_1}{r^2}\\ $$

One Body Problem

The project description can be found here: Project Description. When one object is much more massive than another, the large body will feel an insignificant force as compared to the smaller obejct. When a mathematical model chooses to ignore the pull the smaller mass object creates on the larger object it is known as the one body problem. A simulation of the one body problem is shown below. Notice how the center object does not move, as our mathematical model treats the object with no force and no initial velocity.

Two Body Problem

The two body problem, ran from the Project Code, gets rid of the assumption that the more massive object feels no gravataional attraction. Once this assumption is eliminated, the new differential equations are shown below. The $\Delta$ represnts the difference in distance in the X or Y components of the two objects. R represents the total distance between the two objects. Since the accelration is the second derivative of position, order reduction was used to change two second order differential equations into four first order differnetial equations. Breaking each vector into its components doubles the amount of equations giving us a total of eight. The trapezoid iterative method was used to solve the differential equations for the object's position which is shown in the simulations. \begin{align} \dot{x}_{1x}&=V_{1x}\\ \dot{V}_{1x}&=- \frac{mg_2 \cdot \Delta X}{r^2}\\ \dot{x}_{1y}&=V_{1y}\\ \dot{V}_{1y}&=- \frac{mg_2 \cdot \Delta Y}{r^2}\\ \dot{x}_{2x}&=V_{2x}\\ \dot{V}_{2x}&=- \frac{mg_1 \cdot \Delta X}{r^2}\\ \dot{x}_{2y}&=V_{2y}\\ \dot{V}_{2y}&=- \frac{mg_1 \cdot \Delta Y}{r^2}\\ \end{align}

Adjusting to the Three Body Problem

The instructions for the project and the relevant background information can be found in the reality check. In order to adopt the two body problem to a three body problem the differential equations of motion were changed to account for the new system of three objects. The code can be found here. The updated accelerations are shown in the equations below.

$$ A_{1x}=-\left(\frac{gm_2 \cdot \Delta X_{12}}{R^3_{12}}\right)-\left(\frac{gm_3 \cdot \Delta X_{13}}{R^3_{13}}\right)\\ A_{1y}=-\left(\frac{gm_2 \cdot \Delta Y_{12}}{R^3_{12}}\right)-\left(\frac{gm_3 \cdot \Delta Y_{13}}{R^3_{13}}\right)\\ A_{2x}=\left(\frac{gm_1 \cdot \Delta X_{12}}{R^3_{12}}\right)-\left(\frac{gm_3 \cdot \Delta X_{23}}{R^3_{23}}\right)\\ A_{2y}=\left(\frac{gm_1 \cdot \Delta Y_{12}}{R^3_{12}}\right)-\left(\frac{gm_3 \cdot \Delta Y_{23}}{R^3_{23}}\right)\\ A_{3x}=\left(\frac{gm_2 \cdot \Delta X_{23}}{R^3_{23}}\right)+\left(\frac{gm_1 \cdot \Delta X_{13}}{R^3_{13}}\right)\\ A_{3y}=\left(\frac{gm_2 \cdot \Delta X_{23}}{R^3_{23}}\right)+\left(\frac{gm_1 \cdot \Delta Y_{13}}{R^3_{13}}\right) $$

The variable 'A' represents the acceleration of each planet and is broken into its x and y components. $ \Delta X_{ij}$ and $ \Delta Y_{ij}$ represents the X and Y distance between two planets. $ R_{ij}$ represents the total distance between two planets which can be represented as $R^2_{ij}=\Delta X^2_{ij} + \Delta Y^2_{ij} $.

Three Body Problem Simulations

In this part of the project, the code was run with the following initial parameters represented by the following set of values below. Here m represnts the mass, x and y are initial spatial components, and $\dot{x}$ and $\dot{y}$ are initial velocities.

$$ \{m \,,\, x,\, \dot{x},\, y,\, \dot{y}\} $$ \begin{align} Planet\, 1 \,(Green)&: \{.03 \,,\, 2,\, 0.2,\, 2,\, -0.2\} \\ Planet\, 2 \,(Red)&:\{.3 \,,\, 0,\, 0,\, 0,\, 0\} \\ Planet\, 3 \,(Blue)&:\{.03 \,,\, -2,\, -0.2,\, -2,\, 0.2\}\\ \end{align}

Question 10A:

The resulting three body simulation under these conditons are shown below. As you can see, the green and blue planets orbit around the central red planet. This result is expected as the red planet has the largest mass and has no initial velocity.

Question 10B:

However when $\dot{x_1}$, the initial x velocity of planet 1, is changed by $\frac{1}{100000}$ from 2 to 2.00001 the planets no longer remain in a stable orbit. As you can see below, after about two complete orbits the blue planet begins to drift further away and then suddenly diverges in the positive y-direction. The small changes in initial conditions shows that the system can drastically change due to small disturbances.

Equal Mass Three Body Problem

For this section of the projects, the masses of all three planets are set to be equal. In addition, the new initial conditions are listed below. This figure eight orbit is a fascinating example or gravataional laws.

$$ \{m \,,\, x,\, \dot{x},\, y,\, \dot{y}\} $$ \begin{align} Planet\, 1 \,(Green)&: \{1 \,,\, -0.970 ,\, -0.466,\, 0.243,\, -0.433\} \\ Planet\, 2 \,(Red)&:\{1 \,,\, 0.970 ,\, -0.466,\, -0.243,\, -0.433\} \\ Planet\, 3 \,(Blue)&:\{1 \,,\, 0,\, .932,\, 0,\, .866\}\\ \end{align}

Question 11A:

Below is the simulation given by the three body problems based on the initial conditions listed above.

Question 11B:

In this section, the effect of changing $\dot{x_3}$ is investigated. $\dot{x_3}$ was changed by $10^{-k}$ for $1 \le k \le 5 $. Adding $10^{-k}$ caused the figure eight to shift to the right, where as subtracting $10^{-k}$ resulted in the figure eight shifting to the left. Since the behavior of the orbits is symmetric despite the sign of the displacement, the diagrams below all show subtraction of the displacement from the initial x-velocity of planet 3. When $k\ \le 2 $, the orbit changes and is shown below. For $k\ \ge 3 $ the orbit is unchanged and is the same as shown in 11A.

Resulting orbit for $k=1$ is shown below. As seen by the orbit, this small displacement casues drastic changes in the orbit to happen quickly.

Resulting orbit for $k=2$ is shown below. In this orbit, the smaller displacement in the hundreths place causes the figure eight to slowly drift off center and not be in equilibrium.