The Verlet Integration algorithm is simple, stable and … not ideal yet.
Simple & Nice
The code implementation of the algorithm is pretty simple and straightforward
prev = obj.position df = obj.force_sum * DT * DT / obj.mass obj.position += obj.position - obj.position_prev + df obj.position_prev = prev
Note that this code assumes the time step to be constant, although the variable time step could be used also with some code modification.
Beside the simplicity of implementation the algorithm has another advantage. The constrains on the particle positions are simpler to solve and it works very stable. Usually the distance constrains are modeling with the springs of infinite stiffness. We may code it that way
r = o2.position - o1.position delta = r.magnitude - @original_distance delta *= STIFFNESS k = o1_mass / o2_mass d1 = delta * k / (1 + k) d2 = delta - d1 o1.position += r.normalize * d2 o2.position -= r.normalize * d1
This is a so called relaxation function. It is called for each particle after new positions calculation is finished. Practically it is enough to call it once per step but for complex constrain configurations may be called multiple times to met error threshold.
Also we use masses of the particles to preserve the summary impulse.
I wrote a simple Ruby example ))