Given the centers of gravity of $n$ spherical bodies of unit mass, $p_1$, $p_2$, ...$p_n$, and assuming perfectly inelastic collisions, how does one find the location of the bodies after time $t$?
Note: If $t$ is long enough, I think all bodies will agglomerate to a single body of mass $n$ (but I'm not sure whether the location of this single body will be the same as the center of gravity of the initial state).
Just to clarify, everything has zero velocity at time $t = 0$.
There very likely isn't an analytic way to show where the $n$ bodies will be after time $t$, depending how large $n$ is. The best you can hope for is doing a numerical simulation.
Essentially you are evolving two equations (one to get the new position & one to get the new velocity) for each of the particles in the system. That is, for each body you solve,
$$\mathbf{x}(t+\mathrm{d}t)\simeq\mathbf{x}(t)+\mathbf{v}(t)\,\mathrm{d}t\\ \mathbf{v}(t+\mathrm{d}t)\simeq\mathbf{v}(t)+\mathbf{a}(t)\,\mathrm{d}t$$ where $\mathbf{a}(t)=\mathbf{F}(t)/m$ with $\mathbf{F}(t)$ the gravitational force, $\mathbf{x}$ and $\mathbf{v}$ the position & velocity vectors and $\mathrm{d}t$ an increment in time. More advanced techniques are called symplectic integration techniques, which are generally recommended for $n$-body simulations because they conserve energy whereas other common methods (e.g., Euler scheme) do not.
Assuming you create a class Body that stores the mass, position & velocity of each particle, your time-evolution algorithm, using a symplectic method, would be something like,
// initializes bodies, t_end, dt, etc.
while (t <= t_end) {
// Compute new accel.
for (auto body : bodies) {
accel[body] = calc_force(body, bodies) / body.mass;
}
// Compute new position
for (auto body : bodies) {
body.position += (body.velocity + 0.5 * accel[body] * dt) * dt;
}
// Compute new accel.
for (auto body : bodies) {
accel2[body] = calc_force(body, bodies) / body.mass;
}
// Compute new velocity
for (auto body : bodies) {
body.velocity += 0.5 * (accel[body] + accel2[body]) * dt;
}
t += dt;
}
(I have other details here and here; see also search: verlet or search: symplectic integrator for more)
The function calc_force can be optimized a bit since $F_{i,j}=-F_{j,i}$, but you get the idea. Now note that this is still pretty slow since you have a double sum, which is $\mathcal{O}(n^2)$ operation. There are faster techniques that can drop this to $\mathcal{O}(n\log n)$ (e.g., Barnes-Hut), but this may be a bit more than what you want, depending on the size of $n$.
If you want the particles to aggregate, then you'd need to account for this during the integration. For instance, if the distance between any two bodies is below some threshold, then you need to figure out how the merge algorithm (e.g., do collision analysis to get new velocity & do simple addition for masses) should be inserted in the evolution algorithm (probably check after each time step, but could also do every few steps). You may want to see if there are any questions on collision detection on GameDev.SE to help you on that, as I don't think there are any here on Physics.SE.
