Or: Implementing a damped spring model in 16 lines of code.
Governing Equation:
m x ¨ + c x ˙ + k x = f ( t ) m\ddot{x}+c\dot{x}+kx=f(t) m x ¨ + c x ˙ + k x = f ( t ) Solve
0 = m x ¨ + c x ˙ + k x x o = x ( t o ) x ˙ o = x ˙ ( t o ) \begin{aligned}
0 & = m\ddot{x} + c\dot{x} + kx \\
x_o & = x(t_o) \\
\dot{x}_o & = \dot{x}(t_o)
\end{aligned} 0 x o x ˙ o = m x ¨ + c x ˙ + k x = x ( t o ) = x ˙ ( t o ) Exponential ansatz:
x ( t ) = A e λ t x ˙ ( t ) = A e λ t λ x ¨ ( t ) = A e λ t λ 2 \begin{aligned}
x(t) &= A e^{\lambda t} \\
\dot{x}(t) &= A e^{\lambda t} \lambda \\
\ddot{x}(t) &= A e^{\lambda t} \lambda^2
\end{aligned} x ( t ) x ˙ ( t ) x ¨ ( t ) = A e λ t = A e λ t λ = A e λ t λ 2 Plugging in gives
0 = m λ 2 ( A e λ t ) + c λ ( A e λ t ) + k ( A e λ t ) = m λ 2 + c λ + k \begin{aligned}
0 &= m\lambda^2 (Ae^{\lambda t}) + c\lambda(Ae^{\lambda t}) + k(Ae^{\lambda t}) \\
&= m\lambda^2 + c\lambda + k
\end{aligned} 0 = m λ 2 ( A e λ t ) + c λ ( A e λ t ) + k ( A e λ t ) = m λ 2 + c λ + k ⇒ λ = − c / m 2 ± ( c / m 2 ) 2 − k / m = − μ 2 ± ( μ 2 ) 2 − ω 2 = − μ 2 ± μ 2 1 − ( 2 ω μ ) 2 = − λ ˉ ± σ = λ ˉ ( − 1 ± 1 − ζ 2 ) \begin{aligned}
\Rightarrow
\lambda &= -\frac{c/m}{2} \pm \sqrt{ \left(\frac{c/m}{2}\right)^2 - k/m } \\
&= -\frac{\mu}{2} \pm \sqrt{ \left(\frac{\mu}{2}\right)^2 - \omega^2} \\
&= -\frac{\mu}{2} \pm \frac{\mu}{2}\sqrt{ 1 - \left(\frac{2\omega}{\mu}\right)^2 } \\
&= -\bar{\lambda} \pm \sigma \\
&= \bar{\lambda} \left( -1 \pm \sqrt{1-\zeta^2} \right)
\end{aligned} ⇒ λ = − 2 c / m ± ( 2 c / m ) 2 − k / m = − 2 μ ± ( 2 μ ) 2 − ω 2 = − 2 μ ± 2 μ 1 − ( μ 2 ω ) 2 = − λ ˉ ± σ = λ ˉ ( − 1 ± 1 − ζ 2 ) Two solutions $x_1(t) = \left(A_1 e^{\sigma t}\right) e^{-\bar{\lambda}t}$ and $x_2(t) = \left(A_2e^{-\sigma t}\right) e^{-\bar{\lambda}t}$
To solve for $A_1$ and $A_2$ we need two initial condition $x(t_o)=x_o$ and $\dot{x}(t_o)=\dot{x}_o$ , i.e.
x ( t ) = e − λ ˉ t ( A 1 e σ t + A 2 e − σ t ) x ˙ ( t ) = e − λ ˉ t ( ( − λ ˉ + σ ) A 1 e σ t + ( − λ ˉ − σ ) A 2 e − σ t ) → ( x x ˙ ) ( t ) = ( 1 1 − λ ˉ + σ − λ ˉ − σ ) ( e σ t 0 0 e − σ t ) ( A 1 A 2 ) e − λ ˉ t \begin{aligned}
x (t) &= e^{-\bar{\lambda}t}
\left(
A_1 e^{\sigma t}
+ A_2 e^{-\sigma t}
\right) \\
\dot{x}(t) &= e^{-\bar{\lambda}t}
\left(
(-\bar{\lambda}+\sigma) A_1 e^{ \sigma t}
+ (-\bar{\lambda}-\sigma) A_2 e^{-\sigma t}
\right) \\
\rightarrow
\begin{pmatrix}
x \\
\dot{x}
\end{pmatrix}(t)
&=
\begin{pmatrix}
1 & 1 \\
-\bar{\lambda}+\sigma & -\bar{\lambda}-\sigma
\end{pmatrix}
\begin{pmatrix}
e^{\sigma t} & 0 \\
0 & e^{-\sigma t}
\end{pmatrix}
\begin{pmatrix}
A_1 \\
A_2
\end{pmatrix}
e^{-\bar{\lambda}t}
\end{aligned} x ( t ) x ˙ ( t ) → ( x x ˙ ) ( t ) = e − λ ˉ t ( A 1 e σ t + A 2 e − σ t ) = e − λ ˉ t ( ( − λ ˉ + σ ) A 1 e σ t + ( − λ ˉ − σ ) A 2 e − σ t ) = ( 1 − λ ˉ + σ 1 − λ ˉ − σ ) ( e σ t 0 0 e − σ t ) ( A 1 A 2 ) e − λ ˉ t Let $t_o=0$, then all the exponentials are 1, so it simplifies to
( x o x o ˙ ) = ( 1 1 − λ ˉ + σ − λ ˉ − σ ) ( A 1 A 2 ) ⇒ ( A 1 A 2 ) = 1 − 2 σ ( − λ ˉ − σ − 1 λ ˉ − σ 1 ) ( x o x o ˙ ) \begin{aligned}
\begin{pmatrix}
x_o \\
\dot{x_o}
\end{pmatrix}
&=
\begin{pmatrix}
1 & 1 \\
-\bar{\lambda}+\sigma & -\bar{\lambda}-\sigma
\end{pmatrix}
\begin{pmatrix}
A_1 \\
A_2
\end{pmatrix} \\
\Rightarrow
\begin{pmatrix}
A_1 \\
A_2
\end{pmatrix}
&=
\frac{ 1 }
{-2\sigma}
\begin{pmatrix}
-\bar{\lambda}-\sigma & -1 \\
\bar{\lambda}-\sigma & 1
\end{pmatrix}
\begin{pmatrix}
x_o \\
\dot{x_o}
\end{pmatrix}
\end{aligned} ( x o x o ˙ ) ⇒ ( A 1 A 2 ) = ( 1 − λ ˉ + σ 1 − λ ˉ − σ ) ( A 1 A 2 ) = − 2 σ 1 ( − λ ˉ − σ λ ˉ − σ − 1 1 ) ( x o x o ˙ ) Altogether it becomes
( x x ˙ ) ( t ) = ( 1 1 − λ ˉ + σ − λ ˉ − σ ) ( e σ t 0 0 e − σ t ) ( − λ ˉ − σ − 1 λ ˉ − σ 1 ) ( x o x o ˙ ) e − λ ˉ t − 2 σ \begin{pmatrix}
x \\
\dot{x}
\end{pmatrix}(t)
=
\begin{pmatrix}
1 & 1 \\
-\bar{\lambda}+\sigma & -\bar{\lambda}-\sigma
\end{pmatrix}
\begin{pmatrix}
e^{\sigma t} & 0 \\
0 & e^{-\sigma t}
\end{pmatrix}
\begin{pmatrix}
-\bar{\lambda}-\sigma & -1 \\
\bar{\lambda}-\sigma & 1
\end{pmatrix}
\begin{pmatrix}
x_o \\
\dot{x_o}
\end{pmatrix}
\frac{ e^{-\bar{\lambda}t} }
{ -2\sigma } ( x x ˙ ) ( t ) = ( 1 − λ ˉ + σ 1 − λ ˉ − σ ) ( e σ t 0 0 e − σ t ) ( − λ ˉ − σ λ ˉ − σ − 1 1 ) ( x o x o ˙ ) − 2 σ e − λ ˉ t Which can be rearranged to
( x x ˙ ) ( t ) = ( e σ t + e − σ t 2 + λ ˉ e σ t − e − σ t 2 σ e σ t − e − σ t 2 σ − ( λ ˉ 2 − σ 2 ) e σ t − e − σ t 2 σ e σ t + e − σ t 2 − λ ˉ e σ t − e − σ t 2 σ ) ( x o x o ˙ ) e − λ ˉ t \begin{pmatrix}
x \\
\dot{x}
\end{pmatrix}(t)
=
\begin{pmatrix}
\frac{e^{\sigma t}+e^{-\sigma t}}
{ 2 }
+
\bar{\lambda}
\frac{e^{\sigma t}-e^{-\sigma t}}
{ 2\sigma }
&
\frac{e^{\sigma t}-e^{-\sigma t}}{2\sigma}
\\
-(\bar{\lambda}^2-\sigma^2)
\frac{ e^{\sigma t} - e^{-\sigma t} }
{ 2\sigma }
&
\frac{e^{\sigma t}+e^{-\sigma t}}
{ 2 }
-\bar{\lambda}
\frac{e^{\sigma t}-e^{-\sigma t}}
{ 2\sigma }
\end{pmatrix}
\begin{pmatrix}
x_o \\
\dot{x_o}
\end{pmatrix}
e^{-\bar{\lambda}t} ( x x ˙ ) ( t ) = ( 2 e σ t + e − σ t + λ ˉ 2 σ e σ t − e − σ t − ( λ ˉ 2 − σ 2 ) 2 σ e σ t − e − σ t 2 σ e σ t − e − σ t 2 e σ t + e − σ t − λ ˉ 2 σ e σ t − e − σ t ) ( x o x o ˙ ) e − λ ˉ t
Click to see derivation
$$
\begin{aligned}
\begin{pmatrix}
x \\
\dot{x}
\end{pmatrix}
(t)
&=
\begin{pmatrix}
1 & 1 \\
-\bar{\lambda}+\sigma & -\bar{\lambda}-\sigma
\end{pmatrix}
\begin{pmatrix}
e^{\sigma t} & 0 \\
0 & e^{-\sigma t}
\end{pmatrix}
\begin{pmatrix}
-\bar{\lambda}-\sigma & -1 \\
\bar{\lambda}-\sigma & 1
\end{pmatrix}
\begin{pmatrix}
x_o \\
\dot{x_o}
\end{pmatrix}
\frac{e^{-\bar{\lambda}t}}
{ -2\sigma } \\
&=
\begin{pmatrix}
1 & 1 \\
-\bar{\lambda}+\sigma & -\bar{\lambda}-\sigma
\end{pmatrix}
\begin{pmatrix}
(-\bar{\lambda}-\sigma)e^{ \sigma t} & -e^{ \sigma t} \\
( \bar{\lambda}-\sigma)e^{-\sigma t} & e^{-\sigma t}
\end{pmatrix}
\begin{pmatrix}
x_o \\
\dot{x_o}
\end{pmatrix}
\frac{e^{-\bar{\lambda}t}}
{ -2\sigma } \\
&=
\begin{pmatrix}
(-\bar{\lambda}-\sigma)e^{\sigma t}
+ e^{-\sigma t}(\bar{\lambda}-\sigma) &
(-e^{\sigma t}) + (e^{-\sigma t}) \\
(-\bar{\lambda}+\sigma)(-\bar{\lambda}-\sigma)e^{\sigma t}
+ e^{-\sigma t}(\bar{\lambda}-\sigma)(-\bar{\lambda}-\sigma) &
(-\bar{\lambda}+\sigma)(-e^{\sigma t})
+ (e^{-\sigma t})(-\bar{\lambda}-\sigma)
\end{pmatrix}
\begin{pmatrix}
x_o \\
\dot{x_o}
\end{pmatrix}
\frac{e^{-\bar{\lambda}t}}
{ -2\sigma } \\
&=
\begin{pmatrix}
\frac{ e^{\sigma t} + e^{-\sigma t} }
{ 2 }
+ \bar{\lambda}
\frac{ e^{\sigma t} - e^{-\sigma t} }
{ 2\sigma }
&
\frac{ e^{\sigma t} - e^{-\sigma t} }
{ 2\sigma }
\\
-(\bar{\lambda}^2-\sigma^2)
\frac{ e^{\sigma t} - e^{-\sigma t} }
{ 2\sigma }
&
\frac{ e^{\sigma t} + e^{-\sigma t} }
{ 2 }
- \bar{\lambda}
\frac{ e^{\sigma t} - e^{-\sigma t} }
{ 2\sigma }
\end{pmatrix}
\begin{pmatrix}
x_o \\
\dot{x_o}
\end{pmatrix}
e^{-\bar{\lambda}t}
\end{aligned}
$$
One can see that $e^{-\bar{\lambda}t}$ contributes uniform decay and anchors the characteristic timescale of the system.
Now define the odd and even exponential factors as
T 1 ( x , k ) = lim s → k e s x − e − s x 2 s T 0 ( x , k ) = lim s → k e s x + e − s x 2 = ∂ ∂ x T 1 ( x , k ) \begin{aligned}
T_1(x,k) & = \lim_{ s \rightarrow k }{ \frac{e^{sx}-e^{-sx}}{2s} } \\
T_0(x,k) & = \lim_{ s \rightarrow k }{ \frac{e^{sx}+e^{-sx}}{2 } }
= \frac{\partial}{\partial x} T_1(x,k)
\end{aligned} T 1 ( x , k ) T 0 ( x , k ) = s → k lim 2 s e s x − e − s x = s → k lim 2 e s x + e − s x = ∂ x ∂ T 1 ( x , k ) Thus, the most general form for all cases can be expressed as
( x x ˙ ) ( t ) = ( τ 0 + λ ˉ τ 1 τ 1 − ω 2 τ 1 τ 0 − λ ˉ τ 1 ) ( x o x o ˙ ) e − λ ˉ t = ( ( 1 0 0 1 ) τ 0 + ( λ ˉ 1 − ω 2 − λ ˉ ) τ 1 ) ( x o x o ˙ ) e − λ ˉ t where τ n = T n ( t , σ ) σ = λ ˉ 2 − ω 2 λ ˉ = c / ( 2 m ) ω 2 = k / m \begin{aligned}
\begin{pmatrix}
x \\
\dot{x}
\end{pmatrix}
(t)
=&
\begin{pmatrix}
\tau_0 + \bar{\lambda} \tau_1 & \tau_1 \\
-\omega^2 \tau_1 & \tau_0 - \bar{\lambda} \tau_1
\end{pmatrix}
\begin{pmatrix}
x_o \\
\dot{x_o}
\end{pmatrix}
e^{-\bar{\lambda}t} \\
=&
\left(
\begin{pmatrix}
1&0\\0&1
\end{pmatrix}
\tau_0
+
\begin{pmatrix}
\bar{\lambda}&1\\-\omega^2&-\bar{\lambda}
\end{pmatrix}\tau_1
\right)
\begin{pmatrix}
x_o \\
\dot{x_o}
\end{pmatrix}e^{-\bar{\lambda}t}\\
&\begin{aligned}
\text{where }
&\tau_n = T_n(t,\sigma)\\
&\sigma = \sqrt{\bar{\lambda}^2-\omega^2} \\
&\bar{\lambda} = c/(2m) \\
&\omega^2 = k/m
\end{aligned}
\end{aligned} ( x x ˙ ) ( t ) = = ( τ 0 + λ ˉ τ 1 − ω 2 τ 1 τ 1 τ 0 − λ ˉ τ 1 ) ( x o x o ˙ ) e − λ ˉ t ( ( 1 0 0 1 ) τ 0 + ( λ ˉ − ω 2 1 − λ ˉ ) τ 1 ) ( x o x o ˙ ) e − λ ˉ t where τ n = T n ( t , σ ) σ = λ ˉ 2 − ω 2 λ ˉ = c / ( 2 m ) ω 2 = k / m Sample Lua Implementation
function sprung_response ( t , d , v , k , c , m )
local decay = c / 2 / m
local omega = math.sqrt ( k / m )
local resid = decay * decay - omega * omega
local scale = math.sqrt ( math.abs ( resid ))
local T1 , T0 = t , 1
if resid < 0 then
T1 , T0 = math.sin ( scale * t ) / scale , math.cos ( scale * t )
elseif resid > 0 then
T1 , T0 = math . sinh ( scale * t ) / scale , math . cosh ( scale * t )
end
local dissipation = math.exp ( - decay * t )
local evolved_pos = dissipation * ( d * ( T0 + T1 * decay ) + v * ( T1 ) )
local evolved_vel = dissipation * ( d * ( - T1 * omega ^ 2 ) + v * ( T0 - T1 * decay ) )
return evolved_pos , evolved_vel
end
Click to Open Interactive Demo
Damping Cases
The $e^{-\bar{\lambda}t}$ factor is the global dissipator , and $\tau_1$, $\tau_0$ are the first- and zeroth-order time-integrators .
If all system parameters (m,c,k) are real, then $\sigma$ can be either real, imaginary, or zero, giving rise to the three well-known over/under/critically damped classifications.
Substituting each cases’ $\sigma$ into $T_n$ evaluates to the following:
let σ = α = i κ τ 1 = lim s → σ e s t − e − s t 2 s = { t σ = 0 sin κ t κ σ ∈ ℑ sinh α t α σ ∈ ℜ τ 0 = lim s → σ e s t + e − s t 2 = { 1 σ = 0 cos κ t σ ∈ ℑ cosh α t σ ∈ ℜ \begin{aligned}
&\text{let } \sigma=\alpha=i\kappa \\
\tau_1 =& \lim_{ s \rightarrow \sigma }{ \frac{e^{st}-e^{-st}}{2s} }
= \begin{cases}
t & \sigma =0 \\
\frac{\sin{\kappa t}}{\kappa} & \sigma \in \Im \\
\frac{\sinh{\alpha t}}{\alpha} & \sigma \in \Re
\end{cases} \\
\tau_0 =& \lim_{ s \rightarrow \sigma }{ \frac{e^{st}+e^{-st}}{2} }
=\begin{cases}
1 & \sigma =0 \\
\cos{\kappa t} & \sigma \in \Im \\
\cosh{\alpha t} & \sigma \in \Re
\end{cases}
\end{aligned} τ 1 = τ 0 = let σ = α = iκ s → σ lim 2 s e s t − e − s t = ⎩ ⎨ ⎧ t κ s i n κ t α s i n h α t σ = 0 σ ∈ ℑ σ ∈ ℜ s → σ lim 2 e s t + e − s t = ⎩ ⎨ ⎧ 1 cos κ t cosh α t σ = 0 σ ∈ ℑ σ ∈ ℜ Critically Damped $(\sigma=0)$:
τ 1 = t τ 0 = 1 \begin{aligned}
\tau_1 &= t \\
\tau_0 &= 1
\end{aligned} τ 1 τ 0 = t = 1 ( x x ˙ ) ( t ) = ( 1 + λ ˉ t t − λ ˉ 2 t 1 − λ ˉ t ) ( x o x o ˙ ) e − λ ˉ t \begin{pmatrix}
x \\
\dot{x}
\end{pmatrix}(t)
=\begin{pmatrix}
1+\bar{\lambda}t & t\\
-\bar{\lambda}^2 t & 1-\bar{\lambda}t
\end{pmatrix}
\begin{pmatrix}
x_o \\
\dot{x_o}
\end{pmatrix}
e^{-\bar{\lambda}t} ( x x ˙ ) ( t ) = ( 1 + λ ˉ t − λ ˉ 2 t t 1 − λ ˉ t ) ( x o x o ˙ ) e − λ ˉ t Special Case: Unimpeded $\mu = \omega = 0 \Rightarrow \bar{\lambda}=0$
( x x ˙ ) ( t ) = ( 1 t 0 1 ) ( x o x o ˙ ) \begin{pmatrix}
x \\
\dot{x}
\end{pmatrix}(t)
=
\begin{pmatrix}
1 & t \\
0 & 1
\end{pmatrix}
\begin{pmatrix}
x_o \\
\dot{x_o}
\end{pmatrix} ( x x ˙ ) ( t ) = ( 1 0 t 1 ) ( x o x o ˙ ) Underdamped $(\sigma = i\kappa \textnormal{ , } \kappa \in \mathbb{R})$:
τ 1 = sin κ t κ τ 0 = cos κ t where κ = ω 2 − λ ˉ 2 \begin{aligned}
\tau_1 &= \frac{\sin{\kappa t}}
{ \kappa } \\
\tau_0 &= \cos{\kappa t} \\
\textnormal{where } \kappa &= \sqrt{\omega^2-\bar{\lambda}^2}
\end{aligned} τ 1 τ 0 where κ = κ sin κ t = cos κ t = ω 2 − λ ˉ 2 ( x x ˙ ) ( t ) = ( cos κ t + λ ˉ κ sin κ t 1 κ sin κ t − ω 2 κ sin κ t cos κ t − λ ˉ κ sin κ t ) ( x o x o ˙ ) e − λ ˉ t \begin{pmatrix}
x \\
\dot{x}
\end{pmatrix}(t)
=
\begin{pmatrix}
\cos{\kappa t}+\frac{\bar{\lambda}}{\kappa}\sin{\kappa t} & \frac{1}{\kappa} \sin{\kappa t}\\
-\frac{\omega^2}{\kappa}\sin{\kappa t} & \cos{\kappa t}-\frac{\bar{\lambda}}{\kappa}\sin{\kappa t}
\end{pmatrix}
\begin{pmatrix}
x_o \\
\dot{x_o}
\end{pmatrix}
e^{-\bar{\lambda}t} ( x x ˙ ) ( t ) = ( cos κ t + κ λ ˉ sin κ t − κ ω 2 sin κ t κ 1 sin κ t cos κ t − κ λ ˉ sin κ t ) ( x o x o ˙ ) e − λ ˉ t Special Case: Undamped $\mu = 0 \Rightarrow \bar{\lambda}=0, \kappa = \omega$
( x x ˙ ) ( t ) = ( cos ω t 1 ω sin ω t − ω sin ω t cos ω t ) ( x o x o ˙ ) \begin{pmatrix}
x \\
\dot{x}
\end{pmatrix}(t)
=
\begin{pmatrix}
\cos{\omega t} & \frac{1}{\omega}\sin{\omega t} \\
-\omega\sin{\omega t} & \cos{\omega t}
\end{pmatrix}
\begin{pmatrix}
x_o \\
\dot{x_o}
\end{pmatrix} ( x x ˙ ) ( t ) = ( cos ω t − ω sin ω t ω 1 sin ω t cos ω t ) ( x o x o ˙ ) Overdamped $(\sigma = \alpha \textnormal{ , } \alpha \in \mathbb{R})$:
τ 1 = sinh α t α τ 0 = cosh α t where α = λ ˉ 2 − ω 2 \begin{aligned}
\tau_1 &= \frac{\sinh{\alpha t}}{\alpha} \\
\tau_0 &= \cosh{\alpha t} \\
\textnormal{where } \alpha &= \sqrt{ \bar{\lambda}^2 - \omega^2 }
\end{aligned} τ 1 τ 0 where α = α sinh α t = cosh α t = λ ˉ 2 − ω 2 ( x x ˙ ) ( t ) = ( cosh α t + λ ˉ α sinh α t 1 α sinh α t − ω 2 α sinh α t cosh α t − λ ˉ α sinh α t ) ( x o x o ˙ ) e − λ ˉ t \begin{pmatrix}
x \\
\dot{x}
\end{pmatrix}(t)
=
\begin{pmatrix}
\cosh{\alpha t}+\frac{\bar{\lambda}}{\alpha}\sinh{\alpha t} & \frac{1}{\alpha} \sinh{\alpha t}\\
-\frac{\omega^2}{\alpha}\sinh{\alpha t} & \cosh{\alpha t}-\frac{\bar{\lambda}}{\alpha}\sinh{\alpha t}
\end{pmatrix}
\begin{pmatrix}
x_o \\
\dot{x_o}
\end{pmatrix}
e^{-\bar{\lambda}t} ( x x ˙ ) ( t ) = ( cosh α t + α λ ˉ sinh α t − α ω 2 sinh α t α 1 sinh α t cosh α t − α λ ˉ sinh α t ) ( x o x o ˙ ) e − λ ˉ t Special Case: Unsprung $\omega = 0 \Rightarrow \bar{\lambda} = \alpha = \frac{\mu}{2}$
( x x ˙ ) ( t ) = ( 1 1 − e − μ t μ 0 e − μ t ) ( x o x o ˙ ) \begin{pmatrix}
x \\
\dot{x}
\end{pmatrix}(t)
=
\begin{pmatrix}
1 & \frac{ 1 - e^{-\mu t} }
{ \mu } \\
0 & e^{-\mu t}
\end{pmatrix}
\begin{pmatrix}
x_o \\
\dot{x_o}
\end{pmatrix} ( x x ˙ ) ( t ) = ( 1 0 μ 1 − e − μ t e − μ t ) ( x o x o ˙ )