Matteo & Nadav: Week One

Nadav Hendel & Matteo Bjornson

Week One

Our first model, Planets, assumes the sun’s motion due to orbiting planets is negligible and we can keep it fixed at the center of our system. This approximation is a good starting point because the sun is so massive in comparison to any of the planets (Jupiter, the most massive planet, is still no more than 0.1% of the solar mass). The planet motion results from the gravitational force of the sun. The size of this force is described by

$F_G=\frac{GM_SM_p}{r^2}$

where $G$ is the gravitational constant, $M_S$ and $M_p$ are the Solar and planet mass, respectively. $r$ is the planet-sun distance.

Due to Newton’s Second Law we know the acceleration of the planet is equal to the force experienced by the planet divided by its mass. This results in a second order differential equation that can be written componentwise as

$\frac{d^2x}{dt}=\frac{F_{G,x}}{M_E}$

$\frac{d^2y}{dt}=\frac{F_{G,y}}{M_E}$

These can be rewritten as first order differential equations

$\frac{dv_x}{dt}=-\frac{GM_Sx}{r^3}$

$\frac{dv_y}{dt}=-\frac{GM_Sy}{r^3}$

where $v_x=\frac{dx}{dt}$ and $v_y=\frac{dy}{dt}$ . This program assumes the planet follows uniform circular motion. This is a fairly reasonable approximation to make, as all the planets except Mercury and Pluto have very circular orbits. For example, the eccentricity of Earth’s orbit is under 0.017 [1]. Making this assumption we know that the $GM_S$ product is equal to $v^2r$ . The average speed of the planet is equal to the distance traveled in one orbit divided by the orbital period, thus $GM_S=(\frac{2\pi r}{T})^2r$ . Using units of AU for distance and years for time, the $GM_S$ product for Earth (where T=1yr and r=1AU), for example is simply $4\pi^2$ . Therefore we can rewrite the first order differentials as

$\frac{dv_x}{dt}=-\frac{4\pi^2x}{rT}$

$\frac{dv_x}{dt}=-\frac{4\pi^2y}{rT}$

The solution to these first order differential equations, x and y, can be approximated using numerical methods. We chose to use the Euler-Cromer Method (Runge-Kutta would work as well), outlined by Giordano, Equations 4.7:

$v_{x,i+1}=v_{x,i}-\frac{4\pi^2x_i}{r_i^3}\Delta t$

$x_{i+1}=x_i+v_{x,i+1}\Delta t$

$v_{x,i+1}=v_{y,i}-\frac{4\pi^2y_i}{r_i^3}\Delta t$

$y_{i+1}=x_i+v_{y,i+1}\Delta t$

Thus, the only constants necessary to calculate the position of a planet is its orbital period and radius (in this case the semi-major axis). We will be using Astronomical Units (AU) for this section since 1 AU is defined as the average distance between the earth and sun, which has a nearly circular orbit (eccentricity of under 0.017) [1]. Using the Euler-Cromer method, we can create a plot in Planets that traces the earth’s position around the sun due to gravity for a set amount of time.

As long as the average velocity and semi-major axis are known, our program can currently be used to model any other real or hypothetical orbit of a planet around the sun, given the assumption of circular orbits. In fact, we have altered our final subroutine Planets to include the assumed circular orbit for every planet in our solar system excluding Mercury and Pluto. Both of these two planets have eccentricities that are relatively greater than the rest, and are too large to be considered even close to circular. While this is a somewhat subjective judgment made by observation, it seems logical given the relative magnitude of the other planets eccentricities. Given this, the user is prompted at the beginning of the code to pick which planet they wish to observe from the remaining seven.

Now that we have completed a working model for the two-body systems of each planet in our solar system, sans Mercury and Pluto, we wished to do a comparison between computational methods. Since energy must be conserved in our system, similar to the pendulum problems from earlier this semester, we know that difficulties might arise when using the Euler method. To validate this, we wrote another routine implementing Euler instead of Euler Cromer. A comparison of the methods in GIF form is shown below:

Pictured topis the earth’s orbit calculated using the Euler-Cromer method, and bottom is the same code, only run solely using the Euler Method

The radius and velocity of the orbits that we have been modeling are already known to be stable, but we wanted to note what an unstable system using our model might look like. By increasing or reducing the initial velocity, we can observe such a system.

earthunstable earthlessv0

Pictured top is the earth’s orbit calculated using a larger initial velocity by a factor of 2, and bottom is the same orbit using less than ½ original velocity.

Given our goals for week one, we believe we now have a good foundation to expand our Planet routine into a multi-body system. We have become confident with the Euler-Cromer method over the Euler Method, and we have explored the boundaries of our system given different initial conditions. Also, through exploration of chapter four in the Giordano text, we have expanded our understanding of the computational methods used to solve general celestial mechanical problems.

Please click the Matlab Icon to download our Planet routine code

Modeling and Experimental Tools with Prof. Magnes

Ideal for getting ideas and feedback on upper level physics tools.

1 thought on “Matteo & Nadav: Week One”