Overview and set-up
In this homework, you will numerically evolve a two-body problem forward in time. The particles in this system evolve according to Newton's law of gravitation, and their motion is governed by a pair of coupled ordinary differential equations (ODEs):
m1d2x1/dt2 = F1,2
m2d2x2/dt2 = F2,1
Here the masses are given by m1 and m2, the position vectors are given by x1 and x2, and the force vectors are given by F1,2 and F2,1.
Note that the forces are equal and opposite, F1,2 = -F2,1. We will simplify our notation by dropping the force subscripts and setting F = F1,2. We can write the gravitational force as follows:
F = F1,2 = Gm1m2/||x2 - x1||2.(x2 -x1)/||x2-x1||= Gm1m2(x2-x1)/||x2-x1||3
We can also simplify our lives by taking the gravitational constant to be one, G = 1. This is equivalent to making a change of variables, and hence doesn't effect the dynamics in a meaningful way.
It turns out that this system is much easier to solve if we split it into a set of first order equations. We can do this by defining the momentum, Pi, of the ith particle as:
Pi = mi dxi/dt = mivi
With this we can split our previous second order ODE system into this equivalent first order system:
dx1/dt = p1/m1
dx2/dt = p2/m3
dp1/dt = F
dp2/dt = -F
which could be solved by hand, but it would not be very fun. We will solve these equations for motion in the plane. This means that each position and momentum vector has 2 components, and this is really a system of 8 equations!
This may seem a little daunting, and might seem like jargon soup. Don't worry though! We'll go through all of the relevant details here.
Discretizing the problem
At present, this is a continuous problem which the computer can not solve for us. To make this tractable we will consider only a discrete set of times. The details of this discretization process are important, but we'll delay them for now. We'll consider the discretization of ODEs in more detail in the final projects at the end of the quarter. In this discrete setting we will not consider t ∈ R, but rather a discrete set of times given by:
tn = nΔt
where the superscript is just a label (not an exponent), and Δt is the time step size. As a result the positions and momenta of the particles will only be known at these discrete times.
If you are interested in where the update formulas that we'll use come from you can read more on the Verlet integration wikipedia page. We are using the Velocity Verlet method which is particularly clean for this problem. This method is part of a larger class of methods called symplectic integrators. These methods are intimately tied to the underlying dynamics of Hamiltonian systems, and get used very frequently.
Solving the problem computationally
Ultimately, our goal is to evolve equation (7), paired with some initial conditions, forward to some final time T. To accomplish this task you will fill in two Fortran 90 source files: orbits.f90, and timestep.f90. Before proceeding, go to your git repository and create a hw2 directory with a code subdirectory. Then download the following 5 files to the newly created code directory:
Makefile
utility.f90
timestep.f90
orbits.f90
plotter.py
The orbits.f90 file will contain your main program. This program will define the problem data, call the time stepper to walk forward in time, and write the numerical solution out to a file.
To aid this process and to keep everything organized, the timestep.f90 file will hold a subroutine for producing the next state of the system given the current one. It will also hold an internal private function to help calculate intermediate quantities, which helps avoid duplicating code.
The utility.f90 file currently holds only the kind parameter for double precision reals. Later we'll have more items to put in there. The Makefile will help you compile everything without needing to call gfortran repeatedly. Finally, the file plotter.py contains a small plotting routine to visualize the data generated by the Fortran code. You may use this if you already have Python 3 installed, along with NumPy and Matplotlib. If you don't have those then don't worry. You can generate plots using excel easily enough for now.
The orbits.f90 file
This is the source file that defines the main program. This file is responsible for the following items:
Defining number of time steps and the final time
From these it will set the time step size
Declaring arrays to hold the solution for all time
Set the initial conditions for the system
Evolve the system by calling into the timestep module
Write the solution out to a file
A few of these have already been done for you. The parts you need to fill in are indicated by comments using !!!, with numbers matching the steps that follow.
Declaring variables
We need variables to store the final time and size of the time step to use, as well as a few arrays to hold the solution to the system. To this end declare the following variables:
tFinal and dt as scalars of type real
mass as a rank 1 real array with 2 elements
pos and mom as rank 3 real arrays of shape
After that, set the final time to be tFinal = 50 and the time step size to be dt = tFinal/nSteps. Make sure that everything is using the fp kind parameter imported from the utility module.
We should also discuss the storage convention for the pos and mom arrays. These are rank 3 arrays, and the index for each rank indicates a different thing.
Consider the position array. The first index corresponds to the x and y coordi- nates of the particles, and since we are working in the plane this rank must have a length of 2. The second index labels which particle is which, and since our system consists of 2 particles this rank must also have length 2. The third and final index corresponds to the time step where position data is known. We could read the fol- lowing indices as such:
pos(1,2,10) is the x position of the second particle at the tenth time step
pos(:,1,1) is the position vector for the first particle at the first time step (note the slicing)
The momentum array is stored similarly. The only difference is that the first index refers to the x and y components of each momentum vector.
Setting initial conditions
Navigate to the bottom of the file to find the set_ics subroutine. Inside here you will set up the initial configuration of the system. Let's set the first particle to have
mass m1 = 1 and the second to have mass m2 = 0.01. We could do this in one line by writing:
mass = (/1.0_fp,0.01_fp/)
Next set the first particle to lie at the origin with zero momentum. From our dis- cussion about storage order above, you could set these as:
pos(:,1,1) = (/0.0_fp,0.0_fp/)
mom(:,1,1) = mass(1)*(/0.0_fp,0.0_fp/)
Note: you will later change the momentum of this first particle, so you should go ahead and write it this way even though it's really just a fancy way to multiply by zero for now.
Finally, set the position of the second particle to be x2 = (0, -1). Give it an initial velocity of v2 = (1, 0) recalling that the momentum is p2 = m2v2.
Compile the code and try it out
You can now compile and run the code, though it won't do very much yet. Call make inside your code directory. This should compile without errors and produce orbits.ex. Try running the executable by calling ./orbits.ex. Inspect sol.dat and verify that the first line holds your initial conditions.
Question: What order are these columns in? Can you figure it out from the write_data subroutine?
The timestep.f90 file
This file declares the public subroutine take_step and the private function compute_acceleration. Your tasks here are to fill in the bodies of each of these.
Computing acceleration
The force that each particle exerts on the other depends only on their positions, and not on their momenta. You can see the definition for it in equation (5) above.
Navigate to the bottom of the file and start filling in the compute_acceleration function. Note that the input arguments and the output have already been declared. We will break this into three sub-steps:
Declare local variables
In addition to the input/output variables that have already been defined, we will need a couple local variables to make the calculation cleaner. Declare a real scalar called dist and a rank 1 real array called force with a length of 2.
Find distance between particles
Take note of the shape of the input pos variable. This time it is only a rank 2 array, where each rank has a length of 2. Here we are only taking in the position at a fixed time, so the final rank is no longer present.
The first particle has position pos(:,1) and the second has position pos(:,2). Calculate the distance between these positions and store it in the dist local variable that you declared previously.
Convert force to acceleration
Go ahead and uncomment the line that has the force calculation on it.
The take_step subroutine that we are going to write later needs the acceleration of each particle. Recall that Newton's second law says that the acceleration is related to the force as . Since our particles have different masses we need to be a little careful how we fill the acc variable.
Fill acc(:,1) with the acceleration of the first particle as given by its mass and the force that was found previously. Fill acc(:,2) with the acceleration of the second particle. Be careful with the direction of the second acceleration!
Updating position and momentum
This is the moment we've all been waiting for. Navigate to the take_step subroutine. Here is where we will figure out what the next state of the system will be given its current state. This time all input and output variables have already been declared, as have all of the necessary local variables. Again we will split this into 3 sub-steps.
Get velocity from momentum
The position update is simpler if we already know the velocity of the particles. Set vel(:,1) to the velocity of the first particle by using mass(1) and mom_old(:,1). Do the same for the second particle, but note that because of the different mass you'll need a second line to accomplish this.
Update the position
We now have the velocity and acceleration of all particles at the current time step. We have also been careful to pull off all of the mass-dependent terms, so we can use the nice array arithmetic of Fortran to accomplish this in one line.
The new position is set by assuming that the velocity and acceleration are constant up to the next time step. Then from regular calculus we find that the new positions are given by:
xi(new) = xi(old) + Δtvi(old) + t2/2ai(old)
Use the pos_old, vel and acc_old variables to set the new position, and store it in pos_new.
Recalculate the acceleration and velocity
We are halfway done with the time update. So far we have set new positions using the old positions and the old particle velocities. We now need to find updated velocities. We could try to do this trivially by using just the old acceleration value and the time step size, but this turns out to be a bad idea.
Instead, start by calculating acc_new by calling the compute_acceleration function again, but this time passing in pos_new. Then update the velocity like so:
vi(new) - vi(old) + Δt/2(ai(old) + ai(new))
This is slightly different in that we are using the average acceleration between the old and new time steps to update the velocity. This seemingly simple change actually has a pretty drastic impact on the final result.
Compile and run the code
You have now finished all the hard parts of the code! At this point you should be able to compile without errors (or warnings :0) and run the code to produce more interesting solution files. The next section will help you verify that you are getting the right outputs
Your report should discuss briefly your actions in implementing the code. Were there any unexpected issues? What was your debugging strategy when an issue arose?
Your report should comment on the following items:
1. Include the commit hash for your final code submission.
2. Comment on why real (fp) can be used throughout the files orbits.f90 and
timestep.f90.
3. Comment on the separation of duties between o rbits.f90 and timestep.f90. Does the distinction make sense?
4. Right now the code is limited to handling two particles. What would you need to do to increase this to three particles? What about N particles?
Don't actually do this! Just comment on what you would do.
5. The code also evolves all of the dynamics in the plane. What would you need to do to raise this to dynamics in full 3D space? Would this be easier or harder than in- creasing the number of particles?
6. Present each of the three configurations discussed above. Present the initial condi- tions that give that solution and comment on what you see.
For the first and third configurations find the total momentum of the system. Comment on how this relates to the trajectories you observe.
If you adjusted the final time and/or number of time steps in the third con- figuration, be sure to comment on that as well.
7. Inside the Makefile flags are specified that control the compilation. What are these flags and what do they do?
8. Write a conclusion to wrap up the whole report.
Attachment:- two-body problem.rar