This exercise is set up like a mini research project. You should work through the theory, and make sure you understand what you will be modelling with the program. Next, familiarise yourself with the n-body code provided, making sure you can get it to work. Then, do the exercises examining the motion of massive objects interacting via gravity.

Write down the basic equations governing motion – those for force, acceleration and so on. Write these in terms of vectors (using \(\vec{r_{1}}\) and \(\vec{r_{2}}\) as the position vectors of two particles), and separate out into \(x\), \(y\), and \(z\) components. Now rewrite them using discrete time steps \(\Delta t\). These are the equations that will be used in a computer program.

Draw up a flow chart of how an \(N\)-body algorithm will work. This does not mean write a program! Rather, write the steps required to calculate the necessary forces, positions etc. at each time-step. (This should be more detailed than the basic algorithm given in Section \ref{ssec-alg}) Also, think about what units to use – for instance, what would you measure the distance in? Why? You will need to calculate the gravitational constant \(G\) in your new units. Wait for a demonstrator to check what you have done before continuing.

Now look at the template code that is provided on the computers (this code is written in the Python programming language). Match up what the code does with your flow chart. How well do they agree? Make your own copy of the template, so that you can alter it as you go along.
Please don’t edit the template file!

Run the code using the iPython terminal:

> run nbody.py

Take a look at the outputs of your code. What is the motion of the Earth and Sun? Is this what you expect, given the initial conditions?

For your first exercise, find the speed at which the Earth rotates around the Sun, and set its initial velocity to this speed. Is the orbit stable? If not, what change to the initial conditions should you make?

Try to determine how accurate your code is. Give some thought as to how you can work out the errors in your results, and how you can make your program more accurate. In doing this, explore the parameter space a little and get a feel for the range of parameters for which things still work well. Things to try would be :

• •

  Change \(\Delta t\). (What is the physical interpretation of a large or small \(\Delta t\)?)

• •

  Try increasing/decreasing the initial velocity by a small amount and/or a significant amount (order of magnitude perhaps?) – what happens? Calculate the escape velocity of the Earth, what happens if your initial velocity is near (above and below) this value?

• •

  Find a way to work out how many orbits the planet is doing. – how many orbits can you get in a sensible time? (Make this independent of the initial parameters – it needs to be adaptable in case the orbit changes e.g. due to the influence of a third body.)

• •

  Change the number of orbits your Earth makes. Does it remain stable?

• •

  Check whether energy is being conserved by your program – what is the expression for the total energy in the system? If conservation is not taking place, explain why.

Finally, convert your program to a 3-body problem by adding Jupiter to the simulation. Repeat your investigation of 2-body dynamics for the Sun-Earth-Jupiter system.

See appendix for information on calculating the temperature at Earth (as this is required for both projects).