# MATH 1003, Lecture 21&22¶

## Lecturer: Eduardo G. Altmann, http://www.maths.usyd.edu.au/u/ega/¶

In [1]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

from matplotlib import animation, rc
from IPython.display import HTML


## Damped Harmonic Oscillator $$\frac{d^2 y}{dx^2} + 2 \gamma \frac{dy}{dx} + \omega^2 y =0$$¶

with $\gamma,\omega >0$

## Case 1 (supercritical damping): $\gamma > \omega$¶

In [2]:
gamma=1.8
omega=1
kappa = np.sqrt(np.absolute(gamma**2-omega**2))

In [3]:
dt=0.05
N=200
yh=[]
xh=[]
th=[]
ones=np.power(1,np.zeros(N))
for i in range(N):
t=i*dt
th.append(t)
yh.append(0)
#    xh.append(1+0.8*np.cos(2*np.pi*t))
#    xh.append(1+0.8*np.exp(-gamma*t)*np.cos(2*np.pi*kappa*t))
xh.append(1+np.exp(-gamma*t)*(0.5*np.exp(kappa*t)+0.5*np.exp(-kappa*t)))
#print(np.exp(-(gamma-kappa)*t))

In [4]:
#Plot
%matplotlib inline

plt.plot(th,xh-ones)
plt.xlabel("time t")
plt.ylabel("position x")
plt.show()

# Animation

fig = plt.figure()
ax = fig.add_subplot(111, xlim=(-0.2, 2), ylim=(-1, 1),yticks=[],xticks=[])
ax.grid()

line, = ax.plot([], [], 'o-', lw=25,c="orange",ls=":",ms=50,marker="s",mfc="gray",fillstyle="none",mec="black",markevery=2)
line2,=ax.plot([],[],lw=20,c="black")
line3,=ax.plot([],[],lw=2,c="black",marker="^")
time_template = 'Damped Harmonic Oscillator,\n  K='+str(kappa)[:2]+', Damping $\gamma$='+str(gamma)[:2]+', Freq. w^2='+str(omega*omega)[:4]+'\n time = %.1fs'
time_text = ax.text(0.25, 0.9, '', transform=ax.transAxes)

def init():
line.set_data([], [])
line2.set_data([],[])
time_text.set_text('')
return line, time_text

def animate(i):
thisx = [xh[i],0]
thisy = [yh[i],0]

line.set_data(thisx, thisy)
line2.set_data([-0.05,-0.05],[-1,1])
line3.set_data([1,1],[-2,-0.4])
time_text.set_text(time_template % (i*dt))
return line, time_text

ani = animation.FuncAnimation(fig, animate, np.arange(1, N),
interval=50, blit=True, init_func=init,repeat=False)

#ani.save('double_pendulum.mp4', fps=15)
#plt.show()

HTML(ani.to_html5_video())

Out[4]:

## Case 2 (subcritical damping): $\gamma < \omega$

In [5]:
gamma=0.3
omega=2
kappa = np.sqrt(np.absolute(gamma**2-omega**2))

In [6]:
dt=0.05
N=200
yh=[]
xh=[]
th=[]
ones=np.power(1,np.zeros(N))

for i in range(N):
t=i*dt
th.append(t)
yh.append(0)
#    xh.append(1+0.8*np.cos(2*np.pi*t))
xh.append(1+0.8*np.exp(-gamma*t)*np.cos(kappa*t))
#    xh.append(1+np.exp(-gamma*t)*(0.5*np.exp(kappa*t)+0.5*np.exp(-kappa*t)))
#print(np.exp(-(gamma-kappa)*t))

In [7]:
#Plot
%matplotlib inline

plt.plot(th,xh-ones)
plt.xlabel("time t")
plt.ylabel("position x")
plt.show()

# Animation

fig = plt.figure()
ax = fig.add_subplot(111, xlim=(-0.2, 2), ylim=(-1, 1),yticks=[],xticks=[])
ax.grid()

line, = ax.plot([], [], 'o-', lw=25,c="orange",ls=":",ms=50,marker="s",mfc="gray",fillstyle="none",mec="black",markevery=2)
line2,=ax.plot([],[],lw=20,c="black")
line3,=ax.plot([],[],lw=2,c="black",marker="^")
time_template = 'Damped Harmonic Oscillator,\n  K='+str(kappa)[:2]+', Damping $\gamma$='+str(gamma)[:2]+', Freq. w^2='+str(omega*omega)[:4]+'\n time = %.1fs'
time_text = ax.text(0.25, 0.9, '', transform=ax.transAxes)

def init():
line.set_data([], [])
line2.set_data([],[])
time_text.set_text('')
return line, time_text

def animate(i):
thisx = [xh[i],0]
thisy = [yh[i],0]

line.set_data(thisx, thisy)
line2.set_data([-0.05,-0.05],[-1,1])
line3.set_data([1,1],[-2,-0.4])
time_text.set_text(time_template % (i*dt))
return line, time_text

ani = animation.FuncAnimation(fig, animate, np.arange(1, N),
interval=50, blit=True, init_func=init,repeat=False)

#ani.save('double_pendulum.mp4', fps=15)
#plt.show()

HTML(ani.to_html5_video())

Out[7]:

## Case 3: $\gamma = \omega$¶

In [8]:
gamma=0.5
omega=gamma
kappa = np.sqrt(np.absolute(gamma**2-omega**2))

In [9]:
dt=0.05
N=200
yh=[]
xh=[]
th=[]
ones=np.power(1,np.zeros(N))

for i in range(N):
t=i*dt
th.append(t)
yh.append(0)
#    xh.append(1+0.8*np.cos(2*np.pi*t))
#    xh.append(1+0.8*np.exp(-gamma*t)*np.cos(kappa*t))
xh.append(1+np.exp(-gamma*t)*(1+0.2*t))
#print(np.exp(-(gamma-kappa)*t))

In [10]:
#Plot
%matplotlib inline

plt.plot(th,xh-ones)
plt.xlabel("time t")
plt.ylabel("position x")
plt.show()

# Animation

fig = plt.figure()
ax = fig.add_subplot(111, xlim=(-0.2, 2), ylim=(-1, 1),yticks=[],xticks=[])
ax.grid()

line, = ax.plot([], [], 'o-', lw=25,c="orange",ls=":",ms=50,marker="s",mfc="gray",fillstyle="none",mec="black",markevery=2)
line2,=ax.plot([],[],lw=20,c="black")
line3,=ax.plot([],[],lw=2,c="black",marker="^")
time_template = 'Damped Harmonic Oscillator,\n  K='+str(kappa)[:2]+', Damping $\gamma$='+str(gamma)[:2]+', Freq. w^2='+str(omega*omega)[:4]+'\n time = %.1fs'
time_text = ax.text(0.25, 0.9, '', transform=ax.transAxes)

def init():
line.set_data([], [])
line2.set_data([],[])
time_text.set_text('')
return line, time_text

def animate(i):
thisx = [xh[i],0]
thisy = [yh[i],0]

line.set_data(thisx, thisy)
line2.set_data([-0.05,-0.05],[-1,1])
line3.set_data([1,1],[-2,-0.4])
time_text.set_text(time_template % (i*dt))
return line, time_text

ani = animation.FuncAnimation(fig, animate, np.arange(1, N),
interval=50, blit=True, init_func=init,repeat=False)

#ani.save('double_pendulum.mp4', fps=15)
#plt.show()

HTML(ani.to_html5_video())

Out[10]:

# Varying parameters and initial conditions¶

## Restart Kernel to Work¶

In [2]:
%matplotlib notebook
from pylab import *
import numpy as np

fig, ax = plt.subplots()
t = np.arange(0.0, 10.0, 0.001)
#a0 = 5
#f0 = 0.2
phi=0

x0=1
v0=0
w=0.2
gamma=0

A=np.sqrt(x0**2+(v0/w)**2)
phi=np.arccos(x0/A)

s = A*np.exp(-gamma*t)*np.cos(2*np.pi*np.sqrt(np.absolute(gamma**2-w**2))*t+phi)
l, = plt.plot(t, s, lw=2, color='red')
plt.axis([0, 10, -3, 3])

axcolor = 'lightgoldenrodyellow'
axfreq = plt.axes([0.25, 0.15, 0.65, 0.03], axisbg=axcolor)
axgamma = plt.axes([0.25, 0.1, 0.65, 0.03], axisbg=axcolor)

#axamp = plt.axes([0.25, 0.15, 0.65, 0.03], axisbg=axcolor)
axx0 = plt.axes([0.25, 0.05, 0.65, 0.03], axisbg=axcolor)
axv0 = plt.axes([0.25, 0.0, 0.65, 0.03], axisbg=axcolor)

sfreq = Slider(axfreq, 'Freq $\omega$', 0.01, 1.0, valinit=w)
sgamma = Slider(axgamma, 'Damp $\gamma$', -1, 3.0, valinit=gamma)
sx0 = Slider(axx0, 'Init Pos $x(0)$', -1.0, 1.0, valinit=x0)
sv0 = Slider(axv0, 'Init Vel $v(0)$', -1.0, 1.0, valinit=v0)

def update(val):
x0=sx0.val
v0=sv0.val
gamma=sgamma.val
w = sfreq.val
A=np.sqrt(x0**2+(v0/w)**2)
phi=-np.sign(v0)*np.arccos(x0/A)
#    l.set_ydata(amp*np.sin(2*np.pi*freq*t))
l.set_ydata(A*np.exp(-gamma*t)*np.cos(2*np.pi*np.sqrt(np.absolute(gamma**2-w**2))*t+phi))
fig.canvas.draw_idle()

sfreq.on_changed(update)
sx0.on_changed(update)
sv0.on_changed(update)

plt.show()


## Forced Damped Harmonic Oscillator -- Resoance $$\frac{d^2 y}{dx^2} + 2 \gamma \frac{dy}{dx} + \omega^2 y = F cos(w_f t)$$¶

In [3]:
fig, ax = plt.subplots()
t = np.arange(0.0, 20.0, 0.001)
#a0 = 5
#f0 = 0.2
phi=0

x0=1
v0=0
w=0.2
gamma=0.2

f=1
wf=1

f=1/np.sqrt((2*w*gamma)**2+(w**2-wf**2)**2)

A=np.sqrt(x0**2+(v0/w)**2)
phi=np.arccos(x0/A)

s = A*np.exp(-gamma*t)*np.cos(2*np.pi*np.sqrt(gamma**2-w**2)*t+phi) + f*np.sin(2*np.pi*wf*t)
l, = plt.plot(t, s, lw=2, color='red')
plt.axis([0, 20, -5, 5])

axcolor = 'lightgoldenrodyellow'
axfreq = plt.axes([0.25, 0.15, 0.65, 0.03], axisbg=axcolor)
axgamma = plt.axes([0.25, 0.1, 0.65, 0.03], axisbg=axcolor)

axx0 = plt.axes([0.25, 0.05, 0.65, 0.03], axisbg=axcolor)
axv0 = plt.axes([0.25, 0.0, 0.65, 0.03], axisbg=axcolor)

sfreq = Slider(axfreq, 'Freq $\omega$', 0.01, 2.0, valinit=w)
sgamma = Slider(axgamma, 'Damp $\gamma$', -1, 3.0, valinit=gamma)
sx0 = Slider(axx0, 'Init Pos $x(0)$', -1.0, 1.0, valinit=x0)
sv0 = Slider(axv0, 'Init Vel $v(0)$', -1.0, 1.0, valinit=v0)

def update(val):
x0=sx0.val
v0=sv0.val
gamma=sgamma.val
w = sfreq.val
A=np.sqrt(x0**2+(v0/w)**2)
phi=-np.sign(v0)*np.arccos(x0/A)
f=1/np.sqrt((2*w*gamma)**2+(w**2-wf**2)**2)
#    l.set_ydata(amp*np.sin(2*np.pi*freq*t))
l.set_ydata(A*np.exp(-gamma*t)*np.cos(2*np.pi*np.sqrt(np.absolute(gamma**2-w**2))*t+phi)+f*np.sin(2*np.pi*wf*t))
fig.canvas.draw_idle()

sfreq.on_changed(update)
sx0.on_changed(update)
sv0.on_changed(update)

plt.show()


## Amplitude of oscillation as a function of ratio of frequencies: Resonance Curve¶

In [11]:
w = np.arange(0.1, 3.0, 0.001)
phi=0

x0=1
v0=0

gamma=0.1
wf=1
#a0 = 5
#f0 = 0.2
phi=0

f=np.divide(1,np.sqrt((2*w*gamma)**2+(w**2-wf**2)**2))

%matplotlib inline
plt.plot(w,f)
plt.ylabel("Amplitude of the Oscillation")
plt.xlabel("natural frequency / forcing frequency = $\omega_f / \omega_0$")
plt.show()

In [ ]: