# Initializing positions of $n$-body simulations

I'm working on an $n$-body simulation project, and I have a very basic question: How does an $n$-body simulation start?

In the script I'm working on, there is a range of forces defined, but they all depend on the distance between objects, and that's where it breaks for me, how does the simulation begin? To define a force between two objects, you need to know the distance between them, but their position is exactly what you're trying to simulate--that is what you don't know. My understanding is off somewhere, anyone care to elaborate?

• You will have to give an initial state to the system, i.e., define the initial positions. The final position should ideally be an equilibrium position (in absence of external force field). – Kalpak Gupta Mar 8 '17 at 10:44

Generally you initialise the simulation using a set of particles (or bodies) with some predefined positions and velocities (and masses).

Depending on what system you'd like to simulate you could use a specific setup, for example: put a more massive, stationary particle in the centre with less massive particles further out and with velocities that are perpendicular to their position vectors (like a planetary system). You could also generate (random) positions and velocities.

Given the gravitational potential distribution $\Phi(r)$ you are trying to model (e.g., any of these for a galaxy), you can generate the positions by getting the cumulative distribution function (CDF) of this: $$\Psi(r)=\int_0^r\Phi(\rho)\,\mathrm d\rho$$ You can then pick your positions by inverting the relation and noting that $\Psi(r)\in(0,1)$ (i.e., $\Psi$ is your random number). Then when you draw from $\Psi(r)$, the resulting distribution will follow the potential $\Phi(r)$.

Sadly, distribution functions are not always so nice as to be easily invertible, so you have to use other means to draw from the potential (e.g., the Newton-Raphson root finding I mention here).

To get velocities from these distributions, see this other answer of mine.

I am assuming you are trying to model a galaxy here. For a solar system like simulation, it wouldn't be that much different from the two body case I discuss here.

When you approch to a numeric simulation of N-body there are some basic concept to take into account:

1. Choose the more suitable integrating algoritm for your problem;
2. Ask yourself if the problem has a kind of correlation length;
3. Start from phenomenon observatorions;

In deatail:

1. The choise of the suitable algorith is the first step in this kind of problem. For an evolution on small timescale is better to use algorith like Runge-Kutta (Not symplectic, NOT conserves energy). For large timescale, and large number of evolution steps, is better to use algorith like "Leap Frog" (symplectic, so conserves the system energy during all the simulation).

2. If the interation between particles has a kind of cutoff that permit you to indroduce a kind of correlation length, so you can provide your algorith with a multithreading structure (using MPI or CUDA) to cut the simulation's time.

3. Usually the main goal of a simulation is to predict the evolution (understood evolution as positions of particles) of a initial state almost always obtained from a observation. If you want to predict the position of a particular asteroid which will have among 30years, you must use a model of solar system where the initial position of the planets (usually only the major planets) is obtained from an observation. The time evolution of this configure (a kind of snapshot of our solar system today) give as result the position of your asteroid in 30years.

Finally is obvious that the evolution's algorithm depend to the position of your particles. But the initial positions are defined by your state that you are interested in investigating.