Analytical solution

Let’s concentrate on Eq. (1) , and divide the equation by \(m\). The simplified differential equation for forced harmonic oscillator with damping reads: \[\begin{eqnarray} \ddot{x}+ 2\zeta \omega_0 \dot{x}+ \omega_0^2 x&=& \frac{ f(t)}{m} ,\, x(0)=x_0, \,\text{and}\,\dot{x}(0)=\dot{x}_0,\tag{3} \end{eqnarray}\] where \(\dot x\equiv\frac{dx}{dt}\), \(\omega_0\equiv\sqrt{\frac{k}{m}}\) is the natural frequency of the oscillation, and \(\zeta\equiv\frac{c}{2\sqrt{mk}}\) is the damping ratio. We also included the initial conditions. We are dealing with an in-homogeneous linear differential equation with constant coefficients. One of the best tools to solve such equations is the Laplace transformation:

\[\begin{eqnarray} X(s)=\mathcal{L}\big[x(t)\big]=\int_0^\infty dt \, e^{-s \,t} x(t) \tag{4}. \end{eqnarray}\]

The nice feature of the Laplace transformation is that it converts differential equations to algebraic equations. It follows from the transformation property of the derivatives:

\[\begin{eqnarray} \mathcal{L}\big[\dot{x}(t)\big]&=&\int_0^\infty dt \, e^{-s \,t} \frac{dx}{dt}=\int_0^\infty dt \,\frac{d}{dt}\big(e^{-s \,t} x\big)-\int_0^\infty dt \big(\frac{d}{dt}e^{-s \,t}\big) x\\\nonumber &=& \big(e^{-s \,t} x\big)\bigg\rvert_0^\infty+s\int_0^\infty dt e^{-s \,t} x= s X(s)- x_0 \tag{5}. \end{eqnarray}\]

Similarly the second order derivative transforms as \[\begin{eqnarray} \mathcal{L}\big[\ddot{x}(t)\big]&=& s \mathcal{L}\big[\dot{x}(t)\big]- \dot{x}_0 =s^2 X(s) -s x_0 -\dot{x}_0 \tag{6}. \end{eqnarray}\]

Laplace transforming Eq. (3) we get \[\begin{eqnarray} s^2 X -s x_0 -\dot{x}_0+ 2\zeta \omega_0(s X- x_0)+\omega_0^2 X=\frac{ 1}{m} F(s) \tag{7}. \end{eqnarray}\] Solving Eq. (7) for \(X\), we get

\[\begin{eqnarray} X&=&\frac{s x_0+ 2\zeta \omega_0x_0+\dot{x}_0}{s^2 +2 \zeta \omega_0 s+\omega_0^2 }+\frac{1}{m}\frac{F(s)}{s^2 +2 \zeta \omega_0 s+\omega_0^2 }\nonumber\\ &=&\frac{(s+ \zeta \omega_0)x_0+ \zeta \omega_0 x_0+\dot{x}_0}{(s + \zeta \omega_0)^2 +\omega_0^2(1 -\zeta^2)}+\frac{1}{m}\frac{F(s)}{(s + \zeta \omega_0)^2 +\omega_0^2(1 -\zeta^2)} \tag{8}. \end{eqnarray}\] In order to evaluate the inverse Laplace transform, we need to know the functional form of the driving force. Let’s assume that \(f(t)\) is of the following form:

\[\begin{eqnarray} f(t)&=&f_0\sin(\omega t). \tag{9} \end{eqnarray}\] Its Laplace transform is given by: \[\begin{eqnarray} F(s)\equiv \mathcal{L}\big[f(t)\big]=\frac{f_0 \omega}{s^2+\omega^2}. \tag{10} \end{eqnarray}\] We will have to do some partial fraction expansion: \[\begin{eqnarray} \frac{1}{\left((s + \zeta \omega_0)^2 +\omega_0^2(1 -\zeta^2)\right)(s^2+\omega^2)}=\frac{A(s+\zeta \omega_0)+B}{(s + \zeta \omega_0)^2 +\omega_0^2(1 -\zeta^2)}+\frac{Cs+D}{s^2+\omega^2}, \tag{11} \end{eqnarray}\] which will be easy to convert back to time domain since they will correspond to sines and cosines with exponential functions in front. We now need to figure out \(A\), \(B\), \(C\) and \(D\). If we were to equate the denominators and sum up the resulting numerators, we will see that, in order to set the coefficient of the \(s^3\) term in the numerator to zero we will need \(A=-C.\) To relate \(C\) and \(D\) we can multiply (11) by \(s-i\omega\) and then set \(s=i\omega.\) This will remove the first term on the right hand side and yield:

\[\begin{eqnarray} i\omega C+D=\frac{1}{(i\omega + \zeta \omega_0)^2 +\omega_0^2(1 -\zeta^2)} \tag{12} \end{eqnarray}\] This is a complex equation, and splitting it into the real and imaginary part, we get: \[\begin{eqnarray} D(\omega_0^2-\omega^2)-2C \omega^2 \omega_0\zeta&=&1\nonumber\\ C\omega(\omega_0^2-\omega^2)+2D \omega \omega_0\zeta&=&0. \tag{13} \end{eqnarray}\] Inverting it, we get: \[\begin{eqnarray} C&=&\frac{-2 \omega_0\zeta}{(\omega_0^2-\omega^2)^2+4\omega^2\omega_0^2\zeta^2}\nonumber\\ D&=&\frac{\omega_0^2-\omega^2}{(\omega_0^2-\omega^2)^2+4\omega^2\omega_0^2\zeta^2} \tag{14} \end{eqnarray}\] Finally, setting \(s=-\zeta \omega_0\), and going trough some algebra we get. \[\begin{eqnarray} B&=&\frac{\omega^2 -\omega_0^2+2\omega_0^2\zeta^2}{(\omega_0^2-\omega^2)^2+4\omega^2\omega_0^2\zeta^2}. \tag{15} \end{eqnarray}\]

We can now inverse transform Eq. (8) using elementary properties of the transformation:

\[\begin{eqnarray} x(t)&=&\mathcal{L}^{-1}\big[X(s)\big]. \tag{16} \end{eqnarray}\]

Inverse Laplace transformation yields. \[\begin{eqnarray} x(t)&=&\left[x_0 \cos(\omega_0\sqrt{1 -\zeta^2} \,t)+\frac{ \zeta \omega_0 x_0+\dot{x}_0}{\omega_0\sqrt{1 -\zeta^2}}\sin(\omega_0\sqrt{1 -\zeta^2} \,t)\right]e^{-\zeta \omega_0 t}\nonumber\\ &&+\frac{2f_0 \omega}{m}\frac{ \omega_0\zeta}{(\omega_0^2-\omega^2)^2+4\omega^2\omega_0^2\zeta^2}e^{-\zeta \omega_0 t}\cos(\omega_0\sqrt{1 -\zeta^2} \,t)\nonumber\\ &&+\frac{f_0 \omega}{m}\frac{1}{\omega_0\sqrt{1 -\zeta^2}}\frac{\omega^2 -\omega_0^2+2\omega_0^2\zeta^2}{(\omega_0^2-\omega^2)^2+4\omega^2\omega_0^2\zeta^2}e^{-\zeta \omega_0 t}\sin(\omega_0\sqrt{1 -\zeta^2} \,t)\nonumber\\ &&-\frac{2f_0 \omega}{m}\frac{ \omega_0\zeta}{(\omega_0^2-\omega^2)^2+4\omega^2\omega_0^2\zeta^2}\cos(\omega \,t)+\frac{f_0 }{m}\frac{\omega^2 -\omega_0^2}{(\omega_0^2-\omega^2)^2+4\omega^2\omega_0^2\zeta^2}\sin(\omega \,t). \tag{17} \end{eqnarray}\] We can do one last touch and combine the last two terms into as single function with a phase shift.

The full solution with damping with \(f(t)=f_0\sin(\omega\, t),\quad\) can be written as:

\[\begin{eqnarray} x(t)&=&\left[x_0 \cos(\omega_0\sqrt{1 -\zeta^2} \,t)+\frac{ \zeta \omega_0 x_0+\dot{x}_0}{\omega_0\sqrt{1 -\zeta^2}}\sin(\omega_0\sqrt{1 -\zeta^2} \,t)\right]e^{-\zeta \omega_0 t}\nonumber\\ &+&\frac{f_0 \omega e^{-\zeta \omega_0 t}}{m[(\omega_0^2-\omega^2)^2+4\omega^2\omega_0^2\zeta^2]}\left[2\omega_0\zeta\cos(\omega_0\sqrt{1 -\zeta^2} \,t)+\frac{\omega^2 -\omega_0^2+2\omega_0^2\zeta^2}{\omega_0\sqrt{1 -\zeta^2}}\sin(\omega_0\sqrt{1 -\zeta^2} \,t)\right]\nonumber\\ &+&\frac{ g}{m\sqrt{(\omega_0^2-\omega^2)^2+4\omega^2\omega_0^2\zeta^2}}\sin(\omega \,t-\delta) \tag{18} \end{eqnarray}\] where \(\delta\equiv \arctan\left[\frac{2\omega\,\omega_0\zeta}{\omega_0^2-\omega^2}\right]\quad\), the first line is related to the initial conditions, the second and third lines are the transient response, and finally the last line is the steady state solution.

Resonances at zero damping: The final solution runs into problems when we consider \(\zeta=0\) and \(\omega=\omega_0\): the coefficient of the steady state solution diverges. This is because of the assumptions we made when we were inverting \(X(s)\). At \(\zeta=0\) and \(\omega=\omega_0\), two poles will merge and create a second order pole. Let’s take a closer look:

\[\begin{eqnarray} \lim_{\zeta\rightarrow 0,\,\omega\rightarrow\omega_0}\frac{1}{\left((s + \zeta \omega_0)^2 +\omega_0^2(1 -\zeta^2)\right)(s^2+\omega_0^2)}=\frac{1}{(s^2+\omega_0^2)^2}. \tag{19} \end{eqnarray}\] We can figure out how to inverse transform it by exploiting few features of the Laplace transforms as follows:

\[\begin{eqnarray} \mathcal{L}^{-1}\left[\frac{1}{(s^2+\omega_0^2)^2}\right]&=&\mathcal{L}^{-1}\left[-\frac{1}{2s}\frac{d}{ds}\left(\frac{1}{s^2+\omega_0^2}\right)\right]=-\frac{1}{2}\int_0^td\tau \tau \frac{\sin (\omega_0\tau)}{\omega_0}\nonumber\\ &=&-\frac{1}{2 \omega_0}\frac{d}{d\omega_0}\int_0^td\tau \cos (\omega_0\tau)=-\frac{1}{\omega_0}\frac{d}{d\omega_0}\left[\frac{\sin(\omega_0 t)}{\omega_0}\right]=\frac{\sin(\omega_0 t)-\omega_0 t\cos(\omega_0 t)}{2\omega_0^3} \tag{20} \end{eqnarray}\]

The full solution at the resonance frequency(\(\omega=\omega_0\,\) ) with no damping (\(\zeta=0\)) is:

\[\begin{eqnarray} x(t)&=&\left[x_0 \cos(\omega_0\,t)+\frac{\dot{x}_0}{\omega_0}\sin(\omega_0 \,t)\right]+\left[\frac{f_0}{2 m\omega_0^2}\left( \sin(\omega_0 t)-\omega_0 t\cos(\omega_0 t)\right)\right].\tag{21} \end{eqnarray}\] This shows that the amplitute will grow with time. In reality the model will break at some point since the ampliture of oscillations cannot grow indefinitely. (For example, the spring will literaly break if it is stretched too far.)