I have written an 2d N-body simulation in Python which allows collisions between the bodies.
A body is modeled as a circle whose area is proportional to its mass. Each time the sim advances by one time-step, I identify which bodies have overlapping radii and I replace them by one particle (respecting conservation of momentum and treating the collision as inelastic).
This worked well for N<100 but for higher N I encountered (not systematically, but depending on the initial conditions of the sim) a weird error: the number of particles N actually increased (and quickly exceeded the maximum integer value Python can handle).
I believe what is occuring is, after the new positions of the bodies have been calculated, it happens that (at least) 4 bodies mutually overlap. Each of the 6 pairs creates a new body and but we also delete the 4 bodies,so N increased by os we have 6-4=2. But the 6 new bodies are spawned in roughly the same positions as the original 4 bodies so now we have even more overlap (N can now increase by potentially (6 choose 2)-6=9 bodies). Hence why N can grow very fast.
How can I avoid this error?
The only thing I can think of is to reduce the radii of the bodies to reduce the change of 4 particles mutually overlapping. But this does not definitively fix the problem.