Damped Harmonic Oscillator

F_spring = -kx; and F_resistance = -bv .

Note: We use F_ar = bv instead of bv² 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 2^nd 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

md²x/dt² + 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 x₁(t) is a solution, then C₁x₁(t) is also a solution (where C₁ is some constant). Since we need two constants, we would expect a full solution to include the addition of two solutions: x(t) = C₁x₁(t) + C₂x₂(t) .

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

mp²Ae^pt + bpAe^pt + kAe^pt = 0.

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

mp² + bp + k = 0.

This is the condition that must be satisfied if Ae^pt is a solution.

Solving for p gives:

p = [-b +/- (b²-4mk)^1/2] / 2m = -b/2m +/- [(b/2m)²-(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 w_o = (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)²-(k/m)]^1/2 = [g²- w_o²]^1/2 = d , so that p = -g +/- d .

Note that there are two values of p (unless b²=4mk). This would indicate that there are two solutions of the form Ae^pt, 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)² > k/m

In this case the quantity

[(b/2m)²-(k/m)]^1/2 = [g²- w_o²]^1/2 = d

is real, and so the solution becomes:

x(t) = C₁e^(-g+d)t + C₂e^(-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 x_o from the equilibrium position. Do this for the OVERDAMPED case. This problem amounts to finding C₁ and C₂ in terms of the initial conditions: x(t=0) = x_o and v(t=0) = v_o .

b) (b/2m)² = 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)²-(k/m)]^1/2 = [g²- w_o²]^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) = C₁ t e^pt in addition to the normal x(t) = C₂ e^pt :

x(t) = C₁ t e^pt + C₂ e^pt ,

dx/dt = C₁ e^pt + C₁ t p e^pt + C₂ p e^pt , and

d²x/dt² = C₁ p e^pt + C₁ p e^pt + C₁ t p² e^pt + C₂ p² e^pt .

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

2mC₁p e^pt + mC₁tp² e^pt + mC₂p² e^pt + bC₁ e^pt + bC₁tp e^pt + bC₂p e^pt + kC₁t e^pt +kC₂ e^pt = 0 .

We can cancel out the e^pt in every term to get:

2mC₁p + mC₁tp² + mC₂p² + bC₁ + bC₁tp + bC₂p + kC₁t + kC₂ = 0 .

For this to be true for all time, all terms with t¹ must cancel and all terms with t⁰ must also cancel:

t¹ terms: mC₁p² + bC₁p + kC₁ = 0, or mp² + bp + k = 0 (but this is our original condition);

t⁰ terms: 2mC₁p + mC₂p² + bC₁ + bC₂p + kC₂ = 0 .

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

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

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

x(t) = C₁ t e^-bt/2m + C₂ 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)² = 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)² < k/m

In this case the quantity [(b/2m)²-(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)²]^1/2 = [w_o² - g²]^1/2 which is a real number

so that the solution in this case is:

x(t) = C₁ e^-gt e^+iwt + C₂ e^-gt e^-iwt = e^-gt [C₁ e^+iwt + C₂ e^-iwt].

To see what this is, we recall that: e^iq = 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 [C₁ e^+iwt + C₂ e^-iwt] must also be a real quantity. But C₁ and C₂ may be complex quantities as well! Let's call C₁ = C_1R + iC_1i where both C_1R and C_1i are real, and C₂ = C_2R + iC_2i where both C_2R and C_2i are real.

Working this out gives:

[C₁ e^+iwt + C₂ e^-iwt] = (C_1R+iC_1i)*[cos(wt) + i sin(wt)] + (C_2R+iC_2i)*[cos(wt) - i sin(wt)] = x_R + ix_i .

x_R = C_1R cos(wt) - C_1i sin(wt) + C_2R cos(wt) + C_2i sin(wt) = real number, and

x_i = C_1i cos(wt) + C_1R sin(wt) + C_2i cos(wt) - C_2R sin(wt) = 0 . 

For the x_i = 0 equation to be true, we need: C_1i = -C_2i and C_1R = C_2R . [Note that these conditions require that C₂=C₁^* where C₁^* is the complex conjugate of C₁.]

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

x(t) = e^-gt [C₁ e^+iwt + C₂ 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 (x_o and v_o).

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) = re^iq = r cos(q) + i r sin(q). This implies that

a = r cos(q) and b = r sin(q), or inversely, r = [a²+b²]^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)²]^1/2 = [w_o² - g²]^1/2 , and that this w is less than the w_o which is the natural frequency of oscillation without damping [w_o = (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 = [w_o² - g²]^1/2 = [(k/m) - (b²/4m²)]^1/2 ?

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

d²x/dt² = +g² A e^-gt cos(wt+q) + g w A e^-gt sin(wt+q) + g w A e^-gt sin(wt+q) - w² A e^-gt cos(wt+q)

= (g² - w²) A e^-gt cos(wt+q) + 2g w A e^-gt sin(wt+q)

Substituting these into Newton's Second Law: md²x/dt² + b dx/dt + kx = 0 gives:

m(g² - w²) 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(g² - w²) + -bg +k = 0, and

2mgw -bw = 0.

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

The first equation gives

w² = (k/m) + g² - (b/m)g = (k/m) + g² - 2g² = (k/m) - g ² = w_o² - g² .

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)² > k/m and x(t) = C₁e^(-g+d)t + C₂e^(-g-d)t .

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

d = [(b/2m)²-(k/m)]^1/2 = (b/2m)[1-e]^1/2 , where we define e = {(k/m)/(b/2m)}² so that (k/m)^1/2 = e^1/2 * (b/2m) or (b/2m) = (k/m)^1/2 / e^1/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)²}^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)² = k/m and x(t) = C₁ t e^-bt/2m + C₂ 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)² < 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.

Return to PHYS 380