Stable Rigid Springs with Sequential Impulses

This is a follow-up on my previous article, "A stiffer, simpler, and more stable damped spring" where I derived a damped spring equation that displays maximum stiffness and stability, and is easier to tune than the original Hooke’s law of elasticity equation. You can read an interactive version here in the blog or download it as pdf.

In this article we will take things a step further and show how to implement an iterative constraint solver with warmstart. This method will greatly improve both spring stability and stiffness compared to the naïve approach used in the previous article. For the sake of brevity and simplicity, we will only cover the simplest case of particles connected with springs and save rigid body constraints for a later article.

Yes, I know - You can't be bothered to read the whole thing and just want to play with the simulation. This sim contains the classic wrecking ball and suspension bridge. Red lines represent naïve springs, and green lines represent sequentially solved springs. Play around with the settings to see how they affect spring behavior. Grab and pull particles with mouse.

Before diving into the algorithm we need to set up some basics. First we need a vector library. Brevity and readability takes precedence over efficiency.

class Vec2 {
    constructor(x = 0.0, y = 0.0) { this.x = x; this.y = y; }

    add(v) { return new Vec2(this.x + v.x, this.y + v.y); }
    sub(v) { return new Vec2(this.x - v.x, this.y - v.y); }
    mul(s) { return new Vec2(this.x * s, this.y * s); }
    div(s) { return new Vec2(this.x / s, this.y / s); }
    dot(v) { return this.x * v.x + this.y * v.y; }
    unit() {
        const mag = Math.sqrt(this.dot(this));
        return mag > 0 ? this.div(mag) : new Vec2();
    }
    zero() { this.x = 0.0; this.y = 0.0; ; return this; }
}

Next we need a particle class. In the update method we use the symplectic or semi-explicit Euler integration method. This is important, since the spring equation is tailored to work with this exact method. Notice that the integrator simply treats impulse as change in velocity.

class Particle {
    constructor(position = new Vec2(), velocity = new Vec2(), mass = 1.0) {
        this.position = position;
        this.velocity = velocity;
        this.impulse = new Vec2();
        this.inverseMass = mass > 0 ? 1 / mass : 0;
        this.radius = 0.0
    }

    init() {
        const mass = this.inverseMass > 0 ? 1 / this.inverseMass : 4;
        this.radius = 2 + Math.sqrt((mass * 16) / Math.PI);
    }
    
    update(dt) {
        if(this.inverseMass > 0.0) {
            this.velocity = this.velocity.add(this.impulse);
            this.position = this.position.add(this.velocity.mul(dt));
        }
        this.impulse.zero();
    }
}

Before diving into the iterative approach, let’s have a look at the non-iterative approach and identify its limitations. When implementing a simple, naïve update loop in a physics engine, resolving spring constraints might look something like this:

1. Calculate and apply the rest impulses that will satisfy the springs
2. Calculate new particle velocities and positions from the applied impulses
3. Repeat step 1

In code it might look like this:

class Spring {
    constructor(particleA = new Particle(), particleB = new Particle(), stiffness = 1.0, damping = 1.0) {
        this.particleA = particleA;
        this.particleB = particleB;
        this.cStiffness = stiffness;
        this.cDamping = damping;
        this.restLength = 0.0;
        this.reducedMass = 0.0;
        this.unitVector = new Vec2();
    }

    init() {
        const distance = this.particleB.position.sub(this.particleA.position);
        this.unitVector = distance.unit();
        this.restLength = this.unitVector.dot(distance);
        const invMasses = this.particleA.inverseMass + this.particleB.inverseMass;
        this.reducedMass = invMasses > 0.0 ? 1.0 / invMasses : 0.0;
    }

    update(invDt) {
        const distance = this.particleB.position.sub(this.particleA.position);
        const velocity = this.particleB.velocity.sub(this.particleA.velocity);
        this.unitVector = distance.unit();
        const distanceError = this.unitVector.dot(distance) - this.restLength;
        const velocityError = this.unitVector.dot(velocity);
        const distanceImpulse = this.cStiffness * distanceError * invDt;
        const velocityImpulse = this.cDamping * velocityError;
        const impulse = -(distanceImpulse + velocityImpulse) * this.reducedMass;
        const impulseVector = this.unitVector.mul(impulse);
        this.particleA.impulse = this.particleA.impulse.sub(impulseVector.mul(this.particleA.inverseMass));
        this.particleB.impulse = this.particleB.impulse.add(impulseVector.mul(this.particleB.inverseMass));
    }
}

This setup is very easy to implement, but it comes with a number of issues and limitations. In most cases you wouldn't just implement a single spring between two particles but rather connect tens or hundreds of particles into larger structures, like girders, ropes, or squishy objects. When an impulse is applied to a pair of particles in order to satisfy the spring that connects them, it inevitably pulls all other springs connected to the same particles out of their rest states.

This happens because the physics engine only handles one spring at a time and doesn't take the state of neighbouring springs into consideration when applying impulses. This neverending Tug-O-War between the springs result in a jittery, unstable simulation that may ride the Neptune Express at any moment, for no apparent reason.

Another disadvantage with this approach is a lack of rigidity in the simulation. Even with spring coefficients set to the largest stable values, ropes and chains behave like rubber bands, and stacked objects will bounce and sway. The underlying reason is a low speed of sound in the simulation. When a spring applies an impulse, it will only influence the two connected particles within that simulation loop, and the impulse propagates outwards through neighbouring springs and particles with a speed of one layer per full simulation loop.

The only efficient way to address this problem is by reducing the timestep. While this may seem naïve, it is used by game physics heavyweight Matthias Müller, who have worked for Nvidias physics department for two decades. You can see his implementations in action on his awesome open-sourced website Ten Minute Physics which is stuffed with demos. The drawback is a potential steep increase in computation. As timestep goes towards zero, CPU usage goes towards infinity. We'll be looking into a different solution here.

The iterative algorithm presented in this paper solves the issues with instability and softness in a very elegant way. The improved physics engine update loop can be expressed like this:

1. Calculate - but don't apply - the rest impulses that will satisfy the springs
2. Apply the difference between rest impulses and already applied impulses
3. Repeat step 2 either a fixed number of times or until the system reaches equilibrium
4. Calculate new particle velocities and positions from the applied impulses
5. Repeat step 1

First of all, the corrective impulse applied to a pair of particles by one spring will affect all neighbouring springs connected to the same particles during the next iteration and will influence the calculation of their corrective impulses. In other words, the impulse applied to a pair of particles will propagate outwards through the neighbouring springs, one “layer” of springs and particles further per iteration.

When a spring applies an impulse it will not only influence its two connected particles within that simulation loop, but particles as many “layers” away as there are iterations. The impulse propagates outwards through neighbouring springs and particles with a speed of one layer per simulation iteration. In other words, we have increased the speed of sound in our simulation.

This solution is only slightly more CPU intensive than the naïve approach. All the expensive calculations, primarily distance and unit vector calculations, happen only once per simulation loop in step 1. In step 2 we only do cheap calculations such as vector projection and simple algebra based on the more expensive calculations in the first step. Here is a bare-bones implementation of the described algorithm.

class SequentialSpring {
    constructor(particleA = new Particle(), particleB = new Particle(), cStiffness = 1.0, cDamping = 1.0) {
        this.particleA = particleA;
        this.particleB = particleB;
        this.cStiffness = cStiffness;
        this.cDamping = cDamping;
        this.reducedMass = 0.0;
        this.restLength = 0.0;
        this.restImpulse = 0.0;
        this.unitVector = new Vec2();
    }

    init() {
        const distance = this.particleB.position.sub(this.particleA.position);
        this.unitVector = distance.unit();
        this.restLength = this.unitVector.dot(distance);
        const invMasses = this.particleA.inverseMass + this.particleB.inverseMass;
        this.reducedMass = invMasses > 0.0 ? 1.0 / invMasses : 0.0;
    }

    applyCorrectiveImpulse(){
        const deltaImpulse = this.particleB.impulse.sub(this.particleA.impulse);
        const impulseError = this.unitVector.dot(deltaImpulse) - this.restImpulse;
        const correctiveImpulse = this.unitVector.mul(-impulseError * this.reducedMass);
        this.particleA.impulse = this.particleA.impulse.sub(correctiveImpulse.mul(this.particleA.inverseMass));
        this.particleB.impulse = this.particleB.impulse.add(correctiveImpulse.mul(this.particleB.inverseMass));
    }

    calculateData(invDt){
        const distance = this.particleB.position.sub(this.particleA.position);
        const velocity = this.particleB.velocity.sub(this.particleA.velocity);
        this.unitVector = distance.unit();
        const distanceError = this.unitVector.dot(distance) - this.restLength;
        const velocityError = this.unitVector.dot(velocity);
        const distanceImpulse = this.cStiffness * distanceError * invDt;
        const velocityImpulse = this.cDamping * velocityError;
        this.restImpulse = -(distanceImpulse + velocityImpulse) * this.reducedMass;
    }
}

This wrecking ball simulation clearly shows the huge difference between the simple spring constraint (red) and the iterative version (green). Notice how the sequential spring performs better than the naïve spring at just one iteration. This is because it applies only the difference between the rest impulse and the aleady applied impulse, while the simpler spring implementation always adds the entire rest impulse.

The idea behind warmstarting is very simple: When iterating through a system of springs, the physics engine usually begins with an initial default impulse value of zero. Unless we’re simulating a very boring setup assuming no forces or movement, this initial guess is always wrong.

1. Calculate - but don't apply - the rest impulses that will satisfy the springs
2. Apply warmstart impulse saved from previous update loop
3. Apply the difference between rest impulses and applied impulses. Save impulse for warmstart
4. Repeat step 2 either a fixed number of times or until the system reaches equilibrium
5. Calculate new particle velocities and positions from the applied impulses
6. Repeat step 1

There are different ways to approach warmstarting, but I simply use the corrective impulse vector from the previous loop projected onto the unit vector of the current loop. The smaller difference there is between the previous and current unit vectors, the larger the warmstart impulse. If the projected impulse vector is perpendicular to or points in the opposite direction of the new unit vector, the warmstart impulse is simply set to zero.

This means that in a static, non-moving setup, like a steel girder bridge that just sits there, warmstart plays a huge role because particle state, spring lengths, and unit vectors do not change. In a more dynamic, fast moving setup warmstart only plays a small role. As a rule of thumb, this warmstart method is roughly equivalent to doubling the number of iterations. Here is the updated constraint. Notice the accumulatedImpulse vector that we have added to store the sum of applied impulses.

class SequentialSpring {
    constructor(particleA = new Particle(), particleB = new Particle(), cStiffness = 1.0, cDamping = 1.0, cWarmstart = 1.0) {
        this.particleA = particleA;
        this.particleB = particleB;
        this.cStiffness = cStiffness;
        this.cDamping = cDamping;
        this.cWarmstart = cWarmstart;
        this.reducedMass = 0.0;
        this.restLength = 0.0;
        this.restImpulse = 0.0;
        this.unitVector = new Vec2();
        this.accumulatedImpulse = new Vec2();
    }

    init() {
        const distance = this.particleB.position.sub(this.particleA.position);
        this.unitVector = distance.unit();
        this.restLength = this.unitVector.dot(distance);
        const invMasses = this.particleA.inverseMass + this.particleB.inverseMass;
        this.reducedMass = invMasses > 0.0 ? 1.0 / invMasses : 0.0;
    }

    applyWarmstart(){
        const projectedImpulse = this.unitVector.dot(this.accumulatedImpulse);
        this.accumulatedImpulse.zero();
        if( projectedImpulse > 0.0 ) { return; }
        const correctiveImpulse = this.unitVector.mul(this.cWarmstart * projectedImpulse);
        this.particleA.impulse = this.particleA.impulse.sub(correctiveImpulse.mul(this.particleA.inverseMass));
        this.particleB.impulse = this.particleB.impulse.add(correctiveImpulse.mul(this.particleB.inverseMass));
    }

    applyCorrectiveImpulse(){
        const deltaImpulse = this.particleB.impulse.sub(this.particleA.impulse);
        const impulseError = this.unitVector.dot(deltaImpulse) - this.restImpulse;
        const correctiveImpulse = this.unitVector.mul(-impulseError * this.reducedMass);
        this.particleA.impulse = this.particleA.impulse.sub(correctiveImpulse.mul(this.particleA.inverseMass));
        this.particleB.impulse = this.particleB.impulse.add(correctiveImpulse.mul(this.particleB.inverseMass));
        this.accumulatedImpulse = this.accumulatedImpulse.add(correctiveImpulse);
    }

    calculateData(invDt){
        const distance = this.particleB.position.sub(this.particleA.position);
        const velocity = this.particleB.velocity.sub(this.particleA.velocity);
        this.unitVector = distance.unit();
        const distanceError = this.unitVector.dot(distance) - this.restLength;
        const velocityError = this.unitVector.dot(velocity);
        const distanceImpulse = this.cStiffness * distanceError * invDt;
        const velocityImpulse = this.cDamping * velocityError;
        this.restImpulse = -(distanceImpulse + velocityImpulse) * this.reducedMass;
    }
}

When implementing the updated constraint, we see that warmstart really has a big impact on spring rigidity. Notice how the simulation behaves almost identical at 10 iterations and no warmstart, and at just 5 iterations with max warmstart. This little trick allows us to either improve simulation quality at almost no computational cost or save us a substantial amount of calculations by halving the number of iterations without quality loss.

So far so good! However, there are still some issues. As you might have noticed, the sequential spring is not 100% stable at all configurations, especially when coeffiecients are high and iterations are low. One reason for this is that the solver always applies the entire difference between desired and actual impulse, potentially countering the previously applied impulses during earlier iterations. The result is that the last applied impulse has a much greater influence over simulation behavior than all the others.

One of several ways to handle this is by introducing a correction coefficient. It controls what fraction of the differencebetween desired and current impulse should be applied.

class SequentialSpring {
    
    constructor(particleA = new Particle(), particleB = new Particle(), cStiffness = 1.0, cDamping = 1.0, cWarmstart = 1.0, cCorrection = 1.0) {
        this.particleA = particleA;
        this.particleB = particleB;
        this.cStiffness = cStiffness;
        this.cDamping = cDamping;
        this.cWarmstart = cWarmstart;
        this.cCorrection = cCorrection;
        this.reducedMass = 0.0;
        this.restLength = 0.0;
        this.restImpulse = 0.0;
        this.unitVector = new Vec2();
        this.accumulatedImpulse = new Vec2();
    }

    init() {
        const distance = this.particleB.position.sub(this.particleA.position);
        this.unitVector = distance.unit();
        this.restLength = this.unitVector.dot(distance);
        const invMasses = this.particleA.inverseMass + this.particleB.inverseMass;
        this.reducedMass = invMasses > 0.0 ? 1.0 / invMasses : 0.0;
    }

    applyWarmstart(){
        const projectedImpulse = this.unitVector.dot(this.accumulatedImpulse);
        this.accumulatedImpulse.zero();
        if( projectedImpulse > 0.0 ) { return; }
        const correctiveImpulse = this.unitVector.mul(this.cWarmstart * projectedImpulse);
        this.particleA.impulse = this.particleA.impulse.sub(correctiveImpulse.mul(this.particleA.inverseMass));
        this.particleB.impulse = this.particleB.impulse.add(correctiveImpulse.mul(this.particleB.inverseMass));
    }

    applyCorrectiveImpulse(){
        const deltaImpulse = this.particleB.impulse.sub(this.particleA.impulse);
        const impulseError = this.unitVector.dot(deltaImpulse) - this.restImpulse;
        const correctiveImpulse = this.unitVector.mul(-impulseError * this.reducedMass * this.cCorrection);
        this.particleA.impulse = this.particleA.impulse.sub(correctiveImpulse.mul(this.particleA.inverseMass));
        this.particleB.impulse = this.particleB.impulse.add(correctiveImpulse.mul(this.particleB.inverseMass));
        this.accumulatedImpulse = this.accumulatedImpulse.add(correctiveImpulse);
    }

    calculateData(invDt){
        const distance = this.particleB.position.sub(this.particleA.position);
        const velocity = this.particleB.velocity.sub(this.particleA.velocity);
        this.unitVector = distance.unit();
        const distanceError = this.unitVector.dot(distance) - this.restLength;
        const velocityError = this.unitVector.dot(velocity);
        const distanceImpulse = this.cStiffness * distanceError * invDt;
        const velocityImpulse = this.cDamping * velocityError;
        this.restImpulse = -(distanceImpulse + velocityImpulse) * this.reducedMass;
    }
}

And finally, the complete simulation in action! Play around with the settings to see how they affect spring behavior. Grab and pull particles with mouse.

The advantage of adding a correction factor may not be as apparent as the other improvements, it's more subtle. But by setting it to value lower than one greatly improves stability at the cost of rigidity.

As promised we have derived an iterative solver that gratly improves spring stiffness and stability. It serves as an excellent foundation for a soft body physics engine, which we will explore in detail in later articles.