Damped Harmonic Oscillator

Alexei Gilchrist

2 Oscillators

The solutions to the harmonic oscillator with a velocity dependant drag force are derived. The relaxation time for the oscillator is defined.

❦

1 Introduction

We will introduce a particular type of damping that is proportional to the velocity (simply because it’s easier to deal with). The faster an object is moving the more drag it experiences and this frictional force opposes the objects velocity. That is, \(F=-b v\). This type of damping dominates when an object moves slowly through a viscous medium.

2 Equation of motion and solution

Including the damping, the total force on the object is

\[\begin{equation*}F = m \ddot{x} = -k x -b \dot{x}.\end{equation*}\]

With a little rearranging we get the equation of motion in a familiar form with just an additional term proportional to the velocity:

\[\begin{equation*}\ddot{x} + \gamma \dot{x} + \omega_0^2 x = 0,\end{equation*}\]

where \(\gamma = b/m\) is a constant that determines the amount of damping, and \(\omega_0^2=k/m\) is the angular frequency of the oscillator when there is no damping. The solution to this equation is

\[\begin{equation*}x(t) = e^{-\gamma/2 t}\left(C_1 e^{\sqrt{(\gamma/2)^2-\omega_0^2} t} + C_2 e^{-\sqrt{(\gamma/2)^2-\omega_0^2} t} \right).\end{equation*}\]

The constants \(C_1\) and \(C_2\) are determined by the initial conditions.

3 Solution regimes

The behaviour of the solution depends on what is happening under the square-root in the exponents, and naturally splits into three types depending on whether \((\gamma/2)^2-\omega_0^2\) is positive, negative, or zero. (Actually, the form of the general solution will change slightly for each of these regimes due to technical reasons).

3.1 Over-damped oscillator: \((\gamma/2)^2-\omega_0^2>0\)

Defining \(\Omega = \sqrt{(\gamma/2)^2-\omega_0^2}\) (which will be a real number), the solution becomes

\[\begin{equation*}x(t) = e^{-\gamma/2 t}\left(C_1 e^{\Omega t} + C_2 e^{-\Omega t} \right).\end{equation*}\]

Despite the overall exponential damping term \(e^{-\gamma/2 t}\), you shouldn’t just expect a decaying amplitude—it all depends on the initial conditions. For instance if the oscillator is given an initial position and a finite velocity in the same direction the amplitude will first rise before decaying away. Eventually these initial transients do die away and the position just decays away to zero.

In the plot below the total amplitude is the solid line, the dashed lines are the two terms in the equation above that add to give the total amplitude.

Here is the solution when begun at rest with some initial amplitude. The smaller the damping the faster the decay to zero. If this appears counter intuitive remember that the damping we have included is the sort of drag forces that a marble would experience falling through oil—the thicker the oil the slower it would move.

3.2 Critically-damped oscillator: \((\gamma/2)^2-\omega_0^2=0\)

In this case, naively substituting in the above expression will yield only a single free parameter. The general method of solving such equations (in the case where the complementary equation has repeated roots) will yield an extra term proportional to \(t\), so the general solution in this case is

\[\begin{equation*}x(t) = (C_1 e^{\Omega t} + C_2 t )e^{-\frac{\gamma}{2} t}.\end{equation*}\]

This equation now has two free parameters that are set by the initial position and velocity.

3.3 Under-damped oscillator: \((\gamma/2)^2-\omega_0^2<0\)

Now the exponent is pure imaginary, and the solution reduces to

\[\begin{equation*}x(t) = A e^{-\frac{\gamma}{2} t} \cos(\omega_d t+\phi)\end{equation*}\]

where \(\omega_d = \sqrt{\omega_0^2-(\gamma/2)^2}\). You can see that this is very similar to the undamped case \(x(t) = A \cos(\omega_0 t + \phi)\). The key differences are that the frequency is smaller (\(\omega_d\le \omega_0\)) and there is an overall exponential damping term.

Below is an example where the dashed lines are the undamped case, and the two exponential decays \(\pm e^{-\gamma/2\; t}\) that bracket the amplitude. If you look carefully, you will notice that the frequency of the damped oscillator is slightly smaller than the undamped case.

3.4 Solution in general

It would be nice if we had a single closed form general solution that was valid in all the parameter ranges and initial conditions. This is indeed possible. Assuming the oscillator has initial position \(x_0\) and velocity \(v_0\) then the motion is given by:

\[\begin{equation*}x(t) = e^{-\frac{\gamma}{2} t} \left[x_0 \cosh(i\omega_d t)-\frac{i(2v_0+\gamma x_0)}{2 \omega_d} \sinh(i \omega_d t) \right].\end{equation*}\]

It would be tempting to simplify the expression further using \(\sinh(i x)=i\sin(x)\) and \(\cosh(i x)=\cos(x)\) but remember that \(\omega_d = \sqrt{\omega_0^2-\gamma^2/4}\) and if \(\gamma>2\omega_0\) this will be pure imaginary. Normally, the argument to \(\sin\) and \(\cos\) is real so keeping the solution in the above form is less confusing.

Note: in the case of \(\omega_d = 0\) the second term requires that the limit \(\omega_d\rightarrow 0\) be taken in order to evaluate it (using L’H\^{o}pital’s rule).

Equipped with this solution we can create a simulation for all parameter values and initial conditions:

4 Relaxation time

For an under-damped oscillator we’ve found the position goes as

\[\begin{equation*}x(t) = A e^{-\frac{\gamma}{2} t} \cos(\omega_d t+\phi)\end{equation*}\]

We can think of this as decaying away with the exponential term \(e^{-\gamma/2\; t}\). Instead of using the damping rate \(\gamma\) though a much more accessible parameter is the relaxation time \(\tau\). This is defined as the time for the amplitude to decay to \(e^{-1}\approx 0.368\) of its original value. Clearly this will happen when \(t=2/\gamma\) and so,

\[\begin{equation*}\tau=2/\gamma.\end{equation*}\]

As we’ve defined it the relaxation time is not very precise. For example, in the diagram the oscillator attains the amplitude \(A e^{-1}\) much earlier! Instead, what we mean is either the smooth curve that connects the peaks of the motion, or the curve obtained by maximising the amplitude over all phases at each time.

The peaks of the motion occur when \(\cos(\omega_d t+\phi)=1\). That is, we can imagine a strobe light that goes off at the times \(t_n=n 2\pi/\omega_d - \phi/\omega_d\), and at these times the oscillator’s amplitude will be \(x_n = Ae^{-\gamma/2\; t_n}\). The relaxation time is then when this interpolated curve intersects \(Ae^{-1}\).

Alternatively if at each instance in time we maximise the amplitude by varying the phase to that \(\cos(\omega_d t+\phi)=1\), then the decay of the oscillator will be the smooth curve \(\tilde{x}(t) = Ae^{-\gamma/2\; t}\).