Developing a rigid, stable spring

A few years ago I wrote a paper on how to implement an improved version of Hooke's law springs. You can read it here. Combined with an iterative solver it is an excellent building block for rigid and stable soft body simulations. In this post I will go through the very simple maths behind it and provide examples in JavaScript.

In this article I will derive a very simple damped spring constraint equation that doesn't display any of the problems most game-developers associate with springs. From a game programming and physics engine perspective, the spring equation combines the highest possible stiffness and damping with great numerical stability. From a game design and level editing perspective, the equation is easy to tune and does not need to be re-tuned if parameters like particle mass or time step are changed. This makes the improved spring equation a great tool for game designers and engine programmers alike.

Okay, you're in a hurry - Aren't we all at times... You just want the conclusion and a fun demo to fiddle around with, but you can't be bothered to follow the links. Here's what you need. The spring equation is tailored to work with the very popular symplectic Euler or semi-explicit Euler method. You MUST use this, or the constraint won't work properly.

Spring reduced mass, used for scaling forces and impulses between two particles with different mass:

Force based spring equation, used with the symplectic Euler integration algorithm:

Impulse based spring equation, used with the impulse based integration algorithm:

Spring stiffness and damping coefficient $C_{s}$ and $C_{d}$ are in the range $[0, 1]$. The spring is unaffected by particle mass or time step and will reach its rest position in as little as one loop, depending on the coefficients. Here is an interactive demo showing the equations in action:

In the game development community, the damped spring based on Hooke's law of elasticity is often regarded as the black sheep of game physics, and a number of notable books, papers and websites warn against implementing them. There are at least two good reasons for this.

The most straight-forward reason is that springs are notoriously difficult to tune. It can be very time consuming to adjust the stiffness and damping coefficients in order to make a simulation behave just right. For large spring systems the task can be daunting. If a single variable is changed, for instance a spring coefficient, the mass of a particle, or the physics engine time step, you might have to re-tune everything from scratch.

A deeper reason for avoiding springs is that they are prone to instability. Every once in a while, a system of connected springs and particles will blow up in your face without warning, and for no apparent reason. This happens because either simulation time step or spring coefficients have been set too high, or particle mass has been set too low, placing the simulation outside it's range of stability.

The problem is that there’s no simple or death-proof way to determine what combination of values adds up to a stable system. Thus, the game designer's ability to balance a spring system often relies on experience and “gut feeling”, which makes it look more like black magic than actual physics. This lack of transparency is why the game development community looks upon springs with a lot of suspicion.

In the following pages we will do a complete breakdown of the math behind damped springs and use it to reformulate the old equation in a way that makes sense for game developers. We will only be focusing on how springs work in the discrete time-stepping world of game physics, not how they work in real-life, which is a slightly different thing. First, let's take a look at the equation at hand:

The equation is divided up into two parts. The first part calculates a spring force that is linearly proportional to displacement x, and the second part calculates a damping force that is linearly proportional to velocity v. The value of the stiffness and damping coefficients k and c decide the magnitude of the forces. Game designers usually want a spring to be as rigid as possible. That is, they want it to behave like a steel rod, and not like a soggy old rubber band. But what k and c values do we plug in to achieve this behaviour? 0.001? 1.337E42? For someone who doesn't understand the underlying math it is basically trial-and-error: Plug in a random value. If the spring behaves too sloppy, increase the value, and if the simulation becomes jittery or explodes, decrease it.

In the following section we will derive a spring equation that eliminates this challenge and allows the game designer to implement a spring that will always stay inside the safe region of stability and will reach its rest state in as little as one loop. To explain how this is done we first need to take a look at our integrator – the part of the program that computes new particle velocities and positions based on the forces applied to them. The most commonly used integrator is the so-called symplectic Euler or semi-explicit Euler method:

This integrator is popular because it is very easy to implement, and it displays great energy conservation and stability compared to similar algorithms. Let's move on and look at a damped spring attached to a fixed point in one end and a particle in the other. Assume that the particle is displaced distance x from its rest position. Now let us try to find the force required to move it the opposite distance in one loop with the symplectic Euler algorithm. From Newton's 2nd law of motion we know that \( a = F / m \), which give us the following relationship between distance and force:

The minus sign simply means that we are going in the opposite direction of the displacement. Now we can isolate force:

And since the spring force in Hooke's law of elasticity is defined as:

It becomes immediately clear that:

Which is the exact stiffness coefficient value required to move the particle back to its rest position in one simulation loop. However, since we are not doing anything to stop the particle, it will keep oscillating back and forth through the rest position. We need to add damping, which is done with the second part of the spring equation. Let's assume the particle has a velocity v. By calculating the amount of force needed to exactly counter this velocity, we can make the particle stop at the rest distance. In the Symplectic Euler integration algorithm we have the following relationship between velocity and force:

When isolating force we get:

And since the damping force is defined:

We immediately see by inspection that the damping coefficient is:

This is the exact damping coefficient value needed to make the spring stop the particle in one simulation loop. Now we can join the two parts and write out the complete damped spring equation in its most basic form:

This spring equation has some very interesting properties. The first thing we notice is the lack of coefficients. We've simply replaced \(k\) with \(m / \Delta t ^{2}\) and \(c\) with \(m / \Delta t \). When implementing the equation we see that it really does work! The spring reaches rest length and stops dead without any oscillation in one loop. This is the stiffest possible spring we can imagine, and it displays behaviour more similar to a rigid constraint than a soft, bouncy spring. We have found the steel rod spring we were looking for in the beginning of the chapter.

This is a big deal! In game-physics papers it is considered a universal truth that to get infinite spring stiffness you need an infinitely high k value, which is not possible. This is obviously not true; we get infinite stiffness from finite k and c values. The equation also has another interesting property. It simply cannot blow up no matter what values we throw at it. We can regard the spring as being unconditionally stable. We can refactor the equation a bit so it's easier to use in a programming context. Now time step, state, and mass are separated, which makes it easier to write re-useable code. Since physics engines normally run with a fixed time step, we can also pre-compute the \( 1 / \Delta t \) and \(1 / \Delta t ^{2} \) values and save a few precious CPU-cycles.

Now we have a really powerful spring equation. It is easy to implement, very stable, and it reaches its rest-state in just one loop. But in our quest towards a better spring it has lost its springiness. We need to get the softness and bounciness back again, and for this purpose we will re-introduce stiffness and damping coefficients. To avoid confusion, the new coefficients are named Cs and Cd. While the original coefficients could represent any positive numerical value, these both lie in the interval between zero and one.

As we can see, the new coefficients are simply the fraction of completely rigid behaviour we would like the spring to display. Soft, bouncy springs would usually have values in the range of 0.00001 – 0.001. In the other end of the scale, values of just 0.01 – 0.1 is enough to display rigid behaviour. Setting both values to 1.0 would of course still satisfy the spring in one loop. Please notice that spring behaviour is determined exclusively by these two coefficients. Particle mass or simulation time step has literally no influence on how the spring behaves, and changing them – even during runtime – wouldn't make it necessary to re-tune the system. This makes the spring a much more powerful and developer-friendly tool that can be used safely by designers without any worries about possible future changes in game physics.

Interestingly, the spring will get less rigid and less stable if we increase the Cs or Cd values above one. This happens because the particles overshoot their target rest values, leading to undefined behaviour. If we keep increasing either or both of the coefficients, the system will start to jitter, and at some point it will ride the Neptune Express. In other words, we have – almost by accident – determined the exact upper limit for the two spring coefficients, which we can define:

This also allows us to simplify the spring a bit, making it almost as simple as the equation we started out with:

Now it also becomes clear why so many people dislike the Hooke's law spring equation. Not because it is inherently bad, but because it is so easy to accidentally plug in a coefficient value that doesn't make any sense. The equation derived in this paper always scales spring force to fit time step and mass, but the original equation does not.

This makes it very easy for the less physics-minded game designer to inadvertently plug in a coefficient value above the upper limits we have just derived. In the original equation it isn't even apparent that there is such a thing as an upper limit. I would like to stress the point that this lack of transparency in the original Hooke's law spring is the main reason why so many people have frustrating experiences with it.

It is only slightly more complicated to constrain two free-moving particles with the improved spring. To do this, we need to introduce the concept of reduced mass. This is a quantity that can be used to compute interactions between two bodies as-if one body were stationary, which allows us to reuse the equation we've already derived. The reduced mass for two particles with mass ma and mb is defined:

Since the inverse mass quantity is often pre-computed for other purposes, it can also be useful to define reduced mass like this:

Unsurprisingly, the maximum coefficient values for two connected particles are:

Which gives us the following equation to work with:

This spring equation will move two connected particles of any given mass to rest distance and make them stand still relative to each other in one loop. However, since angular momentum is conserved, the particles may rotate around each other, which will make the bodies come to rest at a larger distance, depending on how fast they spin.

Today, most physics engines are based on impulses, or direct changes in velocities, rather than forces and acceleration. The spring equation described above works just as well if we redesign it to work at impulse-level. Since impulse is defined \(J = f \Delta t\), the symplectic Euler integration algorithm needs a little update:

The maximum coefficient values for the impulse based spring are:

This gives us the following impulse based spring to work with. When used with the integration algorithm shown above, this equation returns the numerically exact same result as the force based spring we derived earlier in the paper.

And that pretty much sums it up. Now I think it's about time we got our hands dirty and turned all the math into a working piece of code.

First we need a vector class. Here, simplicity and brevity takes precedence over optimization.

class Vec2 {
    constructor(x = 0.0, y = 0.0) {
        this.x = x;
        this.y = y;
    }
    static add(a, b) {
        return new Vec2(a.x + b.x, a.y + b.y);
    }
    static sub(a, b) {
        return new Vec2(a.x - b.x, a.y - b.y);
    }
    static mul(a, scalar) {
        return new Vec2(a.x * scalar, a.y * scalar);
    }
    static div(a, scalar) {
        return new Vec2(a.x / scalar, a.y / scalar);
    }
    static dot(a, b) {
        return a.x * b.x + a.y * b.y;
    }
    static normalize(a) {
        const length = Math.sqrt(a.x * a.x + a.y * a.y);
        return length > 0 ? new Vec2(a.x / length, a.y / length) : new Vec2(0, 0);
    }
}

Next we need the particle class. This part of the code implements the symplectic Euler integration algorithm, and it is important to understand that the spring equation is tailored to be used with this update method. If you swap this out with another algorithm, it will break the spring equation.

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;
    }
    update(dt) {
        this.velocity = Vec2.add(this.velocity, Vec2.mul(this.impulse, this.inverseMass));
        this.position = Vec2.add(this.position, Vec2.mul(this.velocity, dt));
        this.impulse = new Vec2(0, 0);
    }
}

Finally we have the spring class. The update function computes and applies the spring impulse to its two endpoint particles. This is the core of our physics engine, and this is where the spring equation dwells.

class Spring {
    constructor(particleA = new Particle(), particleB = new Particle()) {
        this.particleA = particleA;
        this.particleB = particleB;
        this.unitVector = new Vec2();
        this.restDistance = 0.0;
        this.reducedMass = 0.0; 
        this.stiffnessCoefficient = 1.0;
        this.dampingCoefficient = 1.0;
    }
    init() {
        // Calculate the initial rest distance and reduced mass.
        const distance = Vec2.sub(this.particleB.position, this.particleA.position);
        this.restDistance = Math.sqrt(Vec2.dot(distance, distance));
        this.unitVector = Vec2.normalize(distance);
        const sumInverseMasses = this.particleA.inverseMass + this.particleB.inverseMass;
        this.reducedMass = sumInverseMasses > 0 ? 1 / sumInverseMasses : 0;
    }
    update(invDt) {
        // First we calculate the distance and velocity vectors.
        const distance = Vec2.sub(this.particleB.position, this.particleA.position);
        const velocity = Vec2.sub(this.particleB.velocity, this.particleA.velocity);

        // We normalize the distance vector to get the unit vector.
        this.unitVector = Vec2.normalize(distance);

        // Now we calculate the distance and velocity errors.
        const distanceError = Vec2.dot(this.unitVector, distance) - this.restDistance;
        const velocityError = Vec2.dot(this.unitVector, velocity);

        // Now we use the spring equation to calculate the impulse.
        const distanceImpulse = this.stiffnessCoefficient * distanceError * invDt;
        const velocityImpulse = this.dampingCoefficient * velocityError;
        const impulse = -( distanceImpulse + velocityImpulse ) * this.reducedMass;
        const impulseVector = Vec2.mul(this.unitVector, impulse);

        // Finally, we apply opposite equal impulses to the endpoint particles.
        this.particleA.impulse = Vec2.sub(this.particleA.impulse, impulseVector);
        this.particleB.impulse = Vec2.add(this.particleB.impulse, impulseVector);
    }
}

And here we can see the code in action. Pick up masses with left mouse. Use the sliders to experiment with settings in order to get a feel for the relation between stiffness, damping, timestep, gravity, and mass.

When connecting multiple springs and particles into larger bodies, we run into the same trouble as with any other type of constraint. Rather than cooperating, they tend to compete against each other, and this spells trouble. When a spring moves two particles to satisfy distance and velocity, it usually means dissatisfying one or several other springs. It is outside the scope of this article to dig deeply into this problem, but I would like to provide a bit of advice on how to prevent the worst disasters.

When two or more springs are connected to the same particle, which is the case in any kind of rope, mesh, or squishy body, setting the coefficients to the maximum value of 1.0 will lead to stability problems. Although the spring equation is stable when particles are connected to just one spring, this is sadly not the case for higher number of springs. After some lengthy tampering I have worked out that a safe upper limit for both the stiffness and damping coefficient is:

Where n denotes the highest number of springs attached to any of the two particles connected by the spring. So for example, in a rope where any particle is connected by at most two springs, Cs and Cd can both safely be set to 0.33, and in a square mesh, where any particle is at most connected by four springs, they can be set to 0.2.

As promised, I have managed to derive a very stable, rigid, and user-friendly spring equation over the course of a few pages. The math behind the equation is surprisingly simple, which makes me wonder why it hasn't been described or implemented decades ago. In the next article, we will look further into the mentioned limitations when implementing large systems of interconnected springs and improve it even further using an iterative approach with warm starting.

Michael Schmidt Nissen lives in Denmark and currently makes a living as a back-end developer. In his earlier life he worked as a school teacher and historical blacksmith. He belongs to a very small group of people who knows how to extract iron and steel from iron ores using the same techniques as his Viking ancestors. Programming, especially small games and physics simulations, has been a beloved hobby ever since Michael learned to code in the golden era of the Commodore 64 and Amiga 500.

If you enjoyed reading this article, please consider giving it some constructive feedback or sharing it with your friends or colleagues. If you use the spring equation in a game or demo, please do give the author credit for it.

• Wikipedia: Hooke's Law
• Wikipedia: Newton's laws of motion
• Wikipedia: Reduced mass
• Wikipedia: Implicit Euler
• Gaffer On Games: Spring physics