This notebook contains animations of 2D wavefunctions evolving through time, generated using the Crank-Nicolson solver discussed in the report.
This is an IPython notebook, containing Python code (which is evaluated on my machine) and output (which is embedded into the document, for all to see). That said, some modules need to be loaded before we may begin, so please excuse the following bit of setup code.
In [1]:
# Adds interactive figures
%pylab inline
# External modules
import numpy as np
import matplotlib.pyplot as plt
# Local modules
from crank import *
from plotutil import *
from util import *
from animate import *
# http://github.com/jakevdp/JSAnimation
from JSAnimation.IPython_display import display_animation
Let's obtain Crank-Nicolson solvers for the 2D well, and for a free particle with periodic boundary conditions.
In [2]:
n = 100 # Number of spatial positions along each axis
dt = 0.001 # Time progressed by a single iteration of Crank-Nicolson
# Seed random number gen so that same results are produced
# otherwise I imagine the git repo is gonna balloon
random.seed(100)
# An empty square well potential
potEmpty = np.zeros((n,n))
# Parameters for animations later in the document
NFRAMES = 100
FRAMESKIP = 3
# Obtain a Crank Nicolson evolution functor
evoEmptyWell = getCrankNicolEvo(potEmpty, BoundaryType.REFLECTING, dt)
evoEmptySpace = getCrankNicolEvo(potEmpty, BoundaryType.PERIODIC, dt)
First, let's try a couple of test cases, to verify that the Crank-Nicolson solver is working as intended.
The ground state of the 2D square well is
$$ \Psi = \frac{2}{L}\sin\left(\frac{\pi x}{L}\right)\sin\left(\frac{\pi y}{L}\right). $$Because any eigenfunction of the hamiltonian $\hat{H}$ will also be an eigenfunction of the time evolution operator $\hat{U}(t)\equiv e^{-i\hat{H}t/\hbar}$, one would expect that the probability distribution of the wavefunction will remain constant over time.
So here's a simple test: Does our Crank Nicolson evolution function leave the probability distribution of the ground state unchanged?
In [3]:
# Construct the 2D well ground state
def twoDeeGroundState(n,m):
ys = np.linspace(0., np.pi, n+2)[1:-1] # Shifted because the boundary conditions
xs = np.linspace(0., np.pi, m+2)[1:-1] # are not a part of our solution vector
psigrid = np.outer(np.sin(ys), np.sin(xs))
return normalize(psigrid.flatten())
psi = twoDeeGroundState(n,n)
evo = evoEmptyWell
anim = timeEvoAnimation((psi,), evo, (ProbPlotter, PhasePlotter), 100,5)
display_animation(anim, default_mode='once')
Out[3]:
Indeed it does! Only the phase changes over time.
We can test this against other eigenvectors as well. Here are two eigenvectors produced by an eigenvector solver, $\psi_2$ and $\psi_4$, together with their sum, $\frac{1}{\sqrt{2}}\left(\psi_2+\psi_4\right)$. I chose the indices 2 and 4 to ensure they have different energies; for the 2D well, there are at most two degenerate eigenvectors for any given energy level (and getEigens returns vectors sorted by energy).
In [4]:
epsis, evals = getEigens(potEmpty, BoundaryType.REFLECTING, 8)
# Take two eigenfunctions of different energy
psis = (epsis[2], epsis[4], normalize(epsis[2]+epsis[4]))
evo = evoEmptyWell
# Compute how long a single period should take
period = float(2 * np.pi * evo.hbar / (evals[4] - evals[2]))
frames = int(period / dt)
anim = timeEvoAnimation(psis, evo, (ProbPlotter, PhasePlotter), frames//5, frameStep=5, flip=True)
display_animation(anim, default_mode="loop")
Out[4]:
As seen above, the two eigenfunctions (displayed in the first two columns) remain constant. Their sum in the third column, on the other hand, changes over time.
The animation also loops smoothly. Actually, this verifies another property of our Crank-Nicolson evolution function—namely, that it correctly evolves by the specified amount of time, $dt$. If you look again at the fragment of code above, I have computed how many frames a single loop of the animation should take, and passed it in as an argument to the animation. This calculation is based on the condition that, after a time equal to the period $T$, $\psi_2$ and $\psi_4$ should have the same relative phase: \begin{eqnarray} e^{iE_2 T/\hbar}=e^{iE_4 T/\hbar} \\ \implies T=\frac{2 \pi \hbar}{(E_4-E_2)} \end{eqnarray}
For periodic boundary conditions, the simplest test of our Crank-Nicolson solver is likely the plane wave:
$$ \psi(x,y)=e^{iky}, $$where $k$ is quantized such that $e^{ik/L}=1$ to meet the boundary condition of periodicity. This is an easy test to use because the plane wave is an eigenfunction of the Hamiltonian for free space. Furthermore, as with the animation produced above, we can easily compute the period $T$ of a plane wave's evolution from its energy:
\begin{align} E&=\frac{\hbar^2 k^2}{2m}\\ T&=\frac{2\pi\hbar}{E}=\frac{4 \pi m}{\hbar k^2} \end{align}
In [5]:
# Produce a plane wave exp(i*k*y)
def twoDeeVerticalPlaneWave(n,m,k):
psigrid = np.outer(np.exp(1j*k*np.linspace(0.,1.,n)), np.ones(m))
return normalize(psigrid.flatten())
k = -2*np.pi
psi = twoDeeVerticalPlaneWave(n,n,k)
evo = evoEmptySpace
# Compute how long a single period should take
period = float(4*np.pi * evo.mass / evo.hbar / k**2)
frames = int(period / dt)
anim = timeEvoAnimation((psi,), evo, (PhasePlotter,), frames//5, frameStep=5)
display_animation(anim, default_mode='loop')
Out[5]:
I have left out the plot of probability above, as a plane wave has constant amplitude at all points in space.
Now I will try adding boxes of non-zero potential. Here are some potentials containing a square of non-zero potential (either positive or negative) at the center of our 2D space.
In [6]:
halfwidth = 0.1
box = Box((0.5-halfwidth,0.5-halfwidth), (0.5+halfwidth, 0.5+halfwidth))
potCenterNeg = potentialFromBoxes((n,n), (box,), (-100.,))
potCenterPos = potentialFromBoxes((n,n), (box,), (+100.,))
Let's generate Crank-Nicolson solvers for these potentials:
In [7]:
evoCenterNegWell = getCrankNicolEvo(potCenterNeg, BoundaryType.REFLECTING, dt)
evoCenterPosWell = getCrankNicolEvo(potCenterPos, BoundaryType.REFLECTING, dt)
evoCenterNegPBC = getCrankNicolEvo(potCenterNeg, BoundaryType.PERIODIC, dt)
evoCenterPosPBC = getCrankNicolEvo(potCenterPos, BoundaryType.PERIODIC, dt)
How does this tiny square of potential disturb the ground state?
2D Well ground state, with positive square
In [8]:
psi = twoDeeGroundState(n,n)
evo = evoCenterPosWell
anim = timeEvoAnimation((psi,), evo, (ProbPlotter, PhasePlotter), NFRAMES, frameStep=FRAMESKIP)
display_animation(anim, default_mode='once')
Out[8]:
Here, we can see the little potential square expels the particle out of it.
2D Well ground state, with negative square
In [9]:
psi = twoDeeGroundState(n,n)
evo = evoCenterNegWell
anim = timeEvoAnimation((psi,), evo, (ProbPlotter, PhasePlotter), NFRAMES, frameStep=FRAMESKIP)
display_animation(anim, default_mode='once')
Out[9]:
As one might expect, the negative potential square pulls the particle in.
The effects of these potentials on the plane wave in PBC are unusual.
Plane wave, with positive square
In [10]:
psi = twoDeeVerticalPlaneWave(n,n,-2*np.pi)
evo = evoCenterPosPBC
anim = timeEvoAnimation((psi,), evo, (ProbPlotter, PhasePlotter), NFRAMES, frameStep=FRAMESKIP)
display_animation(anim, default_mode='once')
Out[10]:
Plane wave, with negative square
In [11]:
psi = twoDeeVerticalPlaneWave(n,n,-2*np.pi)
evo = evoCenterNegPBC
anim = timeEvoAnimation((psi,), evo, (ProbPlotter, PhasePlotter), NFRAMES, frameStep=FRAMESKIP)
display_animation(anim, default_mode='once')
Out[11]:
And now for something entirely different
(no, not really, just more boxes)
In [12]:
# Random boxes
nboxes = 5 # number of boxes to make
maxpot = 200 # max potential value
boxes = makeRandomBoxes(nboxes,(0.1,0.3))
potvals = [(random.random()-0.5)*2*maxpot for _ in range(nboxes)]
psi = psis[2] # The sum from before
pot = potentialFromBoxes((n,n), boxes, potvals)
evo = getCrankNicolEvo(pot, BoundaryType.REFLECTING, dt)
anim = timeEvoAnimation((psi,), evo, (ProbPlotter, PhasePlotter), NFRAMES, frameStep=FRAMESKIP)
display_animation(anim, default_mode='once')
Out[12]:
In [12]: