 CoCalc Public FilesSkills 2020 / Lecture14-combined.ipynb
Author: Alyssa Gadsby
Compute Environment: Ubuntu 20.04 (Default)

# Approximations worksheet, Part 5

Solving the Hamiltonian

$\hat{H} = \frac12 \hat{p}^2 + \frac12 \hat{x}^2 + \frac12 \alpha^2 \hat{x}^2$

using the truncated basis approximation, using the harmonic oscillator energy basis.

The size of the basis is set by changing the variable nmax below.

The value of $\alpha$ is set by changing the variable alpha below.

In [ ]:
# import some needed libraries
import numpy as np
import matplotlib.pyplot as plt
import scipy.optimize as opt
%matplotlib inline

In [ ]:
nmax = 10 # max energy level to include (counting starts at zero)

In [ ]:
alpha = 0.3 # perturbation parameter value

In [ ]:
# The exact eigenvalues for later comparison...
exactEn = np.sqrt(1+alpha**2)*(np.array(range(nmax+1))+0.5) # exact energies, calculated earlier by solving H by hand.

In [ ]:
print(exactEn)


Below, the lowering operator $\hat{a}$ is defined in matrix form. The code is a shortcut to generate the matrix elements

$a_{mn} = \langle m | \hat{a} | n \rangle = \sqrt{n}\, \langle m | n-1 \rangle = \sqrt{n}\, \delta_{m,n-1}.$

In [ ]:
a = np.matrix(np.diag(np.sqrt(np.array(range(nmax))+1.),k=1)) # lowering operator in matrix form

In [ ]:
print(a) # show the matrix form of a

In [ ]:
print(a.H) # the .H method of a numpy matrix is the Hermitian conjugate, what we call "dagger".

In [ ]:
a*a.H-a.H*a # just checking if the commutator rule works: [a,a.H]=1
#Should yield the identity matrix.  (Last row/column will be wrong because we are approximating.)

In [ ]:
x = (a + a.H)/np.sqrt(2.) # define the position operator x

In [ ]:
print(x)

In [ ]:
p = -1.j/np.sqrt(2)*(a-a.H) # define the momentum operator p (j = sqrt(-1))

In [ ]:
H0 = p**2/2 + x**2/2 # Unperturbed Hamiltonian ( ** means "power" in python).
# (Note * is matrix multiplication and ** is matrix power for x and p, which are np.matrix objects.)

In [ ]:
Hprime = alpha**2/2*x**2 # perturbation to the Hamiltonian

In [ ]:
H = np.real(H0 + Hprime) # full Hamiltonian (We know H is real, but python doesn't.)

In [ ]:
print(H)

In [ ]:
energies, states = np.linalg.eigh(H) # calculate eigenvalues and eigenvectors

In [ ]:
print(energies)

In [ ]:
# calculate errors in the eigenvalues
errors = energies-exactEn

In [ ]:
# show an example eigenstate vector
print(states[:,1])

In [ ]:
print(errors)

In [ ]:
import matplotlib.pyplot as plt

In [ ]:
plt.bar(range(nmax+1),exactEn,label="Exact")
plt.bar(range(nmax+1),energies,label="Approx.")
plt.xlabel("State")
plt.ylabel("Energy")
plt.legend()
plt.title("Energies")


Notice that the errors get larger with larger n, particularly for n > nmax/2

We'll define the error in the states by $\text{State error}= 1-| \langle true | approx \rangle |^2$

In [ ]:
stateErrors = 1.-np.abs(np.diag(states))**2

In [ ]:
plt.bar(range(nmax+1),stateErrors)
plt.xlabel("State")
plt.ylabel("Error")
plt.title("State error")


# Part 6, Numerical solution to the 1-dimensional Time Independent Schroedinger Equation

Based on the paper "Matrix Numerov method for solving Schroedinger's equation" by Mohandas Pillai, Joshua Goglio, and Thad G. Walker, American Journal of Physics 80 (11), 1017 (2012). doi:10.1119/1.4748813

In [ ]:
autoscale = False # set this equal to true to use Pillai's recommended step sizes

In [ ]:
# values of constants
hbar = 1.0
mass = 1.0 # changing the mass will also change the energy scale
omega = 1.0

In [ ]:
# bounds (These are overwritten if autoscale=True)
xmin = -5.0 # lower bound of position
xmax = 5.0 # upper bound of position
n = 100 # number of steps (may be overwritten if autoscale == True)
dx = (xmax-xmin)/(n-1)

In [ ]:
# the function V is the potential energy function
def V(x):
# make sure there is no division by zero
# this also needs to be a "vectorizable" function
# uncomment one of the examples below, or write your own.
return 0.5*mass*omega**2*x*x # harmonic oscillator


In [ ]:
if (autoscale):
#Emax is the maximum energy for which to check for eigenvalues
Emax = 20.0
#The next lines make some reasonable choices for the position grid size and spacing
xt = opt.brentq(lambda x: V(x)-Emax ,0,5*Emax) #classical turning point
dx = 1.0/np.sqrt(2*Emax) #step size
# bounds and number of steps
n = np.int(0.5+2*(xt/dx + 4.0*np.pi)) #number of steps
xmin = -dx*(n+1)/2
xmax = dx*(n+1)/2

In [ ]:
xmin, xmax, n #show the limits and number of steps

In [ ]:
#define the x coordinates
x = np.linspace(xmin,xmax,n)

In [ ]:
#define the numerov matrices
B = np.matrix((np.eye(n,k=-1)+10.0*np.eye(n,k=0)+np.eye(n,k=1))/12.0)
A = np.matrix((np.eye(n,k=-1)-2.0*np.eye(n,k=0)+np.eye(n,k=1))/(dx**2))

In [ ]:
#calculate kinetic energy operator using Numerov's approximation
KE = -0.5*hbar**2/mass*B.I*A

In [ ]:
#calculate hamiltonian operator approximation
H = KE + np.diag(V(x))

In [ ]:
#Calculate eigenvalues and eigenvectors of H
energies, wavefunctions = np.linalg.eigh(H) # "wavefunctions" is a matrix with one eigenvector in each column.

In [ ]:
energies #display the lowest four energies

In [ ]:
number = [0,1,2,3,4] #which wavefunctions to plot, starting counting from zero
zoom = 3.0 # zoom factor for plotting the wavefunctions to make them more visible
plt.plot(x,V(x),'-k',label="V(x)") # plot the potential
for num in number:
plt.plot(x,zoom*wavefunctions[:,num]+energies[num],label=num) #plot the num-th wavefunction
plt.hlines(energies[num],-5,5,color="black",linewidth=0.5)
plt.ylim(-1,10); # set limits of vertical axis for plot
plt.xlim(-5,5); # set limits of horizontal axis for plot
#plt.legend(loc="lower center");
plt.xlabel("x");
plt.ylabel("Energy, ϕ");

In [ ]: