Solving the Schrödinger equation for bound states

Computer Physics Communications 34 (1985) 287 293 North-Holland, Amsterdam

287

SOLVING THE SCHRODINGER EQUATION FOR BOUND STATES P. FALKENSTEINER, H. GROSSE, F. SCHOBERL Institut für Theoreiische Physik, Universität Wien, 1090 Vienna, Boltzmanng. 5, Austria

and P. HERTEL Fachbereich Physik, Universitdt Osnabriwk, Osnabruck, Fed. Rep. Germany Received 6 July 1984

PROGRAM SUMMARY Title of program; SCR2

ments, scale transformation, numerical integration, convexity arguments, node theorem

Catalogue number; ACDQ Nature of the physical problem Program obtainable from; CPC Program Library, Queen’s University of Belfast, N. Ireland (see application form in this issue) Computer. VAX 11/750; Installation; Prozessrechenanlage Physik, Universität Wien

SCR2 calculates energy levels and wave functions of the rotational symmetric Schrodinger equation for a given potential using a simple and accurate method. Method ofsolution An iterative procedure for getting upper and lower bounds to

Operating system; VMS 35

energy values of the radial Schrodinger equation is given. A numerical integration procedure together with convexity argu-

Note; SCR2 runs also on the CDC CYBER 170/720 without modifications

ments and the nodal theorem for wave functions is used.

Programming language u.sed. USANSI FORTRAN 77

Typical running time 5 s for one bound state. Note; The running time depends strongly on the desired accuracy.

Program size; 2 Kbyte References

No. of bits in a word; 32 No. oflines in combinedprogram and test deck; 274 Keywords. Schrodinger equation, rotational symmetric poten-

tial, reduced wave function, bound states, eigenvalues, mo-

[1) For applications see: H. Grosse and A. Martin, Phys. Rep. 60 (1980) 341. [2] M. Abramowitz and l.A. Stegun, Handbook of Mathematical Functions (Dover, New York, 1968). [3] N. Séto, PubI. RIMS, Kyoto University 9 (1974) 429. [4] H. von Eicken, CERN Library D 100.

OO1O-4655/85/$03.30 © Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)

P Falkensteiner et al

288

/

Schrodinger equation for bound stales

LONG WRITE-UP I. Introduction

2. Main idea

The Schrodinger equation is one of the fundamental equations of physics; nevertheless, exact analytic solution are known only for exceptional cases. On the other hand the two body with rotational symmetric potential V( r),problem where r —xl and x denotes the relative coordinate, reduces to an ordinary differential equation for the reduced wave function y~ 1(r): 2

1(1+1)

(v~ (x)

y,, 1(x) 2, ~~~(x) V(x)+/(I+1) v V(x) 2p.a2V(ax). E,,, 2p.a2E,,

y,’, ( x)

=

E,,

,)

r/u

(3)

-

—

s,

,.

In the following we shall consider (3) for values of E,, even if they are not eigenvalues, the differential equation

and integrate

d

[h2(

r~

)

+ V(r)]y,

7 ,(r)

E,,,j~,,1(r), 2 I

y”(x)

(1) —

0 Jdr (y~1(r))

...

—

The energies and normalized functions are to be determined. With the help of the latter the values of 1 I ______ d’~ dr’~1 ~‘

which enter into decay widths and the moments

—f

dr rm53,~,(r)

(~ff(x)

~v(x)

(4)

from the origin. In order to be able to treat ally start from potentials that xare —6singular (which at is close the origin to thewe origin); actuwe integrate with the boundary conditions

and numerical integration should be possible. Here / 0. 1. denotes the angular momentum quanturn number, ~sis the reduced mass, h Planck’s constants and E~ 1denote eigenvalues, which are labelled for fixed / by the number of nodes n 0, 1, 2, the bound state wave function has within (0, oc). The2potential is assumed to be less singular at the origin, locally integrable and than 1/r should have a finite number of zeros only.

I,,,

The scale transformatton from r to allows us to transform eq. (1) into

(2)

0

entering into transition amplitudes can be calculated [I]. The main problem is the boundary condition at infinity such that normalization is possible. Our main idea is to combine two basic facts in order to solve this problem: We use convexity arguments and the nodal theorem.

v( 8) 6’~ v’( 6) (1 + 1)6’. (5) In order to incorporate the integrability condition of y~,(x),we observe that there is a ~ such that V(x)> E,,1 for all x

>

x5j where E,,1 is the

appropriate eigenvalue which we would like to determine. This follows from the assumption that V does not oscillate infinity; x~1 denotes the classical turning pointatwhich is most to the right (largest x solving the equation i’~(x) E~ 1).Our first simple observation is that sign y”(x) sign v(x) if x > (6) which means thaty(x) is convex (concave) if y(x) is positive (negative). Therefore integrating out to some point x> x~ 1for whichy(x)> 0 andy’(x)> 0 implies that the solution of eq. (4) goes to plus infinity for x oo; if y(x) <0 and y’(x) <0 it goes to minus infinity. This will be the reason why —~

we are allowed later on to stop the integration procedure in this case. The next fact concerns the nodal structure of y~1(x). It is well-known that the ground state wave function (n 0) has no nodes within (0, oc); the radial excitations have a finite number of nodes (n 1, 2, ...). Therefore our procedure runs as follows: We fix an interval for the energy variable ~ (EL. E~), —

—

P. Falkensteiner el al. / Schrixiinger equation for bound stales

289

where E[ and E~are lower and upper bounds to a particular E~,which is to be determined. (A bad

The real EL and EU are lower and upper bounds EL and E0 to the energy En,, respectively. The

choice of EL or E~will be discussed below.) Take the arithmetic mean CM = (EL + E~)/2 and integrate eq. (4) with = CM starting from x = 6 with boundary conditions (5) counting nodes. Two possibilities may occur: a) The number of nodes exceeds n before we reach the point I defined above; then we stop further integration since we know E~, x~1,since within the classically forbidden region we know that

real FEH represents a kind of error: SCR2 proceeds by improving the bounds EL and EU until the length of the interval (EL, EU) is less than FEH. The real H is the stepsize for the numerical integration [2] of eq. (4). The real DEL is the starting point 6 for the integration of eq. (4) and should be nonzero but small. For an effective potential having no minimum the real XWMIN has to fulfil 0 XWMIN DEL; otherwise XWMIN represents the minimum of J’~most to the right; for practical purposes one takes an upper bound to that minimum which is within a distance H to it.

(signy”(x))(signy(x)) >0.

The line size of the two-dimensional array FELD is given by NSR. The potential V(x) has to be defined in a Function subprogram FUNCTION vL(x). (9)

(7)

Then we stop integrating further since we know that the true energy lies above EM. Therefore we replace the old EL by M and continue the integration procedure from x = 8. In order to integrate eq. (4) numerically we use the Runge—Kutta method described in ref. [2]. If the bounds EL and E0 are chosen badly, such that the true energy E~1~ (EL, E0), e converges to the bound which lies next to the true En,. 3. Using the program SCR2 SCR2 is a subroutine subprogram which solves the eigenvalue problem of eq. (3) numerically by the method described in section 2. For a given rotational symmetric potential V(x) it calculates the energy and the reduced radial wave function for a bound state with certain number of nodes and certain angular momentum. The calling sequence is: CALL SCR2(NO, RL, EL, EU, FEH, H, DEL, XWMIN, EP, NSR, NSX, FELD).

(8)

The variables NO, RL, EL, EU, FEH, H, DEL, XWMIN and NSR provide the input data for SCR2: The integer No is the number of nodes n(NO = 0, 1, 2, ...). RL is the angular momentum quantum number 1 and usually RL = 0., 1., 2.

It is necessary to take the name VL for that subroutine. The variables EP, NSX and the array FELD provide the output data from SCR2: The real EP is the energy E~,of the bound state. More precisely: the true energy E~,lies within the interval (EP FEH/2, EP + FEH/2), if the integration procedure is accurate enough. For a bad choice of EU or EL such that E~,>EU or E~,< EL the iteration procedure gives an EP being less than FEH/2 away from the bound. One must —

repeat the calculation with new bounds such that finally E~1E (EL, Ed). The integer NSX counts the number of integration steps made by integrating eq. (4) with energy EP until the criteria explained in section 2 apply and integration is stopped. In the real array FELD the reduced radial wave function y (integrating eq. (4) with EP) and its derivative y’ are stored. The dimension of FELD has to be declared in the calling program of SCR2: DIMENSION FELD (2, value of NSR),

(10)

FELD (1, J) (J = 1,... ,NSX 1) stores the value of y(x) at x = DEL + (J 1)H; FELD (2, J) (J = 1,... ,NSX 1) the value of y’(x). One should stress that the wave function fulfils boundary con—

—

—

290

P. Falkensteiner el al.

Schrodinger equation for bound states

ditions eq. (5) but is not (in general) normalized like: 1.

dxy~(x)

(11)

Normalization may be done in the calling program of SCR2 (see section 5). NSX has to be less than NSR. If NSX is equal to NSR one must repeat the calculation using a higher NSR or a larger stepsize H. Otherwise essential parts of the wave function y(x) may not be stored in the array FELD. if the determination of the wave function is sufficiently accurate (which is the case if H and FEH are sufficiently small) the values of FELD (1, NSX-1) and FELD (2, NSX-1) are negligible (since y(x) and y’(x) goes to zero for x oc). (For subtleties see remark 2.) After calling SCR2 the input data EL and EU have been changed to the value of the last improved bounds for the energy En,. At the end of this section let us make two final remarks: I) The angular momentum variable RL has been taken to be real, since non-integer values may be of interest as well: If one starts from the d-dimensional Schrodinger operator with rotationally symmetric potential —*

V(x))~ti(E)=E4i(E), x=~j,

(—~d+

(12)

where ~ denotes the d-dimensional Laplacian, the introduction of the reduced wave function y(x) yields 2 (l+~(d 3))(l+~(d 1)) / dx2 d + + V(x)) —

Xy(x)

—

1(~) x~ d)

{x~fr~11(x) e <0 or x ~

<

XWMIN} (15)

<~5.(~)}.

where x~1(c)denotes the zero of Veii(X) e which is most to the right. For a monotonic effective potential (or one which has only one extremum) XWMIN fulfils condition (15). If V,~has at least one maximum to the left of XWMIN it may occur that there is an such that ~ ~ ~ e <0 or x
}

D

{x~x
(16)

This does not spoil a correct calculation of E,,, and the wave function y. But it might happen that XWMIN is much larger than x~1(E,5,)and the cut off for stopping integration is too large: As a consequence the calculated v (stored in FELD (I, J)) has a minimum (or maximum) for y(x)> 0 (or y(x) < 0) far before the cut off occurs. In that case one should take the last extremum of ~‘(x) as the cut off. This can be done either by a simple criterium in the calling program of SCR2. Or, alternatively, one may determine x~1(E,,,) by determining the zero of V~~1(x) E~,lying most to the right using a suitable subroutine. Finally, one calls SCR2 once more and takes for XWMIN = x~1(E,,,)and for EL, EU the last improved bounds on

(13)

Ey(x),

q

XWMIN. Note that XWMIN was defined to be the minimum of V(x) which is most to the right ~or 0 ~ XWMIN ~ DEL). For any fixed value e (EL, E~)we would like to have equality of the two sets

2y(x)Y 1(n),

n —~/x,

4. Program structure, common blocks

where Y1(n) denotes a harmonic polynomial of degree / in d dimensions [31 Now one may choose RL such that

The subroutine subprogram SCR2 uses one subroutine subprogram named RK2L. RK2L integrates eq. (4) one step H further with the

RL(RL

Runge Kutta method [21.RK2L uses one function subprogram named DIFFL. DIFFL calculates J/~11(x) at the value of x (compare eq. (4)). DIFFL uses one function subprogram named VL which has to be provided by the user (section 3). There exists one common block labelled by

-i-i)

=

(i-i- ~(d— 3))(I+ ~(d— 1))

(14)

and proceed as we did before in the 3-dimensional case. 2) Here we point to a problem which may occur if f’~~~(x) has a local maximum to the left of

P. Falkensteiner et al.

/

Schrodinger equation for bound states

SC93L. The common variables are: 1) U2 calculated in RK2L as DIFFL(X) (in order to avoid a second call of DIFFL in SCR2). 2) L (declared as REAL) is the angular momentum quantum number / 3) EPS has the value of C (compare eq. (4)).

and the square of the reduced radial wave function at the origin together with the square of the derivative which may be expressed as y2(0)

=

_2f

y’2(0)= _2f 5. Test run

291

dxy(x)y’(x),

(19)

dxy’(x)y”(x).

(20)

.

y”(x) is easily obtained from y by eq. (4). Quite

For the test run SCR2 is called by a main program TEST which serves the following purposes: 1) Reading and printing the input data. 2) Obtaining bound state energies and wave functions by calling SCR2. 3) Calculating various matrix elements from the wave function. 4) Printing the output data in a suitable manner. Clearly SCR2 is especially of interest for potentials for which the Schrodinger equation is no: analytically solvable. Nevertheless, for the test run we take the harmonic oscillator potential V(x) = x2 which can be analytically treated and allows to compare the numerical results to the exact ones. V (x) = x2 + 1(1+ 1)/x2 (17) eff

generally the numerical accuracy is improved by expressing local quantities by global ones. Since y(O) has to vanish for all states, calculation of eq. (19) provides a good consistency check. y’(O) is nonvanishing only for 1= 0 states. In order to normalize matrix elements according to (11), one simply divides them by I~= f~ dx y2(x).

For the numerical integration the Simpson rule integration program ARSIMP [4] has been used. Note that the integration from zero to infinity has been replaced by a summation from 6 to 6 + (NSX 1)H, since 8 is close to zero and y and y’ are negligible at x = 8 + (NSX 1)H. In order to test the accuracy of SCR2 one can 1) check with potentials which are exactly solvable. Some results for V(x) = 1/x for variable —

—

—

has its unique minimum at x

=

[1(1+

1)]i ~ For more complicated potentials the minimum lying most to the right has to be determined numerically before calling SCR2. After calling SCR2 the program TEST calculates the moments ~ ~ 2 defined by 2(x) (18) I~ref dx xey

stepsize H and starting points DEL are shown in table 1, 2) examine the stability of results for the one particular problem while increasing the accuracy (decreasing H, DEL and FEH). It should be remarked that the running time of SCR2 dependsproportional strongly on totheI/H. desired accuracy and is almost Furthermore,

Table 1 We give the energy E, the moments I ~ 1 2 and the value for y’(0)2 for the IS, 2S and 2P states as a function of the stepsize H. The starting point is DEL = H/lU., FEH is 0.00001. Exact results are given in brackets State H —

E

1 /

IS 0.2 0.249882 0.4996

2

(y’(0))2

0.4886 0.4791

0.1 0.249973 (0.25) 0.4999 (0.5) 0.4946 (0.5) 0.4898 (0.5)

2S 0.05 0.249996

0.2 0.062489

0.4999

0.1249

0.4974

0.06106

0.4949

0.05987

0.1 0.062496 (0.0625) 0.1250 (0.125) 0.06184 (0.0625) 0.06123 (0.0625)

2P 0.05 0.062496

0.2 0.062504

0.1250

0.1250

0.06221

0.02084

0.06190

2.0x10

0.1 0.062504 (0.0625) 0.1250 (0.125) 0.02082 (0.02083) 6.1x10 ~ (0.)

0.05 0.06504 0.12500.1250 0.02082 l.7><1O

6

292

P. Falkensleiner et al

Schrbdinger equation for bound states

it depends on the complexity of the potential V(x), on the magnitude of n and / and on the difference between EL and EU. Finally, it is worth mentioning that there is no significant difference between the results from a CDC CYBER (60 bits in a word) and those from a VAX (32 bits in a word). Thus, one may conclude that a word length of 32 bits is sufficient for solving Schrodinger problems with SCR2.

References Ill

For applications see:

El.

Urosse and A. Martin, Phys. Rep

60(1980) 341. [2] M. Abramowitz and l.A. Stegun. Handbook of Mathematical Functions (Dover, New York, 1968). [3] N Séto, PubI. RIMS, Kyoto University 9 (1974) 429. [4] H. von Eicken, CERN Library D 100.

P. Falkensteiner et a!.

/

Schrodinger equation for bound states

293

TEST RUN OUTPUT STEPSIZE H = ERROR FEEl = STARTINOPOINT DEL EL EU RI NO

.

.

.

.

-

.

-

.

-

.

.

.

=

0.0100 0.000010 0.00100

LOWER BOUND OF ENERGY UPPER BOUND OF ENERGY ANGULARMOMENTUM NUMBER OF NODES

ENERGY . . . BOUM)STATE ENERGY NSX . . . NUMBER OF INTEGRATION STEPS 1-N . . . MOMENTS YlO)ft*P SQUARE OF THE WAVE FUNCTION Y AT TIlE ORIGIN YP(O)*~2 . . . SQUARE OF THE DERIVATIVE Of: v AT iRE ORIGIN EL 0.00 0.00 5.00 0.00 0.00 3.00

EU

I RL

I 10.00 I I 15.00 1 I 15.00 I I 10.00 1 I 15.00 I

0. 0. 0. 1. 1.

NO I I 0 I I 1 1 I 2 I I 0 I I 1 1

ENERGY

NSX

1-1

1-2

Y(O)**2

YPIOI**2

2.99999’?

453

1.1284

1.9977

0.0000

2.2567

6.999997

501

0.9403

1.9966

0.0000

3.3851

10.999990

557

0.8369

1.9958

0.0000

4.231A

4.999995

467

0.7523

0.666/

0.0000

-0.0001

8.999995

531

0.6770

0.666’?

0.0000

-0.0002

12.999990

609

0.6260

0.666’?

0.0000

-0.0004

6.999996

493

0.6018

0.4000

0.0000

0.0000

1

1.

3.00

20.00 1 I 11.00 I

2.

2 i I 0 I

5.00

20.00 1

2.

1 1

10.999998

547

0.5588

0.4000

0.0000

0.0000

‘1.00

20.00 I I

2.

2 I I

14.999990

603

0.5266

0.4000

0.0000

0.0000