Lesson 1: Simulating orbits
This project/homework assignment is intended to introduce the students to the basic techniques of numerical simulations. The students will learn how to turn a system of continuous ODEs into a set of discrete algebraic equations that a computer can handle and then writing a computer program capable of solving these equations. This class and project will focus on a system of ODEs using orbital mechanics as our physical example.
The equations of motion for self-gravitating body masses
We are interested in a system of point masses moving under their mutual gravity. The diagram below shows the configuration of the system in two dimensions.
The time-evolution of the system is driven by the gravitational force only and can be described by the following four equations of motion:
where r_1,r_2 and v_1,v_2 are the positions and velocities of the two bodies with masses m_1 and m_2 respectively. If we set up our coordinate system on the x-y plane, the equations can be rewritten as:
This is a system of eight ODEs. Given the initial positions x1(0), x2(0), y1(0), y2(0) and velocities vx1(0), vx2(0), vy1(0), vy2(0) the solution of these equations x1(t), x2(t), y1(t), y2(t) and vx1(t), vx2(t), vy1(t), vy2(t) describes the time-evolution of the system in the position and velocity space respectively.
The solution to the Kepler problem
Johannes Kepler proved that in the limit when one body is much more massive that the other (e.g., the Sun and a planet) the orbits are ellipses with the more massive body at one focus. This is known as Kepler's first law.
This ellipse shown in the diagram above has a semi-major axis a and eccentricity e and is described by the following equation:
having an orbital period of:
Solving for Orbital motion on a Computer: The Leapfrog method
In the leapfrog method we discretize time as
This means that we use the velocity at time
t_(n+1/2) to find the position at time
t_(n+1) and the velocity at time
t_(n+3/2) as is shown in the following diagram.
If we know the position at time t_n, x1(t_n) we compute the position at time t_(n+1), x1(t_(n+1)) using the following equation :
However, we will use a simplified version of it which is
(hereafter equation (1))
where vx1(t_(n+1/2)) is the velocity at time t_(n+1/2). Similarly we do with
1(t_(n+1)) and y2
If we know the velocity at time t_(n+1/2), vx1(t_(n+1/2)) we compute the velocity at time t_(n+3/2), vx1(t_(n+3/2))
sing the following equation:
Similarly we do with vx2(t_(n+3/2)), vy1(t_(n+3/2)) and vy2(t_(n+3/2)) as follows:
In a simplified version this equation becomes
(hereafter equation (2))
Numerical modeling in Python
The leapfrog method is implemented in the function leapfrog_step() defined in hw1.py.
1) The initial positions and velocities are stored in an 8-element array called rv.
Organize this array such that
rv[0:2] contains the x and y position of body 1,
rv[2:4] the x and y components of the velocity of body 1 and similarly for body 2
(i.e., correct the first four #FIX ME! statements)
2) Update the positions to their values at time t1 (next two
#FIX ME! statements) using equation (1) (see figure below) but the velocity at time t0. the correct way would be to use the velocity at time n+1/2 but we ignore this subtlety in the current homework.
3) Similarly Compute velocities at time t1 using equation (2) (see figure below)
(next two #FIX ME! statements)
4) Store the new values back into the array rv (last
four #FIX ME! statements)
Solving for Orbital motion on a Computer: Choosing Initial conditions
where m is the total mass of the system . This is implemented in the the function r0v0
A convenient choice of initial positions and velocities is starting the system while at closest approach, i.e., at periapsis. The initial positions and velocities at periapsis relative to center of mass of the system are then simplified as:
where q is the mass ratio defined as q=m1/m2. This is implemented in the main function of our program hw1() that we explain below.
Solving for Orbital motion on a Computer: Running the program
The main function of our program hw1() :
takes as arguments the initial conditions of our problem, i.e., the masses of the bodies m1 and m2, the semi-major axis a, the eccentricity e as well as the maximum time until which the system is evolved. The time step for the function is set to a default value tstep=1.0e-3 so you can use tmax=nsteps*tstep where nsteps is the number of steps you want to evolve the system. Using different values for nsteps and trying different initial conditions run the program clicking on the green play button and calling in the terminal the function hw1(). For example, typing in the terminal
sets m1=1,m2=3,a=1,e=0.5 and runs the program for tmax=1=1000*1.0e-3. i.e., for 1000 time steps which is equivalent to one orbit.
The program shows in the plot produced the positions of bodies 1 and 2 in blue and green respectively as well as their relative position in red. The relative motion of the two bodies is an ellipse with the semi-major axis and eccentricity defined in the initial conditions verifying
Kepler's first law.