``````

In [3]:

%pylab inline
from IPython.display import HTML

``````
``````

Populating the interactive namespace from numpy and matplotlib

``````
``````

In [19]:

# defining our own kronecker delta
kronecker_delta = lambda i,j: 1 if i == j else 0
potential = lambda x: -1 if absolute(x) < 20 else 0

``````

Next we generate the lattice for the simulation. For now we consider 2000 lattice points.

``````

In [20]:

delta = 1 # The spacing between neighboring lattice points
L = 1000 # The ends of the lattice
length = 10
momentum = 1
lattice = arange(-L,L,delta)
v = vectorize(potential)
plot(v(lattice))

``````
``````

Out[20]:

[<matplotlib.lines.Line2D at 0x284937532b0>]

``````
``````

In [23]:

hamiltonian = np.zeros((lattice.shape[0],lattice.shape[0]))
for row in range(lattice.shape[0]):
for col in range(lattice.shape[0]):
hamiltonian[row,col] = (kronecker_delta(col+1, row)\
+ kronecker_delta(col-1, row)\
- 2*kronecker_delta(col,row))/-(delta)**2\
+ potential(row)*kronecker_delta(col,row)
hamiltonian

``````
``````

Out[23]:

array([[ 1., -1.,  0., ...,  0.,  0.,  0.],
[-1.,  1., -1., ...,  0.,  0.,  0.],
[ 0., -1.,  1., ...,  0.,  0.,  0.],
...,
[ 0.,  0.,  0., ...,  2., -1.,  0.],
[ 0.,  0.,  0., ..., -1.,  2., -1.],
[ 0.,  0.,  0., ...,  0., -1.,  2.]])

``````

Now that we've computed the hamiltonian matrix, let's diagonalize it

``````

In [26]:

eigenvalues, eigenvectors = linalg.eigh(hamiltonian)
plot(lattice,absolute(eigenvectors[0]))

``````
``````

Out[26]:

[<matplotlib.lines.Line2D at 0x2849575ff28>]

``````
``````

In [60]:

# eigenvectors[0].shape
wavefunction.shape

``````
``````

Out[60]:

(2000,)

``````

Let's define the initial state of our system now.

``````

In [56]:

wavefunction = (1/sqrt(2)*length)*exp(-((lattice)**2)/2000)*exp(1j*momentum*lattice)
normalize = sum(absolute(wavefunction))*delta
wavefunction /= normalize
print(sum(absolute(wavefunction))*delta)
plot(lattice,absolute(wavefunction))
wavefunction

``````
``````

1.0

Out[56]:

array([  5.05473356e-220 -7.43209685e-220j,
2.44114938e-219 +6.46172756e-221j,
3.43234625e-219 +5.67016730e-219j, ...,
-7.90882423e-219 -1.61383394e-218j,
3.43234625e-219 -5.67016730e-219j,
2.44114938e-219 -6.46172756e-221j])

``````

Now we decompose the state into the projections on the basis vectors

``````

In [48]:

coeff = empty_like(eigenvalues)
for index in range(eigenvalues.shape[0]):
coeff[index] = vdot(eigenvectors[index],wavefunction)

``````
``````

C:\Anaconda\lib\site-packages\ipykernel\__main__.py:3: ComplexWarning: Casting complex values to real discards the imaginary part
app.launch_new_instance()

``````

Now that we know the wavefunction as a linear combination of basis vectors, the time evolution of the complete wavefuction is just the linear combination of the time evolutions of the basis states. Lets do this for 10 time steps.

``````

In [40]:

func = zeros_like(eigenvalues,dtype=complex128)
for index in range(eigenvalues.shape[0]):
func += coeff[index]*eigenvectors[index]
plot(lattice,absolute(func))

``````
``````

Out[40]:

[<matplotlib.lines.Line2D at 0x29ea1225710>]

``````
``````

In [41]:

time_steps = arange(0,200,2)
evolution = []
for time in time_steps:
wavefunc = zeros_like(eigenvalues,dtype=complex128)
for index in range(eigenvalues.shape[0]):
wavefunc += coeff[index]*eigenvectors[index]*exp(-1j*eigenvalues[index]*time)
evolution.append(wavefunc)

``````
``````

In [44]:

fig = plt.figure()
ax = plt.axes(xlim=(-L, L), ylim=(0, 0.2))
line, = ax.plot([], [], lw=2)
def init():
line.set_data([], [])
return line,

def animate(i):
line.set_data(lattice,absolute(evolution[i]))
return line,

``````
``````

``````
``````

In [45]:

from matplotlib import animation
anim = animation.FuncAnimation(fig,animate,init_func=init,frames=100,interval=50,blit=True)
HTML(anim.to_html5_video())

``````
``````

Out[45]:

Your browser does not support the video tag.

``````
``````

In [ ]:

``````