Author Archives: nahendel

Matteo & Nadav: Final Conclusions

For our project we decided to model the perturbations of Earth’s orbit caused by Jupiter. We chose to study the effects of Jupiter because it is the only mass nearly comparable to the sun (about 0.1%). We started with a routine that modeled a single planet orbiting a fixed Sun, and built up to a true three body simulation with the Earth, Sun, and Jupiter bound gravitationally.

Our end goal was to write a three body simulation for the planets because we planned to observe the effects of an increased Jupiter mass on Earth’s orbit. At a certain point the “fixed-Sun” model becomes an unphysical approximation. Three body problems are nearly impossible to solve besides restricted problems, but quite possible to observe computationally.

Routine I: Planet Orbit

Our first model just contained a single planet orbiting around the sun. Mathematically, this model only relied on the inverse square law of gravity to relate the two bodies. Since we know from the physical system that sun is far more massive than any of the planets, we know that they have a negligible effect on the sun’s orbit. Because of this we simply fixed the sun at the origin of our system, and only modeled the movement of the planet around it.

We chose the Euler Cromer method as our primary tool because the orbit of a planet is an oscillatory motion and energy needs to be conserved over each period. This is comparable to the pendulum problems from earlier in the semester where we ran into the same issue. With the Euler method energy is added each period, so we demonstrated here that the Euler-Cromer method conserves energy with each orbit.

Euler Cromer Euler

Thus, the only thing affecting each orbit is actually the initial conditions, we are essentially approximating all of the orbits except mercury as circular because their eccentricities are very close to zero. Due to this nearly circular motion, the gravitational force is equivalent to the centrifugal force of a body in uniform circular motion and we are able to calculate the initial velocity for the planet when it starts at a distance of its semi-major axis from the sun. We will observe circular motion unless the initial velocity is wrong.

Routine II: Jupiter-Earth

Already having completed a foundational program for a gravitationally bound singular planet simulation, it was not much of a reach to add in a second planet for our second routine. The exact same physical equations were used and we continued using the Euler-Cromer approximation method. Initially, once Earth and Jupiter are related gravitationally, we observe no noticeable difference in the Earth’s orbit, which matches our physical observations. This program is very useful however because it allows us to adjust Jupiter’s mass in order to observe how the Earth’s orbit perturbs without affecting the sun’s motion. While this was very useful as an initial exploration into the perturbations of Earth’s orbit caused by Jupiter, it was not complete because it still didn’t take into account the motion of the Sun when Jupiter’s mass was very large. This was the goal of our next routine.

Routine III: Three-Body

The difference between the three-body routine and the earth-jupiter routine is that we no longer fix the sun in place and instead calculate the effects of Jupiter and Earth on its motion. We also made several changes to the initial conditions to make the model work, such as changing the origin to be the center of mass of the system and giving the sun an initial velocity so that the center of mass does not drift. We also set it so that the initial velocity of sun changes with increase in the mass of Jupiter in order to maintain these conditions.

3body 100xjupiter mass

Observations:

These simulations demonstrated that the Earth’s orbit is very stable even up to many times the current mass of Jupiter (as would be expected from observational evidence).

Upon observing the effects of changes in Jupiter’s mass for this system, we were able to conclude that the final three-body system becomes chaotic at sufficient Jupiter masses. Small changes in Jupiter’s mass caused entirely distinct orbits to form. In many of these systems, the motion appeared completely random, sometimes resulting in the Earth’s ejection from the system. Just like the pendulum, below a certain driving force the pendulum’s motion was still predictable, and above a certain force the motion became chaotic.

This threshold for chaotic motion was clearly demonstrated when the mass of Jupiter approaches the sun’s mass (x1000). The system starts to look very much like a binary star system around a center of mass. In this scenario earth’s orbit is VERY unstable and chaotic, completely unpredictable, and extremely sensitive to initial conditions. When Jupiter’s mass is sufficiently small to never directly interact with earth (AKA never crosses or comes near the path of the sun’s motion) then the motion of Earth is not chaotic.

To download our code, please click the following image:

Matteo & Nadav: Week Two

Nadav Hendel & Matteo Bjornson

Week Two

The celestial two body problem can be solved exactly, allowing us to observe Kepler’s Laws in action, but the same is not always true of the three-body case. Thus, solving the three body problem exactly has been a central problem of celestial mechanics. In this section, we model the effect of adding a third planet into our system. We will use Jupiter, since it has the largest mass of the planets and therefore would have the most significant effect on the Earth and Sun. The planetary motion will still be related using the inverse square law for gravitational force, and the sun will still be fixed to the origin for our first routine Jupiter-Earth.

Jupiter-Earth is very similar to our first routine Planets, in terms of methodology. For each time step, the planets positions are calculated using the Euler-Cromer method. This is completed by keeping track of the distances between both planets and the sun, as well as the interplanetary distance. From these distances, we can compute the acceleration of each planet as related by the inverse square law, and from there, the new positions. The code for this program is quite dense, and as noted in the textbook, it can be more efficient to define some recurring variables at the beginning of the code such that they need not be recalculated during each iterative loop. Despite this we opted not to do this in our code for the sake of clarity since this added efficiency is quite small.

From our initial results running this program, we can see that using normal masses for Earth and Jupiter, stable planetary orbits can be achieved. At this point, we will continue our experimentation by observing how a change in Jupiter’s mass would affect the Earth’s orbit. As an initial point of reference, a small increase of ten fold of Jupiter’s mass yields negligible results, as Jupiter’s mass is still dwarfed by that of the sun.

In the textbook, a mass increase of 1000 fold is supposed to yield much larger changes in the earths orbit. Unfortunately, we were not able to achieve this result exactly, and were only able to note significant changes using a mass of closer to times 9000 Jupiter original mass. The results of this are shown below:

The fact that we had to increase the mass of Jupiter to 9000 to observe noticeable changes in earth’s orbit indicates that we may still have some bugs in our code that need correcting. According to Giordano, these effects should have been visible at 1/10th that mass.

Now that we have established this model, we also wished to explore what a true-three body system would look like with our second routine. As the mass of Jupiter approaches that of the Sun, or at the very least significant in comparison, it is unreasonable to simulate this system with the Sun fixed at the origin. The difference between the 3-Body Orbit routine and the Jupiter-Earth routine is that the center of mass is now set as the origin, and we also incorporate the movement of the sun as it is affected by the enlarged mass of Jupiter. The Sun is given an initial velocity so the net initial momentum is zero, keeping the center of mass from leaving the origin. A few interesting things can happen to this system when the initial conditions are adjusted.

The first change we attempted to make with the model was to establish a base case, simply by multiplying the mass of Jupiter by 100. An increase of mass of this magnitude was intended to perturb earth and the sun, but not by much. As you can see below, this change actually had significant effects on both the sun and Earth’s motion. Despite this, the Earth’s motion with respect to the sun still appears very regular.

3body 100xjupiter mass

One of the first things we noticed was extreme sensitivity to initial conditions. In the following two gifs, the mass of Jupiter is multiplied first by 710x, and second by 710.2x. As you can see, this difference of .2 yields a largely different path for the Earth. In fact, it seems that the results from this section are much closer to the chaotic regime than expected due to their sensitivity. This observation might have some insight into why these systems are so hard to calculate in general.

In exploring the effects of different initial conditions, we set the initial velocity of the Sun to zero. We observed the Sun to follow a cusped path, which seemed almost unnatural. Is this even possible? This also made us realize we need to establish a deeper understanding of the implications of giving the Sun an initial velocity.

something wrong with 3body sim

These investigations have revealed another discussion we were not expecting to find: how do we know these simulations are correct? Or, how do we know the effects we are observing are not artifacts of the numerical approximation or specific lines of code introducing error? Looking forward to next week we need to investigate ways we can can validate our code and results in order to find them meaningful. It is important to have a fundamental understanding of the underlying physics relating to a theoretical model otherwise there will be no intuition regarding the correctness of the outcome. Our initial exploration of these models yielded powerful results, yet we hope to finalize this project by reaching a greater conclusion about the nature of computational approximations of multi-body systems.

Click the two Matlab Icons below to download the Jupiter-Earth and 3-Body Orbit routines:

3-Body Orbit

Jupiter-Earth

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

3-Body Celestial Mechanics Project Plan (Nadav and Matteo)

Purpose and Goals

One purpose of this course is to explore the use of computational methods to find solutions to otherwise difficult analytical problems. The 3-body system is a prominent example. It is relatively straightforward to solve the equations of motion of a 2-body problem analytically, but far more difficult to do so for 3 bodies. The introduction of computational methods have allowed for fairly straightforward ways of analyzing 3-body motion.

A prominent example is that of the Sun, the Earth, and Jupiter, bound by the force of gravity. Since Jupiter is about 0.1% of the solar mass, it is the next biggest influence on Earth’s orbit in the Solar system. Using computational methods to analyze Jupiter’s effect would be far more efficient. We seek to create a model of the orbits of the Earth and Jupiter bound to the sun–and interacting with each other–via the inverse square law of gravity.

For this project we plan to make use of sections 4.1 and 4.4 in Giordano’s book. 4.1 outlines a method for calculating the motion of a planet orbiting the sun, whereas 4.4 opens the model to include Jupiter. The program will generate position data for each planet given their acceleration due to gravity, using the Euler-Cromer Method to approximate each new position based on this acceleration. We will start with a simplified model, with the Sun fixed to the origin, then graduate to a true 3-body model, where the origin would be fixed at the system’s center of mass.

Physics

We start with Newton’s Law of Gravitation, where Ms is the solar mass, Mp is the planet Mass, Fp is the gravitational force experience by the planet due to the Sun, and G is the gravitational constant:

Fp=G*Ms*Mp/r^2

Solving for acceleration (given Newton’s 2nd Law: Fp=Mp*a), or d^2x/dt^2:

a=G*Ms/r^2

This solution for acceleration is a second order differential equation, appropriate for use with the Euler Cromer Method.

Foundational Code (Routines provided in Giordano)

Planets: A routine to model the orbit of a planet around the sun as described by the inverse square law of gravitational force, with the sun at the origin. The Euler-Cromer method is used to approximate position. The Euler Method by itself would not be sufficient, as with all oscillatory motion it would introduce energy into the system and the planet would move away from the Sun. The Euler-Cromer Method guarantees a conservation of energy over one orbit.

Jupiter-Earth: Expands the above system to model Jupiter and Earth orbiting the sun. This method will also use the Euler-Cromer method for the same reasons. The suns mass will affect the orbits of the other two planets but the sun itself will sit at the origin of the system and not move.

True 3-body: This model is similar to the Jupiter-Earth model but instead uses the center of mass of the system as the origin. This model will be used when we are altering the mass of jupiter to be close to that of the sun, therefore causing changing orbits for each planet in the system. In this case, we want the sun to be able to move to better describe the full motion.

Initial Resources

http://astro.geo.tu-dresden.de/~klioner/celmech.pdf

http://arxiv.org/pdf/astro-ph/0411043v1.pdf

http://www.iumj.indiana.edu/docs/27030/27030.asp

Computational Physics, Giordano and Nakanishi. Sections 4.1 and 4.4.

Timeline

Week 1 (4/6-4/12)

During the first week we hope to work on organizing sources and solidifying understanding of background materials and physics. We will also use this week to read the Giordano text to begin understanding the setup of the MATLAB routines. We will begin writing the code for the Planets routine. This routine does not directly inform the 3-body problem but is the foundational code for the project. For the entire first week we will split the workload together since this work constitutes the base understanding of our project.

Week 2 (4/13-4/19)

Since the Planets routine should be well defined and operational by this week, we will use this time to finish off the other routines. This should be manageable since both the Jupiter-Earth and True 3-body routines are derivations of the Planets routine. Nadav will work on the Jupiter-Earth model while Matteo tackles the True 3-Body routine. This week will mostly be coding based, as we hope have all models complete by this time.

Week 3 (4/20-4/26)

By now, all routines should be complete. We will use this week to explore the implications of our models. This exploration will consist of multiple demonstrations of the effects of the changed mass of Jupiter on the system and will also include qualitative conclusions. We will also begin to work on writing up our results, with each of us taking charge of specific simulations.

Week 4 (4/27-5/3)

This final week will consist of finishing up our written results and finalizing the report. We will also leave space here for potential further exploration of our system, given appropriate time. Possible extensions include adding in another body to the system or analyzing in more detail the stability of orbits (Section 4.2, Giordano).

Planetary Dynamics: The Jupiter, Sun, Earth 3-body Problem

The advantage of using a computer to find a numerical solution rather than an analytical solution is quite clear when considering the relationship between Jupiter, the Sun, and the Earth. Rather, it is impossible to achieve an analytical solution for most celestial mechanics problems involving 3 or more bodies, thus numerical approaches are the only alternative.

Our project begins with modeling the Sun/Earth/Jupiter system with the Sun fixed at the origin. The three bodies are related via the inverse square law, and the positions of Jupiter and the Earth are calculated using their respective equations of motion (the sum of the forces from the other two bodies). Once a working model of the planetary orbits is established, it is worth investigating the effects of Jupiter on the Earth’s orbit when Jupiter’s mass is increased. We know that the current orbit is stable, but at what point does it become unstable?

If the mass of Jupiter is increased by roughly a thousand fold, it actually becomes comparable to the Sun’s mass. In this case, the model can be expanded to approximate a true 3-body system, with the origin at the center of mass of the system instead of fixed at the sun. The advantage of this model is that it allows us to observe how all three planets are effected by this change in mass. How large an effect, and at what point the orbits becomes unstable, remain open to investigation.

Modeling and Experimental Tools with Prof. Magnes

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