>The only force at play is newtonian gravity, which we modify by adding a softening length, ε.

> This way we avoid numerical instabilities due to divergent forces when two particles get to close together.

So it appears he doesn't handle close encounters that form binaries. Probably reasonable if he's using heavier discrete particles as a proxy for dark matter.

I wrote my own toy gravity simulation. 2 stars that come very close together at a clock tick will experience a massive force that throws them both right out of the system. This is an artefact of artificially dividing time into clock ticks to make the calculation tractable. As I understand it, 'softening length' is a fudge to stop that happening.

You'd have to use a more sophisticated algorithm with variable time steps for closely interacting particles, I think. And then some way to deal with particles that become bound.

That sounds like a bit of a nightmare. Especially when you have a large number of bodies.

yup. molecular dynamics simulations are hard-core. and the parallelization is highly non-trivial so you need fast fabrics and complex message passing to scale. one of the raisons d'être for supercomputers

Did your simulation use the O(n) gravitational potential algorithm? You don't need to spend O(n^2) time per time step.

IIRC I was calculating forces for every pair, which is ok for small N (e.g. a solar system, rather than a galaxy).

Integrating physics simulations like that is very inefficient and inaccurate, as it does not conserve momentum.

I get it, it was a toy, but if you're really interested look up adaptive step symplectic integrators.

I was just messing doing a bit of recreational programming, with my son. It was good enough for my purposes!

Totally get it ! Just thought of suggesting it in case you wanted to dig in deeper :)