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}