The program you will be working with is one of a class of programs called \(N\)-body codes. These codes simulate the motion of compact, massive particles, moving under the influence of their mutual gravitational attraction. \(N\)-body codes are big business in astrophysics, where they are used for diverse purposes, perhaps most notably cosmological simulations of galaxy formation. The Illustris simulation project (http://www.illustris-project.org/) is a recent suite of simulations that modelled the motion of billions of particles from the beginning of the Universe until close to the present day. A similar project currently ongoing here at the University of Melbourne is the DRAGONS project, led by Professor Stuart Wyithe (http://dragons.ph.unimelb.edu.au/index.html).
The program you will be using is somewhat cruder than the state of the art, but the principle is the same as for all \(N\)-body codes.
Warning! Do not attempt to solve the equations on paper – there is no general solution to the three-body problem, and even limited, special-case solutions are hair-raisingly difficult. If you find yourself writing down pages of equations, you’re doing this the wrong way. All the maths you need can be written in about three lines!
How can we simulate particles moving in orbit around each other? The physics is simple; we use Newton’s inverse square law of Gravity,
where \(\bf{F}\) is the force, \(G=6.67\times 10^{-11}\) Nm\({}^{2}\)kg\({}^{-2}\), and \({\bf r}=r\times\hat{\bf{r}}\) is the vectorial separation of the two objects. Assume throughout that all the objects can be treated as points; for example if the Earth is in orbit around the Sun, assume that all the mass of the Sun is at the centre of the Sun for calculation purposes. We proceed as follows:
We choose initial positions and velocities for all our objects.
Using Newton’s law (above) we work out the force on each object due to the gravitational attraction of the others. The force is a vector quantity, and as gravity is linear, the total force on each object is the vector sum of the forces on it due to each of the other bodies.
Using the force we just calculated, we work out the acceleration of the body.
We now take a small step forward in time, and compute the new position and velocity of all the objects. The new position is simply the old position plus the product of the velocity and the time step. The new velocity is the old velocity plus the product of the acceleration and the time step.
Repeat steps 2—4.
Simple! The key to computer programming is to remember that computers are pretty stupid. They can only do simple things, but they do them very fast. Break your problem into tiny little pieces, and get the computer to do them one after another very fast.
“How big should our time step be?”. “How long should we run our code for?”. These are the most common questions demonstrators get asked, and there is no easy answer. In the algorithm above, you’ve approximated the true, curved trajectories of the bodies with a series of straight lines. Every time you take a step, your final position will be slightly in error, and if these errors are too large they can compound one another until your program gives you complete garbage.
It is difficult to know in advance what time step to use for a numerical problem. The only way to know is to see how changing your time step affects your results. For example, in section 4 below you will simulate the Earth going round the Sun. If your program doesn’t start with the Earth in orbit with a period of about a year then you may have set the code up incorrectly, but if your orbit starts out fine and then starts to go haywire you most likely need to reduce your time step. If you decrease the time step too far, however, you will never be able to run enough orbits to see if your simulation is working properly. By adjusting the time step in a careful and systematic fashion for the Earth-Sun system, you should be able to find a happy medium between accuracy and computation time. Unfortunately, a time step that works for the Earth and the Sun won’t necessarily work for all the other scenarios in this exercise, you will need to perform the same analysis in each case and include them in your lab report.
What will happen if your objects get close to each other? If the distance they travel in a time step is comparable to the distance between two objects, you’ve got trouble (why?). There are ways to fix this in the code, but the simplest way out is not to let the objects get this close; if you have to have close encounters in your simulation, use very small time steps.
For advanced students: there is a simple modification you can make to the above program to greatly reduce systematic errors…
Step 4 above is the most obvious way to step positions/velocities forward in time, and is known as Euler’s method. Leonhard Euler was a famous Swiss mathematician who died in 1783; in the subsequent 250 years more sophisticated techniques have been developed. One such technique which is very popular in simple applications is the fourth-order Runge-Kutta method (RK4) – if you have taken the Computational Physics course you should be familiar with Runge-Kutta methods. If you have some experience with programming, try to adapt the nbody.py code described below to use RK4.
In all the projects below, you have a wide choice of free parameters. What are the starting masses? What are the starting positions and speeds? You may feel that you have been given far too little guidance about these things. You are right, but this is a problem professional astronomers face all the time. There are no definite answers, what you have to do is explore parameter space. There will be some plausible range of input parameters. Explore it! Try out different values; home in on the most interesting ones. This is an important and difficult skill, and a lot of marks will depend on how well you do (and document) it.
You are probably used to working in S.I. units. These are fine for things of order 1m long (like you and me), but for your average star (mass \(\sim 10^{30}\) kg), you’ve got problems. Most computers can handle numbers as large as \(10^{\pm 34}\), but all you would need to do is multiply the mass of the Earth and the Sun together, and the program would crash.
The simplest solution is to invent new units for mass, length and time. If you choose these units correctly, your program will be dealing with numbers like 1, rather than \(10^{30}\).