# 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

%matplotlib inline


## Harmonic Oscillator¶

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

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)

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
^
SyntaxError: invalid syntax

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

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))

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))

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¶

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

In [55]:
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*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()
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)

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

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 [ ]: