Lecture14-combined.ipynb

Kernel: Python 3 (system-wide)

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.

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

nmax = 10 # max energy level to include (counting starts at zero)

alpha = 0.3 # perturbation parameter value

# 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.

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}.$$

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

print(a) # show the matrix form of a

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

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

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

print(x)

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

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

Hprime = alpha**2/2*x**2 # perturbation to the Hamiltonian

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

print(H)

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

print(energies)

# calculate errors in the eigenvalues
errors = energies-exactEn

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

print(errors)

import matplotlib.pyplot as plt

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

stateErrors = 1.-np.abs(np.diag(states))**2

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

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

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

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

# 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
    

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

xmin, xmax, n #show the limits and number of steps

#define the x coordinates
x = np.linspace(xmin,xmax,n)

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

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

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

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

energies #display the lowest four energies

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, ϕ");