__ __ BB
16 October 1995
2%
PHYSICS
EISEVIER
LETTERS A
Physics Letters A 206 (1995) 279-282
Stable forward shooting for eigenvalues and expectation values J.P. Killingbeck, N.A. Gordon, M.R.M. Witwit Department of Applied Mathematics, University of Hull, Hull HlJ6 7RX, UK
Received 24 July 1995; accepted for publication 16 August 1995 Communicated by P.R. Holland
Abstract Internal differentiation techniques are used to produce a simple but highly accurate forwards shooting method for eigenvalues and expectation values of the Schriidinger equation. A multi-well potential is used as a test case.
The most simple finite difference Schriidinger equation -D2’P+
approach to the
W=W
(1)
is to use the lowest order finite difference representation of the second derivative term to produce the three-term recurrence relation !P(x+h)
= [2+G(x)]?P((x)
- !.P((n-h),
(2)
where G(x)
= [V(x)
-E]h2
(3)
and where a constant step length h is used in the finite difference calculation. Eq. (2) can be treated as a tridiagonal matrix eigenvalue problem [ 1,2] or it can be handled using shooting methods [3-51. Two main traditions have developed in the search for improvements to Eq. (2). The first uses higher order finite difference representations of the second derivative and leads to matrix eigenvalue problems of increasing bandwidth [6,7]; the second sticks to the three-term recurrence relation but uses more complicated functions of G(x) in the square bracket [8,9]. In the second tradition shooting methods are preferred, since the resulting tridiagonal matrix problem
would involve E in a non-linear manner; nevertheless, matrix formulations of the popular Numerov method have been produced by several workers [lo121, usually leading to generalized rather than ordinary matrix eigenvalue problems. The present work presents an efficient synthesis and improvement of ideas developed in several previous works [3,4,8,9,13,14] leading to a simple finite difference method with the following principal features: (1) It uses an incremental form, which has been found by computational experiment to give more stable results than the direct use of the recurrence relation (2) relating the ?P values. (2) It uses internal differentiation [ 13,141 to compute expectation values. This avoids explicit storage of the wavefunction and makes possible a scaling operation during the shooting process. This scaling enables forward shooting to be used even in extreme circumstances; an example is given in which accurate calculations were carried out for three potential wells separated by thick barriers, with virtually no tunnelling taking place through the barriers. In general, forward shooting is not reliable for computing the wavefunction accurately, particularly near the
0375-9601/95/$09.50 0 1995 Elsevier Science B.V. All rights reserved SSDI 0375-9601(95)00632-X
J. P. Killingbeck
280
et al. / Physics Letters A 206 (I 99.5) 279-282
end of the shooting region; however, the approach outlined here gives accurate energies and expectation values without storing a wavefunction. The incremental form of the finite difference three-term recurrence relation which we use takes the form A(x)
=P(x)?P(x)
‘P((x+h)
= ?P(x)
+A(x-h), +A(x).
(4) (5)
Re-writing some of the results of Ref. [14] in this incremental form gives for the simple method of Eq. (2) the result
rithms in which all calculated quantities depend only on ratios which are not affected by a joint scaling of all the relevant variables. The internal differentiation technique [9,13,14] is appropriate to produce this kind of algorithm. For example, to apply Newton’s root finding method to find eigenvalues from the shooting method would require the evaluation of the quantity ps, the partial derivation of T with respect to E. Using Eqs. (3)-(5) leads to the result As(x)
=P(x)qE(x) +&(x--h),
‘&(x+h) P(X)
=G(x)
and for the Numerov
(6)
= ‘J’s(x)
where the derivative
(9) +&(x)9 P,(x)
(10)
is given by
method the result PE(X)
P(X)
+PE(x)q(x)
= 12G(x)[12-G(x)]-‘.
= 42;
(11)
(7)
Ref. [4] developed a finite difference method which, while giving an eigenvalue error of leading order h4, like the Numerov method, is much more accurate than that method, particularly for excited states. The method is based on an angle-doubling formula which is equally valid for both the cos and cash functions appropriate to the classical and non-classical regions traversed during the shooting process. Algebraic manipulation of the multi-step algorithm originally presented in Ref. [4] enables the incorporation of that calculation in the present incremental formalism. The relevant P function is found to be (in nested form)
and is written as an explicit function of G for the particular P which is being used in Eq. (4). Clearly, to work out eigenvalues using Newton’s method requires knowledge of the ratio q/qE at the endpoint (x = L) at which we set q(L) = 0. Propagation of ‘ly and !?E throughout the shooting is accomplished using the recurrence relations given above, scaling them both together whenever the wavefunction threatens to diverge. To work out the expectation value of a given function U(x) we introduce a dummy parameter p and suppose that the potential takes the form V(x) + pU(x). The expectation value of V( x> is then given exactly by the expression [ 131
(12)
q1+$y1+Jy)]).
(8)
For the record, we note that a direct four-term Taylor series expansion approach to the original S&r& dinger equation (1) would lead to a P function with the coefficients 32 and 96 instead of 30 and 56, provided that all derivatives of V(x) are set equal to zero. Calculations indicate that (8) gives the more accurate results, presumably because it gives a more uniformly accurate approximation to the exact local constant potential solutions over a region involving both classical and non-classical portions. To permit scaling of the wavefunction ?? and the increment A during the shooting process, we must devise algo-
(13) if we suppose that our eigenvalue criterion is F(E) = 0. The derivative *E has already been treated above. A little algebra using (3)-(5) gives the results for the derivatives with respect to F, A,( x> = P( x) qp:( x) + P,( x> *( x> + A,( x - h), (14) P$x+h)
= ‘$(.)
+A+),
(15)
with P,(x)
=hW(x)~.
(16)
J. P. Killingbeck
et al. /Physics
The equations established so far thus enable the propagation of ?P,tY, and PK, with joint scaling throughout the shooting process, so that E and (U) can be calculated using the ratios of these quantities at the endpoint x = L. To start off the calculation (at x = 0, for example) assuming that ?P= 0 is required, we set the initial values of 9, ‘1yE, Y;, A, A,, A,, at x = 0 and x = h equal to zero, with the sole exceptions P(h) = A(0) = h. The incremental form of the shooting equations means that, provided the quantities at each step are worked out in the order qfi, !PE, ?P, only one computer storage element is required for each of the six quantities involved, with the new values overwriting the old ones. The calculation involves two specified functions; V(x) is the potential function and U(x) is the function for which the expectation value is required. By using appropriate forms of U(X) it is possible to calculate the properly normalized 9’ value at a point or the integrated probability over an extended region, even though no quadratures are performed and the wavefunction is not stored [14]. Counting the number of sign changes in !P during the shooting process indicates the number of the excited state to which the calculated eigenvalue belongs. The formulae developed above refer to the use of a fixed step length h. This leads to particularly simple algorithms and programs, with the added benefit that the calculated quantities can be improved systematically by the use of Richardson extrapolation [1,3]. Variable step methods [15,16] often lose this capability by using local criteria for adjusting the step length throughout the shooting process. As a test calculation we treated the region 0 G x < 36rr* for the potential V = cod&>, which is a special case of the class of slowly oscillating potentials studied by Stolz [17]. We did this in order to provide a severe test for our method; it is clear that the method itself will be generally applicable to the one-dimensional Schriidinger equation. Of finite difference methods which have previously been used for potentials with two or more wells (see, e.g., Refs. [18-201) the renormalized Numerov method of Johnson [ 181 is perhaps the closest in spirit to our approach, although it is not as accurate and does not handle expectation values directly. Our test potential has three wells. At the bottom of each well (i.e. at x = r2, 9rr2, 25~~) the potential can bc expanded
Letters A 206 (1995) 279-282
281
Table 1 The lowest 20 levels for the cos(&) potential, showing the truncated ( x> value and asociated well number; the well centres are at 9.870, 88.826 and 246.740 Level
-E
0.977554578327 0.962655726436 0.932793527194 0.8884 I 1244570 0.888339927083 0.888289070323 0.844043472771 0.8 1474326595 1 0.800059067405 0.756338258322 0.741876660414 0.712883524502 0.669818342815 0.66975 1597346 0.669697423737 0.626782596842 0.598380183821 0.584141772192 0.54 1777770605 0.527775204315
0
1 2 3 4 5 6 7 8 9 10 I1 12 13 14 15 16 17 18 19
lx>
N
246.808 88.939 246.944 10.214 89.169 247.082 247.222 89.402 247.364 247.507 89.641 247.652 10.931 89.885 247.799 247.948 90. I34 248.098 248.25 1 90.389
3 2 3
I 2 3 3 2 3 3 2 3
I 2 3 3 2 3 3 2
locally as a harmonic oscillator potential, leading to a set of energy levels. A little algebra shows that within this crude local harmonic oscillator approximation the energy level number n ( = 0, 1, 2, . . . > in the well number N (= 1, 2, 3) is given by the formula I E(n,
N)=
-1
+--
2n+l
(17)
nJ82N-1’
This model predicts the existence of groups of degenerate levels, e.g., E(0, 11, E(2, 21, E(3, 3) should all equal roughly -0.8875. The calculated energy levels (Table 1) show a tight group of three levels centred around - 0.8883 and in general formula (17) is a good guide to the location of the levels. The expectation values (x) in Table 1, obtained by setting U(x) = x in the finite difference formalism, indicate clearly the well within which the wavefunction for each level is predominantly located. Highly accurate results were obtained by using the strip numbers N,, kN,, k2No,. . . with k = 3/2, followed by Richardson extrapolation based on an error law of the form E(h)=E+Ah4+Bh6+Chs...
.
(18)
282
J.P. Killingbeck et al. / Physics Letters A 206 (199.5) 279-282
for the eigenvalue computed using step length h. A similar formula holds for the (x) values also. For the test potential used an initial step length around 1 was found to be adequate to put the calculation in the asymptotic region where (18) applies. For three wells, with x = 367r2, an N, value of 352 was used; in general, choosing N, as a multiple of 32 allows up to six calculations to be performed (if necessary) before Richardson extrapolation. The remarkable feature of the calculation is that forward shooting starting from x = 0 manages to find E and (x) even for eigenfunctions separated from the origin by one or two thick potential barriers. This indicates the effectiveness of combining wavefunction scaling with the internal differentiation technique. As a check on the three-well calculations we also performed some calculation using only one well at a time. The results confirmed that tunnelling effects for the low energy levels are extremely small. For example, the energies of levels 12, 13, and 14 were unchanged (to the digits quoted) even when the calculation of each level was restricted to the appropriate single well (indicated by the (x) value from the three-well calculation). A further checking calculation was performed for the case of the first well on its own. The potential was expanded in a power series about the minimum at x = yr2 and the energy levels and ( x) values were then computing using the hypervirial perturbation method [21,22]. Although the accuracy of the hypervirial method is decreased by the presence of odd powers of x in the potential which perturbs the harmonic oscillator potential, the calculation serves to give a good estimate of the corrections to the crude formula (17) for the special case N = 1. Table 2 shows some results which compare perturbation and finite difference results and also show the magnitude of the corrections to formula (17) for the first well. Although the present calculations were carried out using the P function of Eq. (81, it will be clear that the formalism described can be used with the Numerov P function of Eq. (7) or with the P function appropriate to other finite difference approximations.
Table 2 Hypervirial perturbation results for the first well of the potential cos(&/;, (AE= E, - E, where E is given by Eq. (17). (x) is determined bv the hvuervirial method) -E
finite difference
hypervirial
0.88841124 0.6698 I834 0.45795018 0.253 10674 0.05544925
0.88841202 0.66982865 0.45803321 0.2533933 0.05635 I
I
IAEI
(x)
0.00095078 0.00743694 0.026927 18 0.04883539 0.0683051 I
10.214 10.930 11.695 12.517 13.408
The present work was carried out as part of a research project investigating Sturm-Liouville problems; the authors wish to thank EPSR for financial support under research grant No. GR/K17750.
References (11 D.G. Truhlar. J. Comput. Phys. 10 (1972) 123. [2] P.J. Cooney, E.P. Kanter and 2. Vager, Am. J. Phys. 49 (1981) 76. [3] J. Killingbeck, Microcomputer algorithms (Hilger, Bristol, 1991). [4] J. Killingbeck and G. Jolicard, Phys. Lett. A 172 (1993) 313. [5] F.Y. Hajj, J. Phys. B 13 (1980) 4521. [6] V. Fack and G. Vanden Berghe, J. Phys. A 18 (1985) 3355. [7] G.C. Groenenboom and H.M. Buck, J. Chem. Phys. 92 ( 1990) 4374. [8] J. Killingbeck, Phys. Lett. A 1 I5 (1986) 301. [9] M.R.M. Witwit, J. Phys. A 25 (1992) 503. [IO] M.M. Chawla and C.P. Katti, Nord. Tidskr. Inform. Beh. 20 (1980) 107. [l I] B. Lindberg, J. Chem. Phys. 88 (1988) 3805. (121 B. Grieves and D. Dunn, J. Phys. A 23 (1990) 5479. [I31 J. Killingbeck, J. Phys. A 18 (1985) 245. 1141 J. Killingbeck, J. Phys. A 21 (1988) 3399. [I51 A.D. Raptis and J.R. Cash, Comput Phys. Commun. 36 (1985) 113. [16] T.E. Simos, J. Comput. Phys. 108 (1993) 175. [17] G. Stolz, Manus. Mathematics 84 (1994) 245. 1181 B.R. Johnson, J. Chem. Phys. 67 (1977) 4086. [19] L. Wolniewicz and T. Orlikowski, J. Comput. Phys. 27 (1978) 16. (201 A. Wierzbicki and J.M. Bowman, Int. J. Quantum Chem. 35 (1989) 297. [21] J. Killingbeck, J. Phys. A 14 (1981) 1005. [22] M.R.M. Witwit, Indian J. Phys. B 68 (1994) 515.