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.

Advantages

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.

Examples

I wrote a simple Ruby example ))