Numerical Methods Used

N-body simulations require two main numerical methods. One for time advancement, i.e. integrating the equations of motion, and another for calculating the accelerations of particles from their interactions.

1. Calculating the accelerations of particles

The simplest, most accurate but most computationally expensive algorithm to calculate the accelerations of particles is to calculate their forces on each other directly (direct PP method). This algorithm is ~O(N^2).

A fairly simple alternative is the Particle-Mesh (PM) method which relies on constructing a potential from all particles on a mesh and calculating their approximate accelerations from the potential. This method can be made ~ O (N ln(N)) by using FFT. But since it involves approximating the location of a particle to within a cell in the mesh, it is not accurate for simulating systems where close encounters occur frequently.

Our main goal is to try and simulate a globular cluster, so we will have many close encounters especially in the core of the cluster. Hence, our choice for the numerical method is direct PP.

Prototyping the Main Engine in MATLAB

A basic N-body code was developed with the particle-particle (PP) method and simple Euler time stepping. The code's time complexity was analyzed as below:

Being a PP code, its time complexity is O(N^2) and based on the runtimes recorded this code would take about 65 years to run for a million particles (10^6) and about 230 days to run for 10^5 particles which is our required threshold for good impression. So, a lot of work has to be done.

The ways to reduce runtime are:
1. Use a more accurate time-stepping method which allows for larger time steps with tolerable numerical diffusion
2. Change the time complexity to Nlog(N) using Particle-Mesh method instead of particle-particle. This method uses FFT and iFFT. With Nlog(N) complexity, assuming that it does not increase the time for each calculation very significantly, the current MATLAB code could complete a simulation for a million particles in 9.5 days! Using FFT might increase the time for every calculation upto 5-10 times, so may be 90 days.
3. Rewriting the Code in FORTRAN. Unknown speed up. I would expect it to be at least 5 times? May be.
4. Parallelization: Speed up of at least 2-3 times.

If we are able to implement all of these steps  (which is unlikely since we only have about 45 days and we also have to factor in time for studying for finals, homeworks and doing science with this code) we could reduce the runtime to 90/15=6 days for a million particles! This corresponds to 12 hours for 10^5 particles.

These numbers might be revised once we get the hang of the actual amount of speed up FORTRAN can provide.

Oh yes, and how can I forget some other stuff we have to include:
1. Variable time-stepping depending on the velocities and separations of particles
2. Collisions - Inelastic?

Lets go about implementing these improvements step by step.

Project: N-body Simulations of Globular Clusters or Dwarf Spheroidal Galaxies

Class: Galactic Astronomy
Semester: Fall, 2015


Project Goals: 

1) Build you own N-body simulator. 
2) Run a simulation with a large number of point masses, or “stars”, in an initially uniform distribution over a spherical volume, let the simulation evolve to convergence. 
3) Compare simulations carried out with different initial parameters for the stars, this might include stars of different masses - perhaps following some predefined mass function, or some initial velocities for the stars. 
4) Compare simulations where the stars initially have no net angular momentum, to simulations where the system of stars, as a whole, has significant angular momentum.

A quote from Prof. Lepine: "I would be really really impressed if you can do this for more than 10^5 particles" which means that this is our REAL goal!