   Verlet/Gauss-Seidel-Relaxation stability problem
  I've been experimenting with physics like described in this commonly referenced article article many times. I started building another simulator with the algorith but to my surprise I encountered some strange, instability-like behaviour.Here I have a "rope" consisting of 30 constraints and 31 particles (1 attached, 30 freely-moving). The rope is hanging from the mouse pointer. When I gently deflect the equilibrium it starts to oscillate chaotically:I'm using a Finnish-made BASIC-dialect called CoolBasic just because of it's simplicity. Unfortunately, I couldn't find the sources of my old engines which means I have nothing to compare this with. That's why I decided to post it here if someone had any ideas of what's the cause. Maybe floating-point number rounding errors cumulate somehow, but where?PARTICLE and CONSTRAINT are C-struct-like objects which carry the data for each particle/constraint in the simulation:Type CONSTRAINTField ID%Field p1_ID%Field p2_ID%Field restLength#EndTypeType PARTICLEField ID%Field posX#Field posY#Field oldPosX#Field oldPosY#Field accX#Field accY#Field invMass#EndTypeThe two essential functions of the engine://Verlet Integration schemeFunction particle.Verlet(parID%, dt#, damp# = 0)Dim p.PARTICLEp = ConvertToType(parID)//Save current position for laterDim tmpPosX, tmpPosY As FloattmpPosX = p\posXtmpPosY = p\posY//Verlet Integration step//x(t + dt) = 2x(t) - x(t - dt) + a(t) dt^2p\posX = (1+damp)*p\posX - damp*p\oldPosX + p\accX * dt*dtp\posY = (1+damp)*p\posY - damp*p\oldPosY +//Update old positionp\oldPosX = tmpPosXp\oldPosY = tmpPosY//Reset accelerationp\accX = 0.0p\accY = 0.0EndFunction//Constraint RelaxationFunction constraint.Satisfy(cID%)Dim c.CONSTRAINTc = ConvertToType(cID)//Get the particlesDim p1.PARTICLE, p2.PARTICLEp1 = ConvertToType(c\p1_ID)p2 = ConvertToType(c\p2_ID)Dim dx#, dy# As FloatDim d#, p# As FloatDim sumInvMass# As Float//Distance componentsdx = p2\posX - p1\posXdy = p2\posY - p1\posY//Actual distanced = Distance(0.0, 0.0, dx, dy)//Distance errorp = (d - c\restLength)/(d*(p1\invMass + p2\invMass))//Error correctionp1\posX = p1\posX + p1\invMass*p*dxp1\posY = p1\posY + p1\invMass*p*dyp2\posX = p2\posX - p2\invMass*p*dxp2\posY = p2\posY - p2\invMass*p*dyEndFunctionAnd finally the Simulation part of the main loop:...For p = Each PARTICLEIf gravity Thenp\accY = -GElsep\accY = 0EndIfparticle.Verlet(p\ID, DELTA_TIME, F)Next pFor itr = 1 To RELAXATION_ITERATIONS//Glues the first particle to mousemp\posX = MouseWX()mp\posY = MouseWY()For c = Each CONSTRAINTconstraint.Satisfy(c\ID)Next cNext itr...Simulation runs with a fixed timestep and about 20-30 iterations of relaxation per frame. I have tried a wide range of damping values and time step lengths. No matter what I do, the effect remains. It actually happens with just 2 particles and 1 constraints O.o   