CoCalc Public FilesSkills 2020 / Lecture14-combined.ipynbOpen with one click!
Author: Alyssa Gadsby
Compute Environment: Ubuntu 20.04 (Default)

Approximations worksheet, Part 5

Solving the Hamiltonian

H^=12p^2+12x^2+12α2x^2 \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 a^ \hat{a} is defined in matrix form. The code is a shortcut to generate the matrix elements

amn=ma^n=nmn1=nδm,n1. 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 State error=1trueapprox2 \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 [ ]: