Reminder: PIB energies given by $$ E_n = n^2\frac{h^2}{8mL^2}. $$ I work with dimensionless quantities: \begin{eqnarray*} x^\prime &\equiv& \frac{x}{L},\\ E^\prime &\equiv& \frac{E}{\frac{h^2}{8mL^2}},\\ t^\prime &\equiv& \frac{t}{\frac{2mL^2}{h}}. \end{eqnarray*} The time is scaled to be the time for a classical particle with momentum $p = h/2L$ to travel a distance $L$.
In these units $$e^{-iEt/\hbar} \rightarrow e^{-i2\pi E^\prime t^\prime},$$ effectively setting $h \rightarrow 1$.
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import matplotlib.animation as animation
# Following is an Ipython magic command that puts figures in notebook.
%matplotlib notebook
# M.L. modification of matplotlib defaults
# Changes can also be put in matplotlibrc file,
# or effected using mpl.rcParams[]
plt.style.use('classic')
plt.rc('figure', figsize = (6, 4.5)) # Reduces overall size of figures
plt.rc('axes', labelsize=16, titlesize=14)
plt.rc('xtick', labelsize=12)
plt.rc('ytick', labelsize=12)
plt.rc('figure', autolayout = True) # Adjusts supblot params for new size
def phiPIB(n,x):
'''Particle-in-a-box wavefunction'''
return np.sqrt(2.)*np.sin(n*np.pi*x)
def psiTotal(x,t):
'''Linear combination of PIB wave functions'''
return np.exp(-1j*2*np.pi*1*t)*phiPIB(1,x)/np.sqrt(2.) \
+ np.exp(-1j*2*np.pi*4.*t)*phiPIB(2,x)/np.sqrt(2.)
def animate(i): # Update time, title (containing time), and data
dt = 0.01
t0 = 0
t = t0 + i*dt
plt.title("t = {0:.2f}".format(t))
P.set_ydata(abs(psiTotal(x,t))**2) # update the data
return P,
fig = plt.figure()
N = 50
x = np.linspace(0, 1, N, endpoint=True)
plt.grid()
plt.ylim(-.5,4)
plt.xlabel('$x$') # Label for horizontal axis
plt.ylabel('$|\psi|^2$') # Label for vertical axis
plt.axhline(0,color='green') # Makes solid green x-axis
y = psiTotal(x,0)
P, = plt.plot(x, abs(y)**2)
ani = animation.FuncAnimation(fig, animate, interval = 200) # interval --> time between frames
For the linear combination of $n=1$ and $n=2$ states, I get 3 periods in a time of 1. Compare this with "prediction":
The frequency is \begin{eqnarray*} f = \frac{\Delta E}{h} &=& \frac{2^2 - 1^1}{1}\\ &=& 3 \end{eqnarray*} which gives the period $$ T = \frac{1}{3}. $$
%version_information
is an IPython magic extension for showing version information for dependency modules in a notebook;
%version_information
is available on Bucknell computers on the linux network. You can easily install it on any computer.
%load_ext version_information
version_information numpy, matplotlib