• # Towards a Simpler, Stiffer, and more Stable Spring

Math and Physics

• Posted By h4tt3n

# 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$ In this particular case there is no difference between the force-based and impulse-based spring. They should return the exact same results.

## 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

Report Article

## User Feedback

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.

##### Share on other sites

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 :P

@OP: Very nice article, btw! Lots of thought went into this, I can tell :)

##### Share on other sites
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?

##### Share on other sites

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.

##### Share on other sites

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

##### Share on other sites

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!

##### Share on other sites

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!

##### Share on other sites

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.

##### Share on other sites

@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 :-)

##### Share on other sites

@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.

Cheers Mike

##### Share on other sites

Well, that's an awful lot of downvotes someone just gave me... But since he just wants to discredit me anonymously rather than give constructive criticism as the site managers encourage readers to do, I guess I'll never know what he doesn't like about the article :-)

##### Share on other sites

Nice article - I've derived similar stiffness limits myself, but I didn't find a discussion of this when I looked. So now the internet has one :)

I find that when connecting multiple springs together, as long as each spring adds mass to each point, then using maximum stiffness for each spring is still fine.

I think I might have defined the maximum force for the spring and damper as half the values you calculated, as I think the spring and damper forces can work together to reinforce rather than cancel each other, depending on the starting conditions.

##### Share on other sites

I find that when connecting multiple springs together, as long as each spring adds mass to each point, then using maximum stiffness for each spring is still fine.

That's a nice idea. In ny approach I have always assumed constant mass, but of course you could do it the other way around and assume constant coefficients and sum up the spring endpoint masses instead. I'll take a closer look at this and maybe implement it in a later version of the article. :-)

##### Share on other sites

h4tt3n, you would define the global timestep based on the limiting stiffness. That is, you would compute the eigen-frequency for all springs in your system, and base your global timestep on the maximum. If the stiffnesses are constant for all springs, you only ever have to do this once so its hardly impractical.

What you have done basically relies on the critical timestep. Look at the equations in the link again and you will see the critical stiffness is,

K_crit = 4m/(dt^2)

Which is suspiciously similar to the form of your coefficients in the equation that follows the "reintroducing coefficients" section, ie Cm/(dt^2) for the first term. Your formulation and guidance therefore sets the stiffness well below the critical stiffness for the chosen timestep.

I feel this article would benefit from a discussion of the critical timestep and how it relates to your work. Despite my comments I like the conceptual explanation of the coefficients relating to rigid and springy springs.

PS. Your quote from the YADE docs's refers to the case where timestep leads to very large displacements relative to the system size, your formulation would suffer problems in this case also.

##### Share on other sites

brucedjones, I'm a bit intrigued by the similarities you point out and I'm taking a serious closer look at it. But one thing is clear at this point already: A critical stiffnes of 4m / dt^2 does *not* apply to my work, that's for sure. K_max in my implementation is m / dt^2, that's been both derived mathematically and tested and tried a hundred times over. Max allowed stiffnes appears to be very integrator dependant, as buckeye pointed out in the first post of this thread, where he derived a k_max = 2m / dt^2 for second order integration. Maybe there's a simple connection, and it turns out 4m / dt^2 is max stiffnes for a 4th order integrator :-) I will thest this tonight (european time) with a number of different integration algorithms.

##### Share on other sites

As expected, it turned out that indeed for the 2nd order velocity Verlet integration algorithm K_max = 2m / dt^2 and D_max = 2m / dt. With these coefficients any pair of particles reach rest state in one loop. That doesn't make it better than the one described in the article, though - they appear to return the exact same results. It displays good stability, so this might be a path for further improvement. As for Higher order integration, I've tested the therwise excellent 4th and 6th order algorithms by David Whysong (link below) with different coefficient pairs without much sucess. They handle gravitational interaction *very* well, but show less stability when simulating springs than 1st and 2nd order integration.

##### Share on other sites

Was the example source code ever released for this?

##### Share on other sites

@Spookycat

No, I didn't get around to release any source for this. Currently I am working on a follow-up article that explains how to implement an iterative spring with really awesome properties. It is very easy to implement and - as opposed to matrix based constraints - it is very intuitive and easy to understand. There will be a demo app with source released along with the article.

Cheers, Mike

##### Share on other sites

I am working on a vehicle physics implementation of a simple suspension that will use this concept.

?Thanks for posting this.

##### Share on other sites

Great! It's nice to hear someone found the article useful :-)

##### Share on other sites

Well I got my vehicle physics implementation run works nice.

##### Share on other sites

I am trying too add this to box2d prismaticjoints (and once that works I'd also like to use it on revolute joints), but I can't get it to work. There are at least two questions that are not completely clear to me, probably because I am not very deep into physics.

At this spot in the velocity constrain calculations of box2d https://github.com/ColaColin/liquidfun/blob/master/liquidfun/Box2D/Box2D/Dynamics/Joints/b2PrismaticJoint.cpp#L293 I add this code:

	// if my spring is enabled
if (m_enableSpring && m_limitState != e_equalLimits)
{
// calculate the reduced mass
float32 redM = (1/(mA+mB));
// take the inverse timestep
float32 idt = data.step.inv_dt;
// calculate x. I am not quite sure if this is correct.
float32 x = m_springRestTranslation - GetJointTranslation();
// the impulse based formula
float32 springImpulse = -redM * idt * m_ck * x - m_cd * GetJointSpeed();
// calculate the impulse in the direction of the axis of the prismatic joint
b2Vec2 p = springImpulse * m_axis;
// apply the movement change on the bodies
// is this correct? Just apply the impulse to both bodies?
vA += mA * p;
vB -= mB * p;
}


My problem is that I am not getting it to be stable at all. Quite the opposite, it tends to massively explode. In the few moments it doesn't completely explode it does show some spring-like behavior though.

Is there anything obviously wrong with my code (especially with my assumptions as pointed out in the comments?)

It might also be that this just doesn't play nice with all the other things box2d does, I don't have a full overview of those :s

Hoping for help, thanks in advance.

##### Share on other sites
Hello ColaColin,

At a first quick glance it looks like you are not using the right equation for reduced mass:
Assuming mA and mB are not the inverse masses, 1/(mA+mB) should instead be 1/(1/mA+1/mB)

What happens if you set the stiffness and damping coeffs at a very low value?
What happens if either stiffness or damping is zero?

I am studying Erin's awesome Box2D code, but atm I can't say for sure wether the spring equation is implemented right,
or wether it's going to play nicely with the rest of box2d.

Could you please explain why you are adding the spring to the prismatic joint,
rather than implement it as a stand-alone alternative to the different constraints?

Cheers,
Mike

##### Share on other sites
Nice article and with an interesting premise! While I have yet to implement this system and play around with it, my first concerns are regarding reduced mass. I'm probably overlooking something obvious, but using that equation with a particle of zero or infinite mass, connected to a particle with mass, neither will move, right?

## Create an account

Register a new account

• 0
• 0
• 1
• 1
• 3

• 9
• 12
• 9
• 33
• 13