- Email: [email protected]
A numerical method for the nonlinear oscillator problem
6 March 1998
Chemical Physics Letters 284 Ž1998. 359–362
A numerical method for the nonlinear oscillator problem J. Killingbeck, G. Jolicard UniÕersite´ de Franche Comte, ´ ObserÕatoire de Besanc¸on, 41 bis AÕenue de l’ObserÕatoire, BP 1615, 25010Besanc¸on, France Received 24 November 1997
Abstract The treatment of a nonlinear Schrodinger equation with power law potential terms by means of hypervirial perturbation ¨ theory ŽHVPT. is considered. Previous workers have tried to handle the nonlinearity by constructing an HVPT which also contains a nonlinear term. We show that higher numerical accuracy can be obtained by reverting to the usual linear HVPT in combination with a simple numerical procedure. The procedure also works with finite-difference shooting calculations, which would permit the calculations to be extended to handle more general potentials. q 1998 Published by Elsevier Science B.V.
The interaction of a quantum mechanical system with its surroundings is often represented by a mathematical model in which the effects of action and reaction are simulated by introducing a nonlinear term into the Schrodinger equation of the system ¨ w1–3x. The nonlinearity is comparitively weak within the context of nonlinear theory; the Schrodinger ¨ equation is modified to take the form H0 q l² B : A c s Ec ,
Ž 1.
where H0 is the Hamiltonian of the isolated system, A and B are functions of the system coordinates and ² B : is the expectation value of B for the eigenfunction c . This problem apparently calls for some kind of self-consistent approach to its solution. Several authors have proposed perturbation formalisms for the nonlinear eigenvalue problem Ž1.; Surjan and Angyan w4x proposed a modified form of traditional Rayleigh–Schrodinger theory and Cioslowski w5x de¨ scribed a steepest descent form of perturbation theory Žfor the ground state only.. Two works w6,7x have applied hypervirial perturbation theory ŽHVPT. to
numerical calculations for the case in which H describes a harmonic oscillator and A and B are even powers of the coordinate; yd 2crd x 2 q x 2c q l² x 2 P : x 2 Qc s Ec .
Ž 2.
Surprisingly, the results obtained in w6,7x were of much lower accuracy than would be expected for the case Q s 2, for which the standard renormalized HVPT is known to give highly accurate results in the case of the linear Schrodinger equation w8x. Inspec¨ tion of w6,7x shows that in both works the expectation value perturbation series ² x N : s Ý X Ž N, J . l J
Ž 3.
J
was introduced into both factors of the product ² x 2 P : x 2 Q , so that the usual hypervirial recurrence relations acquired an extra nonlinear term involving sums of products of the X coefficients. It is this apparently ‘obvious’ way of handling the nonlinearity which we believe gives the reduction in numerical accuracy. We show here that the use of the usual
0009-2614r98r$19.00 q 1998 Published by Elsevier Science B.V. All rights reserved. PII S 0 0 0 9 - 2 6 1 4 Ž 9 7 . 0 1 4 4 1 - 3
J. Killingbeck, G. Jolicardr Chemical Physics Letters 284 (1998) 359–362
360
linear HVPT together with a simple numerical procedure leads to highly accurate results. Further, our numerical approach will work in conjunction with non-perturbative methods; we give some illustrative results for the finite-difference method described previously w9,10x. The method is straightforward; we simply solve the linear eigenvalue problem yd 2crd x 2 q x 2c q m x 2 Qc s Ec
Ž 4.
for a selected state and form the ratio mr² x 2 P :, varying m until the ratio equals l. Both linear HVPT and the finite-difference method of w9,10x can compute not only E but also ² x 2 P : and mr² x 2 P : quickly and accurately, and only ten or so repetitions of the calculation are needed to home in on the required m using any simple root finding procedure. HVPT has the advantage that it can pick out any desired state number n; it has been described many times in the literature Že.g., w8,11x. and so we give only a brief summary here. The Schrodinger Eq. Ž4. ¨ can be written in the equivalent form yd 2crd x 2 q b x 2c q m Ž x 2 Q y Kx 2 . c s Ec ,
Ž 5.
with b s 1 q K m. As K Žand b . vary the total Hamiltonian remains invariant but the potential is redistributed between an unperturbed part b x 2 and an effective perturbation m Ž x 2 Q y Kx 2 .. This serves to improve the behaviour of the Rayleigh–Schrodi¨ nger perturbation series for the energy and the expectation values ² x N :. To calculate these series we use the hypervirial relations, which can be derived either from operator commutation relations or by first principles using integration by parts w8,11x. For a Schrodinger equation of the form ¨ yd 2crd x 2 q Ý VJ x Jc s Ec
Ž 6.
J
these relations become Žfor N s 1, 2, 3, etc.. 1
Ž 2 N q 2 . E² x N : q N Ž N 2 y 1 . ² x Ny2 : 2
s Ý Ž 2 N q 2 q J . VJ ² x Nq J :
Ž 7.
J
The potential in Eq. Ž5. has even parity, so only even integers N are relevant in Eq. Ž7., all ² x N : for
odd N being zero. To apply HVPT we introduce the perturbation series ² x 2 N : s Ý m M X Ž N, M . ,
Ž 8.
E s Ý mMEŽ M .
Ž 9.
into Eq. Ž7., taking the coefficients V2 and V2 Q from Eq. Ž5.. Comparing coefficients of the m M term gives the recurrence relation
b Ž N q 1 . X Ž N q 1, M . sN N 2y
ž
1 4
/
X Ž N y 1, M .
q K Ž N q 1 . X Ž N q 1, M y 1 .
ž ž
q Nq
y Nq
1 2
M
/
Ý E Ž J . X Ž N, M y J . 0
1
Q q
2
2
/
X Ž N q Q, M y 1 . .
Ž 10 .
To find the energy coefficients EŽ M . we use the Hellman–Feynman equation appropriate to Eq. Ž5. EE Em
s ² x 2 Q y Kx 2 : ,
Ž 11 .
which leads to
Ž M q 1 . E Ž M q 1 . s X Ž Q, M . y KX Ž 1, M . . Ž 12 . Starting from the initial values X Ž 0, 0 . s 1;
'
E Ž 0. s Ž 2 n q 1. b
Ž 13 .
which are appropriate to the state number n for the unperturbed potential b x 2 , the recurrence relations Ž10. and Ž11. lead to the perturbation series for the energy and the ² x 2 N :. In the non-linear HVPT w6,7x the quantity m , which is a number here, was given the full form l² x 2 P : taken directly from Eq. Ž2.. This led to a term ² x 2 P :² x 2 QqJ : in Eq. Ž7. and thus to a double sum of X Ž J, M . X Ž K, N . type products when the perturbation series were used. In the linear HVPT used here we replace this attempt to incorporate the non-linearity into the perturbation equations by using the more simple and Žas it turns out. more accurate linear HVPT together with a
J. Killingbeck, G. Jolicardr Chemical Physics Letters 284 (1998) 359–362
search for the m value which makes m equal to l² x 2 P : numerically. The adjustable parameter K is the essential ingredient which ensures high accuracy for the case Q s 2 w8x. Rather than giving a plethora of tables of not very informative data we shall concentrate here on the most difficult case studied in w6,7x, since this suffices to demonstrate the relative efficiency of our alternative method. We take the case of the potential term x 2 q 20000² x 2 P :x 4 . For such an enormous perturbation it might appear that perturbation theory is quite useless; in fact, the use of a K value of around 1.5 and terms up to order 50 leads to very good results for Q s 2. The nonlinear HVPT of w6,7x is much less accurate, even when aided by Pade´ extrapolation. Our HVPT series sums for E and ² x 2 P : are already quite accurate Žat Q s 2. and can be further improved using the flexible Wynn algorithm described previously w12x. In the search process the ratio mr² x 2 P : can be formed directly as the E and ² x 2 P : sums are being computed, thus producing a sequence of mr² x 2 P : estimates which itself can be put through the Wynn analysis to yield an accurate ratio at each m value. m is then varied to make this ratio equal to 20000 for the particular calculations reported here. As an example of the accuracy attainable we note that for the state n s 10 the HVPY calculation for P s 1 and Q s 2 allows the m value to be obtained as 4060.173841578; the m value can be reconstructed from our tabulated results for that case, of course, by forming the product 20000² x 2 :. Table 1 compares our results with those of w6x for the ground state for P s 1 to 3 and then continues up
Table 1 Ground-state HVPT energies for the nonlinear case y D 2 q x 2 q l² x 2 P : x 4 with l s 20000 P
E
Ref. w7x
1 2 3 4 5 6
9.8150273496 6.2911985466 5.0298320135 4.4568416090 4.1705911121 4.0274609206
9.815026"4 6.29118"6 5.028"6
The nonlinear HVPT results of Ref. w7x are also shown, with the quoted error in the last digit.
361
Table 2 Results for the case y D 2 q x 2 q l² x 2 : x 4 with l s 20000 n
E
² x2:
0 1 2 3 4 5 6 7 8 9 10
9.815027349Ž5721. 44.0982711Ž899269. 93.73913798Ž77320. 154.81565370Ž07275. 225.399808377Ž7721. 304.35363705Ž33539. 390.864881506Ž3403. 484.32158638Ž87631. 584.240441533Ž0399. 690.22711322Ž41641. 801.95183779Ž4505.
0.03918003388Ž356. 0.07774917652Ž911. 0.0990581218Ž9978. 0.1172296580Ž7140. 0.1329089159Ž6697. 0.1469307047Ž0319. 0.15972787253Ž205. 0.17157406067Ž384. 0.18265406380Ž837. 0.19309965380Ž694. 0.20300869207Ž884.
The extra digits in brackets illustrate the typical extra accuracy obtained Žusing 16 decimal digit precision. on switching from linear HVPT to the fourth-order finite-difference method of Refs. w9,10x.
to P s 6 to show the extra range of our method. Table 2 gives results for the first 11 states for the case P s 1, Q s 2, showing also the extra accuracy which is possible if the method used to compute E and ² x 2 : is the finite-difference method of w9,10x. Only the ground state was treated in the nonlinear HVPT calculations of w6,7x. The finite-difference recurrence relation used to propagate the wavefunction along the x axis with steplength h is a simple three-term one of the form
c Ž x q h. s 2 q F Ž x , E . c Ž x . y c Ž x y h. ,
Ž 14 . where F is a polynomial in h 2 Ž V Ž x . y E . for the potential V and the trial energy E. The choice of F determines which of various standard integrators is obtained and Eq. Ž14. can be combined with an internal differentiation technique which produces accurate expectation values ²U : for any desired function UŽ x . w9,10x. The three-term form of Eq. Ž14. means that in a shooting process only three c values need to be stored at any step, so that scaling can be used to overcome any problems of wavefunction divergence; the resulting method is very stable, being able to integrate even through a sequence of thick potential barriers w10x. For the case of our test potential m x 4 the range of integration needed is small, for example, y5 - x - 5 is quite sufficient to include the region of the x axis where the wavefunction is at
362
J. Killingbeck, G. Jolicardr Chemical Physics Letters 284 (1998) 359–362
all appreciable in magnitude. For the case Q s 3 the HVPT is less accurate, giving only a few digits of E and ² x 2 :; nevertheless the HVPT’s ability to pick out a specified state makes it useful in giving a good initial value to start off the more accurate finite-difference calculation. For more general cases, of course, the finite-difference approach comes into its own; for a one-dimensional problem it will work for any choice of the functions A and B in the general nonlinear Eq. Ž1.. In our calculations the F polynomial used in Eq. Ž14. was that of Killingbeck and Jolicard w9x, which gives a fourth-order method more accurate than the traditional Numerov method and produces highly accurate E and ²UŽ x .: values when used with Richardson extrapolation. If the finite-difference method has to be used without preliminary HVPT results then the state number is picked out by
counting the number of nodes in the wavefunction as the shooting process proceeds along the x axis. References w1x w2x w3x w4x w5x w6x w7x w8x w9x w10x
S. Yomosa, J. Phys. Soc. Jpn. 35 Ž1972. 1738. M.D. Newton, J. Chem. Phys. 58 Ž1973. 5833. J.G. Kirkwood, J. Chem. Phys. 2 Ž1934. 351. P.R. Surjan, J. Angyan, Phys. Rev. A 28 Ž1983. 45. J. Cioslowski, Chem. Phys. Lett. 137 Ž1987. 78. Manda, Vishwamittar, Chem. Phys. Lett. 232 Ž1995. 35. E.R. Vrscay, J. Math. Phys. 29 Ž1988. 901. J. Killingbeck, J. Phys. A 20 Ž1987. 601. J. Killingbeck, G. Jolicard, Phys. Lett. A 172 Ž1993. 313. J. Killingbeck, N.A. Gordon, M.R.M. Witwit, Phys. Lett. A 206 Ž1995. 279. w11x J. Killingbeck, Microcomputer Algorithms, Hilger, Bristol, 1991. w12x J. Killingbeck, Phys. Lett. A 132 Ž1988. 223.
Chemical Physics Letters 284 Ž1998. 359–362
A numerical method for the nonlinear oscillator problem J. Killingbeck, G. Jolicard UniÕersite´ de Franche Comte, ´ ObserÕatoire de Besanc¸on, 41 bis AÕenue de l’ObserÕatoire, BP 1615, 25010Besanc¸on, France Received 24 November 1997
Abstract The treatment of a nonlinear Schrodinger equation with power law potential terms by means of hypervirial perturbation ¨ theory ŽHVPT. is considered. Previous workers have tried to handle the nonlinearity by constructing an HVPT which also contains a nonlinear term. We show that higher numerical accuracy can be obtained by reverting to the usual linear HVPT in combination with a simple numerical procedure. The procedure also works with finite-difference shooting calculations, which would permit the calculations to be extended to handle more general potentials. q 1998 Published by Elsevier Science B.V.
The interaction of a quantum mechanical system with its surroundings is often represented by a mathematical model in which the effects of action and reaction are simulated by introducing a nonlinear term into the Schrodinger equation of the system ¨ w1–3x. The nonlinearity is comparitively weak within the context of nonlinear theory; the Schrodinger ¨ equation is modified to take the form H0 q l² B : A c s Ec ,
Ž 1.
where H0 is the Hamiltonian of the isolated system, A and B are functions of the system coordinates and ² B : is the expectation value of B for the eigenfunction c . This problem apparently calls for some kind of self-consistent approach to its solution. Several authors have proposed perturbation formalisms for the nonlinear eigenvalue problem Ž1.; Surjan and Angyan w4x proposed a modified form of traditional Rayleigh–Schrodinger theory and Cioslowski w5x de¨ scribed a steepest descent form of perturbation theory Žfor the ground state only.. Two works w6,7x have applied hypervirial perturbation theory ŽHVPT. to
numerical calculations for the case in which H describes a harmonic oscillator and A and B are even powers of the coordinate; yd 2crd x 2 q x 2c q l² x 2 P : x 2 Qc s Ec .
Ž 2.
Surprisingly, the results obtained in w6,7x were of much lower accuracy than would be expected for the case Q s 2, for which the standard renormalized HVPT is known to give highly accurate results in the case of the linear Schrodinger equation w8x. Inspec¨ tion of w6,7x shows that in both works the expectation value perturbation series ² x N : s Ý X Ž N, J . l J
Ž 3.
J
was introduced into both factors of the product ² x 2 P : x 2 Q , so that the usual hypervirial recurrence relations acquired an extra nonlinear term involving sums of products of the X coefficients. It is this apparently ‘obvious’ way of handling the nonlinearity which we believe gives the reduction in numerical accuracy. We show here that the use of the usual
0009-2614r98r$19.00 q 1998 Published by Elsevier Science B.V. All rights reserved. PII S 0 0 0 9 - 2 6 1 4 Ž 9 7 . 0 1 4 4 1 - 3
J. Killingbeck, G. Jolicardr Chemical Physics Letters 284 (1998) 359–362
360
linear HVPT together with a simple numerical procedure leads to highly accurate results. Further, our numerical approach will work in conjunction with non-perturbative methods; we give some illustrative results for the finite-difference method described previously w9,10x. The method is straightforward; we simply solve the linear eigenvalue problem yd 2crd x 2 q x 2c q m x 2 Qc s Ec
Ž 4.
for a selected state and form the ratio mr² x 2 P :, varying m until the ratio equals l. Both linear HVPT and the finite-difference method of w9,10x can compute not only E but also ² x 2 P : and mr² x 2 P : quickly and accurately, and only ten or so repetitions of the calculation are needed to home in on the required m using any simple root finding procedure. HVPT has the advantage that it can pick out any desired state number n; it has been described many times in the literature Že.g., w8,11x. and so we give only a brief summary here. The Schrodinger Eq. Ž4. ¨ can be written in the equivalent form yd 2crd x 2 q b x 2c q m Ž x 2 Q y Kx 2 . c s Ec ,
Ž 5.
with b s 1 q K m. As K Žand b . vary the total Hamiltonian remains invariant but the potential is redistributed between an unperturbed part b x 2 and an effective perturbation m Ž x 2 Q y Kx 2 .. This serves to improve the behaviour of the Rayleigh–Schrodi¨ nger perturbation series for the energy and the expectation values ² x N :. To calculate these series we use the hypervirial relations, which can be derived either from operator commutation relations or by first principles using integration by parts w8,11x. For a Schrodinger equation of the form ¨ yd 2crd x 2 q Ý VJ x Jc s Ec
Ž 6.
J
these relations become Žfor N s 1, 2, 3, etc.. 1
Ž 2 N q 2 . E² x N : q N Ž N 2 y 1 . ² x Ny2 : 2
s Ý Ž 2 N q 2 q J . VJ ² x Nq J :
Ž 7.
J
The potential in Eq. Ž5. has even parity, so only even integers N are relevant in Eq. Ž7., all ² x N : for
odd N being zero. To apply HVPT we introduce the perturbation series ² x 2 N : s Ý m M X Ž N, M . ,
Ž 8.
E s Ý mMEŽ M .
Ž 9.
into Eq. Ž7., taking the coefficients V2 and V2 Q from Eq. Ž5.. Comparing coefficients of the m M term gives the recurrence relation
b Ž N q 1 . X Ž N q 1, M . sN N 2y
ž
1 4
/
X Ž N y 1, M .
q K Ž N q 1 . X Ž N q 1, M y 1 .
ž ž
q Nq
y Nq
1 2
M
/
Ý E Ž J . X Ž N, M y J . 0
1
Q q
2
2
/
X Ž N q Q, M y 1 . .
Ž 10 .
To find the energy coefficients EŽ M . we use the Hellman–Feynman equation appropriate to Eq. Ž5. EE Em
s ² x 2 Q y Kx 2 : ,
Ž 11 .
which leads to
Ž M q 1 . E Ž M q 1 . s X Ž Q, M . y KX Ž 1, M . . Ž 12 . Starting from the initial values X Ž 0, 0 . s 1;
'
E Ž 0. s Ž 2 n q 1. b
Ž 13 .
which are appropriate to the state number n for the unperturbed potential b x 2 , the recurrence relations Ž10. and Ž11. lead to the perturbation series for the energy and the ² x 2 N :. In the non-linear HVPT w6,7x the quantity m , which is a number here, was given the full form l² x 2 P : taken directly from Eq. Ž2.. This led to a term ² x 2 P :² x 2 QqJ : in Eq. Ž7. and thus to a double sum of X Ž J, M . X Ž K, N . type products when the perturbation series were used. In the linear HVPT used here we replace this attempt to incorporate the non-linearity into the perturbation equations by using the more simple and Žas it turns out. more accurate linear HVPT together with a
J. Killingbeck, G. Jolicardr Chemical Physics Letters 284 (1998) 359–362
search for the m value which makes m equal to l² x 2 P : numerically. The adjustable parameter K is the essential ingredient which ensures high accuracy for the case Q s 2 w8x. Rather than giving a plethora of tables of not very informative data we shall concentrate here on the most difficult case studied in w6,7x, since this suffices to demonstrate the relative efficiency of our alternative method. We take the case of the potential term x 2 q 20000² x 2 P :x 4 . For such an enormous perturbation it might appear that perturbation theory is quite useless; in fact, the use of a K value of around 1.5 and terms up to order 50 leads to very good results for Q s 2. The nonlinear HVPT of w6,7x is much less accurate, even when aided by Pade´ extrapolation. Our HVPT series sums for E and ² x 2 P : are already quite accurate Žat Q s 2. and can be further improved using the flexible Wynn algorithm described previously w12x. In the search process the ratio mr² x 2 P : can be formed directly as the E and ² x 2 P : sums are being computed, thus producing a sequence of mr² x 2 P : estimates which itself can be put through the Wynn analysis to yield an accurate ratio at each m value. m is then varied to make this ratio equal to 20000 for the particular calculations reported here. As an example of the accuracy attainable we note that for the state n s 10 the HVPY calculation for P s 1 and Q s 2 allows the m value to be obtained as 4060.173841578; the m value can be reconstructed from our tabulated results for that case, of course, by forming the product 20000² x 2 :. Table 1 compares our results with those of w6x for the ground state for P s 1 to 3 and then continues up
Table 1 Ground-state HVPT energies for the nonlinear case y D 2 q x 2 q l² x 2 P : x 4 with l s 20000 P
E
Ref. w7x
1 2 3 4 5 6
9.8150273496 6.2911985466 5.0298320135 4.4568416090 4.1705911121 4.0274609206
9.815026"4 6.29118"6 5.028"6
The nonlinear HVPT results of Ref. w7x are also shown, with the quoted error in the last digit.
361
Table 2 Results for the case y D 2 q x 2 q l² x 2 : x 4 with l s 20000 n
E
² x2:
0 1 2 3 4 5 6 7 8 9 10
9.815027349Ž5721. 44.0982711Ž899269. 93.73913798Ž77320. 154.81565370Ž07275. 225.399808377Ž7721. 304.35363705Ž33539. 390.864881506Ž3403. 484.32158638Ž87631. 584.240441533Ž0399. 690.22711322Ž41641. 801.95183779Ž4505.
0.03918003388Ž356. 0.07774917652Ž911. 0.0990581218Ž9978. 0.1172296580Ž7140. 0.1329089159Ž6697. 0.1469307047Ž0319. 0.15972787253Ž205. 0.17157406067Ž384. 0.18265406380Ž837. 0.19309965380Ž694. 0.20300869207Ž884.
The extra digits in brackets illustrate the typical extra accuracy obtained Žusing 16 decimal digit precision. on switching from linear HVPT to the fourth-order finite-difference method of Refs. w9,10x.
to P s 6 to show the extra range of our method. Table 2 gives results for the first 11 states for the case P s 1, Q s 2, showing also the extra accuracy which is possible if the method used to compute E and ² x 2 : is the finite-difference method of w9,10x. Only the ground state was treated in the nonlinear HVPT calculations of w6,7x. The finite-difference recurrence relation used to propagate the wavefunction along the x axis with steplength h is a simple three-term one of the form
c Ž x q h. s 2 q F Ž x , E . c Ž x . y c Ž x y h. ,
Ž 14 . where F is a polynomial in h 2 Ž V Ž x . y E . for the potential V and the trial energy E. The choice of F determines which of various standard integrators is obtained and Eq. Ž14. can be combined with an internal differentiation technique which produces accurate expectation values ²U : for any desired function UŽ x . w9,10x. The three-term form of Eq. Ž14. means that in a shooting process only three c values need to be stored at any step, so that scaling can be used to overcome any problems of wavefunction divergence; the resulting method is very stable, being able to integrate even through a sequence of thick potential barriers w10x. For the case of our test potential m x 4 the range of integration needed is small, for example, y5 - x - 5 is quite sufficient to include the region of the x axis where the wavefunction is at
362
J. Killingbeck, G. Jolicardr Chemical Physics Letters 284 (1998) 359–362
all appreciable in magnitude. For the case Q s 3 the HVPT is less accurate, giving only a few digits of E and ² x 2 :; nevertheless the HVPT’s ability to pick out a specified state makes it useful in giving a good initial value to start off the more accurate finite-difference calculation. For more general cases, of course, the finite-difference approach comes into its own; for a one-dimensional problem it will work for any choice of the functions A and B in the general nonlinear Eq. Ž1.. In our calculations the F polynomial used in Eq. Ž14. was that of Killingbeck and Jolicard w9x, which gives a fourth-order method more accurate than the traditional Numerov method and produces highly accurate E and ²UŽ x .: values when used with Richardson extrapolation. If the finite-difference method has to be used without preliminary HVPT results then the state number is picked out by
counting the number of nodes in the wavefunction as the shooting process proceeds along the x axis. References w1x w2x w3x w4x w5x w6x w7x w8x w9x w10x
S. Yomosa, J. Phys. Soc. Jpn. 35 Ž1972. 1738. M.D. Newton, J. Chem. Phys. 58 Ž1973. 5833. J.G. Kirkwood, J. Chem. Phys. 2 Ž1934. 351. P.R. Surjan, J. Angyan, Phys. Rev. A 28 Ž1983. 45. J. Cioslowski, Chem. Phys. Lett. 137 Ž1987. 78. Manda, Vishwamittar, Chem. Phys. Lett. 232 Ž1995. 35. E.R. Vrscay, J. Math. Phys. 29 Ž1988. 901. J. Killingbeck, J. Phys. A 20 Ž1987. 601. J. Killingbeck, G. Jolicard, Phys. Lett. A 172 Ž1993. 313. J. Killingbeck, N.A. Gordon, M.R.M. Witwit, Phys. Lett. A 206 Ž1995. 279. w11x J. Killingbeck, Microcomputer Algorithms, Hilger, Bristol, 1991. w12x J. Killingbeck, Phys. Lett. A 132 Ž1988. 223.