Damped Harmonic Oscillator

Fspring = -kx;   and   Fresistance = -bv  .

Note:  We use Far = bv instead of bv2 for two reasons:
1.  usually, because the motion will go back and forth, the speeds don’t reach very high values, and so F = -bv is more appropriate; and

2.  if we use F = -bv, Newton’s 2nd law gives a linear differential equation which will probably be easier than a non-linear one.

If there are no other forces, then Newton's Second Law becomes:

-kx -bv  =  ma, or

md2x/dt2 + b dx/dt + kx  =  0 .

Note that this is a linear, homogenous, second order differential equation.  Because it is second order, it involves x, v, and t and can’t be separated like we were able to do before.  We need to use a different fundamental technique.

Because it is second order, there should be two constants of integration (usually initial conditions on x and v).

Because it is linear, if x1(t) is a solution, then C1x1(t) is also a solution (where C1 is some constant). Since we need two constants, we would expect a full solution to include the addition of two solutions: x(t)  =  C1x1(t) + C2x2(t) .

Since the equation is linear (that is, we only have x1 in any term), and since we know that d(Aeat)/dt = aAeat (that is, the differential of the exponential function gives itself back again), we might guess that x(t) = Aept.
To see if and when this works, we simply substitute it back into our differential equations to get:

mp2Aept + bpAept + kAept  =  0.

We can cancel out the Aept in each term to get an algebraic equation for p:

mp2 + bp + k  =  0.

This is the condition that must be satisfied if Aept is a solution.

Solving for p gives:

p  =  [-b +/- (b2-4mk)1/2] / 2m  =  -b/2m +/- [(b/2m)2-(k/m)]1/2 .  Note that we usually get two solutions since this is a quadratic equation.

To save writing time, let's define g = (b/2m) and let's use wo = (k/m)1/2 since (k/m)1/2 is the frequency for an undamped harmonic oscillator.

For ease in what follows, let’s further define d to be:   [(b/2m)2-(k/m)]1/2  =  [g2- wo2]1/2  =  d , so that p  =  -g +/- d .

Note that there are two values of p (unless b2=4mk).  This would indicate that there are two solutions of the form Aept, one for each of the two values of p.

Also note that g is a real number, but d could be either real or imaginary.

There are three distinct cases:

a) (b/2m)2 > k/m

In this case the quantity

[(b/2m)2-(k/m)]1/2  =  [g2- wo2]1/2  =  d

is real, and so the solution becomes:

x(t)  =  C1e(-g+d)t + C2e(-g-d)t .

Note that g > d so this is obviously a combination of two dying exponentials.  This is called the overdamped case.

Homework Problem: Problem #11:   Find x(t) and v(t) for a particle of mass, m, subject to a spring force of –kx and a damping force of –bv if it is displaced an initial distance of xo from the equilibrium position.  Do this for the OVERDAMPED case.  This problem amounts to finding C1 and C2 in terms of the initial conditions: x(t=0) = xo and v(t=0) = vo .

b) (b/2m)2 = k/m

This is the special case between case a) above: the dying exponential, and case c) below: which we will soon see to be the dying but oscillating while it dies case.

In this case, since p  =  -g +/- d and d  =  [(b/2m)2-(k/m)]1/2  =  [g2- wo2]1/2 = 0, the two p’s and hence the two exponentials become the same. But we must have two, rather than just one, functions.

Let's try the function x(t) = C1 t ept in addition to the normal x(t) = C2 ept :

x(t)  =  C1 t ept + C2 ept ,

dx/dt  =  C1 ept + C1 t p ept + C2 p ept ,  and

d2x/dt2  =  C1 p ept + C1 p ept + C1 t p2 ept + C2 p2 ept .

Plugging them into the differential equation (Newton's Second Law:  md2x/dt2 + b dx/dt + kx  =  0) gives:

2mC1p ept + mC1tp2 ept + mC2p2 ept + bC1 ept + bC1tp ept + bC2p ept + kC1t ept +kC2 ept  =  0 .

We can cancel out the ept in every term to get:

2mC1p + mC1tp2 + mC2p2 + bC1 + bC1tp + bC2p + kC1t + kC2  =  0 .

For this to be true for all time, all terms with t1 must cancel and all terms with t0 must also cancel:

t1 terms: mC1p2 + bC1p + kC1  =  0, or mp2 + bp + k  =  0 (but this is our original condition);

t0 terms: 2mC1p + mC2p2 + bC1 + bC2p + kC2  =  0 .

This last equation says that: C1(2mp + b) + C2(mp2 + bp + k)  =  0. Note that the coefficient of the C2 term is zero; this means that the coefficient of the C1 term must also be zero:

2mp + b  =  0, or p = -b/2m.

But this is true for the C2 coefficient only if (b/2m)2 = k/m, which was our condition for this case b!  Hence our solution for this case is:

x(t)  =  C1 t e-bt/2m + C2 e-bt/2m .

This solution also dies out, since the exponential will die out quicker than the t will go to infinity. At first it looks like there is no k in the solution, but when you realize that the condition for this solution is: (b/2m)2 = k/m, you see that by knowing b and m, you also determine k!  This is called the critically damped case, since x will approach zero quicker than either case a above or c below.  This is detailed in the section on comparisons below.

c) (b/2m)2 < k/m

In this case the quantity [(b/2m)2-(k/m)]1/2  =  d is imaginary, so let's take out the negative sign from the square root and so define d = iw where

w  =  [(k/m)-(b/2m)2]1/2  =  [wo2 - g2]1/2 which is a real number

so that the solution in this case is:

x(t)  =  C1 e-gt e+iwt + C2 e-gt e-iwt  =  e-gt [C1 e+iwt + C2 e-iwt].

To see what this is, we recall that: eiq  =  cos(q) + i sin(q) and e-iq  =  cos(q) - i sin(q) .

[You can show that this is true by expanding each term in a power series, and seeing that you get the same series on both sides of the equality sign.]

Thus the above solution indicates some type of exponentially dying oscillation because of the factor e-gt.

But how do we work with imaginary solutions to real physical problems? We start by saying that x(t) is a real quantity, so [C1 e+iwt + C2 e-iwt] must also be a real quantity. But C1 and C2 may be complex quantities as well!  Let's call C1  =  C1R + iC1i where both C1R and C1i are real, and C2  =  C2R + iC2i where both C2R and C2i are real.

Working this out gives:

[C1 e+iwt + C2 e-iwt]  =  (C1R+iC1i)*[cos(wt) + i sin(wt)] + (C2R+iC2i)*[cos(wt) - i sin(wt)]  =  xR + ixi .

xR  =  C1R cos(wt) - C1i sin(wt) + C2R cos(wt) + C2i sin(wt)  =  real number, and

xi  =  C1i cos(wt) + C1R sin(wt) + C2i cos(wt) - C2R sin(wt)  =  0

For the xi = 0 equation to be true, we need: C1i = -C2i and C1R = C2R . [Note that these conditions require that C2=C1* where C1* is the complex conjugate of C1.]

Note that we are back to having two (real) constants for our second order differential equation: C1R and C1i. To simplify the notation, let's rename them C1R = a, and C1i = b. Our solution is now:

x(t)  =  e-gt [C1 e+iwt + C2 e-iwt]  =  e-gt [(a+ib) e+iwt + (a-ib) e-iwt] , where a and b are real constants normally determined from the initial conditions (xo and vo).

We can make our answer look even "nicer" if we recognize that a complex number (a+ib) can be expressed in "polar" form: (a+ib) = reiq = r cos(q) + i r sin(q). This implies that

a = r cos(q) and b = r sin(q), or inversely, r = [a2+b2]1/2 and q = tan-1(b/a) . Let's now apply this to our solution for x(t):

xt)  =  e-gt [(a+ib) e+iwt + (a-ib) e-iwt]   =  e-gt [r e+iq e+iwt + r e-iq e-iwt]  =

r e-gt [ e+i(wt+q) + e-i(wt+q) ]   =

r e-gt [cos(wt+q) + i sin(wt+q) + cos(wt+q) - i sin(wt+q)]   =

2 r e-gt cos(wt+q)  =  A e-gt cos(wt+q),

where (A=2r) and q are the two constants that are normally determined from initial conditions. It is easy to see from this last expression that the motion is an oscillation that has an amplitude that dies off exponentially.

Note that the frequency of this oscillation is  w where  w  =  [(k/m)-(b/2m)2]1/2  =  [wo2 - g2]1/2  , and that this w is less than the wo which is the natural frequency of oscillation without damping [wo = (k/m)1/2].

If we started from a physically inspired guess of x(t)  =  A e-gt cos(wt+q ) , can we show that this is indeed a solution of the differential equation with the conditions that g = b/2m and that w  =  [wo2 - g2]1/2  =  [(k/m) - (b2/4m2)]1/2 ?

dx/dt  =  -g A e-gt cos(wt+q) + -w A e-gt sin(wt+q)

d2x/dt2  =  +g2 A e-gt cos(wt+q) + g w A e-gt sin(wt+q) + g w A e-gt sin(wt+q) - w2 A e-gt cos(wt+q)

=  (g2 - w2) A e-gt cos(wt+q) + 2g w A e-gt sin(wt+q)

Substituting these into Newton's Second Law: md2x/dt2 + b dx/dt + kx  =  0  gives:

m(g2 - w2) A e-gt cos(wt+q) + 2mgw A e-gt sin(wt+q)

-bg A e-gt cos(wt+q) + -bw A e-gt sin(wt+q)

+k A e-gt cos(wt+q)  =  0 .

There is an A e-gt in every term that can cancel out. But we have two different functions of time in this equation: cos(wt+q) and sin(wt+q) . For the equation to be satisfied for all times, then the coefficient of each function must be zero. This leads to two equations:

m(g2 - w2) + -bg +k  =  0,  and

2mgw -bw  =  0.

The second equation is easy and gives: g = b/2m.

The first equation gives

w2  =  (k/m) + g2 - (b/m)g  =  (k/m) + g2 - 2g2  =  (k/m) - g 2  =  wo2 - g2 .

This is exactly what we need!  This is called the underdamped case.

Comparisons:

All three cases result in solutions that have a dying exponential, so all three do die out (due to losing energy to air resistance).  Which dies out the quickest?  The longest lived term is the one with the smallest coefficient of t in the exponent.  The quickest to die out has the biggest coefficient of t in the exponential.  In all cases below, remember that  g = (b/2m).

For case a)  (b/2m)2 > k/m   and     x(t)  =  C1e(-g+d)t + C2e(-g-d)t  .

The smallest coefficient of t in the exponent in this case is g-d where

d  =  [(b/2m)2-(k/m)]1/2  =  (b/2m)[1-e]1/2  , where we define  e  =  {(k/m)/(b/2m)}2  so that (k/m)1/2  =  e1/2 * (b/2m)   or   (b/2m)  =  (k/m)1/2 / e1/2 ;

note that  e < 1 ,   so  g-d   =  (b/2m){1 - [1-e]1/2 }

a)  for small e,   [1-e]1/2    1 - (½)e   so   g-d       (b/2m)*(e/2)  =  (1/2)*(k/m)/(b/2m)    =  (½)*{(k/m)/(b/2m)2}1/2 *(k/m)1/2  =  (½)*e * (k/m)1/2  < (k/m)1/2.

b)  for e = .99,  {1 - [1-.99]1/2 }  =  .9,  so  g-d  =  .9 * (b/2m)   =  .9 * (k/m)1/2 /.995  <  (k/m)1/2

For case b)  (b/2m)2 = k/m   and     x(t)   =  C1 t e-bt/2m +  C2 e-bt/2m .

The coefficients of t in the exponent in this case are both equal to (b/2m) = (k/m)1/2.

For case c)  (b/2m)2 < k/m   and     x(t)   =  A e-gt cos(wt+q)  .

The coefficient of t in the exponent in this case is (b/2m), but (b/2m) <  (k/m)1/2.

From this, we see that case b) has the largest coefficient of t in the exponential and so dies out the quickest.

This completes the homogeneous case. The next thing to consider is what happens when we have an applied force on this damped harmonic oscillator.