Numerical integration of the Schrödinger

Numerical integration of the Schrödinger equation

***Please read the file under additional files.

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