Damped and Forced Oscillations.ipynb

Kernel: VPython

ParseError: KaTeX parse error: Undefined control sequence: \require at position 1: \̲r̲e̲q̲u̲i̲r̲e̲{cancel}ParseError: KaTeX parse error: Undefined control sequence: \require at position 1: \̲r̲e̲q̲u̲i̲r̲e̲{enclose}

PY345

Damped and Forced Oscillations

Puzzle of the Day

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

Practice

---

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

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

b) Try the solution $x(t)=e^{rt}$ then solve for $r$ in terms of $\beta$ and $\omega_0$.

c) Interpret the solutions for $\beta<\omega_0$ and $\beta>\omega_0$.

Solution

---

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

$-F_{sc}^s-F_{ac}^d=m\ddot{x}$.

$-kx-bv=m\ddot{x}$

$m\ddot{x}+b\dot{x}+kx=0$

$\ddot{x}+\frac{b}{m}\dot{x}+\frac{k}{m}x=0$

$\ddot{x}+2\beta\dot{x}+\omega_0^2x=0$

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

$\dot{x}(t)=re^{rt}$

$\ddot{x}(t)=r^2e^{rt}$

Therefore the equation from Newton's Second Law becomes

$r^2\cancel{e^{rt}}+2\beta r\cancel{e^{rt}}+\omega_0^2\cancel{e^{rt}}=0$.

We can solve this for $r$.

$r=\frac{-2\beta \pm\sqrt{4\beta^2-4\omega_0^2}}{2}$

$r=-\beta \pm \sqrt{\beta^2-\omega_0^2}$

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

$x(t)=C_1e^{r_1t}+C_2e^{r_2t}$

$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 $\beta<\omega_0$, then the square roots become imaginary.

$\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 $\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 $\beta<\omega_0$ the result is an exponentially decaying oscillator.

$x(t)=e^{-\beta t}\left (C_1e^{i\omega_1t}+C_2e^{-i\omega_1t}\right )$

If $\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.

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, color=color.blue, 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-")

<IPython.core.display.Javascript object>

<IPython.core.display.Javascript object>

<IPython.core.display.Javascript object>

<IPython.core.display.Javascript object>

<IPython.core.display.Javascript object>

<IPython.core.display.Javascript object>

<IPython.core.display.Javascript object>

[<matplotlib.lines.Line2D at 0x7f6651917780>]

Practice

---

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

a) Show that the general solution $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)=0$ with $\dot{x}(0)=v_0\ne 0$.

b) Show that $x(t)=C_1e^{-\beta t}+C_2te^{-\beta t}$ is a solution to $\ddot{x}+2\beta\dot{x}+\omega_0^2x=0$ when $\beta=\omega_0$.

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

Solution

---

a) $x(0)=C_1+C_2=0$

$\dot{x}(0)=-\beta (C_1+C_2)=v_0$

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

b) $\dot{x}(t)=-\beta C_1e^{-\beta t}-\beta C_2te^{-\beta t}+C_2e^{-\beta 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 $\ddot{x}+2\beta\dot{x}+\omega_0^2x=0$ becomes

$\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$

$\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$

$0\stackrel{?}{=}0$

c) $x(0)=C_1=0$

ParseError: KaTeX parse error: Undefined control sequence: \cancelto at position 19: …t{x}(0)=-\beta \̲c̲a̲n̲c̲e̲l̲t̲o̲{0}{C_1}-0+C_2=…

$C_2=v_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 $$\ddot{x}+2\beta\dot{x}+\omega_0^2x=0,$$ where $\beta=\frac{b}{2m}$ is a decay parameter and $\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 $\omega_0$?	decay parameter	$\hspace{1cm}$general solution$\hspace{1cm}$
none	$\beta=0$	0	$x(t)=C_1e^{i\omega_0 t}+C_2e^{-i\omega_0t}$
underdamped	$\beta<\omega_0$	$\beta$	$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	$\beta=\omega_0$	$\beta$	$x(t)=C_1e^{-\beta t}+C_2te^{-\beta t}$
overdamped	$\beta>\omega_0$	$\beta-\sqrt{\beta^2-\omega_2^2}$	$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=\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.$$

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

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

First, if $D$ is a linear operator, then $$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.

$$Dx_p=f$$

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

$$Dx_h=0$$

Now because $D$ is a linear operator, then $$D(x_p+x_h)=Dx_p+Dx_h=f+0=f.$$

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

Practice

---

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

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

Solution

---

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

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

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

Then the forced oscillator equation is $$\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 $e^{i\omega t}$, so we can cancel it out.

$$\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 $A$.

$$A=\frac{f_0e^{i\delta }}{\omega_0^2-\omega^2+2\beta\omega i}$$$$A^2=AA^*=\frac{f_0^2}{(\omega_0^2-\omega^2)^2+4\beta^2\omega^2}$$

We can also solve for $\delta$.

$$e^{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.

$$\sin\delta=\frac{2A\beta\omega}{f_0}$$$$\cos\delta=\frac{A(\omega^2-\omega_0^2)}{f_0}$$

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

$$\tan\delta=\frac{\frac{2\cancel{A}\beta\omega}{\cancel{f_0}}}{\frac{\cancel{A}(\omega^2-\omega_0^2)}{\cancel{f_0}}}$$$$\delta=\arctan\left (\frac{2\beta\omega}{\omega^2-\omega_0^2}\right )$$

Driven Oscillators

The general solution to $\ddot{x}+2\beta\dot{x}+\omega_0^2=f_0\cos(\omega t)$ is $$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 $$A^2=\frac{f_0^2}{(\omega_0^2-\omega^2)^2+4\beta^2\omega^2}$$ and $$\delta =\arctan\left (\frac{2\beta\omega}{\omega_0^2-\omega^2}\right ).$$

Practice

---

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?

Solution

---

a) As time goes on, the $e^{-\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) $A$ is maximized when $\omega=\omega_0$. Recall that $\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.

Practice

---

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

$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 $C_1$ and $C_2$ in terms of $A$ and $\delta$.

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

Solution

---

a) The general solution is

$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

$\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 $\omega_1\equiv i\sqrt{\omega_0^2-\beta^2}$ (notice that $\beta<\omega_0$ in this problem, which is why I switched the order and added $i$).

Substituting initial conditions, we get

$x(0)=A\cos(-\delta)+C_1+C_2=0$

$\dot{x}(0)=-A\sin(-\delta)-\beta(C_1+C_2)+i\omega_1C_1-i\omega_1C_2=0$.

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

Then the first equation becomes

$A\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)=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)$ must be real, we take the real part of the function.

$\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)$. The solution can be simplified further.

$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

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>