CoCalc Public FilesNotes / Damped and Forced Oscillations.ipynbOpen with one click!
Author: Evan Halstead
Views : 71
Compute Environment: Ubuntu 18.04 (Deprecated)


Damped and Forced Oscillations

Puzzle of the Day

Why did this happen to the Millenium Bridge in London on the day that it opened?


When a horizontal, frictionless oscillator is subject to linear drag, the resulting equation of motion is mx¨+bx˙+kx=0.m\ddot{x}+b\dot{x}+kx=0.

a) Verify that the above equation is true, then simplify it using β=b2m\beta=\frac{b}{2m} and ω0=km\omega_0=\sqrt{\frac{k}{m}}.

b) Try the solution x(t)=ertx(t)=e^{rt} then solve for rr in terms of β\beta and ω0\omega_0.

c) Interpret the solutions for β<ω0\beta<\omega_0 and β>ω0\beta>\omega_0.


a) If we suppose that the cart is to the right of the equilibrium point (i.e. in the +x-direction), then according to the free body diagram






b) Let's try the solution x(t)=ertx(t)=e^{rt}.



Therefore the equation from Newton's Second Law becomes

r2ert+2βrert+ω02ert=0r^2\cancel{e^{rt}}+2\beta r\cancel{e^{rt}}+\omega_0^2\cancel{e^{rt}}=0.

We can solve this for rr.

r=2β±4β24ω022r=\frac{-2\beta \pm\sqrt{4\beta^2-4\omega_0^2}}{2}

r=β±β2ω02r=-\beta \pm \sqrt{\beta^2-\omega_0^2}

c) Since rr has two different values, that means that our general solution can be written as a linear sum of those two solutions.


x(t)=eβt(C1eβ2ω02t+C2eβ2ω02t)x(t)=e^{-\beta t}\left (C_1e^{\sqrt{\beta^2-\omega_0^2}t}+C_2e^{-\sqrt{\beta^2-\omega_0^2}t}\right )

If β<ω0\beta<\omega_0, then the square roots become imaginary.

β2ω02=iω02β2ω1\sqrt{\beta^2-\omega_0^2}=i\sqrt{\omega_0^2-\beta^2}\equiv \omega_1

This means the terms in parentheses represent oscillations with frequency ω1=ω02β2\omega_1=\sqrt{\omega_0^2-\beta^2}. The stuff out front of the parentheses, meanwhile, represent an exponential decay. Putting it all together, when β<ω0\beta<\omega_0 the result is an exponentially decaying oscillator.

x(t)=eβt(C1eiω1t+C2eiω1t)x(t)=e^{-\beta t}\left (C_1e^{i\omega_1t}+C_2e^{-i\omega_1t}\right )

If β>ω0\beta>\omega_0, then the square root is positive, which means that no oscillation occurs.

In the simulation below, try playing around with the parameters so that you can see what each of these situations looks like.

In [1]:
from __future__ import division, print_function from vpython import * import matplotlib.pyplot as plt %matplotlib inline L=1 #relaxed length of spring omega0=10 #natural frequency beta=15 #decay parameter x0=L+0.2 #starting displacement t=0 dt=0.01 #time increment xdata=[] #set up data array for plotting x-coordinates tdata=[] #set up data array for plotting t-coordinates #draw box and spring box=box(pos=vec(x0,0.2,0),length=0.4,height=0.4, width=0.4,, v=vec(0,0,0)) spring=helix(pos=vec(0,0.2,0),axis=box.pos-vec(0,0.2,0), radius=0.1) while t<=10: rate(1/dt) box.a=-2*beta*box.v-omega0**2*(box.pos-vec(L,0.2,0)) #calculate acceleration box.v=box.v+box.a*dt #update velocity box.pos=box.pos+box.v*dt #update position spring.axis=box.pos-vec(0,0.2,0) #change spring length t=t+dt #increment time xdata.append(box.pos.x) #add x-coordinate to data array for plotting purposes tdata.append(t) #add t to data array for plotting purposes plt.xlabel("t (s)") plt.ylabel("x (m)") plt.plot(tdata,xdata,"b-")
(not running untrusted Javascript)
(not running untrusted Javascript)
(not running untrusted Javascript)
(not running untrusted Javascript)
(not running untrusted Javascript)
(not running untrusted Javascript)
(not running untrusted Javascript)
[<matplotlib.lines.Line2D at 0x7f6651917780>]


In the previous problem, we didn't look at the case β=ω0\beta=\omega_0.

a) Show that the general solution x(t)=eβt(C1eβ2ω02t+C2eβ2ω02t)x(t)=e^{-\beta t}\left (C_1e^{\sqrt{\beta^2-\omega_0^2}t}+C_2e^{-\sqrt{\beta^2-\omega_0^2}t}\right ) is not able to describe the cart if the initial conditions are such that x(0)=0x(0)=0 with x˙(0)=v00\dot{x}(0)=v_0\ne 0.

b) Show that x(t)=C1eβt+C2teβtx(t)=C_1e^{-\beta t}+C_2te^{-\beta t} is a solution to x¨+2βx˙+ω02x=0\ddot{x}+2\beta\dot{x}+\omega_0^2x=0 when β=ω0\beta=\omega_0.

c) Solve for the unknown coefficients in the special case x(0)=0x(0)=0 and x˙(0)=v0\dot{x}(0)=v_0.


a) x(0)=C1+C2=0x(0)=C_1+C_2=0

x˙(0)=β(C1+C2)=v0\dot{x}(0)=-\beta (C_1+C_2)=v_0

According to the first equation, C1=C2C1=-C_2. But if that is true, then the second equation is immediately violated.

b) x˙(t)=βC1eβtβC2teβt+C2eβt\dot{x}(t)=-\beta C_1e^{-\beta t}-\beta C_2te^{-\beta t}+C_2e^{-\beta t}

x¨(t)=β2C1eβt+β2C2teβtβC2eβtβC2eβt=β2C1eβt+β2C2teβt2βC2eβt\ddot{x}(t)=\beta^2C_1e^{-\beta t}+\beta^2C_2te^{-\beta t}-\beta C_2e^{-\beta t}-\beta C_2e^{-\beta t}=\beta^2C_1e^{-\beta t}+\beta^2C_2te^{-\beta t}-2\beta C_2e^{-\beta t}

Therefore the equation x¨+2βx˙+ω02x=0\ddot{x}+2\beta\dot{x}+\omega_0^2x=0 becomes

β2C1eβt+β2C2teβt2βC2eβt+2β(βC1eβtβC2teβt+C2eβt)+β2(C1eβt+C2teβt)=?0\beta^2C_1e^{-\beta t}+\beta^2C_2te^{-\beta t}-2\beta C_2e^{-\beta t}+2\beta \left (-\beta C_1e^{-\beta t}-\beta C_2te^{-\beta t}+C_2e^{-\beta t}\right )+\beta^2\left (C_1e^{-\beta t}+C_2te^{-\beta t}\right )\stackrel{?}{=}0

β2C1eβt+β2C2teβt2βC2eβt2β2C1eβt2β2C2teβt+2βC2eβt+β2C1eβt+β2C2teβt=?0\cancel{\beta^2C_1e^{-\beta t}}+\cancel{\beta^2C_2te^{-\beta t}}-\cancel{2\beta C_2e^{-\beta t}}-\cancel{2\beta^2C_1e^{-\beta t}}-\cancel{2\beta^2C_2te^{-\beta t}}+\cancel{2\beta C_2e^{-\beta t}}+\cancel{\beta^2C_1e^{-\beta t}}+\cancel{\beta^2C_2te^{-\beta t}}\stackrel{?}{=}0


c) x(0)=C1=0x(0)=C_1=0


Damped Oscillations

  • in damped oscillators, the amplitude of oscillations decreases as energy is dissipated away
  • in the special case of linear drag, the equation of motion is x¨+2βx˙+ω02x=0,\ddot{x}+2\beta\dot{x}+\omega_0^2x=0, where β=b2m\beta=\frac{b}{2m} is a decay parameter and ω0=km\omega_0=\sqrt{\frac{k}{m}} is the natural frequency (i.e. the oscillation frequency in the absence of damping)
damping how does β\beta compare to ω0\omega_0? decay parameter \hspace{1cm}general solution\hspace{1cm}
none β=0\beta=0 0 x(t)=C1eiω0t+C2eiω0tx(t)=C_1e^{i\omega_0 t}+C_2e^{-i\omega_0t}
underdamped β<ω0\beta<\omega_0 β\beta x(t)=eβt(C1eβ2ω02t+C2eβ2ω02t)x(t)=e^{-\beta t} (C_1e^{\sqrt{\beta^2-\omega_0^2}t}+C_2e^{-\sqrt{\beta^2-\omega_0^2}t})
critically damped β=ω0\beta=\omega_0 β\beta x(t)=C1eβt+C2teβtx(t)=C_1e^{-\beta t}+C_2te^{-\beta t}
overdamped β>ω0\beta>\omega_0 ββ2ω22\beta-\sqrt{\beta^2-\omega_2^2} x(t)=eβt(C1eβ2ω02t+C2eβ2ω02t)x(t)=e^{-\beta t} (C_1e^{\sqrt{\beta^2-\omega_0^2}t}+C_2e^{-\sqrt{\beta^2-\omega_0^2}t} )

Linear Differential Operators

Suppose we define D=d2dt2+2βddt+ω02.D=\frac{d^2}{dt^2}+2\beta\frac{d}{dt}+\omega_0^2.

Then the linear damped oscillator equation can be written as Dx=0.Dx=0.

Now suppose that there is an additional external force that doesn't depend on x(t)x(t) at all: x¨+2βx˙+ω02x=f(t),\ddot{x}+2\beta\dot{x}+\omega_0^2x=f(t), where f(t)F(t)mf(t)\equiv \frac{F(t)}{m}. (This is called a forced oscillator.) Then the linear damped oscillator equation could be written as Dx=f.Dx=f.

We will now show that the full solution to Dx=fDx=f consists of the solution to Dx=fDx=f as well as the solution to Dx=0Dx=0.

First, if DD is a linear operator, then D(ax1+bx2)=aDx1+bDx2.D(ax_1+bx_2)=aDx_1+bDx_2. Suppose we somehow found a solution that works in the forced oscillator equation. We will call this the particular solution.


We will call the solution to the unforced oscillator the homogeneous solution.


Now because DD is a linear operator, then D(xp+xh)=Dxp+Dxh=f+0=f.D(x_p+x_h)=Dx_p+Dx_h=f+0=f.

Therefore the most general complete solution to Dx=fDx=f must also include solutions to the corresponding unforced (i.e. homogeneous) oscillator equation.


An oscillator is forced with a forcing function f(t)=f0cos(ωt)f(t)=f_0\cos(\omega t), where ωω0\omega\ne\omega_0. (An example similar to this would be pushing someone who is on a swing.)

Show that x(t)=Acos(ωtδ)x(t)=A\cos(\omega t-\delta) satisfies the equation, where A2=f02(ω02ω2)2+4β2ω2A^2=\frac{f_0^2}{(\omega_0^2-\omega^2)^2+4\beta^2\omega^2} and δ=arctan(2βωω02ω2).\delta =\arctan\left (\frac{2\beta \omega}{\omega_0^2-\omega^2}\right ).


First it may help to write x(t)=[Aei(ωtδ)]x(t)=\Re\left [Ae^{i(\omega t-\delta)}\right ].

x˙=[iωAei(ωtδ)]\dot{x}=\Re\left [i\omega Ae^{i(\omega t-\delta )}\right ]

x¨=[ω2Aei(ωtδ)]\ddot{x}=\Re\left [-\omega^2Ae^{i(\omega t-\delta )}\right ]

Then the forced oscillator equation is [ω2Aei(ωtδ)+2βωiAei(ωtδ)+ω02Aei(ωtδ)=f0eiωt].\Re\left [-\omega^2Ae^{i(\omega t-\delta)}+2\beta\omega iAe^{i(\omega t-\delta)}+\omega_0^2Ae^{i(\omega t-\delta)}=f_0e^{i\omega t}\right ].

Every term contains eiωte^{i\omega t}, so we can cancel it out.

[ω2Aeiδ+2βωiAeiδ+ω02Aeiδ=f0]\Re\left [-\omega^2Ae^{-i\delta}+2\beta\omega iAe^{-i\delta}+\omega_0^2Ae^{-i\delta}=f_0\right ]

Now we can turn it around and solve for AA.

A=f0eiδω02ω2+2βωiA=\frac{f_0e^{i\delta }}{\omega_0^2-\omega^2+2\beta\omega i}


We can also solve for δ\delta.

eiδ=cosδ+isinδ=A(ω2ω02+2βωi)f0=A(ω2ω02)f0+i2Aβωf0e^{i\delta}=\cos\delta +i\sin\delta=\frac{A(\omega^2-\omega_0^2+2\beta\omega i)}{f_0}=\frac{A(\omega^2-\omega_0^2)}{f_0}+i\frac{2A\beta\omega}{f_0}

Then we can match the real terms and the imaginary terms.



To solve for δ\delta we can divide the first equation by the second.


δ=arctan(2βωω2ω02)\delta=\arctan\left (\frac{2\beta\omega}{\omega^2-\omega_0^2}\right )

Driven Oscillators

The general solution to x¨+2βx˙+ω02=f0cos(ωt)\ddot{x}+2\beta\dot{x}+\omega_0^2=f_0\cos(\omega t) is x(t)=Acos(ωtδ)+eβt(C1eβ2ω02t+C2eβ2ω02t),x(t)=A\cos(\omega t-\delta)+e^{-\beta t}\left (C_1e^{\sqrt{\beta^2-\omega_0^2}t}+C_2e^{-\sqrt{\beta^2-\omega_0^2}t}\right ), where A2=f02(ω02ω2)2+4β2ω2A^2=\frac{f_0^2}{(\omega_0^2-\omega^2)^2+4\beta^2\omega^2} and δ=arctan(2βωω02ω2).\delta =\arctan\left (\frac{2\beta\omega}{\omega_0^2-\omega^2}\right ).


a) Describe the long-term behavior of a damped, sinusoidally forced oscillator.

b) For what value of the forcing frequency ω\omega is the amplitude of the long-term oscillations the greatest?


a) As time goes on, the eβte^{-\beta t} term makes the natural oscillation disappear. It is because of this that the natural oscillation terms are referred to as transients. What remains is the oscillation due to the forcing.

b) AA is maximized when ω=ω0\omega=\omega_0. Recall that ω0\omega_0 is the frequency of oscillation without forcing. This means that the amplitude of oscillation is biggest when the driving frequency is equal to the natural frequency. This condition is known as resonance. The idea of resonance makes sense when you think about pushing someone on a swing. You push only when they come back to you. You wouldn't give them a starting push and then push them again while they are on their way back to you mid-swing.

This, by the way, helps us answer what happened to the Millenium Bridge. By an unfortunate accident, the tension cables of the bridge happened to create a system in which the natural frequency matched the frequency of footfalls on the bridge on that day. It didn't help that the swinging made people adjust their gait in such a way that amplified the problem.

So how was the problem solved? The bridge was closed down and extra dampeners were added to the bridge.


The equation of motion of a damped, forced oscillator is x¨+2βx˙+ω02x=f0cos(ωt)\ddot{x}+2\beta\dot{x}+\omega_0^2x=f_0\cos(\omega t) with x(0)=0x(0)=0, x˙(0)=0\dot{x}(0)=0, ω=2π\omega=2\pi, ω0=5ω\omega_0=5\omega, β=ω0/20\beta=\omega_0/20, and f0=1000f_0=1000. Recall that the general solution is

x(t)=Acos(ωtδ)+eβt(C1eβ2ω02t+C2eβ2ω02t)x(t)=A\cos(\omega t-\delta)+e^{-\beta t}\left (C_1e^{\sqrt{\beta^2-\omega_0^2}t}+C_2e^{-\sqrt{\beta^2-\omega_0^2}t}\right ).

a) Use the initial conditions to solve for C1C_1 and C2C_2 in terms of AA and δ\delta.

b) Create a plot of x(t)x(t).


a) The general solution is

x(t)=Acos(ωtδ)+eβt(C1eβ2ω02t+C2eβ2ω02t)x(t)=A\cos(\omega t-\delta)+e^{-\beta t}\left (C_1e^{\sqrt{\beta^2-\omega_0^2}t}+C_2e^{-\sqrt{\beta^2-\omega_0^2}t}\right ).

The first derivative is

x˙(t)=Aωsin(ωtδ)βeβt(C1eiω1t+C2eiω1t)+eβt(iω1C1eiω1tiω1C2eiω1t)\dot{x}(t)=-A\omega \sin(\omega t-\delta)-\beta e^{-\beta t}\left (C_1e^{i\omega_1t}+C_2e^{-i\omega_1t}\right )+e^{-\beta t}\left (i\omega_1C_1e^{i\omega_1 t}-i\omega_1C_2e^{-i\omega_1t}\right ),

where ω1iω02β2\omega_1\equiv i\sqrt{\omega_0^2-\beta^2} (notice that β<ω0\beta<\omega_0 in this problem, which is why I switched the order and added ii).

Substituting initial conditions, we get



Since the right hand side of the second equation is purely real, that means that iω1(C1C2)=0C1=C2i\omega_1(C_1-C_2)=0\Rightarrow C_1=C_2.

Then the first equation becomes

Acos(δ)+2C1=0C1=C2=Acos(δ)2A\cos(-\delta)+2C_1=0\Rightarrow C_1=C_2=-\frac{A\cos(-\delta)}{2}.

Putting this back into the general solution, we get

x(t)=Acos(ωtδ)Acos(δ)2eβt(eiω1t+eiω1t)x(t)=A\cos(\omega t-\delta)-\frac{A\cos(-\delta)}{2}e^{-\beta t}\left(e^{i\omega_1 t}+e^{-i\omega_1t}\right ).

Since x(t)x(t) must be real, we take the real part of the function.

[x(t)]=[eiω1t+eiω1t]=[cos(ω1t)+isin(ω1t)+cos(ω1t)+isin(ω1t)]=cos(ω1t)+cos(ω1t)=2cos(ω1t)\Re[x(t)]=\Re\left [e^{i\omega_1t}+e^{-i\omega_1t}\right ]=\Re\left [\cos(\omega_1t)+i\sin(\omega_1t)+\cos(-\omega_1t)+i\sin(-\omega_1t)\right]=\cos(\omega_1t)+\cos(-\omega_1t)=2\cos(\omega_1t)

In the last step, we used the fact that cos(x)=cos(x)\cos(x)=\cos(-x). The solution can be simplified further.

x(t)=Acos(ωtδ)Acos(δ)2eβt(2cos(ω1t)=Acos(ωtδ)Acos(δ)eβtcos(ω1t)x(t)=A\cos(\omega t-\delta)-\frac{A\cos(-\delta)}{2}e^{-\beta t}(2\cos(\omega_1t)=A\cos(\omega t-\delta)-A\cos(-\delta)e^{-\beta t}\cos(\omega_1 t)

b) see the python code below

In [14]:
from __future__ import division, print_function import matplotlib.pyplot as plt import numpy as np %matplotlib inline x=0 xdot=0 omega=2*pi omega0=5*omega beta=omega0/20 f0=1000 A=sqrt(f0**2/((omega0**2-omega**2)**2-4*beta**2*omega**2)) delta=np.arctan(2*beta*omega/(omega0**2-omega**2)) omega1=sqrt(omega0**2-beta**2) t=0 dt=0.001 xdata=[] xexact=[] tdata=[] while t<10: xddot=f0*cos(omega*t)-2*beta*xdot-omega0**2*x xdot=xdot+xddot*dt x=x+xdot*dt t=t+dt xdata.append(x) xexact.append(A*cos(omega*t-delta)-A*cos(-delta)*exp(-beta*t)*cos(omega1*t)) tdata.append(t) print(A, delta) plt.xlabel("t (s)") plt.ylabel("x (m)") plt.plot(tdata,xdata,"b.", tdata, xexact, "r-") plt.legend(["numerical","exact"])
1.05565811361 0.0208303200362
<matplotlib.legend.Legend at 0x11a076310>
In [ ]: