4.3.1 Shooting methods

The popular shooting methods integrate the second-order di¬erential equa-

tion

d2 ψ(x) 2m

= 2 [V (x) ’ E]ψ(x) (4.4)

dx2

from one end with an estimate of the eigenvalue E and iterate over changes

of E; only when E is equal to an eigenvalue will the wave function ful¬ll

the proper boundary condition at the other end. In fact, one replaces a

boundary value problem with an iterative initial value problem.2 The Nu-

merov method is recommended; it consists of solving (4.4) to fourth-order

precision in the grid spacing, requiring the second derivative of the poten-

tial, which is given by (4.4) and which is discretized using three points. It

is left to the reader as an exercise to derive the corresponding algorithm.

Function numerov(m,E,V) ¬nds the nearest eigenvalue by iteration, starting

from a guessed value for E. It shoots from both sides to a common point,

and after scaling, compares the ¬rst derivatives at the end points (Pang,

1997). Then E is adjusted to equalize both relative derivatives with a linear

inter-/extrapolation root search method. The shooting method is in princi-

ple exact (limited only by discretization errors) and has the advantage that

it can also generate excited states, but it is not very suitable for higher

dimensional cases.

python program 4.1 numerov(m,E,V)

Finds the nearest eigenvalue for the single-particle Schr¨dinger equation.

o

01 def numerov(m,E,V):

02 # m=index for matching point

03 # E=trial energy*mass*delx**2/hbar**2

04 # V=array of pot. energy*mass*delx**2/hbar**2

05 # returns [nr of zerocrossings, difference in relat. first

derivatives]

06 E1=E; E2=E1*1.05

07 F1=shoot(m,E1,V)

08 while (abs(E2-E1)> 1.e-8):

09 nF=shoot(m,E2,V); F2=nF[1]

10 Etemp=E2

11 E2=(F1*E2-F2*E1)/(F1-F2)

12 E1=Etemp

13 F1=F2

14 print ™%3d %13.10f™ %(nF[0], E2)

15 return [nF[0],E2]

16

2 The one-dimensional Schr¨dinger equation is a special case of the class of Sturm“Liouville

o

df (x)

d

problems: dx (p(x) dx ) + q(x)f (x) = s(x). See Pang (1997) or Vesely (2001) for a discussion

of such methods in a more mathematical context. Both books describe Numerov™s method in

detail.

4.3 The few-particle problem 81

17 def shoot(m,E,V):

18 # m=index of matching point, should be near right end

19 # E=trial energy*mass*delx**2/hbar**2

20 # V=array of pot. energy*mass*delx**2/hbar**2

21 # returns [nr of zerocrossings, difference in first derivatives]

22 nx=len(V)

23 ypresent=0.; yafter=0.001

24 i=1

25 sign=1.

26 nz=0

27 while i <= m: # shoot from left to right

28 ybefore=ypresent; ypresent=yafter

29 gplus=1.-(V[i+1]-E)/6.

30 gmed=2.+(5./3.)*(V[i]-E)

31 gmin=1.-(V[i-1]-E)/6.

32 yafter=gmed*ypresent/gplus -gmin*ybefore/gplus

33 if (yafter*sign < 0.):

34 nz=nz+1

35 sign=-sign

36 i=i+1

37 ym=ypresent

38 forwardderiv=yafter-ybefore

39 ypresent=0.; yafter=0.001

40 i=nx-2

41 while i >= m: #shoot from right to left

42 ybefore=ypresent; ypresent=yafter

43 gplus=1.-(V[i-1]-E)/6.

44 gmed=2.+(5./3.)*(V[i]-E)

45 gmin=1.-(V[i+1]-E)/6.

46 yafter=gmed*ypresent/gplus -gmin*ybefore/gplus

47 i=i-1

48 backwardderiv=(ybefore-yafter)*ym/ypresent

49 return [nz,forwardder-backwardder]

Comments

Line 02: m is the point where the forward and backward “shooting” should match. It is best

taken near the right border, say at 80% of the length of the vectors.

Line 06: the ¬rst two guesses are E1 and E1 + 5%.

Line 09: this produces the next guess, using the value of E produced in line 12 of the previous

step, which is based on a linear relation between E and the output of shoot (di¬erence between

derivatives at matching point).

The routine may not always converge to the expected (nearest) eigenvalue, as the output of shoot

is very erratic when E deviates far from any eigenvalue.

Line 15: The routine also returns the number of nodes in the wave function, which is an indication

of the eigenvalue number.

Note that numerov does not produce the wave function itself. In order to generate φ when E is

already known, the function psi(m,E,V) can be called.

python program 4.2 psi(m,E,V)

Constructs the wave function for a given exact eigenvalue from the

single-particle Schr¨dinger equation.

o

01 def psi(m,E,V):

02 # m=index of matching point

03 # E=energy*mass*delx**2/hbar**2 must be converged for same m and V

04 # V=array of pot. energy*mass*delx**2/hbar**2

05 # returns wave function y; sum(y**2)=1

06 nx=len(V)

82 Quantum chemistry: solving the time-independent Schr¨dinger equation

o

07 y=zeros(nx,dtype=float)

08 y[1]=0.001

09 i=1

10 while i < m: # shoot from left to right

11 gplus=1.-(V[i+1]-E)/6.

12 gmed=2.+(5./3.)*(V[i]-E)

13 gmin=1.-(V[i-1]-E)/6.

14 y[i+1]=gmed*y[i]/gplus -gmin*y[i-1]/gplus

15 i=i+1

16 ym=y[m]

17 y[-2]=0.001

18 i=nx-2

19 sign=1

20 while i > m: #shoot from right to left

21 gplus=1.-(V[i-1]-E)/6.

22 gmed=2.+(5./3.)*(V[i]-E)

23 gmin=1.-(V[i+1]-E)/6.

24 y[i-1]=gmed*y[i]/gplus -gmin*y[i+1]/gplus

25 i=i-1

26 scale=ym/y[m]

27 for i in range(m,nx): y[i]=scale*y[i]

28 y=y/sqrt(sum(y**2))

29 return y

Comments

This algorithm is similar to the shoot algorithm, except that an array of y-values is kept and the

number of zero crossings is not monitored. Lines 26“27 scale the backward shot on the forward

one; line 28 normalizes the total sum of squares to 1.

Figure 4.1 illustrates the use of numerov and psi to compute the ¬rst six

eigenvalues and eigenfunctions for a double-well potential, exemplifying the