Damped Harmonic Oscillator

Nov 21, 2022

Or: Implementing a damped spring model in 16 lines of code.

Governing Equation:

$m\ddot{x} + c\dot{x} + kx = f(t)$

Zero-input response f(t) = 0

Solve

$\begin{aligned} 0 & = m\ddot{x} + c\dot{x} + kx \\ x_o & = x(t_o) \\ \dot{x}_o & = \dot{x}(t_o) \end{aligned}$

Exponential ansatz:

$\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}$

Plugging in gives

$\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}$ $\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}$

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.

$\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}$

Let t_o=0 , then all the exponentials are 1, so it simplifies to

$\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}$

Altogether it becomes

$\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 }$

Which can be rearranged to

$\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}$

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

$\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}$

Thus, the most general form for all cases can be expressed as

$\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}$

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

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:

$\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}$

Critically Damped (\sigma = 0):

$\begin{aligned} \tau_1 &= t \\ \tau_0 &= 1 \end{aligned}$ $\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}$

Special Case: Unimpeded \mu = \omega = 0 \Rightarrow \bar{\lambda} = 0

$\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}$

Underdamped (\sigma = i\kappa \textnormal{ , } \kappa \in \mathbb{R}):

$\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}$ $\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}$

Special Case: Undamped \mu = 0 \Rightarrow \bar{\lambda} = 0, \kappa = \omega

$\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}$

Overdamped (\sigma = \alpha \textnormal{ , } \alpha \in \mathbb{R}):

$\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}$ $\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}$

Special Case: Unsprung \omega = 0 \Rightarrow \bar{\lambda} = \alpha = \frac{\mu}{2}

$\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}$