• Create Account

Like
12Likes
Dislike

# Towards a Simpler, Stiffer, and more Stable Spring

By Michael Schmidt Nissen | Published Dec 19 2014 01:53 AM in Math and Physics
Peer Reviewed by (Hodgman, apatriarca, jbadams)

spring damping mass-spring-damper harmonic oscillator reduced mass

# Towards a Simpler, Stiffer, and more Stable Spring

Written by Michael Schmidt Nissen
Version 0.5.2, December 18th 2014

## The trouble with springs

The damped spring based on Hooke's law of elasticity is often regarded as a black sheep in the field of game physics, and some papers and websites even warn the reader against implementing them. There are a couple of 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 spring stiffness and damping coefficients in order to make a simulation behave just right. For large systems of interconnected particles and springs the task can be daunting. If a variable is changed, for instance a spring coefficient, the mass of a particle, or the simulation time-step, you might have to re-tune the entire system from scratch.

A more profound reason for avoiding springs is that they are potentially unstable. When increasing the spring coefficients above some unknown threshold, the simulation will start to oscillate chaotically and explode in all directions. Everyone who ever worked with springs has experienced this frustrating phenomenon. The problem is that there’s no good method to determine in advance if a given spring system will show this behavour or not. Thus, the programmers or game designers ability to balance a spring system often rely on experience and gut feeling, which makes it look more like black magic than actual physics.

In this article I’ll introduce a modified spring equation that is easy to tune, displays maximum stiffness and damping, and is guaranteed to keep the spring system stable under all conditions. The math behind the equation is surprisingly simple, which makes me wonder why it doesn’t appear to have been mentioned or implemented earlier.

## Making a better spring

The improved spring equation stems from trying to answer a simple question: Is it possible to make a spring that reaches its rest state in just one simulation loop? It turns out that the answer is yes! This wouldn’t be possible in the real world, since it would require infinitely high stiffness and damping coefficients or zero mass. But in the not-quite real world of discrete, time-stepping physics simulation, this is achievable.

To explain how this is done, let's look at a one dimensional harmonic oscillator with a damped spring attached to a fixed point in one end and to a moving body in the other. Assume that the body is displaced distance x from its rest position. Now let us try to find the exact amount of force required to move the mass this distance in one simulation loop.

Distance can be expressed as velocity multiplied by time

$\large x = -{v} {\Delta t}$

The minus sign simply means that we're moving in the opposite direction of the displacement. Velocity can be expressed as acceleration multiplied by time

$\large x = -{a} {\Delta t^{2}}$

Newton's 2nd law of motion states that acceleration equals force over mass

$\large x = -\dfrac{f}{m} {\Delta t^{2}}$

Now we can isolate force

$\large f= -\dfrac{m}{\Delta t^{2}} x$

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

$\large f= -k x$

It becomes immediately clear that

$\large k = \dfrac{m}{\Delta t^{2}}$

Which is the exact spring 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.

Velocity is equal to acceleration multiplied by time step

$\large v = -{a} {\Delta t}$

Which is equal to force over mass

$\large v = -\dfrac{f}{m} {\Delta t}$

When isolating force we get

$\large f = -\dfrac{m}{\Delta t} v$

And since the damping force is defined

$\large f = -d v$

We immediately see by inspection that the damping coefficient is

$\large d = \dfrac{m}{\Delta t}$

This is the exact damping coefficient value needed to stop the particle in one simulation loop. Now we can write the complete damped spring equation.

$\large f= -\dfrac{m}{\Delta t^{2}} x -\dfrac{m}{\Delta t} v$

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/Δt2 and d with m/Δt. When implementing the equation we see that it really does work! The spring reaches equilibrium in one loop, completely independent of particle position, velocity, mass, and simulation time-step. This is the stiffest possible spring, and it displays behavior more similar to a rigid constraint than a soft, bouncy spring. The equation also has another interesting property. It simply cannot blow up no matter what values you feed it. Practically speaking, we can regard the spring as being unconditionally stable.

## Re-introducing the spring coefficients

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

$\large f= -\dfrac{m}{\Delta t^{2}} C_{k} x -\dfrac{m}{\Delta t} C_{d} v \qquad [ 0 \leq C_{k}, C_{d} \leq 1 ]$

As we can see, the new coefficients are simply the fraction of completely rigid behavior we would like the spring to display. Soft, bouncy springs would usually have values in the range of 0.001 – 0.00001. In the other end of the scale, values of just 0.1 – 0.01 is enough to display rigid behavior. Setting both values to 1.0 would of course still satisfy the constraint in one loop.

Please notice that spring behavior is determined exclusively by these two coefficients. Particle mass or simulation time-step has no influence on how the spring behaves, and changing them wouldn't make it necessary to tune the system!

Interestingly, the spring will get less rigid and less stable if we increase Ck or Cd above 1. If we keep increasing the coefficient values, the system will start oscillating chaotically, and at some point it will explode. In other words, we have determined the exact upper limit for the two spring coefficients, which we define

$\large k_{max} = \dfrac{m}{\Delta t^{2}}$
$\large d_{max} = \dfrac{m}{\Delta t}$

This allows us to simplify the spring equation

$\large f= -k_{max} C_{k} x -d_{max} C_{d} v$

There is an important conclusion to be drawn from this. It is a misunderstanding to think that spring stiffness and damping can be increased towards infinity by simply increasing k and d. A system of very stiff springs doesn’t necessarily blow up because the integration algorithm can't handle it, but because the coefficients might have been set to a value that simply does not make sense.

## Two connected particles with different mass

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 was stationary, which allows us to use the equation we've already derived. The reduced mass for two particles with mass m1 and m2 is defined as

$\large m_{red} = \dfrac{m_1 m_2}{(m_1 + m_2)}$

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

$\large m_{red}= \dfrac{1}{ ( \dfrac{1}{m_1} + \dfrac{1}{m_2} ) }$

For two connected particles, the maximum coefficient values are

$\large k_{max} = \dfrac{m_{red}}{\Delta t^{2}}$
$\large d_{max} = \dfrac{m_{red}}{\Delta t}$

When replacing mass with reduced mass we get

$\large f= -\dfrac{m_{red}}{\Delta t^{2}} C_{k} x -\dfrac{m_{red}}{\Delta t} C_{d} v$

This spring equation will bring the two connected particles to equilibrium distance in one loop and make them stand still relative to each other. However, since energy and momentum are 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 rotate.

## Angular springs

The spring equation can also be used for angular springs, also commonly named rotational or torsional springs. Rather than keeping two particles at a fixed distance by applying opposite equal forces, the angular spring will try to keep two rotating bodies at a fixed angle by applying opposite equal torques. The equation introduces the concept reduced moment of inertia, which is calculated in a similar way as reduced mass

$\large I_{red} = \dfrac{I_1 I_2}{(I_1 + I_2)}$

Or alternatively

$\large I_{red}= \dfrac{1}{ ( \dfrac{1}{I_1} + \dfrac{1}{I_2} ) }$

The maximum coefficient values for angular springs are

$\large k_{max} = \dfrac{I_{red}}{\Delta t^{2}}$
$\large d_{max} = \dfrac{I_{red}}{\Delta t}$

When replacing the variables of linear motion with those of angular motion we get

$\large \tau = -\dfrac{I_{red}}{\Delta t^{2}} C_{k} \theta -\dfrac{I_{red}}{\Delta t} C_{d} \omega$

This spring equation will keep two rotating bodies at any given rest angle. If both coefficients are set to 1, the constraint will be solved in one loop. The spring allows for angular displacements larger than one full turn, making it possible to “wind up” bodies like the coil spring of a mechanical clock.

## Impulse-based springs

Today a lot of physics engines and simulation methods are based on impulses - direct changes in velocities - rather than forces and acceleration. The linear and angular spring equations described above works equally well if we redesign them to work with impulses. Here are the equations without further ado

$\large J_{linear}= -\dfrac{m_{red}}{\Delta t} C_{k} x - C_{d} v$
$\large J_{angular} = -\dfrac{I_{red}}{\Delta t} C_{k} \theta - C_{d} \omega$

The only notable difference is that when dealing with continuous forces that change gradually over time – like in springs – the impulse based equations are limited to 1st order of magnitude accuracy when it comes to conserving energy and momentum. Force based equations are, on the other hand, only limited by the integration algorithm, and can return results that are much more accurate.

## Limitations and pitfalls

When connecting a multitude of springs and particles into larger bodies, we run into the same trouble as any other type of distance and velocity constraint. Rather than cooperating, the constraints tend 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 the author 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 deformable body, setting the coefficients to the maximum value of 1.0 will lead to stability problems. Although the spring equation is stable when used with a single or two particles, this is sadly not the case for higher number of particles. The author of this article has after some lengthy tampering worked out that a safe upper limit for both the stiffness and damping coefficient is

$\large c_{max} \approx \dfrac{1}{n + 1}$

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, Ck 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.

## Reference

http://en.wikipedia.org/wiki/Reduced_mass

Michael Schmidt Nissen lives in Denmark and currently makes a living as a traditional blacksmith. As opposed to most other smiths, Michael extracts iron and steel from local ores using the same techniques as his iron age and viking era ancestors. Programming - especially physics simulations - has been a beloved hobby ever since Michael learned to code in the golden era of the Commodore 64 and Amiga 500.

Your article provides a useful approach to simulating a spring + damper system. However, you characterize a system of a spring + damping as just a spring, which is quite misleading. The principle of adding a damper to an oscillatory system has been around for a long time. Although your article seems to be directed to particle simulation, a spring + damper (as separate components) is the more common approach to (for example) vehicle simulation.

The math behind the equation is surprisingly simple, which makes me wonder why it doesn’t appear to have been mentioned or implemented earlier.

Although I haven't seen your particular presentation before, the math behind your equations is just a (slightly incorrect**) derivation of the math related to a classical mass-spring-damper system (linked above,) with the goal of determining spring and damping constants to achieve certain conditions.

Hopefully the following is correct, as it's mostly off the top of my head, trying to remember course-work from several decades ago-

If the required conditions are:
- initial displacement from spring rest position = X
- initial velocity = V
- after time T, displacement = -X (returns to rest position)
- after time T, velocity = 0 (at rest)

Given the classical displacement equation:

x(t) = x0 + v0 * t + 0.5 * a0 * t^2

then the required x(t) = 0 (back to rest position)

0 = X + VT + 0.5*a*T^2

Solve for a:
a = -2X/T^2 -2V/T

Classical spring-damper equation:
F = -kx - cv = ma

Substitute the a above which reflects required conditions:

a = -kx/m - cv/m = -2X/T^2 - 2V/T

or

k = 2m/T^2
c = 2m/T

** You assumed x = -a*t^2, rather than -1/2*a*t^2, I believe that's the factor of 2 difference in your coefficients.

You assumed x = -a*t^2, rather than -1/2*a*t^2, I believe that's the factor of 2 difference in your coefficients.

The equation that you cited, the one with the 1/2 term, only holds when the timestep is differential instead of discrete. When working with timesteps that aren't infinitely small, that equation holds. Trying to make the classical equations of motion work in a discrete world, like the motion equation you cited, is one of the reasons that simulations explode while real life doesn't

@OP: Very nice article, btw! Lots of thought went into this, I can tell
The equation that you cited, the one with the 1/2 term, only holds when the timestep is differential instead of discrete. When working with timesteps that aren't infinitely small, that equation holds.

Do you have a reference for that concept? Particularly the concept that the classical displacement equation doesn't apply for timesteps larger than infinitely small?

The equation that you cited, the one with the 1/2 term, only holds when the timestep is differential instead of discrete. When working with timesteps that aren't infinitely small, that equation holds.

Do you have a reference for that concept? Particularly the concept that the classical displacement equation doesn't apply for timesteps larger than infinitely small?

I think the biggest thing here is that the symbols are a little confusing. The OP was saying that

delta_x = a * delta_t^2.

While it's true that

x = a/2*t^2

, it's NOT true that

delta_x =  a/2 * delta_t^2

It should instead be equal to

delta_x = a/2 * (t + delta_t)^2 - a/2*t^2 = a * t * delta_t + a/2 * delta_t^2

The issue with that equation is that it only supports a constant acceleration (because of the a*t term), and because you have to keep track of the t value. A Euler integrator uses

delta_x = a * delta_t^2,

because it too converges to the correct a/2*t^2 formula for a small enough delta_t, but you don't have to be confined to constant acceleration or have to keep track of t.

Buckeye, CulDeVU - thanks for your constructive feedback, and thanks for taking time reading my article.

You've pointed my attention to some unclear parts in my article, which I will fix asap. I've failed to explain clearly that this method does not in the strictest sense apply to real-world physics, but is tailored to exploit the limitations in discrete physics as seen in impulse based methods and the 1st order symplectic Euler differential algorithm, which are to my knowledge the two most common methods in use. In the real world indeed x(t+dt) = v(t) * dt + 0.5 * a(t) * dt^2, as seen expressed in 2nd order algorithms like the leap-frog and verlet integration algorithms. However, this is *not* true when working with impulse based physics or the symplectic Euler algorithm. In the strictest sense, the algorithms return wrong results, since per definition they assume constant force during the entire timestep. However, when using very small timesteps they converge towards the right result. It's a compromise that works satisfactory in the realm of game physics which doesn't need to be too accurate to produce rewarding, believable results.

Cheers Mike

Also, I'd like to add that the equations in the article have been developed across a timespan of several years, and that they have been implemented into a physics engine and *extensively* twisted and tested for some time. There may be some formal errors in the article, mostly probably stemming from the fact that english is not my mother tongue or that I'm self-taught and do not have much education on this topic (I'm actually a traditional blacksmith). But the equations work for sure, wether I have succeeded in explaining them sufficiently or not.

Oh, and there are some fun-to-play-with code samles coming up during the christmas holidays!

Without doubt, your English is better than my Danish, and the math is well expressed. The intent of my comments, perhaps, could've been better expressed. Not to diminish your math work, or the work you put into the article: it appears your approach is intended to determine stable tuning parameters for a spring-damper system for a restricted case. That could, perhaps, be emphasized more.

FYI: your approach is closely related to the general method for tuning a system for critical damping. Note, however, that critical damping is not always the desired response for more general systems. Critical damping is the slowest system response that is not under-damped. Your approach to determining stable tuning parameters may not be in more general use because, in many cases, a faster (though potentially less stable) response is desired. E.g., the goal of industrial control system tuning is often to achieve quarter-wave damping, which is a compromise between stability, and approaching the desired state in a reasonable time. The "need for speed" is related to what you've mentioned: getting to a stable state before external forces cause additional changes.

Thanks for your effort! Overall, a very nice article.

And Merry Christmas to you, too!

When increasing the spring coefficients above some unknown threshold, the simulation will start to oscillate chaotically and explode in all directions

Your motivation for all this is flawed. This "unkown threshold" is actually a known threshold, the term that applies here is the critical timestep. You can compute the critical timestep from known quanitities of mass and stiffness. Further reading is found here (see the section for a general mass-spring system):

Now in a game environment there are likely limitations on timestep. However it is very easy to do a "back of the napkin" calculation to see if a simple mass-spring system will work in your case.

@Buckeye

Thank you for your kind words. As already mentioned, I will ephasize in the article that the solution is restricted to 1st order discrete physics simulations and have no real-world application. As for the "need for speed" damping, my equation brings the object to rest exactly at the equilibriun point in one loop, and thats the best possible result you can imagine. However, this is not the same as critical damping. I use d = m/dt, where real-world critical damping is d = 2*sqr(m*k). Again, this is only possible within the limited field of 1st order discrete physics, not in the real world.

Cheers & a Merry christmas to you too :-)

@brudedjones

The equations seems to apply to real-world cases or cases using higher order integration, which makes it incompatible with my method that deliberatle exploits the limitations of 1st order integration. Also, defining global timestep from a single spring's stiffnes coefficient doesn't seem too practical. But I'll be looking into the math.

I'd like to point out that the author of the website you link to explicitly admits that the calculated critical timestep may lead to problems, and that in some cases lower timesteps are necessary in order to get good results. But he cannot tell how much lower. In other words, he admits that the critical timestep is just a loose guesstimate:

Let us note at this place that not only  assuring numerical stability of motion integration is a constraint. In systems where particles move at relatively high velocities, position change during one timestep can lead to non-elastic irreversible effects such as damage. The  needed for reasonable result can be lower . We have no rigorously derived rules for such cases.