Programming for Physicists: Exercise 5

Exercise 4: Basic program units

Introduction:

This exercise is based on an old problem, the three body problem. How does three bodies move in a mutual gravitational field? This is a complicated problem to solve accurately, as many specific situations have to be taken into account. In some case, when two of the objects have a close encounter part of the calculations can be done as a two body calculations, because the motion of the distant partner happens on a much longer time scale, or you have a collision between two of the objects... but most of the time the distances are such that all three bodies have to be treated equally. We are not going to deal with any of the complicated, but very interesting, items here. By taking the problem in it's simplest approach, it provides a good possibility for making use of the different building blocks in a code: Main, Subroutine, Function and Modules. The exercise is setup to train you in these skills, more than requiring you to solve the problem in all its glorious details, while using what you have already learned. Naturally, you will be expected to deliver the goods and make a fully working program! To do this you first have to understand the problem in detail, then you can start dividing the problem into smaller building blocks, of which some may be reused in different branches of the code. Therefore, read the instructions carefully to understand what you need to do. Start with planing how the puzzle fits together. When all of this is in place, then it is time to start the actual programming -- and not any earlier!

Have fun!

Getting the required files:

Download the tar file to your ~/Dat_F directory and unpack it:

> cd ~/Dat-F
> tar -zxvf fortran4.tgz
> rm fortran4.tgz

This creates a new directory, "Exercise_4", that contains the files needed for the project.

Some find it difficult to use the "initial" files provided for the exercises. You don't need to use them. Feel free to start writing the programs from scratch. In fact it is the best way to learn all the different small parts need for a full program.

The three body problem:

How does the position and velocity change of three bodies that are initially bound by their mutual gravity? This problem is to complicated to have a general analytical solution, and a lot of effort has been invested over the years working on solutions to the problem. Initially by hand for a number of special cases and more recently by the use of computers. To achieve a robust and fast solution to the problem for all possible situations one has to develop a rather complex program. That is far to complicated for this course and we are going to take a simple approach. By doing this, it is relatively easy to write a code that can solve the problem. The trait off, naturally, is that we are not guarantied to obtained the correct solution as time progresses, nor that the method is suited for getting a reliable solution. Instead we setup a problem for which it is possible to write a simple code that can solve the relevant equations which gives an impression of the dynamical behavior of the three objects.

The background for the problem is simple, the governing force is the gravitational attraction between the objects. The force on one object is simply given by the sum of forces from the other objects:

(1)

Here m_i is the mass of the object, r_i the position vector of the i_th object, G the constant of gravity and the variables with index j represents the remaining two objects. In the following we use A_i as the acceleration of the object. The acceleration is obtained from eq.(1) by dividing the right hand side by m_i.
From a more detailed analysis of the problem, it has been shown that to determine the initial condition of the three dimensional three body problem unambiguously, seven parameters are needed for each object. The choice, for the representation we are going to use here, are three values for position vector, r, three values for the velocity vector, v, and the mass, m.

To solve for the time evolution of the objects, we take the Taylor expansion of the position vector, r, in time around the position r. This gives the following expression for the position vector at the time delta t later,

(3)

Here we have written the 2ed and higher order derivatives of r in terms of the acceleration A, and only included terms up to fourth order in time. Notice, that the second term, the first time derivative of r, is equivalent to the velocity and the fourth term, the third time derivative of r, is the first time derivative of the acceleration. Using eq.(1) written in the from as an acceleration with the following substitutions

and

the acceleration, A, can be written as;

(4)

By taking the first and second order time derivatives of eq.(4), we obtain an expression of the third time derivative in eq.(3);

(5)

and the fourth order time derivative from eq.(3);

(6)

Ignoring the higher order terms in the Taylor expansion we could try to solve this equation numerically. The problem is that it only gives a description of the change in the spacial position. To be able to advance the solution more than one time step, we need to include a description of the evolution of the velocity. This can be obtained by deriving a Taylor expansion of the velocity. For this we take the time derivative of eq.(3),

(7)

By solving eq.(3) and eq.(7) simultaneous we can advance the solution in time. Looking at the expressions for eq.(5) and eq.(6), it is clear that calculating these terms requires a lot of computing time, and for very detailed calculations we have to look for an approximation for these terms. One way to save CPU time is to make a simple cheat; first, only include the first order derivative of the acceleration in eq.(3) and secondly by substituting the expression in eq.(5) with a simple backward time difference term;

(8)

where the index i represents the time step number. Then eq.(3) can be written:

(9)

A similar substitution is done for eq. (7) where we are then ignoring the 2ed order time derivative of A. Depending on the relative position of the three objects, they move fast or slow in their orbit. The magnitude of the time step where the third order Taylor expansion provides a good representation of the time evolution therefore depend on the rate of change in the objects orbit. From experience it is found that a simple way to control the time step is to use the following expression;

(10)

where eta is an input parameter. This gives a small time step when the ratio of the amplitude of the acceleration is small compared to the rate of change of the acceleration. Notice that this value has to be calculated for each of the bodies and the smallest of these value of these values has to be used for the given timestep.

This is the end of the section containing all the equations. They are needed to give you the background for the problem. Don't worry, there may be a number of them and they do look complicated to program. The good news is that we only need a few of them and that the calculation can be broken into a number of small tasks. This makes it possible to write a program that calculates the trajectories of the three objects by advancing their positions in time by using eq.(4),(7)-(10) ignoring the second order derivative of the acceleration in eq.(7). Writing the program you must use modules, subroutines, functions and a single main program. Below is a list of key points in the process. Read the list and use it to guide you through the steps needed to write your code.

Look carefully at the equations, eq(4),(7)-(10) and the definitions of the various vector expressions. Find out how to break the expressions into smaller units that can be easily calculated.
The equations apply for an arbitrary number of objects. Make it possible to do the calculations with an arbitrary number of objects and not only 3!
For data IO you need to handle both input and output:
- As input you need to give the number of objects, the start position, velocity and mass for each object, a value for eta, the end time and an initial value for the first time step.
- The positions and velocity of the objects needs to be written to a data file after each N time integrations. The file needs to have a format that can subsequently be used for plotting the trajectories with gnuplot.
As you can see from eq.(10), it is required to know the time derivative of the acceleration for the first timestep to derive the new position and velocity. You could implement eq.(5), which takes nearly as much effort as writing the rest of the code. Initially you should ignore this higher order contribution in the first time step and use only up to the acceleration term in eq.(7) and eq.(9). For the following time steps use the expression given by eq.(8) as an approximation for the first order derivative of the acceleration.
Use subroutines or functions for each of the trivial low level calculations. This makes it easy to reuse the same calculation several times in the code and makes it easier to find and eliminate errors. Examples of this are IO, various algebraic expression, time stepping.....
Before you start coding, make a paper sketch of the code to make sure you have remembered everything, give variables and routines sensible names.
As we are not interested in doing real 3-body calculations, set G=1. This is not cheat, but a simple way to scale the equations. Numerically this can be of interested as numbers between 0 and 1 are much better represented than number with large differences in their exponents.
There is a skeleton Makefile, you have to change the names of the routines to the ones you choose. To compile the code simply type
$ make
If you have put the correct names into the Makefile it will take care of compiling the relevant parts of the code and link them to one final executable.

If you are not sure that the different parts are written correctly, then it is a good idea to make a simple test of each of the code components, This can be done by making a simple main.f90 program that feed appropriate data into the routine. Knowing the input you can, by hand, calculate the expected reply and compare with the result you get from the program. If it is wrong, you need to find out why and correct it. When you believe all parts work fine, then test the code using the following initial conditions, known as the Pythagorean problem:

parameter/object	A	B	C
x	1	-2	1
y	3	-1	-1
z	0	0	0
vx	0	0	0
vy	0	0	0
vz	0	0	0
mass	3	4	5

This gives a case that only evolves in the x-y plane, but as any three points can be used to define a plane, we have the freedom to rotate the coordinate system onto this plane anyway. Further more, with no initial velocities, these objects will remain in this plane for the entire run. In reality, the objects may have velocities in the 3rd direction and this will break the symmetry, making the problem a real 3D problem.

You can use the data.in as input for the run --- you have to be modified the order of the data to fit with your definition of input --- Try different values for especially dt0 - the first time step - and eta. You will notice that the solution is extremely sensitive to these values, showing that the approach we have taken here is not very robust... Even with a much better method, it is found that the behavior is in general chaotic. So even if we do not get the mathematically correct solution, it shows how complicated the orbits can become. One problem is the accuracy when treating very close encounters between two objects. If the time stepping in this phase is only a little to large, then the objects end on different (wrong!) orbits, and the plot is going to look quit different. This is why it is very difficult to make calculations representing very long periods of time and with high accuracy predict where the objects would have been or will be located. For doing this you need to use a method that carefully conserves the energy in the system.

This graph shows the orbits until t=8.75. For the purpose of checking the result, recreate this graph and call it orbit.png. The values for dt0 and eta are given in the data.in file. Make sure to include proper axis and further information for easy identification -- incl. your userid. If your result has large deviations, then there is a problem in the calculations -- this solution has been checked against a totally independent source.

I have reproduced the data plot and produced a png file called orbit.png Credits: 50/-1

This picture shows the orbits until the time t=72.2.

To get the last point for this exercise you have to run a script by typing

> csh get_info.csh

This creates a file that needs to exist before submitting the next question. This file contains information about the files you have written, indicating the structure of your program. It allows me to see that you have split the program into smaller units as you were asked to do!

The result of running get_info in the exercise directory Credits: 50/-1

Now, play with different initial values and see how dependent the result is to the initial conditions. For instance change the value of eta while keeping all other parameters constant...... In fact it is simple to get chaotic orbits. To see this requires lots of computing power, letting the objects orbit MANY times. One can plot the objects positions in a 2D plane each time they cross from one side to the other (same direction each time). This creates a Pointcaree map in which cyclic or chaotic behavior is easily observed.

Reference: Mauri Valtonen & Hannu Karttunen: The Three-Body Problem, Cambridge University Press, 2006

$Id: index.php,v 1.9 2009/09/19 15:36:52 kg Exp $ kg@astro.ku.dk