Computing the wave function of a square potential

[CrosisBH]

Homework Statement

Recreate the given plot.

Relevant Equations

The Shooting Method algorithm.

The book's procedure for the "shooting method"

The point of this program is to compute a wave function and to try and home in on the ground eigenvalue energy, which i should expect pi^2 / 8 = 1.2337...

This is my program (written in python)

import matplotlib.pyplot as plt
import numpy as np

dx = 0.01
V0 = 1000 #idk units
cutoff = 2.0#number of points of size dx from -1 to 1
N = int(2/dx)

#idk the units on these, natural units i believe.
E = 1
dE = 0.1

#x_axis to make our potential
x_axis = np.arange(-3,3,dx)

#x_axis for our wavefunction
x_axis_wave = np.arange(-1,1,dx)

#actual potential
potential = np.zeros(int(6/dx))#setting up the square well potential
for i in range(0,len(potential)):
    if(abs(x_axis[i])<=1):
        potential[i] = 0
    else:
        potential[i] = V0last_diverge = 1
#implementing psi(0) = psi(-1) = 0, which makes the derivitive 0, required for even parity
wavefunction = np.zeros(int(2/dx))
wavefunction[0]=1
wavefunction[int(len(wavefunction)/2)]=1

fig = plt.figure(figsize = (5,5))

#the book has manual dE messing with but i found this to be easier to implement
while(abs(dE)>0.025):

    #for reach point the wave function is defined by
    for i in range(0,N-1):
        #step function for the 2nd derivitive
        wavefunction[i+1] = 2*wavefunction[i] - wavefunction[i-1]-2*(E-potential[i])* dx**2 *wavefunction[i]

        #divergence check
        if(abs(wavefunction[i+1])>cutoff):
            break

    #plot the wave function and log the data points of interest.
    plt.plot(x_axis_wave, wavefunction)
    print("ΔE = ", dE, "E = ",f"{E:0.2f}","Ψ(1) = ", f"{wavefunction[i+1]:0.2f}", end="\r")

    #updating energy
    #updating dE if the signs are different
    if  wavefunction[i+1]*last_diverge<0:
        dE /= -2

    E = E+dE
    last_diverge = np.sign(wavefunction[i+1])
    
print("ΔE = ", dE, "E = ",f"{E:0.2f}","Ψ(1) = ", f"{wavefunction[i+1]:0.2f}")
#potential lines
plt.plot([-1,-1],[-100,100], color='black',linestyle=':')
plt.plot([ 1, 1],[-100,100], color='black',linestyle=':')

plt.xlim(-2,2)
plt.ylim(-2,2)

plt.xlabel("x")
plt.ylabel(r"$\psi$")

plt.xticks(np.linspace(-2,2,5))
plt.yticks(np.linspace(-2,2,5))

plt.title("Even-Parity Wave Function Calculation")

plt.show()

The program with the code is that it doesn't want to converge on an energy, it kind goes all over the place until the dE condition is met. The associated "eigenvalue" (if you can even call it that at this point) is 2339.77. This is the closest I've gotten to it. Also the (incomplete) plot.

Which of course doesn't match the original plot at all.

 

[TSny]

I have little experience with Python. But, it appears to me that there might be a problem with how you are relating values of the index i to values of x. For example, in the code you set wavefunction[0] = 1. For what value of x does this correspond? What is the value of x corresponding to x_axis_wave[0]?

Likewise, if you consider potential[0], for what value of x does this correspond? What is the value of x corresponding to x_axis[0]?

(I hope I'm not just showing my ignorance of Python. )

 

[CrosisBH]

I have little experience with Python. But, it appears to me that there might be a problem with how you are relating values of the index i to values of x. For example, in the code you set wavefunction[0] = 1. For what value of x does this correspond? What is the value of x corresponding to x_axis_wave[0]?

Likewise, if you consider potential[0], for what value of x does this correspond? What is the value of x corresponding to x_axis[0]?

(I hope I'm not just showing my ignorance of Python. )

Every one of those would represent the first point in the list, corresponding to the very left of it's domain. wavefunction[0] and x_axis_wave[0] correspond to the their function at -1, while with potential and it's domain's 0th element is their function at -3. I have constructed the lists so that their elements run from -1 to 1 or -3 to 3 with increments of dx. I'm not even sure if this is the best way to do it. I'm open for any suggestions.

 

[TSny]

Every one of those would represent the first point in the list, corresponding to the very left of it's domain. wavefunction[0] and x_axis_wave[0] correspond to the their function at -1, while with potential and it's domain's 0th element is their function at -3. I have constructed the lists so that their elements run from -1 to 1 or -3 to 3 with increments of dx. I'm not even sure if this is the best way to do it. I'm open for any suggestions.

Ok. But then it seems to me that for a specific value of the index i, say i = 10, wavefunction[10] and potential[10] would evaluate the wavefunction and the potential at different spatial positions x. wavefunction[10] would correspond to the value of the wavefunction at 10*dx to the right of x = -1, while potential[10] would correspond to the value of the potential at 10*dx to the right of x = -3.

So, in the line

wavefunction(i) and potential(i) are not being evaluated at the same spatial position x. But they should be evaluated at the same x.

Also, I think you want to start the "shooting" process at the center of the well where x = 0. But, you're starting the process at i = 0, which is not the center of the well.

 

[CrosisBH]

Ok. But then it seems to me that for a specific value of the index i, say i = 10, wavefunction[10] and potential[10] would evaluate the wavefunction and the potential at different spatial positions x. wavefunction[10] would correspond to the value of the wavefunction at 10*dx to the right of x = -1, while potential[10] would correspond to the value of the potential at 10*dx to the right of x = -3.

So, in the line

View attachment 260050

wavefunction(i) and potential(i) are not being evaluated at the same spatial position x. But they should be evaluated at the same x.

Also, I think you want to start the "shooting" process at the center of the well where x = 0. But, you're starting the process at i = 0, which is not the center of the well.

I changed how the potential is handled with a simple function

def potential(x):
    if (abs(x)<=1):
        return 0
    else:
        return V0

So now my stepping function is

Which surprisingly yielded the same plot. I will look into starting at x = 0. Will come back with updates.UPDATE:
I remembered my professor posted the website to the author's source code if we get stuck. This is the program in question. http://www.physics.purdue.edu/~hisao/book/www/examples/chap10/SHOOT.TRU It's in True BASIC but I was able to parse it. I learned a few things about from this code:

• I have no idea what psi_-1 represents, especially in the true BASIC context.
• It's probably better to evaluate potential using potential(i*dx) (worked better for me anyway)
• I think the author only cared about psi(0) to psi(1) due to symmetry in the problem.