Now let's try a more significant change to the potential. It will be instructive to place a quantum harmonic oscillator potential into our square well, and then build the Hamiltonian, again using the square well eigenvectors as a basis set. We will look at two cases: first, a small magnitude for the potential (a small perturbation, where we would expect to see small changes in the square well energies); second, a large magnitude for the potential (where the first few eigenstates should be very close to the QHO eigenstates).
First we set up the functions that we need (you should be familiar with these by now !).
# Import libraries and set up in-line plotting.
%matplotlib inline
import matplotlib.pyplot as pl
import numpy as np
# This is a new library - linear algebra includes solving for eigenvalues & eigenvectors of matrices
import numpy.linalg as la
# Define the eigenbasis - normalisation needed elsewhere
def square_well_eigenfunctions(n,width,norm,x):
"""The eigenbasis for a square well, running from 0 to a (width), sin(n pi x/a).
N.B. requires a normalisation factor, norm."""
wavevector = np.pi*n/width
return norm*np.sin(wavevector*x)
# We will also define the second derivative for kinetic energy (KE)
def second_derivative_square_well_eigenfunctions(n,width,norm,x):
"""The second derivative of the eigenbasis for a square well, running from 0 to a, sin(n pi x/a)"""
wavevector = np.pi*n/width
return -wavevector*wavevector*norm*np.sin(wavevector*x)
# Define the x-axis
width = 1.0
num_x_points = 101
x = np.linspace(0.0,width,num_x_points)
dx = width/(num_x_points - 1)
# Integrate the product of two functions over the width of the well
# NB this is a VERY simple integration routine: there are much better ways
def integrate_functions(function1,function2,dx):
"""Integrate the product of two functions over defined x range with spacing dx"""
# We use the NumPy dot function here instead of iterating over array elements
integral = dx*np.dot(function1,function2)
return integral
# Now set up the array of basis functions - specify the size of the basis
num_basis = 10
# These arrays will each hold an array of functions
basis_functions_array = np.zeros((num_basis,num_x_points))
second_derivative_basis_functions_array = np.zeros((num_basis,num_x_points))
# Loop over first num_basis basis states, normalise and create an array
# NB the basis_functions_array will start from 0
for i in range(num_basis):
n = i+1
# Calculate A = <phi_n|phi_n>
integral = integrate_functions(square_well_eigenfunctions(n,width,1.0,x),square_well_eigenfunctions(n,width,1.0,x),dx)
# Use 1/sqrt{A} as normalisation constant
normalisation = 1.0/np.sqrt(integral)
basis_functions_array[i,:] = square_well_eigenfunctions(n,width,normalisation,x)
second_derivative_basis_functions_array[i,:] = second_derivative_square_well_eigenfunctions(n,width,normalisation,x)
# Define a function to act on a basis function with the potential
def add_potential_on_basis(Hphi,V,phi):
for i in range(V.size):
Hphi[i] = Hphi[i] + V[i]*phi[i]
//anaconda/lib/python2.7/site-packages/matplotlib/font_manager.py:273: UserWarning: Matplotlib is building the font cache using fc-list. This may take a moment. warnings.warn('Matplotlib is building the font cache using fc-list. This may take a moment.')
As before, we'll start with a small perturbation (though you should think carefully about the following question: what does small mean in this context ?)
# Define the potential in the square well
def square_well_and_QHO_potential(x,V,a):
"""QHO Potential for a particle in a square well, expecting two arrays: x, V(x), and potential height, a"""
# Note that we offset x so that the QHO well is centred on the square well
for i in range(x.size):
V[i] = 0.5*a*(x[i]-0.5)**2
# Plot to ensure that we know what we're getting
pl.plot(x,V)
omega = 1.0
omega_2 = omega**2
Potential_with_QHO = np.linspace(0.0,width,num_x_points)
potential_bump = square_well_and_QHO_potential(x,Potential_with_QHO,omega_2)
Now that we have the potential defined, we will make matrix elements (again, outputting the potential matrix so that we can think about the perturbation this gives), diagonalise and output the first few eigenvalues and eigenvectors. We'll compare these to two limits: the full QHO solution; and the square well solution.
Note that we have defined $\hbar = m_e = 1$, so the QHO eigenenergies are just:
$$ E_n = \omega\left(n+\frac{1}{2}\right) $$The contribution to the eigenenergies from the square well is just:
$$\hat{H} = -\frac{1}{2} \frac{d^{2}}{dx^{2}}$$# Declare space for the matrix elements
H_matrix = np.eye(num_basis)
# Loop over basis functions phi_n (the bra in the matrix element)
print "Full Hamiltonian"
for n in range(num_basis):
# Loop over basis functions phi_m (the ket in the matrix element)
for m in range(num_basis):
# Act with H on phi_m and store in H_phi_m
H_phi_m = -0.5*second_derivative_basis_functions_array[m]
add_potential_on_basis(H_phi_m,Potential_with_QHO,basis_functions_array[m])
# Create matrix element by integrating
H_matrix[m,n] = integrate_functions(basis_functions_array[n],H_phi_m,dx)
# The comma at the end prints without a new line; the %8.3f formats the number
print "%8.3f" % H_matrix[m,n],
# This print puts in a new line when we have finished looping over m
print
print "Perturbation matrix elements:"
# Output the matrix elements of the potential to see how large the perturbation is
# Loop over basis functions phi_n (the bra in the matrix element)
for n in range(num_basis):
# Loop over basis functions phi_m (the ket in the matrix element)
for m in range(num_basis):
# Act with H on phi_m and store in H_phi_m
H_phi_m = np.zeros(num_x_points)
add_potential_on_basis(H_phi_m,Potential_with_QHO,basis_functions_array[m])
# Create matrix element by integrating
H_mn = integrate_functions(basis_functions_array[n],H_phi_m,dx)
# The comma at the end prints without a new line; the %8.3f formats the number
print "%8.3f" % H_mn,
# This print puts in a new line when we have finished looping over m
print
# Solve using linalg module of numpy (which we've imported as la above)
eigenvalues, eigenvectors = la.eigh(H_matrix)
# This call above does the entire solution for the eigenvalues and eigenvectors !
# Print results roughly, though apply precision of 4 to the printing
print
print "Eigenvalues and eigenvector coefficients printed roughly"
np.set_printoptions(precision=4)
print eigenvalues
print eigenvectors[0]
print eigenvectors[1]
print eigenvectors[2]
print
print " QHO Square Perfect QHO Difference Perfect Square Difference"
for i in range(num_basis):
n = i+1
print " %8.3f %8.3f %8.3f %8.3f %8.3f" % (eigenvalues[i],omega*(i+0.5),eigenvalues[i] - omega*(i+0.5),n*n*np.pi*np.pi/2.0,eigenvalues[i] - n*n*np.pi*np.pi/2.0)
Full Hamiltonian 4.951 -0.000 0.019 -0.000 0.004 0.000 0.001 -0.000 0.001 -0.000 0.000 19.775 -0.000 0.023 0.000 0.005 0.000 0.002 -0.000 0.001 0.019 0.000 44.452 -0.000 0.024 0.000 0.005 -0.000 0.002 0.000 0.000 0.023 0.000 78.997 0.000 0.024 -0.000 0.006 -0.000 0.002 0.004 0.000 0.024 0.000 123.411 -0.000 0.025 0.000 0.006 -0.000 0.000 0.005 0.000 0.024 -0.000 177.694 0.000 0.025 0.000 0.006 0.001 0.000 0.005 -0.000 0.025 0.000 241.846 -0.000 0.025 -0.000 0.000 0.002 -0.000 0.006 0.000 0.025 -0.000 315.869 -0.000 0.025 0.001 -0.000 0.002 -0.000 0.006 0.000 0.025 -0.000 399.760 -0.000 -0.000 0.001 0.000 0.002 -0.000 0.006 -0.000 0.025 0.000 493.522 Perturbation matrix elements: 0.016 -0.000 0.019 0.000 0.004 0.000 0.001 -0.000 0.001 0.000 -0.000 0.035 -0.000 0.023 0.000 0.005 -0.000 0.002 -0.000 0.001 0.019 -0.000 0.039 -0.000 0.024 0.000 0.005 0.000 0.002 0.000 0.000 0.023 -0.000 0.040 -0.000 0.024 -0.000 0.006 -0.000 0.002 0.004 0.000 0.024 -0.000 0.041 -0.000 0.025 0.000 0.006 -0.000 0.000 0.005 0.000 0.024 -0.000 0.041 -0.000 0.025 -0.000 0.006 0.001 -0.000 0.005 -0.000 0.025 -0.000 0.041 -0.000 0.025 -0.000 -0.000 0.002 0.000 0.006 0.000 0.025 -0.000 0.041 -0.000 0.025 0.001 -0.000 0.002 -0.000 0.006 -0.000 0.025 -0.000 0.041 -0.000 0.000 0.001 0.000 0.002 -0.000 0.006 -0.000 0.025 -0.000 0.041 Eigenvalues and eigenvector coefficients printed roughly [ 4.9511 19.7745 44.4521 78.9969 123.4107 177.6938 241.8465 315.8686 399.7603 493.5216] [ 1.0000e+00 4.9433e-17 4.8092e-04 -2.2455e-17 2.9745e-05 2.3690e-17 -5.2026e-06 -1.6900e-17 -1.4445e-06 -5.6894e-17] [ 4.9440e-17 -1.0000e+00 0.0000e+00 3.8017e-04 -2.2204e-16 3.0108e-05 5.5511e-17 6.0876e-06 2.3592e-16 1.8574e-06] [ -4.8092e-04 -1.1102e-16 1.0000e+00 9.6711e-17 3.0075e-04 1.2229e-16 -2.6972e-05 -1.1600e-16 -5.9446e-06 2.5646e-16] QHO Square Perfect QHO Difference Perfect Square Difference 4.951 0.500 4.451 4.935 0.016 19.775 1.500 18.275 19.739 0.035 44.452 2.500 41.952 44.413 0.039 78.997 3.500 75.497 78.957 0.040 123.411 4.500 118.911 123.370 0.041 177.694 5.500 172.194 177.653 0.041 241.846 6.500 235.346 241.805 0.041 315.869 7.500 308.369 315.827 0.041 399.760 8.500 391.260 399.719 0.041 493.522 9.500 484.022 493.480 0.041
Notice how the eigenvalues are almost the same as for the square well. More interestingly, the change in the eigenvalues is given by the diagonal elements of the matrix of the potential - which is exactly what we would expect from first order perturbation theory. You could explore how well this holds, by gradually increasing the size of the potential - we will try larger values below (and we could also code second order changes in the energy for greater accuracy, but we won't for now).
Now we will plot the eigenvectors, and the difference between the eigenvectors with this small perturbation, and the unperturbed eigenvectors (the basis functions).
from scipy.special import hermite
from scipy.misc import factorial
from math import pi
root_pi = np.sqrt(pi)
def N(n, alpha):
return np.sqrt(alpha / (root_pi * (2.0**n) * factorial(n)))
def phi(x,n,alpha):
return N(n,alpha) * hermite(n)(alpha * x) * np.exp(-0.5 * alpha**2 * x**2)
x2 = np.linspace(-width/2,width/2,num_x_points)
#pl.plot(x,phi(x2,1,np.sqrt(np.sqrt(omega))))
# Define a figure to take two plots
fig3 = pl.figure(figsize=[12,3])
# Add subplots: number in y, x, index number
axb = fig3.add_subplot(121,autoscale_on=False,xlim=(0,1),ylim=(-2.1,2.1))
axb.set_title("Eigenvectors for perturbed system")
axb2 = fig3.add_subplot(122,autoscale_on=False,xlim=(0,1),ylim=(-0.001,0.001))
axb2.set_title("Difference to perfect square well eigenvectors")
#axb2.set_title("Difference to QHO eigenvectors")
for m in range(3): # Plot the first four states
psi = np.zeros(num_x_points)
for i in range(num_basis):
psi = psi+eigenvectors[i,m]*basis_functions_array[i]
if 2*(m/2)!=m: # This is just to ensure that psi and the basis function have the same phase
psi = -psi
axb.plot(x,psi)
psi = psi - basis_functions_array[m]
axb2.plot(x,psi)
As we would expect, the eigenvectors are almost unchanged.
Now we will go beyond a small perturbation, by making the QHO potential much stronger. We have to be a little careful how we do this: we'll see the effects of making the potential too strong lower down.
We'll also plot the eigenvectors, this time comparing them to the full QHO eigenvectors.
# Make omega larger, so that the QHO energy dominates the square well
omegaLarger = 50.0
omegaLarger2 = omegaLarger**2
Potential_with_QHO2 = np.linspace(0.0,width,num_x_points)
potential_bump = square_well_and_QHO_potential(x,Potential_with_QHO2,omegaLarger2)
# Declare space for the matrix elements
H_matrix3 = np.eye(num_basis)
# Loop over basis functions phi_n (the bra in the matrix element)
for n in range(num_basis):
# Loop over basis functions phi_m (the ket in the matrix element)
for m in range(num_basis):
# Act with H on phi_m and store in H_phi_m
H_phi_m = -0.5*second_derivative_basis_functions_array[m]
add_potential_on_basis(H_phi_m,Potential_with_QHO2,basis_functions_array[m])
# Create matrix element by integrating
H_matrix3[m,n] = integrate_functions(basis_functions_array[n],H_phi_m,dx)
# The comma at the end prints without a new line; the %8.3f formats the number
print "%8.3f" % H_matrix3[m,n],
# This print puts in a new line when we have finished looping over m
print
print "Perturbation matrix elements:"
# Loop over basis functions phi_n (the bra in the matrix element)
for n in range(num_basis):
# Loop over basis functions phi_m (the ket in the matrix element)
for m in range(num_basis):
# Act with H on phi_m and store in H_phi_m
H_phi_m = np.zeros(num_x_points)
add_potential_on_basis(H_phi_m,Potential_with_QHO2,basis_functions_array[m])
# Create matrix element by integrating
H_mn = integrate_functions(basis_functions_array[n],H_phi_m,dx)
# The comma at the end prints without a new line; the %8.3f formats the number
print "%8.3f" % H_mn,
# This print puts in a new line when we have finished looping over m
print
# Solve using linalg module of numpy (which we've imported as la above)
eigenvalues, eigenvectors = la.eigh(H_matrix3)
# This call above does the entire solution for the eigenvalues and eigenvectors !
# Print results roughly, though apply precision of 4 to the printing
print
print "Eigenvalues and eigenvector coefficients printed roughly"
np.set_printoptions(precision=4)
print eigenvalues
print eigenvectors[0]
print eigenvectors[1]
print eigenvectors[2]
print
print " QHO Square Perfect QHO Difference Perfect Square Difference"
for i in range(num_basis):
n = i+1
print " %8.3f %8.3f %8.3f %8.3f %8.3f" % (eigenvalues[i],omegaLarger*(i+0.5),eigenvalues[i] - omegaLarger*(i+0.5),n*n*np.pi*np.pi/2.0,eigenvalues[i] - n*n*np.pi*np.pi/2.0)
45.776 -0.000 47.494 -0.000 8.795 0.000 3.078 -0.000 1.425 -0.000 -0.000 108.074 -0.000 56.290 0.000 11.874 0.000 4.503 -0.000 2.199 47.494 -0.000 141.544 0.000 59.368 0.000 13.298 -0.000 5.277 0.000 0.000 56.290 0.000 179.166 0.000 60.793 -0.000 14.072 -0.000 5.744 8.795 -0.000 59.368 -0.000 225.004 -0.000 61.567 0.000 14.539 -0.000 0.000 11.874 0.000 60.793 -0.000 280.060 -0.000 62.033 0.000 14.842 3.078 -0.000 13.298 -0.000 61.567 -0.000 344.679 -0.000 62.336 -0.000 -0.000 4.503 -0.000 14.072 0.000 62.033 -0.000 419.004 -0.000 62.544 1.425 -0.000 5.277 -0.000 14.539 -0.000 62.336 -0.000 503.104 -0.000 0.000 2.199 0.000 5.744 -0.000 14.842 -0.000 62.544 -0.000 597.013 Perturbation matrix elements: 40.841 -0.000 47.494 0.000 8.795 0.000 3.078 -0.000 1.425 0.000 -0.000 88.335 -0.000 56.290 0.000 11.874 -0.000 4.503 -0.000 2.199 47.494 -0.000 97.130 -0.000 59.368 0.000 13.298 0.000 5.277 0.000 0.000 56.290 -0.000 100.209 -0.000 60.793 -0.000 14.072 -0.000 5.744 8.795 0.000 59.368 -0.000 101.634 0.000 61.567 0.000 14.539 -0.000 0.000 11.874 0.000 60.793 0.000 102.407 -0.000 62.033 -0.000 14.842 3.078 -0.000 13.298 -0.000 61.567 -0.000 102.874 -0.000 62.336 -0.000 -0.000 4.503 -0.000 14.072 0.000 62.033 -0.000 103.177 -0.000 62.544 1.425 -0.000 5.277 -0.000 14.539 -0.000 62.336 -0.000 103.385 -0.000 0.000 2.199 0.000 5.744 -0.000 14.842 -0.000 62.544 -0.000 103.533 Eigenvalues and eigenvector coefficients printed roughly [ 25.0007 75.0161 125.1641 175.9898 228.9861 286.5714 351.5116 424.9498 529.4435 620.7907] [ 9.0718e-01 9.9888e-16 3.8803e-01 -1.5581e-16 1.5570e-01 8.0255e-17 -4.5610e-02 4.1496e-17 -1.1746e-02 -5.1886e-17] [ 6.1062e-16 -8.4746e-01 1.1102e-15 4.9390e-01 -2.2204e-16 1.8744e-01 -6.1062e-16 5.0480e-02 3.7470e-16 1.3696e-02] [ -4.1202e-01 5.5511e-16 7.4170e-01 -1.6653e-16 4.9802e-01 2.7756e-16 -1.7329e-01 6.9389e-18 -4.5479e-02 1.5092e-16] QHO Square Perfect QHO Difference Perfect Square Difference 25.001 25.000 0.001 4.935 20.066 75.016 75.000 0.016 19.739 55.277 125.164 125.000 0.164 44.413 80.751 175.990 175.000 0.990 78.957 97.033 228.986 225.000 3.986 123.370 105.616 286.571 275.000 11.571 177.653 108.919 351.512 325.000 26.512 241.805 109.706 424.950 375.000 49.950 315.827 109.122 529.443 425.000 104.443 399.719 129.724 620.791 475.000 145.791 493.480 127.310
# Define a figure to take two plots
fig4 = pl.figure(figsize=[12,3])
# Add subplots: number in y, x, index number
axc = fig4.add_subplot(121,autoscale_on=False,xlim=(0,1),ylim=(-2.1,2.1))
axc.set_title("Eigenvectors for perturbed system")
axc2 = fig4.add_subplot(122,autoscale_on=False,xlim=(0,1),ylim=(-0.1,0.1))
#axc2.set_title("QHO eigenvectors")
axc2.set_title("Difference to QHO eigenvectors")
for m in range(3): # Plot the first four states
psi = np.zeros(num_x_points)
for i in range(num_basis):
psi = psi+eigenvectors[i,m]*basis_functions_array[i]
#if 2*(m/2)!=m: # This is just to ensure that psi and the basis function have the same phase
# psi = -psi
axc.plot(x,psi)
psi = psi - phi(x2,m,np.sqrt(omegaLarger))
axc2.plot(x,psi)
Now we can see two things:
This second result comes about when the perturbation (the QHO well) becomes comparable to the square well energy. If we went higher in energy, we'd gradually recover the square well eigenvalues.
We'll now see the problem with using a finite basis set (we defined num_basis
, the number of basis functions to include, as 10 above, to make matrices a sensible size). If we make the QHO potential huge compared to the square well eigenvalues, we'll find that the resulting eigenfunctions are very compressed, and that our ten square well eigenstates are not enough to properly represent them.
The calculation takes the usual form: define a potential; form a Hamiltonian matrix (we won't print the potential matrix this time); diagonalise; analyse the results.
# Make omega larger, so that the QHO energy dominates the square well
omegaH = 500.0
omegaH2 = omegaH**2
Potential_with_QHO3 = np.linspace(0.0,width,num_x_points)
potential_bump = square_well_and_QHO_potential(x,Potential_with_QHO3,omegaH2)
# Declare space for the matrix elements
H_matrix4 = np.eye(num_basis)
# Loop over basis functions phi_n (the bra in the matrix element)
for n in range(num_basis):
# Loop over basis functions phi_m (the ket in the matrix element)
for m in range(num_basis):
# Act with H on phi_m and store in H_phi_m
H_phi_m = -0.5*second_derivative_basis_functions_array[m]
add_potential_on_basis(H_phi_m,Potential_with_QHO3,basis_functions_array[m])
# Create matrix element by integrating
H_matrix4[m,n] = integrate_functions(basis_functions_array[n],H_phi_m,dx)
# The comma at the end prints without a new line; the %8.3f formats the number
print "%8.3f" % H_matrix4[m,n],
# This print puts in a new line when we have finished looping over m
print
# Solve using linalg module of numpy (which we've imported as la above)
eigenvalues, eigenvectors = la.eigh(H_matrix4)
print " QHO Square Perfect QHO Difference"
for i in range(num_basis):
n = i+1
print " %8.3f %8.3f %8.3f" % (eigenvalues[i],omegaH*(i+0.5),eigenvalues[i] - omegaH*(i+0.5))
4089.027 -0.000 4749.429 0.000 879.522 0.000 307.831 -0.000 142.479 0.000 -0.000 8853.261 -0.000 5628.951 -0.000 1187.353 -0.000 450.310 -0.000 219.873 4749.429 -0.000 9757.457 -0.000 5936.782 0.000 1329.832 0.000 527.703 0.000 0.000 5628.951 -0.000 10099.831 -0.000 6079.261 -0.000 1407.225 -0.000 574.367 879.522 0.000 5936.782 -0.000 10286.723 0.000 6156.655 0.000 1453.889 -0.000 0.000 1187.353 0.000 6079.261 0.000 10418.400 -0.000 6203.318 -0.000 1484.172 307.831 -0.000 1329.832 -0.000 6156.655 -0.000 10529.215 -0.000 6233.601 -0.000 -0.000 450.310 -0.000 1407.225 0.000 6203.318 -0.000 10633.521 -0.000 6254.361 142.479 -0.000 527.703 -0.000 1453.889 -0.000 6233.601 -0.000 10738.172 -0.000 0.000 219.873 0.000 574.367 -0.000 1484.172 -0.000 6254.361 -0.000 10846.779 QHO Square Perfect QHO Difference 341.200 250.000 91.200 1091.149 750.000 341.149 2607.504 1250.000 1357.504 3836.361 1750.000 2086.361 6937.708 2250.000 4687.708 8364.067 2750.000 5614.067 13428.721 3250.000 10178.721 14703.795 3750.000 10953.795 22085.462 4250.000 17835.462 22856.419 4750.000 18106.419
Here's the first clue that we have a problem: the eigenvalues from the matrix using the limited square well basis set are too high. There is a significant error. Let's plot the resulting eigenfunctions, along with the perfect eigenfunctions, to compare. (As before, don't worry about the minus sign in the functions - this is just a phase.)
# Define a figure to take two plots
fig5 = pl.figure(figsize=[12,3])
# Add subplots: number in y, x, index number
axd = fig5.add_subplot(121,autoscale_on=False,xlim=(0,1),ylim=(-3.1,3.1))
axd.set_title("Eigenvectors for perturbed system")
axd2 = fig5.add_subplot(122,autoscale_on=False,xlim=(0,1),ylim=(-3.1,3.1))
#axc2.set_title("QHO eigenvectors")
axd2.set_title("QHO eigenvectors")
for m in range(3): # Plot the first four states
psi = np.zeros(num_x_points)
for i in range(num_basis):
psi = psi+eigenvectors[i,m]*basis_functions_array[i]
if 2*(m/2)!=m: # This is just to ensure that psi and the basis function have the same phase
psi = -psi
axd.plot(x,psi)
#psi = psi - phi(x2,m,np.sqrt(omegaH))
#axd2.plot(x,psi)
axd2.plot(x,phi(x2,m,np.sqrt(omegaH)))
Here we see the problem: the sum over the basis functions (left graph) contains ringing. We would need to include many more basis functions to correctly represent these rapidly varying functions (if you remember Fourier analysis, you might like to think about what this means in terms of Fourier analysis).
Feel free to play with num_basis
and see how the agreement improves with basis size.