# Damped Springs – RyanJuckett.com

*by*Phil Tadros

I’ve labored on third-person digicam methods for quite a few video games and I typically want a technique for smoothing digicam movement because the participant strikes by way of the world. The participant’s movement might create sharp, discontinuous modifications in course, velocity and even place. The digicam, nonetheless, ought to:

- Keep away from discontinuities in movement. Speed up and decelerate as wanted relatively than snap to a brand new velocity.
- By no means let the participant outrun the digicam. The farther away he will get, the quicker the digicam wants to maneuver.
- Transfer precisely the identical with respect to time no matter body fee.

To resolve this downside, I typically simulate the digicam with a set of damped springs.

If the digicam wants to stay a relentless distance from the participant (e.g. a Mario type digicam), I can use a spring for its radial distance. If the digicam place is mounted in area and it rotates round to trace the participant, I can simulate springs for pitch and yaw.

Now that we have now some context, let’s dive into the main points and work in direction of with the ability to write a perform that may combine the movement of a damped spring. I will attempt to cowl the complete derivation of the damped harmonic movement formulation for these , however be warned that there’s a lot of math. In the event you simply need to seize the code, be at liberty to skip forward to the final web page.

Earlier than we get into damped springs, I’m going to speak about regular springs. We’ll derive the equations of movement for a standard spring with the identical strategy that we’ll use for damped springs. It will assist us be taught a few of math concerned with less complicated equations.

Your on a regular basis spring system strikes in accordance with simple harmonic motion. I’m certain you might be acquainted with how a spring behaves, however let’s spell it out anyhow for the sake of completeness. In the event you maintain a small coiled metallic spring in your hand, it is going to stay at some relaxation size. In the event you attempt to compress the spring, it is going to apply a pressure attempting to develop again to its relaxation size. In the event you stretch the spring, it is going to apply a pressure attempting to shrink again to its relaxation size. We name this relaxation size, the equilibrium place of the spring. In contrast to an actual world coiled spring, our easy harmonic movement spring may have no drag, friction or different complicated forces on it. It solely has the push and pull attributable to combating compression or extension.

One other fascinating property of springs is that the farther you compress or stretch them from their equilibrium place, the extra pressure they may apply. Based on the wants of the digicam instance I began with, that is nice as a result of the digicam will transfer again to the participant quicker if it begins to pull behind.

Let’s faux one finish of our spring is rigidly mounted to the place zero and a degree mass is connected to the opposite finish (e.g. our digicam). At relaxation, the spring has a size, (r). As the present size, (p), deviates farther from (r), we are going to generate a pressure in the other way. We are able to outline this pressure as:

(F=-k(p-r))

The worth, (okay), known as the spring constant and controls how sturdy the spring is. Increased values of (okay) will create a stronger pressure for a similar displacement from the equilibrium place, (r) .

As a way to make our lives simpler, let’s outline the variable, (x), to be the place relative to equilibrium and simplify the pressure equation.

(x = p – r)

(F = -kx)

As said earlier than, we actually need damped harmonic movement, however it is going to be a helpful train to learn the way we derive the formulation for easy harmonic movement first. We will even be taught why a easy harmonic oscillator (the spring) is just not adequate for the wants of a simulated online game digicam.

Our purpose is to discover a perform (x(t)) that solves our place, (x), primarily based on time, (t).

In the event you recall your physics, acceleration is the by-product of velocity and velocity is the by-product of place. I’ll seek advice from the place, velocity and acceleration of our equilibrium relative place, (x) , as follows.

(x = place)

(v = {dx over dt} = velocity)

(a = {dv over dt} = acceleration)

The subsequent prerequisite earlier than we will proceed is Newton’s second law of motion. That’s the one that claims that the acceleration induced by a pressure is inversely proportional to mass (i.e. huge rocks are laborious to push). That is additionally the legislation that provides us the equation (F = ma), the place (m) is the mass of our object connected to the spring.

Earlier we mentioned that pressure was equal to (-kx) as a result of spring. We are able to mix these pressure equations into one equality as follows:

(F =ma)

(F = -kx)

(ma = -kx)

(ma + kx = 0)

As a way to simplify the remaining math, we’re going rearrange the equation and outline a brand new variable, (omega). We’ll later discover that (omega) is the angular frequency of the spring.

(a + {okay over m}x = 0)

(omega = sqrt{okay over m})

(a + omega^2 x = 0)

This equation says that our place depends upon our acceleration and vice versa. It’s known as an ordinary differential equation (ODE). An ODE is an equation containing the derivatives of a single impartial variable. In our case this variable is x and its second by-product a. Our ODE can also be thought of linear and homogenous. It’s a linear ODE as a result of all of its derivatives have an exponent of 1. It’s a homogenous ODE as a result of the sum of all of the phrases containing (x) and its derivates is the same as zero (i.e. there is no such thing as a remaining fixed time period).

As a result of the very best by-product in our linear ODE is the second by-product, it’s thought of to be of second diploma. The diploma of a linear ODE determines what number of linearly independent options it has. These linearly impartial options create a foundation for an answer vector space through which any scaled mixture can also be an answer. For instance, if (s_1) and (s_2) are two linearly impartial options to the system, then (3s_1 + 5 s_2) can also be an answer to the system. Which means that there are an infinite variety of options however they’re all constrained inside the vector area.

I am not going to enter an excessive amount of element about what make a set of options linearly impartial or the way to take a look at their independence, however if you’re acquainted with linear algebra you possibly can imaging an n-dimensional vector the place n is the diploma of the ODE (in our case it’s 2). The weather of the vector are an answer and its derivatives as much as the nth diploma. Our resolution foundation is a set of those vectors – one for every resolution. For the idea to be linearly impartial, you shouldn’t be capable of derive anybody vector as a scaled mixture of the others.

The ODE’s diploma additionally determines what number of preliminary circumstances should be specified to unravel a selected resolution. This is sensible for a degree mass connected to a spring. When simulating the spring, we might want to specify our preliminary spring place and preliminary spring velocity. These are the 2 preliminary circumstances that decide how the system will behave over time.

Earlier than we will enter our preliminary place and velocity, we have to clear up an equation for the final movement of the system.

As a way to clear up for (x) relative to (t) we have to cancel the derivatives out of the equation. Fortunately our buddy, Euler, got here up with a intelligent method for doing simply this for homogenous linear ODEs. We have to clear up for the characteristic equation equivalent to the ODE. In the event you recall from calculus, the by-product of (e^u) with respect to (y) is fairly particular in that it doesn’t change its form.

({d over dy} e^u = {du over dy} e^u)

If we let (x) equal (e^{z t}), we will simplify as follows

(x = e^{z t})

(v = {dx over dt} = ze^{z t})

(a = {dv over dt} = z^2 e^{z t})

(a + omega^2 x = 0)

(z^2 e^{z t} + omega^2 e^{z t} = 0)

((z^2 + omega^2) e^{z t} = 0)

Dividing either side by (e ^ {z t}) provides us the attribute equation:

(z^2 + omega^2 = 0)

As a way to clear up for (z), we will rearrange the equation and issue right into a difference of squares. To do that, we make use of the imaginary unit, (i), which is the same as the (sqrt{-1})

(z^2 – i^2 omega^2 = 0)

((z + i omega) (z – i omega) = 0)

The factored equation exhibits us two attainable options for (z)

(z_1 = -i omega)

(z_2 = i omega)

By substituting (z) again into our authentic equations for (x), we get

(x_1 = e ^ {z_1 t} = e ^ {-i omega t})

(x_2 = e ^ {z_2 t} = e ^ {i omega t})

In the event you recall, we decided that our equation, (ma + kx = 0), was a linear ODE and thus any linear mixture of options can also be an answer. This lets us mix our two options for (x) into one equation representing the final resolution for (x(t)). The variables (c_1) and (c_2) will probably be used as linear scalars for every resolution.

(x(t) = c_1 x_1 + c_2 x_2)

(x(t) = c_1 e ^ {-i omega t} + c_2 e ^ {i omega t})

Utilizing Euler’s formula, (e^{i y} = cos (y) + i sin (y)), we will convert to a trigonometric equation

(x(t) = c_1 Huge( cos(-omega t) + i sin(-omega t) Huge) + c_2 Huge( cos(omega t) + i sin(omega t) Huge))

Utilizing the trigonometric identities (cos(-theta) = cos(theta)) and (sin(-theta) = -sin(theta)), let’s simplify and regroup phrases.

(x(t) = c_1 Huge( cos(omega t) – i sin(omega t) Huge) + c_2 Huge( cos(w t) + i sin(omega t) Huge))

(x(t) = (c_1 + c_2) cos(omega t) + i (c_2 – c_1) sin(omega t))

By taking the by-product with respect to (t), we will clear up for (v(t))

(v(t) = -omega (c1 + c2) sin(omega t) + i omega (c2 – c1) cos(omega t))

The phrases (c1) and (c2) are arbitrary fixed scalars that choose from an infinite variety of options to the unique equation, (a + omega^2 x = 0). We are able to derive their values primarily based on the preliminary state of the spring to generate a selected resolution.

First, we clear up (x) for a time worth of (t=0).

(x(0) = (c_1 + c_2) cos(omega 0) + i (c_2 – c_1) sin(omega 0))

(x(0) = (c_1 + c_2))

Subsequent, we clear up (v(t)) for a time worth of (t=0).

(v(0) = -omega (c_1 + c_2) sin(w 0) + i omega (c_2 – c_1) cos(w 0))

(v(0) = i omega (c_2 – c_1))

Primarily based on the (x(0)) and (v(0)) equations, we will clear up for the ((c_1 + c_2)) and ((c_2 – c_1)) phrases. Let’s additionally let (x_0 = x(0)) and (v_0 = v(0)) to simplify the equations a bit.

((c_1 + c_2) = x_0)

((c_2 – c_1) = {v_0 over i omega})

Now substitute the values for ((c_1 + c_2)) and ((c_2 – c_1)) again into the equation for (x(t)).

(x(t) = x_0 cos(omega t) + i {v_0 over i omega} sin(omega t))

(x(t) = x_0 cos(omega t) + {v_0 over omega} sin(omega t))

Subsequent we do the identical for(v(t)).

(v(t) = -omega x_0 sin(w t) + i omega {v_0 over i omega} cos(omega t))

(v(t) = -w x_0 sin(omega t) + v_0 cos(omega t))

Utilizing our equations for (x(t)) and (v(t)), we will simulate easy harmonic movement over an elapsed time frame. For instance, let’s contemplate the case the place we’re simulating the yaw rotation of a digicam over a single recreation body.

- Let (okay) be our spring fixed.
- Let (m) be our spring level mass.
- Let (x_0) be our present distance from the specified yaw angle (chosen to be the shortest angle across the circle).
- Let (v_0) be our present yaw-velocity.
- Let (t) be the length of the sport body.
- Let (x(t)) be our new distance from the specified yaw angle.
- Let (v(t)) be our new yaw-velocity.

We are able to compute our new yaw-position and yaw-velocity by working the next equations.

(omega = sqrt{okay over m})

(x(t) = x_0 cos(omega t) + {v_0 over omega} sin(omega t))

(v(t) = -w x_0 sin(omega t) + v_0 cos(omega t))

As an alternative of getting the consumer specify (okay) and (m), it may be extra helpful to solely outline (omega) which seems to correspond to the angular frequency of the ocillation. Which means that selecting increased values for (omega) will create the next frequency movement.

It’s value noting that in case you googled easy harmonic movement elsewhere, you’d most likely discover a totally different equation that may consider a price for (x) at time (t) primarily based on angular frequency, amplitude and part shift. Utilizing the suitable trigonometric identities, we should always be capable to convert from one resolution to the opposite.

Let’s check out the ensuing movement from simulating place over time.

(start{array}{lcl} omega & = & 1 x_0 & = & 1 v_0 & = & 0end{array}) |

We begin with none velocity however will not be at equilibrium so we begin to transfer in direction of it after which oscillate endlessly |

(start{array}{lcl} omega & = & 1 x_0 & = & 0 v_0 & = & 1end{array}) |

We begin at equilibrium, however have an preliminary velocity so we transfer away after which oscillate endlessly. |

(start{array}{lcl} omega & = & 1 x_0 & = & 0 v_0 & = & 0end{array}) |

We’re already at equilibrium so we keep there. |

As we will see from these examples, easy harmonic movement alone is just not appropriate for a recreation digicam. The digicam must settle at equilibrium relatively than capturing previous it and oscillating endlessly. It is velocity must dampen over time which brings us to the issue of damped harmonic movement

A damped spring is rather like our easy spring with an extra pressure. Within the easy spring case, we had a pressure proportional to our distance from the equilibrium place. This brought on us to at all times speed up in direction of equilibrium, however by no means come to relaxation at it (aside from the trivial case the place we began there with no movement). As a way to settle at equilibrium, we add a pressure proportional to our velocity within the case of damped springs.

(x = p – r)

(F = -beta v – kx)

As soon as once more, (p) is the present place, (r) is the equilibrium place, (okay) is the spring fixed and (v) is the present velocity. The brand new variable, (beta), is the viscous damping coefficient. Increased values of (beta) will create a stronger pressure towards the present velocity.

Utilizing Newton’s second, we will work in direction of our abnormal differential equation.

(F = ma)

(F= -kx – beta v)

(ma = -beta v – kx)

(ma + beta v + kx = 0)

Rearranging the equation and defining some new variables will simplify our lives within the following math. The (omega) variable is identical angular frequency variable we had in easy harmonic movement. The brand new variable, (zeta), is our damping ratio which I’ll talk about later.

(a + {beta over m}v + {okay over m}x = 0)

(a + {beta over sqrt m sqrt m}v + {okay over m}x = 0)

(a + 2 {beta over 2 sqrt m sqrt m }v + {okay over m}x = 0)

(a + 2 {beta sqrt okay over 2 sqrt m sqrt m sqrt okay}v + {okay over m}x = 0)

(a + 2 {sqrt okay over sqrt m} {beta over 2 sqrt m sqrt okay}v + {okay over m}x = 0)

(a + 2 sqrt{ okay over m} {beta over 2 sqrt{m okay} }v + {okay over m}x = 0)

(omega = sqrt{okay over m})

(zeta = { beta over 2 sqrt{mk} })

(a + 2 omega zeta v + omega^2 x = 0)

We now have our homogeneous linear ODE for damped easy harmonic movement.

## Basic Answer

To resolve for the final resolution of (x) we as soon as once more observe Euler’s lead and let (x = e^{zt}) to simplify.

(x = e^{z t})

(v = {dx over dt} = ze^{z t})

(a = {dv over dt} = z^2 e^{z t})

(a + 2 omega zeta v + omega^2 x = 0)

(z^2 e^{z t} + 2 omega zeta z e^{zt} + omega^2 e^{z t} = 0)

((z^2 + 2 omega zeta z + omega^2) e^{z t} = 0)

Dividing either side by (e ^ {z t}) provides us the attribute equation:

(z^2 + 2 omega zeta z + omega^2 = 0)

As a way to clear up for (z), we have to use the quadratic formula which states that the roots of an equation within the type (Az^2+Bz+C=0) are given by

(z=frac{-B pm sqrt {B^2-4AC}}{2 A })

Plugging in our coefficients we get

(A = 1)

(B = 2 omega zeta)

(C = omega^2)

(z=frac{-2 omega zeta pm sqrt {(2 omega zeta)^2-4 omega^2}}{2})

(z=frac{-2 omega zeta pm sqrt {4 omega^2 zeta^2-4 omega^2}}{2})

(z=frac{-2 omega zeta pm sqrt {4 omega^2 (zeta^2-1)}}{2})

(z=frac{-2 omega zeta pm 2 omega sqrt {zeta^2-1}}{2})

(z=-omega zeta pm omega sqrt {zeta^2-1})

This offers us two options

(z_1=-omega zeta – omega sqrt {zeta^2-1})

(z_2=-omega zeta + omega sqrt {zeta^2-1})

Each of those equations have a sq. root with a (zeta) inside it. (zeta) represents the damping ratio and defines how our movement decays over time. Relying on the worth of (zeta), we are going to find yourself with two actual roots, one actual root or two complicated roots. These correspond to over-damped, critically damped and under-damped harmonic movement respectively. We’ll clear up each in flip.

For over-damped movement, (zeta ; > ; 1). An over-damped spring won’t ever oscillate, however reaches equilibrium at a slower fee than a critically damped spring.

As a result of (zeta > 1), we all know that (z_1) and (z_2) will probably be actual numbers. Relatively than increasing them again out, let’s proceed to make use of these variables in our two linearly impartial options.

(x_1 = e^{z_1 t})

(x_2 = e^{z_2 t})

We are able to now construct our normal resolution to the linear ODE as a liner mixture of our impartial options by including the scalar constants (c_1) and (c_2).

(x(t) = c_1 x_1 + c_2 x_2)

(x(t) = c_1 e^{z_1 t} + c_2 e^{z_2 t})

Taking the by-product provides us the final resolution for velocity.

(v(t) = c_1 z_1 e^{z_1 t} + c_2 z_2 e^{z_2 t})

As a way to get the specified explicit resolution, we will clear up for the constants primarily based on our preliminary state.

First, we clear up (x(t)) for a time worth of (t=0).

(x(0) = c_1 e^{z_1 0} + c_2 e^{z_2 0})

(x(0) = c_1 + c_2)

Subsequent, we clear up (v(t)) for a time worth of (t=0).

(v(0) = c_1 z_1 e^{z_1 0} + c_2 z_2 e^{z_2 0})

(v(0) = c_1 z_1 + c_2 z_2)

Primarily based on the (x(0)) and (v(0)) equations, we will clear up for (c_1) and (c_2). Let’s additionally let (x_0 = x(0)) and (v_0 = v(0)) to simplify the equations a bit.

(c_2 = x_0 – c_1)

(v_0 = c_1 z_1 + (x_0 – c_1) z_2)

(v_0 = c_1 z_1 + x_0 z_2 – c_1 z_2)

(v_0 = x_0 z_2 + c_1 (z_1 – z_2))

(c_1 = {v_0 – x_0 z_2 over z_1 – z_2})

(c_2 = x_0 – {v_0 – x_0 z_2 over z_1 – z_2})

Now substitute the values for (c_1) and (c_2) again into the equation for (x(t)).

(x(t) = {v_0 – x_0 z_2 over z_1 – z_2} e^{z_1 t} ; + ; (x_0 – {v_0 – x_0 z_2 over z_1 – z_2}) e^{z_2 t})

Subsequent we do the identical for(v(t)).

(v(t) = {v_0 – x_0 z_2 over z_1 – z_2} z_1 e^{z_1 t} ; + ; (x_0 – {v_0 – x_0 z_2 over z_1 – z_2}) z_2 e^{z_2 t})

These equations allow us to calculate a brand new place and velocity for an over-damped spring primarily based on an elapsed time from an preliminary place and velocity.

Let’s check out the ensuing movement from simulating place over time.

(start{array}{lcl} omega & = & 1 zeta & = & 2 x_0 & = & 1 v_0 & = & 0end{array}) |

We begin with none velocity however will not be at equilibrium so we begin to transfer in direction of it. On account of being over-damped, it takes us a very long time to succeed in equilibrium. |

(start{array}{lcl} omega & = & 1 zeta & = & 2 x_0 & = & 0.5 v_0 & = & 1end{array}) |

We begin with velocity transferring away from equilibrium, however our movement shortly dampens out. On account of being over-damped, it takes us a very long time to succeed in equilibrium. |

For critically damped movement, (zeta = 1). A critically damped spring will attain equilibrium as quick as attainable with out oscillating.

Let’s clear up for (z_1) and (z_2).

(z_1 = -omega zeta – omega sqrt {zeta^2-1} = -omega 1 – omega sqrt {1-1} = -omega)

(z_2 = -omega zeta + omega sqrt {zeta^2-1} = -omega 1 + omega sqrt {1-1} = -omega)

We now have an issue right here as a result of (z_1) and (z_2) are each equal to the identical worth. These are known as repeated roots of the attribute equation. If we had been to substitute (z_1) and (z_2) again into our authentic equations for (x_1) and (x_2), we’d get

(x_1 = e^{z_1 t} = e^{-omega t})

(x_2 = e^{z_2 t} = e^{-omega t})

That is usually the place we’d mix (x_1) and (x_2) right into a normal resolution, however that will not work right here as a result of they don’t seem to be linearly impartial options. They’re the identical resolution. I mentioned beforehand {that a} second diploma linear ODE should have two linearly impartial options and that also stands. In the event you had been to observe the maths alongside from right here below the belief that we solely have one resolution, you’d finally hit a wall the place your equations grew to become nonsensical.

So, we all know one resolution, however how do we discover our second resolution? We are able to use a technique known as reduction of order, which is credited to Jean le Rond d’Alembert. This is the way it works. Let’s guess that our normal resolution is a few perform of (t) multiplied by our present resolution.

(x(t) = f(t) x_1 = f(t) e^{-omega t})

Now we are going to attempt to clear up for (f(t)) by fixing the rate and acceleration derivatives of (x(t)). To do that, we have to use the product rule. Let (f^prime(t)) and (f^{prime prime}(t)) be the primary and second derivatives of (f(t)) respectively.

(v(t) = {d over dt}x(t))

(v(t) = f^prime(t) e^{-omega t} – f(t) omega e^{-omega t})

(a(t) = {d over dt}v(t))

(a(t) = f^{prime prime}(t) e^{-omega t} – f^{prime}(t) omega e^{-omega t} – f^{prime}(t) omega e^{-omega t} + f(t) omega^2 e^{-omega t})

(a(t) = f^{prime prime}(t) e^{-omega t} – 2 f^{prime}(t) omega e^{-omega t} + f(t) omega^2 e^{-omega t})

Simplify by factoring out (e^{-omega t}).

(v(t) = e^{-omega t} Huge(f^prime(t) – f(t) omega Huge))

(a(t) = e^{-omega t} Huge(f^{prime prime}(t) -2 f^{prime}(t) omega + f(t) omega^2 Huge))

Now substitute again into our authentic equation of (a + 2 omega zeta v + omega^2 x = 0) and simplify

(a(t) + 2 omega zeta v(t) + omega^2 x(t) = 0)

(e^{-omega t} Huge(f^{prime prime}(t) -2 f^{prime}(t) omega + f(t) omega^2 Huge) ; + ; 2 omega zeta e^{-omega t} Huge(f^prime(t) – f(t) omega Huge) ; + ; omega^2 f(t) e^{-omega t} = 0)

(e^{-omega t} Huge(f^{prime prime}(t) ; – ; 2 f^{prime}(t) omega ; + ; f(t) omega^2 ; + ; 2 omega zeta f^prime(t) ; – ; 2 omega^2 zeta f(t) ; + ; omega^2 f(t) Huge) = 0)

(e^{-omega t} Huge(f^{prime prime}(t) ; – ; 2 f^{prime}(t) omega ; + ; 2 f(t) omega^2 ; + ; 2 omega zeta f^prime(t) ; – ; 2 omega^2 zeta f(t) Huge) = 0)

(e^{-omega t} Huge(f^{prime prime}(t) ; + ; 2 f^{prime}(t) omega(zeta – 1) ; + ; 2 f(t) omega^2 ; – ; 2 omega^2 zeta f(t) Huge) = 0)

(e^{-omega t} Huge(f^{prime prime}(t) ; + ; 2 f^{prime}(t) omega(zeta – 1) ; + ; 2 f(t) omega^2 ( 1 – zeta) Huge) = 0)

As a result of we’re within the critically damped case, we all know (zeta = 1) and might simplify additional.

(e^{-omega t} Huge(f^{prime prime}(t) ; + ; 2 f^{prime}(t) omega(1 – 1) ; + ; 2 f(t) omega^2 ( 1 – 1) Huge) = 0)

(e^{-omega t} f^{prime prime}(t) = 0)

The primary time period, (e^{-omega t}), won’t ever be zero whatever the values for (omega) and (t). Thus we will divide it out of the equation to get:

(f^{prime prime}(t) = 0)

So we all know our capabilities second by-product is zero, however what we actually need to know is the perform itself. So, let’s combine again up.

(f^{prime prime}(t) = 0)

(f^{prime}(t) = c_1)

(f(t) = c_1 t ; + ; c_2)

Now, let’s plug our resolution for (f(x)) again into our guess on the normal resolution to the system.

(x(t) = f(t) e^{-omega t})

(x(t) = (c_1 t + c_2) e^{-omega t})

As anticipated, we have now discovered a normal resolution that’s the sum of two linearly impartial options scaled by arbitrary constants.

Taking the by-product provides us the final resolution for velocity.

(v(t) = c_1 e^{-omega t} – (c_1 t + c_2) omega e^{-omega t})

(v(t) = Huge(c_1 – (c_1 t + c_2) omega Huge) e^{-omega t})

(v(t) = Huge( c_1(1 – omega t) – c_2 omega Huge) e^{-omega t})

As a way to get the specified explicit resolution, we will clear up for the constants primarily based on our preliminary state.

First, we clear up (x(t)) for a time worth of (t=0).

(x(0) = (c_1 0 + c_2) e^{-omega 0})

(x(0) = c_2)

Subsequent, we clear up (v(t)) for a time worth of (t=0).

(v(0) = Huge( c_1(1 – omega 0) – c_2 omega Huge) e^{-omega 0})

(v(0) = c_1 – c_2 omega)

Primarily based on the (x(0)) and (v(0)) equations, we will clear up for (c_1) and (c_2). Let’s additionally let (x_0 = x(0)) and (v_0 = v(0)) to simplify the equations a bit.

(c_2 = x_0)

(v_0 = c_1 – x_0 omega)

(c_1 = v_0 + x_0 omega)

Now substitute the values for (c_1) and (c_2) again into the equation for (x(t)).

(x(t) = Huge( (v_0 + x_0 omega) t + x_0 Huge) e^{-omega t})

Subsequent we do the identical for (v(t)).

(v(t) = Huge( (v_0 + x_0 omega)(1 – omega t) – x_0 omega Huge) e^{-omega t})

(v(t) = Huge( v_0 + x_0 omega – v_0 omega t – x_0 omega^2 t – x_0 omega Huge) e^{-omega t})

(v(t) = Huge( v_0 – (v_0 + x_0 omega) omega t Huge) e^{-omega t})

These equations allow us to calculate a brand new place and velocity for a critically damped spring primarily based on an elapsed time from an preliminary place and velocity.

Let’s check out the ensuing movement from simulating place over time.

(start{array}{lcl} omega & = & 1 zeta & = & 1 x_0 & = & 1 v_0 & = & 0end{array}) |

We begin with none velocity however will not be at equilibrium so we begin to transfer in direction of it. On account of being critically-damped, our place and velocity attain equilibrium as quickly as attainable for our given angular frequency. |

(start{array}{lcl} omega & = & 1 zeta & = & 1 x_0 & = & 0.5 v_0 & = & 1end{array}) |

We begin with velocity transferring away from equilibrium, however our movement shortly dampens out. On account of being critically-damped, our place and velocity attain equilibrium as quickly as attainable for our given angular frequency. |

For under-damped movement, (0 ; le ; zeta ; < ; 1). An under-damped spring will attain equilibrium the quickest, but additionally overshoots it and continues to oscillate as its amplitude decays over time.

As a result of (zeta ; < ; 1), we all know that (sqrt{zeta^2 – 1}) is a fancy quantity. By rearranging, we will expose the imaginary time period as follows.

(sqrt{zeta^2 – 1} = sqrt{-1 * -(zeta^2 – 1)} = sqrt{-1} sqrt{-(zeta^2 – 1)} = i sqrt{1 – zeta^2})

Utilizing this equality, we will replace our roots.

(z=-omega zeta pm omega i sqrt{1 – zeta^2)})

Earlier than persevering with, let’s make our equations a bit much less verbose by defining a brand new variable, (alpha).

(alpha = omega sqrt{1 – zeta^2})

(z=-omega zeta pm i alpha)

Plugging again into our equation for (x), we get:

(x = e^{z t})

(x = e^{(-omega zeta pm i alpha) t})

Breaking apart the (pm) into our two options, we have now:

(x_1 = e^{-omega zeta t + i alpha t})

(x_2 = e^{-omega zeta t – i alpha t})

If we break up the exponent, we will apply Euler’s formulation.

(x_1 = e^{-omega zeta t} e^{i alpha t})

(x_1 = e^{-omega zeta t} Huge( cos(alpha t) + i sin(alpha t) Huge))

(x_2 = e^{-omega zeta t} e^{-ialpha t})

(x_2 = e^{-omega zeta t} Huge( cos(-alpha t) + i sin(-alpha t) Huge))

Utilizing the trig identities (cos(-theta)=cos(theta)) and (sin(-theta)=-sin(theta)), we will simplify (x_2).

(x_2 = e^{-omega zeta t} Huge( cos(alpha t) – i sin(alpha t) Huge))

We now have two linearly impartial options to the ODE, however they nonetheless include imaginary phrases. We’re solely excited by actual options so it could be preferable to take away the imaginary components prior to later. To do that we are going to generate two new foundation equations as linear sums of our current equations. We’ll select scaling constants that cancel out the imaginary whereas nonetheless getting linearly impartial outcomes, (x_3) and (x_4).

(x_3 = {1 over 2} x_1 + {1 over 2} x_2)

(x_3 = {1 over 2} e^{-omega zeta t} Huge( cos(alpha t) + i sin(alpha t) Huge) + {1 over 2} e^{-omega zeta t} Huge( cos(alpha t) – i sin(alpha t) Huge))

(x_3 = {1 over 2} e^{-omega zeta t} Huge( cos(alpha t) + i sin(alpha t) + cos(alpha t) – i sin(alpha t) Huge))

(x_3 = e^{-omega zeta t} cos(alpha t))

(x_4 = {1 over {2 i}} x_1 – {1 over {2 i}} x_2)

(x_4 = {1 over {2 i}} e^{-omega zeta t} Huge( cos(alpha t) + i sin(alpha t) Huge) – {1 over {2 i}} e^{-omega zeta t} Huge( cos(alpha t) – i sin(alpha t) Huge))

(x_4 = {1 over {2 i}} e^{-omega zeta t} Huge( cos(alpha t) + i sin(alpha t) – cos(alpha t) + i sin(alpha t) Huge))

(x_4 = e^{-omega zeta t} sin(alpha t))

We are able to now construct our normal resolution to the linear ODE as a liner mixture of our impartial options by including the scalar constants (c_1) and (c_2).

(x(t) = c_1 x_3 + c_2 x_4)

(x(t) = c_1 e^{-omega zeta t} cos(alpha t) + c_2 e^{-omega zeta t} sin(alpha t))

(x(t) = e^{-omega zeta t} Huge( c_1 cos(alpha t) + c_2 sin(alpha t) Huge))

Taking the by-product provides us the final resolution for velocity.

(v(t) = -omega zeta e^{-omega zeta t} Huge( c_1 cos(alpha t) + c_2 sin(alpha t) Huge) + e^{-omega zeta t} Huge( -c_1 alpha sin(alpha t) + c_2 alpha cos(alpha t) Huge))

(v(t) = -e^{-omega zeta t} Huge( omega zeta c_1 cos(alpha t) + omega zeta c_2 sin(alpha t) + c_1 alpha sin(alpha t) – c_2 alpha cos(alpha t) Huge))

(v(t) = -e^{-omega zeta t} Huge( (c_1 omega zeta – c_2 alpha) cos(alpha t) + ( c_1 alpha + c_2 omega zeta) sin(alpha t) Huge))

As a way to get the specified explicit resolution, we will clear up for the constants primarily based on our preliminary state.

First, we clear up (x(t)) for a time worth of (t=0).

(x(0) = e^{-omega zeta 0} Huge( c_1 cos(alpha 0) + c_2 sin(alpha 0) Huge))

(x(0) = c_1)

Subsequent, we clear up (v(t)) for a time worth of (t=0).

(v(0) = -e^{-omega zeta 0} Huge( (c_1 omega zeta – c_2 alpha) cos(alpha 0) + ( c_1 alpha + c_2 omega zeta) sin(alpha 0) Huge))

(v(0) = -(c_1 omega zeta – c_2 alpha))

(v(0) = -c_1 omega zeta + c_2 alpha)

Primarily based on the (x(0)) and (v(0)) equations, we will clear up for (c_1) and (c_2). Let’s additionally let (x_0 = x(0)) and (v_0 = v(0)) to simplify the equations a bit.

(c_1 = x_0)

(v_0 = -x_0 omega zeta + c_2 alpha)

(c_2 = { v_0 + omega zeta x_0 over alpha})

Now substitute the values for (c_1) and (c_2) again into the equation for (x(t)).

(x(t) = e^{-omega zeta t} Huge( x_0 cos(alpha t) ; + ; { v_0 + omega zeta x_0 over alpha} sin(alpha t) Huge))

Subsequent we do the identical for(v(t)).

(v(t) = -e^{-omega zeta t} Huge( (x_0 omega zeta ; – ; { v_0 + omega zeta x_0 over alpha} alpha) cos(alpha t) ; + ; (x_0 alpha ; + ; { v_0 + omega zeta x_0 over alpha} omega zeta) sin(alpha t) Huge))

These equations allow us to calculate a brand new place and velocity for an under-damped spring primarily based on an elapsed time from an preliminary place and velocity.

Let’s check out the ensuing movement from simulating place over time.

(start{array}{lcl} omega & = & 1 zeta & = & 0.5 x_0 & = & 1 v_0 & = & 0end{array}) |

We begin with none velocity however will not be at equilibrium so we begin to transfer in direction of it. On account of being under-damped, we proceed to overshoot equilibrium with every oscillation, however we overshoot it much less and fewer every time. |

(start{array}{lcl} omega & = & 1 zeta & = & 0.5 x_0 & = & 0.5 v_0 & = & 1end{array}) |

We begin with velocity transferring away from equilibrium, however our movement shortly dampens out. On account of being under-damped, we proceed to overshoot equilibrium with every oscillation, however we overshoot it much less and fewer every time. |

We are able to now program a damped easy harmonic oscillator. As a result of simulating damped springs requires calls to probably costly trigonometric and exponential capabilities, I’ve break up the method into two steps. The primary computes a set of coefficients for the place and velocity parameters by increasing the related equations. These coefficients can then be used to shortly replace a number of springs utilizing the identical angular frequency, damping ratio and time step. In case your simulation updates at a locked time step, you possibly can even cache off the coefficients as soon as at initialization time and use them each body!

These code samples are launched below the next license.

```
/******************************************************************************
Copyright (c) 2008-2012 Ryan Juckett
http://www.ryanjuckett.com/
This software program is supplied 'as-is', with none specific or implied
guarantee. In no occasion will the authors be held answerable for any damages
arising from the usage of this software program.
Permission is granted to anybody to make use of this software program for any goal,
together with industrial purposes, and to change it and redistribute it
freely, topic to the next restrictions:
1. The origin of this software program should not be misrepresented; you have to not
declare that you just wrote the unique software program. In the event you use this software program
in a product, an acknowledgment within the product documentation could be
appreciated however is just not required.
2. Altered supply variations should be plainly marked as such, and should not be
misrepresented as being the unique software program.
3. This discover might not be eliminated or altered from any supply
distribution.
******************************************************************************/
```

Name the CalcDampedSpringMotionParams perform as soon as along with your parameters to initialize a tDampedSpringMotionParams construction. Then name UpdateDampedSpringMotion with pointers referencing the place and velocity values for every spring that wants updating.

```
//******************************************************************************
// Cached set of movement parameters that can be utilized to effectively replace
// a number of springs utilizing the identical time step, angular frequency and damping
// ratio.
//******************************************************************************
struct tDampedSpringMotionParams
{
// newPos = posPosCoef*oldPos + posVelCoef*oldVel
float m_posPosCoef, m_posVelCoef;
// newVel = velPosCoef*oldPos + velVelCoef*oldVel
float m_velPosCoef, m_velVelCoef;
};
//******************************************************************************
// This perform will compute the parameters wanted to simulate a damped spring
// over a given time frame.
// - An angular frequency is given to regulate how briskly the spring oscillates.
// - A damping ratio is given to regulate how briskly the movement decays.
// damping ratio > 1: over damped
// damping ratio = 1: critically damped
// damping ratio < 1: below damped
//******************************************************************************
void CalcDampedSpringMotionParams(
tDampedSpringMotionParams* pOutParams, // movement parameters consequence
float deltaTime, // time step to advance
float angularFrequency, // angular frequency of movement
float dampingRatio) // damping ratio of movement
{
const float epsilon = 0.0001f;
// pressure values into authorized vary
if (dampingRatio < 0.0f) dampingRatio = 0.0f;
if (angularFrequency < 0.0f) angularFrequency = 0.0f;
// if there is no such thing as a angular frequency, the spring won't transfer and we will
// return identification
if ( angularFrequency < epsilon )
{
pOutParams->m_posPosCoef = 1.0f; pOutParams->m_posVelCoef = 0.0f;
pOutParams->m_velPosCoef = 0.0f; pOutParams->m_velVelCoef = 1.0f;
return;
}
if (dampingRatio > 1.0f + epsilon)
{
// over-damped
float za = -angularFrequency * dampingRatio;
float zb = angularFrequency * sqrtf(dampingRatio*dampingRatio - 1.0f);
float z1 = za - zb;
float z2 = za + zb;
float e1 = expf( z1 * deltaTime );
float e2 = expf( z2 * deltaTime );
float invTwoZb = 1.0f / (2.0f*zb); // = 1 / (z2 - z1)
float e1_Over_TwoZb = e1*invTwoZb;
float e2_Over_TwoZb = e2*invTwoZb;
float z1e1_Over_TwoZb = z1*e1_Over_TwoZb;
float z2e2_Over_TwoZb = z2*e2_Over_TwoZb;
pOutParams->m_posPosCoef = e1_Over_TwoZb*z2 - z2e2_Over_TwoZb + e2;
pOutParams->m_posVelCoef = -e1_Over_TwoZb + e2_Over_TwoZb;
pOutParams->m_velPosCoef = (z1e1_Over_TwoZb - z2e2_Over_TwoZb + e2)*z2;
pOutParams->m_velVelCoef = -z1e1_Over_TwoZb + z2e2_Over_TwoZb;
}
else if (dampingRatio < 1.0f - epsilon)
{
// under-damped
float omegaZeta = angularFrequency * dampingRatio;
float alpha = angularFrequency * sqrtf(1.0f - dampingRatio*dampingRatio);
float expTerm = expf( -omegaZeta * deltaTime );
float cosTerm = cosf( alpha * deltaTime );
float sinTerm = sinf( alpha * deltaTime );
float invAlpha = 1.0f / alpha;
float expSin = expTerm*sinTerm;
float expCos = expTerm*cosTerm;
float expOmegaZetaSin_Over_Alpha = expTerm*omegaZeta*sinTerm*invAlpha;
pOutParams->m_posPosCoef = expCos + expOmegaZetaSin_Over_Alpha;
pOutParams->m_posVelCoef = expSin*invAlpha;
pOutParams->m_velPosCoef = -expSin*alpha - omegaZeta*expOmegaZetaSin_Over_Alpha;
pOutParams->m_velVelCoef = expCos - expOmegaZetaSin_Over_Alpha;
}
else
{
// critically damped
float expTerm = expf( -angularFrequency*deltaTime );
float timeExp = deltaTime*expTerm;
float timeExpFreq = timeExp*angularFrequency;
pOutParams->m_posPosCoef = timeExpFreq + expTerm;
pOutParams->m_posVelCoef = timeExp;
pOutParams->m_velPosCoef = -angularFrequency*timeExpFreq;
pOutParams->m_velVelCoef = -timeExpFreq + expTerm;
}
}
//******************************************************************************
// This perform will replace the provided place and velocity values over
// in accordance with the movement parameters.
//******************************************************************************
void UpdateDampedSpringMotion(
float* pPos, // place worth to replace
float* pVel, // velocity worth to replace
const float equilibriumPos, // place to strategy
const tDampedSpringMotionParams& params) // movement parameters to make use of
{
const float oldPos = *pPos - equilibriumPos; // replace in equilibrium relative area
const float oldVel = *pVel;
(*pPos) = oldPos*params.m_posPosCoef + oldVel*params.m_posVelCoef + equilibriumPos;
(*pVel) = oldPos*params.m_velPosCoef + oldVel*params.m_velVelCoef;
}
```