Numerical integration of the Schrödinger equation

6.1 The Problem
Solve for the stationary states of an electron in a ramped infinite square well.

6.2 Introduction
A basic problem in quantum mechanics is to find the stationary states (“energy levels”) of a
bound system, using the time-independent Schrödinger equation (TISE)

(-h^2(bar)/2m)(d^2/dx^2)(psi(x)) + V(x)*psi(x) = E*psi(x)

In this activity, we will look at the numerical
integration of the differential equation.
The problem we face here is not that of the initial value problem of mechanics, that we dealt
with earlier. We can’t just start from some initial value, use numerical integration, and end up
somewhere after some time. Instead, solving the TISE is a boundary value problem. There is
an unknown quantity (E) in the differential equation, and we need to solve for this, subject to
constraints on the boundary conditions of the wave function.
How to go about this? How do we tell if we have a correct value for E (called an eigenvalue)?
We have to look at the properties of the required solution. We find that if E is an energy
eigenvalue, the wave function ψ(x) has the expected behaviour at the boundaries. If we choose
the wrong value for E, then the wave function will not have this property.
Let us consider a specific example. Suppose that we have an electron in a ramped infinite square
well with walls at x = 0 and x = a:

V (x) = {bx for 0 < x < a
{∞ otherwise

Since the potential energy is infinite outside the well, the electron cannot be beyond the well’s
walls. Continuity of the wave function then implies that ψ(0) = ψ(a) = 0, exactly as in the case
of the unramped infinite square well considered in lectures. However, owing to the changing
potential energy in the well, this is a more interesting problem than the one done in class.
A way of finding a solution numerically then suggests itself:

1. Assume a value for E.
2. Start at one of the walls (e.g. the left wall at x = 0), setting ψ to zero there.
3. Integrate numerically to the other wall.
4. If the numerical solution is close (?!) to zero there, we have a correct value for E
5. If we don’t have a correct value, choose another and try again.

Ideally, we would use some sort of root-finding algorithm to drive this process, but we will see
that it can be done by hand.

6.3 Numerical Methods
We can do the numerical integration as we have for Newton’s equation, by splitting the second
order ODE into coupled first order equations. But it’s just as easy to use a numerical
approximation to the second derivative

d^2ψ(x)/dx2 ≈ ψ(x − ∆x) − 2ψ(x) + ψ(x + ∆x)/(∆x)^2

Where does this approximation come from? Substitute the Taylor expansion of ψ(x ± ∆x) into! the RHS of Eq. 6.1 and see.

Note that the integration is over x, using a fixed step size ∆x.
Then we obtain for the TISE:

Then we obtain for the TISE:
ψi+1 = 2ψi − ψi−1 + (∆x)^2g(E, xi) ψi

where
g(E, xi) = −2m/h^2(bar)(E − V (x)) = −2mc^2/h^2(bar)c^2(E − V (x))

for i = 1 . . . N − 1.

Note that, to start this off, we need two values of the wave function. One is at the left wall:
ψ0 = 0. What about the next, ψ1? We can take any value for this. The wave function is usually
taken to be normalised, but this is not a consequence of the solution, but is imposed by our
interpretation of the wave function. Any multiple of ψ(x) is still a solution. We can take any
value for ψ1 — this simply affects the scaling of the wave function — and normalise afterwards.

6.4 Model
Use the following: a = 3.0 nm, b = 0.5 eV nm^−1,mc^2 = 5.11×105 eV, h(bar)c = 197.3 eV nm, where
c is the speed of light. Note that these values allow you to use “natural scale” units of eV and
nm for an atom-sized problem. (And you don’t have to bother with laborious conversions to SI
units).

6.5 Computer Solution
Given the discussion above, we can outline an algorithm for searching for the energy levels in
our problem.
1. Integrate from x = 0 nm to x = a = 3.0 nm with a step size of, say, 0.01 nm and a suitable
initial value of E (perhaps close to the ground state of the infinite unramped square well
of the same dimensions). For the initial values choose ψ0 = 0 and ψ1 something small,
say 10^−6.Some experimenting might be necessary. Print out the value ψN of the wave function at the right wall.
——————————————————————
psi = np . zeros ( N+1,float )
psi [0] = 0. 0 # two initial conditions

psi [1] = 1.0 e−6
# i n t e g r a t e
for i in range ( 1 , N ) :
. . .

——————————————————————-
2. You can adjust the value of E in order to home in on the eigenvalue for the ground state.
Find the region where ψN changes sign. This indicates that an eigenvalue has been stepped
over. You will probably need to read in a value of E from the keyboard. This is easy:
——————————————————————-
. . .
while 1 : # l o o p f o r e v e r w hil e t r yi n g new v al u e s o f E
E = input ( ’ Enter a value for E : ’ )
. . .
# i n t e g r a t e
for i in range ( 1 , N ) :
. . .
# print out last value
——————————————————————-
3. Normalize and then plot out the wave function. A simple approximation for an integral is
given by the trapezoidal rule

(Please look at the file under “Additional files”
to see the equation)

Can you interpret this equation?
Using numpy this can be written: …

***Please look at the file under “Additional files” to see the rest of the assignment,