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


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



These can be rewritten as first order differential equations



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



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:


earthEC earthE


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.




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



1 thought on “Matteo & Nadav: Week One

  1. Avatar photoJenny Magnes

    The graph above the last would be easier to follow with fixed axis. I am looking forward to you adding a body to your simulation. Good job in commenting on your code.

Comments are closed.