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


from matplotlib.widgets import Slider,Button, RadioButtons

fig, ax = plt.subplots()
plt.subplots_adjust(left=0.25, bottom=0.25)
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()
plt.subplots_adjust(left=0.25, bottom=0.25)
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 [ ]: