Numerical solution to Schrödinger equation - eigenvalues

Question

This is my first question on here. I'm trying to numerically solve the Schrödinger equation for the Woods-Saxon Potential and find the energy eigenvalues and eigenfunctions but I am confused about how exactly this should be done.

I've solved some initial value problems in the past using iterative methods such as Runge–Kutta. I've read that Numerov's method is the way to solve Schrödinger's equation but Wikipedia also describes it as an iterative method for initial value problems.

How do I use it to solve an eigenvalue problem?

This confuses me for the following reasons:

• Wouldn't iteratively solving the DE require knowledge of the energy eigenvalues to use as input to the calculation? I don't know the eigenvalues yet; they're precisely what I'm trying to calculate.
• If I did that, wouldn't I simply get a unique solution, instead of a family of eigenfunctions and eigenvalues?

I've seen some mention of "tridiagonal matrices" being generated somehow, but am not sure what the elements of that matrix would be or how that applies to the problem. Leandro M. mentioned that "the discretization defines a finite dimensional (matrix) eigenvalue problem". This seems like the correct road I should be going down, but I haven't been able to find anything that explicitly explains this process or how the matrix is constructed. If this is the correct procedure, how is such a matrix constructed?

Answer 1

I'm glad I can finally answer this.

Numerov's method as described on Wikipedia is not how you want to proceed here. To give you an idea of how to proceed, let's start with a simplified version of the method. What I'm going to do is to just naively discretize the differential operator, like so:

$$ \frac{d^2}{dx^2}\psi \approx \frac{\psi(x-d)-2\psi(x) +\psi(x+d)}{d^2} \equiv \frac{\psi_{n-1} - 2 \psi_n + \psi_{n+1}}{d^2}$$

The last equation is just a definition -- I'm treating space as if it's a lattice with points $d$ apart and I'm calling $\psi_n$ the value of the wavefunction on the $n$th point.. Now the Schrödinger equation reads in this notation:

$$-\frac{\hbar^2}{2m}\frac{\psi_{n-1} - 2 \psi_n + \psi_{n+1}}{d^2} + V_n \psi_n = E \psi_n$$

But this is a matrix equation! Let me be more explicit:

\begin{align} -\frac{\hbar^2}{2md^2}&\left(\begin{array}{cccccccc} -2 & 1 & & & \\ 1 & -2 & 1 & & \\ &\ddots&\ddots&\ddots& \\ & & 1 & -2 & 1 \\ & & & 1 & -2 \\ \end{array}\right) \left(\begin{array}{c} \psi_1 \\ \psi_2 \\ \vdots \\ \psi_{N-1} \\ \psi_N \end{array}\right) +\\ &\qquad\qquad\quad\left(\begin{array}{cccccccc} V_1 & & & & \\ & V_2 & & & \\ & & \ddots & & \\ & & & V_{N-1} & \\ & & & & V_N\\ \end{array}\right)\left(\begin{array}{c} \psi_1 \\ \psi_2 \\ \vdots \\ \psi_{N-1} \\ \psi_N \end{array}\right) = E \left(\begin{array}{c} \psi_1 \\ \psi_2 \\ \vdots \\ \psi_{N-1} \\ \psi_N \end{array}\right) \end{align}

Of course, I had to pick some integer $N$ that just corresponds to placing the system in a box that's "big enough" in order to have a finite system.

It is clear now that what you have is a matrix eigenvalue problem of the form

$$A\psi = E \psi$$

and you may proceed to diagonalize it in whatever way you choose. Note that we call $A$ a tridiagonal matrix, for obvious reasons. You might want to take care of the boundary conditions first -- you do that by setting $\psi_1 = \psi_n = 0$ before you construct the matrix, which corresponds to setting the first and last columns to zero. This will give you some spurious eigenfunctions with zero eigenvalue that you can just discard. If you have different boundary conditions, you're out of luck -- I don't know of a way that makes it work in general.

Now you just have to redo the above with the full Numerov's method, which will be a bit more complicated, and you're all set.

This seems to be the way you're looking for but of course it's not the only way to do this. Griffiths describes one called "wag the dog" that consists of the following: you guess an initial value for an energy and compute the wavefunction however you want (RK, for instance). Chances are it won't be an eigenvalue so the wavefunction will blow up at infinity. It might go to $+\infty$ for large $x$ or it might go to $-\infty$. You now slowly vary the energy until the behavior at infinity "flips", that is, until the tail "wags". That will allow you to constrain the value of a single energy eigenvalue and give you the form of the wavefunction out to some maximum size that increases as the chosen $E$ approaches the exact eigenvalue.

Answer 2

You can search for eigenvalues using the bisection method.

Priliminaries: To get the eigenvalues from Numerov method you will need to know the wavefunction at the boundaries. Generally this would mean that you need to set the potential to infinity at the boundaries hence putting the wavefunction to zero at those points. For your potential, modify it as follows:

V = infinity if -50fm>r>50fm V = 0 if |r|>8.5fm V = Woods potential otherwise

The real deal: Now write a program to calculate the wavefunction in the above potential using any arbitrary energy in the domain -50fm

Here is a python code that finds an eigen value in the given interval (E1,E2) if it exists.

---

import numpy as np

rmin = 0 rmax = 0

def vradial(r): global rmin global rmax if r < rmin or r > rmax: #This is important. >= and <= won't work #as they interfere with the bc on psi return np.inf elif r > rmin and r < rmin + 2: return (r - rmin - 2)**2 elif r > rmax - 2 and r < rmax: return (r - rmax + 2)**2 else: return 0

def f(l,E,r): """Calculate the f(r) in the Schrodinger equation of the form D2(Psi(r)) = f(r)Psi(r)""" if r == 0: return np.inf return l*(l+1)/(r**2) + vradial(r) - E

def psitophi(psi, l, E, r, delta): if psi == 0: return 0 #To handle the 0*infinity case of boundary return psi*(1-(f(l, E, r)*delta**2)/12)

def phitopsi(phi, l, E, r, delta): #if phi == 0: return 0 return phi/(1-(f(l, E, r)*delta**2)/12)

def calcpsiradial(rmin, rmax, N, psibcmin, psibcmax, l, E):

psiarray = []
r = np.linspace(rmin, rmax, N)
delta = r[1]-r[0]

#Assume psi(rmin) = psibcmin and psi(rmin+delta) = 1 and then
#calculate phi0 and phi1
psiarray.append(psibcmin)
psiarray.append(1)
phi0 = psitophi(psiarray[0], l, E, r[0], delta) #r[0]=rmin
phi1 = psitophi(psiarray[1], l, E, r[1], delta)

#Now populate the psiarray for each value of r
for i in range(2, len(r)):
    phi2 = 2*phi1 - phi0 + (delta**2)*f(l, E, r[i-1])*psiarray[i-1]
    psi = phitopsi(phi2, l, E, r[i], delta)
    psiarray.append(psi)
    phi0 = phi1
    phi1 = phi2

return r, psiarray

def normalize(delta, psi): area = 0 for i in range(1,len(psi)-1): area = area + abs(psi[i])*delta

for i in range(len(psi)):
    psi[i] = psi[i]/area

return psi

def locateEvalueInBracket(rmin, rmax, N, psibcmin, psibcmax, l, e1, e2,\ tol): #Any value of Psi smaller than psi is while abs(e2-e1) > tol:
r, psi = calcpsiradial(rmin, rmax, N, psibcmin, psibcmax, l, e1) psi1 = psi.pop() r, psi = calcpsiradial(rmin, rmax, N, psibcmin, psibcmax, l, e2) psi2 = psi.pop()

    if psi1*psi2 < 0:
        emid = e1 + (e2 - e1) * 0.5
        r, psi = calcpsiradial(rmin, rmax, N, psibcmin, psibcmax, l, emid)
        psimid = psi.pop()

        if psimid*psi1 < 0:
            e2 = emid

        elif psimid*psi2 < 0:
            e1 = emid
    elif psi1*psi2 > 0:
        print "There are either no eigenvalues or too many of them in"+\
        " the given energy interval. For e1={} psi={} and e2={} psi=".format(e1, psi1, e2)+\
        "{}".format(psi2)
        return e1,e2
return e1,e2

def main(): import sys import matplotlib.pyplot as plt global rmin global rmax

e1 = float(sys.argv[1])
e2 = float(sys.argv[2])
l = int(sys.argv[3])
if l == 0:
    rmin = -1
else:
    rmin = 0
tol =1e-5
rmax = 1
psibcmin = 0
psibcmax = 0
N = 100

e1, e2 = locateEvalueInBracket(rmin, rmax, N, psibcmin, psibcmax,\
        l, e1, e2, tol)

fig = plt.figure()
ax = fig.add_subplot(111)
r, psi = calcpsiradial(rmin, rmax, N, psibcmin, psibcmax, l, e1)
ax.plot(r,psi, 'o')
r, psi = calcpsiradial(rmin, rmax, N, psibcmin, psibcmax, l ,e2)
ax.plot(r, psi, 'g')
ax.set_title("l = {} and E_blue = {}, E_green={}".format(l,e1,e2))
plt.show()

if name=="main": main()