MATH 1003, Lecture 20

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

In [3]:
import matplotlib.pyplot as plt
from scipy import *
from scipy import integrate
from scipy.integrate import ode
import numpy as np
            
%load_ext autoreload
%autoreload 2

%matplotlib inline

Harmonic Oscillator

In [4]:
fig = plt.figure(num=1,figsize=(10,4))
ax=fig.add_subplot(111)

xmin=0
xmax=20

plt.xlim([xmin,xmax])
#plt.ylim([ymin,ymax])
plt.xlabel(r"$t$",fontsize=18)
plt.ylabel(r"$x$",fontsize=18)
plt.tick_params(axis='both', which='major', labelsize=16)

x=np.linspace(xmin,xmax,1000)

v0=-0.2
x0=0.5
w=1

A=v0/w 
B=x0    #=xo


# Ploting solutions derived in classroom
f1= A*sin(w*x)
f2= B*cos(w*x)
f= f1+f2


+-
ax.plot(x,f,label="y")
#ax.plot(x,f1,label=r"$y_1$")
#ax.plot(x,f2,label=r"$y_2$")

plt.legend()
plt.show()

2nd order, linear, homogeneous differential equation $$\frac{d^2 y}{dx^2} + a \frac{dy}{dx} +b y =0$$

Case 1: $\Delta \equiv a^2-4b$>0

In [5]:
fig = plt.figure(num=1,figsize=(10,5)
ax=fig.add_subplot(111)

xmin=0
xmax=5

plt.xlim([xmin,xmax])
#plt.ylim([ymin,ymax])
plt.xlabel(r"$t$",fontsize=18)
plt.ylabel(r"$x$",fontsize=18)
plt.tick_params(axis='both', which='major', labelsize=16)

x=np.linspace(xmin,xmax,1000)


m1=-0.1
m2=-0.2
A=1
B=1

# Ploting solutions derived in classroom
f1= A*exp(m1*x)
f2= B*exp(m2*x)
f= f1+f2
ax.plot(x,f)
ax.plot(x,f1)
ax.plot(x,f2)

plt.show()

  File "<ipython-input-5-876a57804a5d>", line 2
    ax=fig.add_subplot(111)
     ^
SyntaxError: invalid syntax
In [6]:
fig = plt.figure(num=1,figsize=(10,5))
ax=fig.add_subplot(111)

xmin=0
xmax=15

plt.xlim([xmin,xmax])
#plt.ylim([ymin,ymax])
plt.xlabel(r"$x$",fontsize=18)
plt.ylabel(r"$y$",fontsize=18)
plt.tick_params(axis='both', which='major', labelsize=16)

x=np.linspace(xmin,xmax,1000)


m1=-0.1
m2=-0.2
A=1
B=1

# Ploting solutions derived in classroom
f1= A*exp(m1*x)
f2= B*exp(m2*x)
f= f1+f2
ax.plot(x,f)
ax.plot(x,f1,label="m1")
ax.plot(x,f2,label="m2")
plt.legend()
plt.show()

Case 2 $\Delta < 0$

$$ y = A e^{-a t/2}cos(kx+\phi)$$
In [7]:
fig = plt.figure(num=1,figsize=(10,10))
ax=fig.add_subplot(111)

xmin=0
xmax=10

plt.xlim([xmin,xmax])
#plt.ylim([ymin,ymax])
plt.xlabel(r"$t$",fontsize=18)
plt.ylabel(r"$x$",fontsize=18)
plt.tick_params(axis='both', which='major', labelsize=16)

x=np.linspace(xmin,xmax,1000)


k=5
phi=0
A=1
a=1

# Ploting solutions derived in classroom
f=A*exp(-a*x/2)*cos(k*x+phi)
ax.plot(x,f)
ax.plot(x,A*exp(-a*x/2),label=r"$e^{-a t/2}$")
plt.legend()
plt.show()
In [8]:
fig = plt.figure(num=1,figsize=(10,10))
ax=fig.add_subplot(111)

xmin=0
xmax=10

plt.xlim([xmin,xmax])
#plt.ylim([ymin,ymax])
plt.xlabel(r"$t$",fontsize=18)
plt.ylabel(r"$x$",fontsize=18)
plt.tick_params(axis='both', which='major', labelsize=16)

x=np.linspace(xmin,xmax,1000)


k=5
phi=0
A=1
a=-1

# Ploting solutions derived in classroom
f=A*exp(-a*x/2)*cos(k*x+phi)
ax.plot(x,f)
ax.plot(x,A*exp(-a*x/2),label=r"$e^{-a t/2}$")
plt.legend()
plt.show()

Initial conditions in the Harmonic Oscilator case

$\frac{d^2 x}{dt^2} + \omega^2 x =0$, where $\omega^2=k/m$

$x(t) = A cos (\omega t + \phi)$

In [55]:
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*f0*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*w*t+phi))
    fig.canvas.draw_idle()

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


plt.show()
In [10]:
import numpy as np
import matplotlib.pyplot as plt

%matplotlib notebook
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, 1.0, 0.001)
a0 = 5
f0 = 3
s = a0*np.sin(2*np.pi*f0*t)
l, = plt.plot(t, s, lw=2, color='red')
plt.axis([0, 1, -10, 10])

axcolor = 'lightgoldenrodyellow'
axfreq = 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)


sfreq = Slider(axfreq, 'Freq', 0.1, 30.0, valinit=f0)
samp = Slider(axamp, 'Amp', 0.1, 10.0, valinit=a0)


def update(val):
    amp = samp.val
    freq = sfreq.val
    l.set_ydata(amp*np.sin(2*np.pi*freq*t))
    fig.canvas.draw_idle()
sfreq.on_changed(update)
samp.on_changed(update)

resetax = plt.axes([0.8, 0.025, 0.1, 0.04])
button = Button(resetax, 'Reset', color=axcolor, hovercolor='0.975')


def reset(event):
    sfreq.reset()
    samp.reset()
button.on_clicked(reset)

rax = plt.axes([0.025, 0.5, 0.15, 0.15], axisbg=axcolor)
radio = RadioButtons(rax, ('red', 'blue', 'green'), active=0)


def colorfunc(label):
    l.set_color(label)
    fig.canvas.draw_idle()
radio.on_clicked(colorfunc)

plt.show()
In [32]:
np.arccos(-0.1)
Out[32]:
1.6709637479564565
In [33]:
np.arccos(0.1)
Out[33]:
1.4706289056333368
In [35]:
np.sign(-2)
Out[35]:
-1
In [ ]: