Main Content

The Physics of the Damped Harmonic Oscillator

This example explores the physics of the damped harmonic oscillator by solving the equations of motion in the case of no driving forces. This example investigates the cases of under-, over-, and critical-damping.

Contents

  1. Derive Equation of Motion

  2. Solve the Equation of Motion (F = 0)

  3. Underdamped Case (ζ<1)

  4. Overdamped Case (ζ>1)

  5. Critically Damped Case (ζ=1)

  6. Conclusion

1. Derive Equation of Motion

Consider a forced harmonic oscillator with damping shown below. Model the resistance force as proportional to the speed with which the oscillator moves.

Define the equation of motion where

  • m is the mass

  • c is the damping coefficient

  • k is the spring constant

  • F is a driving force

syms x(t) m c k F(t)
eq = m*diff(x,t,t) + c*diff(x,t) + k*x == F
eq(t) = 

m2t2 x(t)+ct x(t)+kx(t)=F(t)

Rewrite the equation using c=mγ and k=mω02.

syms gamma omega_0
eq = subs(eq, [c k], [m*gamma, m*omega_0^2])
eq(t) = 

m2t2 x(t)+γmt x(t)+mx(t)ω02=F(t)

Divide out the mass m. Now we have the equation in a convenient form to analyze.

eq = collect(eq, m)/m
eq(t) = 

2t2 x(t)+γt x(t)+x(t)ω02=F(t)m

2. Solve the Equation of Motion where F = 0

Solve the equation of motion using dsolve in the case of no external forces where F=0. Use the initial conditions of unit displacement and zero velocity.

vel = diff(x,t);
cond = [x(0) == 1, vel(0) == 0];
eq = subs(eq,F,0);
sol = dsolve(eq, cond)
sol = 

e-tγ2-σ12γ+σ12σ1-e-tγ2+σ12γ-σ12σ1where  σ1=γ-2ω0γ+2ω0

Examine how to simplify the solution by expanding it.

sol = expand(sol)
sol = 

σ1σ22+σ1etγ2-4ω0222-γσ1σ22γ2-4ω02+γσ1etγ2-4ω0222γ2-4ω02where  σ1=e-γt2  σ2=e-tγ2-4ω022

Notice that each term has a factor of σ1, or e-γt/2, use collect to gather these terms

sol = collect(sol, exp(-gamma*t/2))
sol = 

e-σ12+eσ12-γe-σ12γ2-4ω02+γeσ12γ2-4ω02e-γt2where  σ1=tγ2-4ω022

The term γ2-4ω02 appears in various parts of the solution. Rewrite it in a simpler form by introducing the damping ratio ζγ2ω0.

Substituting ζ into the term above gives:

γ24ω02=2ω0(γ2ω0)21=2ω0ζ21,

syms zeta;
sol = subs(sol, ...
    sqrt(gamma^2 - 4*omega_0^2), ...
    2*omega_0*sqrt(zeta^2-1))
sol = 

e-γt2σ22+σ12+γσ24ω0ζ2-1-γσ14ω0ζ2-1where  σ1=e-ω0tζ2-1  σ2=eω0tζ2-1

Further simplify the solution by substituting γ in terms of ω0 and ζ,

sol = subs(sol, gamma, 2*zeta*omega_0)
sol = 

e-ω0tζσ22+σ12+ζσ22ζ2-1-ζσ12ζ2-1where  σ1=e-ω0tζ2-1  σ2=eω0tζ2-1

We have derived the general solution for the motion of the damped harmonic oscillator with no driving forces. Next, we'll explore three special cases of the damping ratio ζ where the motion takes on simpler forms. These cases are called

  • underdamped (ζ<1),

  • overdamped (ζ>1), and

  • critically damped (ζ=1).

3. Underdamped Case (ζ<1)

If ζ<1, then ζ2-1=i1-ζ2 is purely imaginary

solUnder = subs(sol, sqrt(zeta^2-1), 1i*sqrt(1-zeta^2))
solUnder = 

e-ω0tζσ12+σ22+ζσ1i21-ζ2-ζσ2i21-ζ2where  σ1=e-ω0t1-ζ2i  σ2=eω0t1-ζ2i

Notice the terms eiω0tζ2-1±e-iω0tζ2-1 in the above equation and recall the identity eix=cos(x)+isin(x).

Rewrite the solution in terms of cos.

solUnder = coeffs(solUnder, zeta);
solUnder = solUnder(1);
c = exp(-omega_0 * zeta * t);
solUnder = c * rewrite(solUnder / c, 'cos')
solUnder = e-ω0tζcos(ω0t1-ζ2)
solUnder(t, omega_0, zeta) = solUnder
solUnder(t, omega_0, zeta) = e-ω0tζcos(ω0t1-ζ2)

The system oscillates at a natural frequency of ω01-ζ2 and decays at an exponential rate of 1/ω0ζ.

Plot the solution with fplot as a function of ω0t and ζ.

z = [0 1/4 1/2 3/4];
w = 1;
T = 4*pi;
lineStyle = {'-','--',':k','-.'};

fplot(@(t)solUnder(t, w, z(1)), [0 T], lineStyle{1});

hold on;
for k = 2:numel(z)
    fplot(@(t)solUnder(t, w, z(k)), [0 T], lineStyle{k});
end
hold off;
grid on;
xticks(T*linspace(0,1,5));
xticklabels({'0','\pi','2\pi','3\pi','4\pi'});
xlabel('t / \omega_0');
ylabel('amplitude');
lgd = legend('0','1/4','1/2','3/4');
title(lgd,'\zeta');
title('Underdamped');

Figure contains an axes object. The axes object with title Underdamped, xlabel t / blank omega indexOf 0 baseline, ylabel amplitude contains 4 objects of type functionline. These objects represent 0, 1/4, 1/2, 3/4.

4. Overdamped Case (ζ>1)

If ζ>1, then ζ2-1 is purely real and the solution can be rewritten as

solOver = sol
solOver = 

e-ω0tζσ22+σ12+ζσ22ζ2-1-ζσ12ζ2-1where  σ1=e-ω0tζ2-1  σ2=eω0tζ2-1

solOver = coeffs(solOver, zeta);
solOver = solOver(1)
solOver = 

e-ω0tζeω0tζ2-12+e-ω0tζ2-12

Notice the terms (eω0tζ2-1+e-ω0tζ2-1)2 and recall the identity cosh(x)=ex+e-x2.

Rewrite the expression in terms of cosh.

c = exp(-omega_0*t*zeta);
solOver = c*rewrite(solOver / c, 'cosh')
solOver = cosh(ω0tζ2-1)e-ω0tζ
solOver(t, omega_0, zeta) = solOver
solOver(t, omega_0, zeta) = cosh(ω0tζ2-1)e-ω0tζ

Plot the solution to see that it decays without oscillating.

z = 1 + [1/4 1/2 3/4 1];
w = 1;
T = 4*pi;
lineStyle = {'-','--',':k','-.'};

fplot(@(t)solOver(t, w, z(1)), [0 T], lineStyle{1});

hold on;
for k = 2:numel(z)
    fplot(@(t)solOver(t, w, z(k)), [0 T], lineStyle{k});
end
hold off;
grid on;
xticks(T*linspace(0,1,5));
xticklabels({'0','\pi','2\pi','3\pi','4\pi'});
xlabel('\omega_0 t');
ylabel('amplitude');
lgd = legend('1+1/4','1+1/2','1+3/4','2');
title(lgd,'\zeta');
title('Overdamped');

Figure contains an axes object. The axes object with title Overdamped, xlabel omega indexOf 0 baseline blank t, ylabel amplitude contains 4 objects of type functionline. These objects represent 1+1/4, 1+1/2, 1+3/4, 2.

5. Critically Damped Case (ζ=1)

If ζ=1, then the solution simplifies to

solCritical(t, omega_0) = limit(sol, zeta, 1)
solCritical(t, omega_0) = e-ω0tω0t+1

Plot the solution for the critically damped case.

w = 1;
T = 4*pi;

fplot(solCritical(t, w), [0 T])
xlabel('\omega_0 t');
ylabel('x');
title('Critically damped, \zeta = 1');
grid on;
xticks(T*linspace(0,1,5));
xticklabels({'0','\pi','2\pi','3\pi','4\pi'});

Figure contains an axes object. The axes object with title Critically damped, zeta blank = blank 1, xlabel omega indexOf 0 baseline blank t, ylabel x contains an object of type functionline.

6. Conclusion

We have examined the different damping states for the harmonic oscillator by solving the ODEs which represents its motion using the damping ratio ζ. Plot all three cases together to compare and contrast them.

zOver  = pi;
zUnder = 1/zOver;
w = 1;
T = 2*pi;
lineStyle = {'-','--',':k'};

fplot(@(t)solOver(t, w, zOver), [0 T], lineStyle{1},'LineWidth',2);
hold on;
fplot(solCritical(t, w), [0 T], lineStyle{2},'LineWidth',2)
fplot(@(t)solUnder(t, w, zUnder), [0 T], lineStyle{3},'LineWidth',2);
hold off;

textColor = lines(3);
text(3*pi/2, 0.3 , 'over-damped'      ,'Color',textColor(1,:));
text(pi*3/4, 0.05, 'critically-damped','Color',textColor(2,:));
text(pi/8  , -0.1, 'under-damped');

grid on;
xlabel('\omega_0 t');
ylabel('amplitude');
xticks(T*linspace(0,1,5));
xticklabels({'0','\pi/2','\pi','3\pi/2','2\pi'});
yticks((1/exp(1))*[-1 0 1 2 exp(1)]);
yticklabels({'-1/e','0','1/e','2/e','1'});
lgd = legend('\pi','1','1/\pi');
title(lgd,'\zeta');
title('Damped Harmonic Oscillator');

Figure contains an axes object. The axes object with title Damped Harmonic Oscillator, xlabel omega indexOf 0 baseline blank t, ylabel amplitude contains 6 objects of type functionline, text. These objects represent \pi, 1, 1/\pi.

See Also

Functions