Computational Astrophysics: N-Body Exercise

Example: Stardard vs Advanced? We should change the structure of this lab.

Background

It is a remarkable and sobering thought that it is impossible to solve the equations for three or more bodies flying around in space under the influence of their mutual gravity. In special cases you can make useful approximations, but in general, the problem is insoluble.

This is a major menace in astronomy, which after all is the study of large numbers of objects flying around in space under the influence of gravity.

In the late 18th century, a partial solution was found. While it is impossible to calculate the orbits exactly, you can calculate a good approximation valid for a short period of time. You start by putting all your bodies in their starting positions. You then work out the gravitational force on each object at this time. Using this force, you work out the acceleration on each body. Using this acceleration, you work out the velocity of each body. Using this velocity, you work out where each body will be a short time in the future. You move each body to its new position. You then go through the whole process again, starting at the new position.

In this slow, tedious way you can calculate trajectories for even the most complex situations, though thousands of calculations are required. In the past, teams of students would sit at rows of desks, each calculating one tiny stage over and over again, these were the original computers. Though it sounds like a monstrous waste of time, the rewards were great  –  it was from calculations such as these that the planet Neptune was discovered.

Luckily for you, computers have evolved from rooms full of hapless students into the silicon chip-based devices we know and love. They are perfectly suited for this job; computers love nothing more than doing simple calculations over and over again. In this exercise, you will use a computer program that will do in seconds the same calculations that took some of the greatest astronomers of the 19th and 20th centuries most of their lives.

\(N\)-body Codes

\label{sec-nbody}

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!

Algorithm

\label{ssec-alg}

How can we simulate particles moving in orbit around each other? The physics is simple; we use Newton’s inverse square law of Gravity,

\begin{equation} \label{eq:newtongrav}{\bf F}=\frac{GM_{1}M_{2}}{r^{2}}\hat{\bf{r}}\\ \end{equation}

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