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.

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

https://yade-dem.org/doc/formulation.html

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.

https://yade-dem.org/doc/formulation.html#non-elastic-constraints

Cheers Mike

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

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.

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

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.

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.

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.

http://read.pudn.com/downloads72/sourcecode/math/261769/symplectic.cpp__.htm

Was the example source code ever released for this?

@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

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

?Thanks for posting this.

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

Well I got my vehicle physics implementation run works nice.

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.

One issue with parameterizing the damping like this is that you'll find that you need to change the damping coefficient every time you change the stiffness. It's easier to define the damping coefficient using the damping ratio. (https://en.wikipedia.org/wiki/Damping_ratio)

i.e.

d = cd * 2 * Sqrt(k * m)

Once you expand out for the definition of k using the dimensionless coefficient ck, you get the following:

d = cd * 2 * m / dt * Sqrt(ck)

Here is the updated and expanded version of the article. It clarifies a few things and includes code samples.

Stable springs 0.7.3.pdf

On 7/28/2019 at 9:03 PM, h4tt3n said:

Here is the updated and expanded version of the article. It clarifies a few things and includes code samples.

Stable springs 0.7.3.pdf

Hi @h4tt3n, thank you for the amazing article!  I implemented your linear spring equation in Unity and it works amazingly. Very helpful to a complete physics beginner like me. But I was sad to see that you removed your angular spring info from your latest PDF version.

I'm now trying to implement an angular spring with your equation, and I'm having trouble as Unity provides the inertia tensor as a Vector3 diagonal along with a Quaternion inertia tensor rotation, instead of a 3x3 matrix. I believe this is the way to apply the inertia tensor to torque using Unity's provided data: 

            var q = Rigidbody.rotation * Rigidbody.inertiaTensorRotation;
            var T = q * Vector3.Scale(Rigidbody.inertiaTensor, (Quaternion.Inverse(q) * torque));

I'm unsure of how to use the inertia tensor diagonal/rotation with your angular spring equation, as I assume your equation expects 'I' to be a 3x3 matrix? Any advice on how I can implement your angular spring equation with the inertia tensor diagonal/rotation?

Thank you!