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?
 A: 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.
A: 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.
A: The $N$-body problem is an initial value problem, i.e., given a system in a given state at time $t=t_0$, what will its state be at an arbitrary time $t$?
For the $N$-body (gravitational) problem, the question would be: given a system with $N$ masses, each with a given position and velocity, where will they be at a time $t$ in the future?
Thus, to be able to determine (numerically) the positions of $N$ bodies interacting under gravity, their initial positions and velocities must be known.
Given these initial positions and velocities, you can then start to time-step the simulation forward using an integrator of your choosing.
As you have stated, the force between two masses interacting under gravity is determined by the distance between them. This is given by Newton's law of gravitation:
$$F=\frac{Gm_1m_2}{r^2}\tag{1}$$
In vector form, the force on body 1 due to body 2 is given by:
$$\textbf{F}_{12}=-\frac{Gm_1m_2}{r_{12}^2}\hat{\textbf{r}}_{21}\tag{2}$$
... where $\hat{\textbf{r}}_{21}$ is the unit vector pointing from body 2 to body 1.
This is trivially generalised to $N$ bodies, where the force on a body $i$ due to all other bodies is:
$$\textbf{F}_i=-Gm_i\sum_{\substack{j=1 \\ j \neq i}}^N{\frac{m_j}{r_{ij}^2}\hat{\textbf{r}}_{ji}}\tag{3}$$
... which can be divided by $m_i$ to give the acceleration experienced by body $i$:
$$\textbf{a}_i=-G\sum_{\substack{j=1 \\ j \neq i}}^N{\frac{m_j}{r_{ij}^2}\hat{\textbf{r}}_{ji}}\tag{4}$$
Given all of the initial positions and velocites of the bodies, the initial acceleration $\textbf{a}_i$ for each body will be known (exactly), and can then be used to approximate the velocity one time-step later. Conversely, the initial velocity of each body (which is also known) can be used to approximate the position one time-step later.
Using the Euler method of integration (probably the worst, but also the simplest to grasp), the time-stepping would look like this:
$$\textbf{r}_{n+1} = \textbf{r}_n + \textbf{v}_n\Delta{t}\tag{5}$$
$$\textbf{v}_{n+1} = \textbf{v}_n + \textbf{a}_n\Delta{t}\tag{6}$$
... where $\textbf{r}_n$, $\textbf{v}_n$ and $\textbf{a}_n$ are the position, velocity and acceleration vectors of a body at the start of the time-step, respectively, $\textbf{r}_{n+1}$ and $\textbf{v}_{n+1}$ are the position and velocity vectors at the end of the time-step, respectively, and $\Delta{t}$ is the size of the time-step itself. $\textbf{a}_n$ is determined by plugging all the positions of the bodies at the start of the step into $(4)$.
The global error per step for the Euler method is $\mathcal{O}(\Delta{t})$ (not so good).
Thus, to elaborate on your question, in the case of the Euler method, you would simply plug your initial positions and velocities of all the particles into equations $(5)$ and $(6)$ (having already calculated their accelerations $\textbf{a}_i$ using $(4)$) to obtain approximations for their positions and velocities at a time $\Delta{t}$, and then you would repeat the procedure for as many iterations as you like to obtain the positions and velocities of the particles at a time $2\Delta{t}$, $3\Delta{t}$, $4\Delta{t}$, ... so on and so forth.
I hope this helped.
A: When you approch to a numeric simulation of N-body there are some basic concept to take into account:


*

*Choose the more suitable integrating algoritm for your problem;

*Ask yourself if the problem has a kind of correlation length;

*Start from phenomenon observatorions;


In deatail:


*

*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).

*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.

*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.
