- Email: [email protected]
A new effective algorithm for the resonant state of a Schrödinger equation
Computer Physics Communications 167 (2005) 1–6 www.elsevier.com/locate/cpc
A new effective algorithm for the resonant state of a Schrödinger equation Zhongcheng Wang Department of Physics, Shanghai University, 99 ShangDa Road, Shanghai 200436, P.R. China Received 9 September 2004; received in revised form 7 December 2004; accepted 13 December 2004
Abstract In this paper we present a new effective algorithm for the Schrödinger equation. This new method differs from the original Numerov method only in one simple coefficient, by which we can extend the interval of periodicity from 6 to infinity and obtain an embedded correct factor to improve the accuracy. We compare the new method with the original Numerov method by the well-known problem of Woods–Saxon potential. The numerical results show that the new method has great advantage in accuracy over the original. Particularly for the resonant state, the accuracy is improved with four orders overall, and even six to seven orders for the highest oscillatory solution. Surely, this method will replace the original Numerov method and be widely used in various area. 2004 Elsevier B.V. All rights reserved. Keywords: Numerov method; Trigonometric fitting; P-stable; High-oscillatory solution; Schrödinger equation PACS: 02.60.Cb; 02.60.Lj; 02.70.Bf
1. Introduction Many physical problems occurred in various scientific areas such as theoretical physics, nuclear physics, physical chemistry, quantum chemistry, and molecular physics [1–4] can be reduced into an one-dimensional Schrödinger equation, which can be written as y (x) = f (x)y(x),
(1)
E-mail address: [email protected] (Z. Wang). 0010-4655/$ – see front matter 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.cpc.2004.12.004
with f (x) = v(x) − Es , where the function v(x) is given and denotes the effective potential and Es is a real number or complex denoting the energy. The earliest method to solve (1) numerically is the Numerov’s method. Although this method is only of algebraic order four, it is still widely used in many area, because of its simplicity and wider interval of periodicity. On the other hand, improving this method and seeking a more accurate, fast and reliable algorithm to achieve the numerical solution have been the subject of great activity for many years [5–18]. From the numerical solution to the well-known Woods–Saxon potential problem by
2
Z. Wang / Computer Physics Communications 167 (2005) 1–6
the available methods [9–18], it is interesting to find that the error in the calculated eigenvalue always increases with its value. Hence, to obtain the highest resonant state with high accuracy is a headache problem and has long been the motivation for developing a new method, such as the exponential fitting [9–12], the embedded methods [14], and the eighth-order formula [15], Obrechkoff method [16–18]. In order to obtain a high accurate solution, we must adapt a method with a high algebraic order which usually involves a complex algorithm. Is it possible to find a simple and effective way to obtain the high accurate solution? We find it is possible. The motivation of this paper is to report a simple and effective algorithm to improve the accuracy of the resonant state. This algorithm is nearly the same as the original Numerov method. We just use trigonometric fitting to extend the interval of periodicity of the traditional Numerov method to infinity. Then we use an embedded parameter to correct the error. This paper is arranged as follows: in the next section, we will give the derivation for the new method. In Section 3 we will analyze the stability and the local truncation error. In Section 4 we will test our new method by the Woods–Saxon potential. In the final section we will give a conclusion.
h2 (1 − 2c1 − c2 ) ω2 y(x) + y (x) h4 (1 − 12c1 ) ω4 y(x) − y (4) (x) + O(h6 ) = 0. 12 (5) From (5), we obtain the coefficients as
−
1 5 , c2 = . (6) 12 6 After substituting (4) and (6) into (2) and changing the notations such as y(x) → yk and y (x) → fk yk , we obtain the following iterative formula 5h2 αs + yk+2 = fk+1 yk+1 6 −1 2 h2 h − 1 − fk yk 1 − fk+2 , (7) 12 12 c1 =
where
5 2 2 1 2 2 αs = h ω + 2 1 + h ω cos(hω). 6 12
(8)
The local truncation error is LET(h) = −
h6 6 ω y(x) + y (6) (x) . 240
(9)
3. Stability and local truncation error 2. Derivation We consider the following structure which are the same for the Numerov method y(x + h) + y(x − h) − αs y(x) = h2 c1 y (x + h) + y (x − h) + c2 y (x) .
(2)
According to the stability theory by Lambert and Watson [6], the stability of (2) can be examined through the following test equation
y (x) = −ω y(x). 2
(3)
We apply (3) and its solutions to (2) and obtain the following equation αs = h2 ω2 c2 + 2 1 + h2 ω2 c1 cos(hω). (4) We use (4) to eliminate αs in (2) and convert the resultant expression to a Taylor series such as
It is easy to see that the traditional Numerov method is just the special case of (7) with αs = 2. We can show the interval of periodicity of (7) to be infinite by making the following substitutions: fk , fk+1 , fk+2 → −ω2 , and yk+2 → λ2 , yk+1 → λ and yk → 1, and hω → H to (7) and (8). The characteristic equation can be obtained as H2 2 λ − cos(H )λ + 1 = 0. 1+ (10) 12 For comparison, we let αs = 2 for (8) and obtain the characteristic equation 5H 2 H2 2 H2 − 2− λ+ 1+ λ 1+ (11) 12 6 12 for the traditional Numerov method. The characteristic roots are λ± = e±iH for (10), which is independent of H . So it is a P-stable and the interval of periodicity is infinite.
Z. Wang / Computer Physics Communications 167 (2005) 1–6
Although for the traditional Numerov method, the characteristic roots of (11) can √ be written as λ± = e±iθ(H ) , but it is only for H 6. So the corresponding interval of periodicity is 6. From (9), we can see that the local truncation error is ω dependent, so that it should be examined with more terms such as 1 6 6 h ω y(x) + y (6) (x) LTE(h) = − 240 11 8 8 h ω y(x) − y (8) (x) + 60480 13 − h10 ω10 y(x) + y (10) (x) 3628800 1 h12 ω12 y(x) − y (12) (x) + 23950080 17 h14 ω14 y(x) + y (14) (x) − 52306974720 + ···. (12) After rearranging, the above expression can be written as LTE(h) =
∞ (−1)n (n − 2)(2n + 3)h2n
3(2n)! 2n × ω y(x) − (−1)n y (2n) (x) . n=3
(13)
Because ω2n y(x) − (−1)n y (2n) (x) ≡ 0 for (3), the total local truncation error (3) is zero for the test equation. Therefore, the trigonometrically-fitted P-stable Numerov method not only extends the interval of periodicity to infinity but also provides an exact formula for the test-like problem. Can we use this property to improve the numerical accuracy for Schrödinger equation with periodic solution? We find the answer is “yes”. If we explain ω as the wave vector, which appears as a parameter in the characteristic root y(x) = e√iωx of the test equation (3), then we can let ω = |v(x) − Es |. Such that ω is position sensitive. If this assumption is correct, we can correct the error for every step by ω. We need to test our idea by the numerical test.
4. Numerical illustrations We will use the Woods–Saxon potential [15] v(x) = c0 z 1 − a(1 − z) (14)
3
Fig. 1. The Woods–Saxon potential defined by (14) is plotted in the range of x ∈ (0, 15).
with z = {exp[a(x − b)] + 1}−1 , c0 = −50, a = 5.0/3 and b = 7.0 to test our new method and compare it with the traditional Numerov method. The plot of the Woods–Saxon potential is given in Fig. 1. The test domain is 0 x 15, which is the same as in the literature [18]. The first problem is to find the eigenvalues in the bound state, for which −50 < Es < 0 and the √ boundary conditions are: y(0) = 0 and y(15) = e− Es Nk h/2 ≈ 0. The Schrödinger equation (1) is integrated forwards and backwards from the two boundary nodes, and an approximate eigenvalue can be found by the fact that the wave function and its derivative from the two sides at the matching point M are equal. − The connecting condition can be written as yM y˜M y˜M yM = 0, where yM and yM are obtained from are from the backward. the forward and y˜M and y˜M The matching√point is chosen to be x = 7.5. y0 = 0 and y1 = sin[ (50 + Es )h] are the initial wave functions for the forward integral; yNk = 0 and yNk −1 = √ − −E s (Nk −1)h/2 are for the backward, where N is e k the total number of steps, h is the step size. Like the other authors [15], the second problem is to find the real eigenvalues in the range 5 Es 1000 for which the phase shift δ(Es ) is just to π/2 in the resonant state with the same domain. The √ boundary conditions are: y(0) = 0 and y(15) = cos[ Es x]. The matching point√is also chosen to be x = 7.5. y0 = 0 and y1 = sin[ ES h] are the initial wave √ functions = cos( Es Nk h) and for the forward integral; y N k √ yNk −1 = cos( ES (Nk − 1)h) for the backwards. − We use the new method to calculate (yM y˜M m y˜M yM )/Es as a function of energy Es for the bound state with h = 1/32 and m = 6, and for the real res-
4
Z. Wang / Computer Physics Communications 167 (2005) 1–6
(a)
(b)
− y˜ y )/E 6 to energy E ∈ (−50, −5) with h = 1/32 and m = 6; (b) is the plot of (y y˜ − y˜ y ) to Fig. 2. (a) is the plot of (yM y˜M s M M M M M M s energy Es ∈ (5, 1005) with h = 1/64 and δ(Es ) = π/2.
Table 1 Error comparison between the new P-stable Numerov and the traditional Numerov method for the seven eigenvalues in the bound states [15] Ei
h
Es1
t (s)
Es2
t (s)
Es2 /Es1
−49.457788728
1/16 1/32 1/64 1/16 1/32 1/64 1/16 1/32 1/64 1/16 1/32 1/64 1/16 1/32 1/64 1/16 1/32 1/64 1/16 1/32 1/64
15 1 0 2.4 × 103 152 10 3.1 × 104 1.9 × 103 120 1.6 × 105 1.0 × 104 634 5.5 × 105 3.4 × 104 2.1 × 103 1.4 × 106 8.8 × 104 5.5 × 103 2.9 × 106 1.8 × 105 1.1 × 104
0.28 0.55 1.11 0.19 0.39 0.73 0.17 0.33 0.66 0.17 0.33 0.66 0.16 0.33 0.64 0.14 0.33 0.66 0.14 0.28 0.56
2 0 0 163 10 1 1.0 × 103 64 4 3.1 × 103 192 13 6.5 × 103 403 26 1.1 × 104 663 41 1.4 × 104 879 55
0.22 0.45 0.92 0.11 0.30 0.61 0.13 0.27 0.53 0.13 0.27 0.53 0.14 0.27 0.53 0.13 0.27 0.53 0.13 0.22 0.45
1 × 10−1 6 × 10−2 5 × 10−1 7 × 10−2 7 × 10−2 1 × 10−1 3 × 10−2 3 × 10−2 3 × 10−2 2 × 10−2 2 × 10−2 2 × 10−2 1 × 10−2 1 × 10−2 1 × 10−2 8 × 10−3 8 × 10−3 8 × 10−3 5 × 10−3 5 × 10−3 5 × 10−3
−46.290753954
−41.232607772
−34.672313205
−26.873448915
−18.094688282
−8.676081670
The true eigenvalues with nine decimal places are shown in the first column, the step-size with 1/16, 1/32 and 1/64 are in the second. The absolute deviations of the computed eigenvalues from the exact values, Es1 for Numerov’s method and Es2 for P-stable Numerov’s method, in 109 units and the consumed CPU time (in seconds) are shown in each sub-column. In the last column, the ratio of the absolute deviations of Es2 /Es1 .
onant states with h = 1/64, m = 0 and δ(Es ) = π/2, which are shown in Fig. 2 (a) and (b), respectively. From the figures, we find that there exists total nineteen eigenstates, thirteen for the bound states and six for resonant states, which are in agreement with those in [18]. In Table 1, we compare the absolute error, |Es | = Es − Ei for seven bound states with the step size of 1/16, 1/32, and 1/64, and the consuming CPU time.
In Table 2, we compare the four resonant states with the size of 1/32, 1/64, and 1/128. Where the true solution Ei to the Woods–Saxon potential problem were obtained by using analytic method which is correct to nine decimal places for the bound states and six decimal places for the resonant states [15]. In order to give an impression how great the improvement has been made to the resonant state, in Fig. 3 we give the plot of the log10 |Es × 106 |
Z. Wang / Computer Physics Communications 167 (2005) 1–6
5
Table 2 Error comparison between the new P-stable Numerov and the traditional Numerov method for the four eigenvalues in the resonant states [15] Ei 53.588872
163.215341
341.495874
989.701916
h 1/32 1/64 1/128 1/32 1/64 1/128 1/32 1/64 1/128 1/32 1/64 1/128
Es1 1.4 × 104 878 55 4.8 × 105 2.9 × 104 1.8 × 103 7.5 × 106 4.4 × 105 2.7 × 104 6.4 × 108 2.9 × 107 1.6 × 106
t (s)
Es2
t (s)
Es2 /Es1
0.23 0.47 0.95 0.19 0.48 0.97 0.19 0.38 0.75 0.19 0.38 0.75
3 1 0 8 1 0 15 2 1 49 6 1
0.14 0.38 0.75 0.16 0.31 0.55 0.27 0.31 0.63 0.41 0.53 0.94
2 × 10−4 6 × 10−4 3 × 10−3 2 × 10−5 5 × 10−5 3 × 10−4 2 × 10−6 5 × 10−6 4 × 10−5 8 × 10−8 2 × 10−7 1 × 10−6
The true eigenvalues, Ei , with six decimal places are shown in the first column, the step-size with 1/32, 1/64 and 1/128 are in the second. The absolute deviations of the computed eigenvalues from the exact values in 106 units and the consumed CPU time (in seconds) are shown in each sub-column. As the same as in Table 1, Es 1 is for Numerov’s method and Es 2 for P-stable Numerov’s method. In the last column, the ratio of the absolute deviations of |Es 2 /Es 1 |.
Fig. 3. Es = Es − Ei as a function of h is plotted for the new method, where Ei = 989.701916. In the figure, 10−6 and 10−7 are denoted for the plots of Es × 10−6 and Es × 10−7 obtained by the traditional Numerov method.
as a function of h for the last eigenvalue, Ei = 989.701916, obtained by the new method, and the plots of log10 |Es | and log10 |Es × 10−1 | by the Numerov method. From the tables and figures, it is easy to see the accuracy has been greatly improved. It is interesting to note (1) Es2 /Es1 decreases with the value of the eigenvalue. It seems that the larger Es is the higher the accuracy is. This result can be easily understood following the discussion of (13). In the resonant state, Es v(x) can be satisfied. Furthermore, if Es is large enough to satisfy
Fig. 4. log10 |Es2 /Es1 | as a function of Es − v(x) and the fitting function are plotted, where v(x) ≈ −50.
√ √ ω = |v(x) − Es | ≈ Es , then Shrödinger equation describes a plane-wave and (1) becomes a test-like equation as (3). In this case the P-stable Numerov method is an exact iterative formula. So it is not surprising to see from Table 2 that Es2 /Es1 can approach 10−8 for h = 1/32 for the last eigenvalue. (2) For the bound state, Es2 /Es1 for a given Es is essentially independent of step size. In Fig. 4, we plot log10 |Es 2 /Es 1 | as a function of Es − v(x). From the figure, we can see Es 2 /Es 1 ≈ −6 2 0.026 × 10−0.011Es +6.62×10 Es . It tells us that it is not necessary to use a method with a high algebraic order to obtain a high accurate oscillatory solution.
6
Z. Wang / Computer Physics Communications 167 (2005) 1–6
All the computations are carried out by the famous algorithm system Mathematica4.2 on an IBM PC-AT with AMD Athlon XP 2600+, 512M memory.
the progress of this work. This work is supported by National Natural Science Foundation of China (Grant no. 60371033).
5. Conclusion
References
In this paper we modify the traditional Numerov method by using trigonometric fitting to extend the interval of period to infinity. We also take advantage of the parameter to improve the accuracy. Although our new method has the same structure and algebra order as the traditional Numerov method, it is more accurate than the original. For the resonant state, it is very powerful and it can surpass the traditional method by 6 orders. The idea behind this method is easy to generalize to the other method. We think that the new method can replace the Numerov method for solving Schrödinger equation.
Acknowledgements Authors would like to thank Department of Physics, Shanghai University, for its continuous support during
[1] J.M. Blatt, J. Comput. Phys. 1 (1967) 382. [2] A.C. Allison, J. Comput. Phys. 6 (1970) 378. [3] G. Herzberg, Spectra of Diatomic Molecules, Van Nostrand, Toronto, 1950. [4] W. Cooley, Math. Comput. 15 (1961) 363. [5] J.R. Cash, A.D. Raptis, Comput. Phys. Comm. 33 (1984) 299. [6] J.D. Lambert, I.A. Watson, J. Inst. Math. Appl. 18 (1976) 189. [7] A.D. Raptis, J.R. Cash, Comput. Phys. Comm. 36 (1985) 113. [8] P.S. Williams, T.E. Simos, Int. J. Mod. Phys. C 11 (2000) 785. [9] L.Gr. Ixaru, M. Rizea, Comput. Phys. Comm. 19 (1980) 23. [10] A.D. Raptis, A.C. Allison, Comput. Phys. Comm. 14 (1978) 1. [11] A.D. Raptis, Computing 28 (1982) 373. [12] A.D. Raptis, J.R. Cash, Comput. Phys. Comm. 44 (1987) 95. [13] T.E. Simos, J. Comput. Phys. 148 (1999) 305. [14] G. Avdelas, T.E. Simos, Comput. Math. Appl. 31 (1996) 85. [15] A.C. Allison, A.D. Raptis, T.E. Simos, J. Comput. Phys. 97 (1991) 240. [16] Z. Wang, Y. Dai, N. Math. J. Chin. Univ. (Suppl.) 12 (2003) 146. [17] Z. Wang, Y. Dai, Int. J. Mod. Phys. C 14 (2003) 1087. [18] Z. Wang, Y. Ge, Y. Dai, D. Zhao, Comput. Phys. Comm. 160 (2004) 23.
A new effective algorithm for the resonant state of a Schrödinger equation Zhongcheng Wang Department of Physics, Shanghai University, 99 ShangDa Road, Shanghai 200436, P.R. China Received 9 September 2004; received in revised form 7 December 2004; accepted 13 December 2004
Abstract In this paper we present a new effective algorithm for the Schrödinger equation. This new method differs from the original Numerov method only in one simple coefficient, by which we can extend the interval of periodicity from 6 to infinity and obtain an embedded correct factor to improve the accuracy. We compare the new method with the original Numerov method by the well-known problem of Woods–Saxon potential. The numerical results show that the new method has great advantage in accuracy over the original. Particularly for the resonant state, the accuracy is improved with four orders overall, and even six to seven orders for the highest oscillatory solution. Surely, this method will replace the original Numerov method and be widely used in various area. 2004 Elsevier B.V. All rights reserved. Keywords: Numerov method; Trigonometric fitting; P-stable; High-oscillatory solution; Schrödinger equation PACS: 02.60.Cb; 02.60.Lj; 02.70.Bf
1. Introduction Many physical problems occurred in various scientific areas such as theoretical physics, nuclear physics, physical chemistry, quantum chemistry, and molecular physics [1–4] can be reduced into an one-dimensional Schrödinger equation, which can be written as y (x) = f (x)y(x),
(1)
E-mail address: [email protected] (Z. Wang). 0010-4655/$ – see front matter 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.cpc.2004.12.004
with f (x) = v(x) − Es , where the function v(x) is given and denotes the effective potential and Es is a real number or complex denoting the energy. The earliest method to solve (1) numerically is the Numerov’s method. Although this method is only of algebraic order four, it is still widely used in many area, because of its simplicity and wider interval of periodicity. On the other hand, improving this method and seeking a more accurate, fast and reliable algorithm to achieve the numerical solution have been the subject of great activity for many years [5–18]. From the numerical solution to the well-known Woods–Saxon potential problem by
2
Z. Wang / Computer Physics Communications 167 (2005) 1–6
the available methods [9–18], it is interesting to find that the error in the calculated eigenvalue always increases with its value. Hence, to obtain the highest resonant state with high accuracy is a headache problem and has long been the motivation for developing a new method, such as the exponential fitting [9–12], the embedded methods [14], and the eighth-order formula [15], Obrechkoff method [16–18]. In order to obtain a high accurate solution, we must adapt a method with a high algebraic order which usually involves a complex algorithm. Is it possible to find a simple and effective way to obtain the high accurate solution? We find it is possible. The motivation of this paper is to report a simple and effective algorithm to improve the accuracy of the resonant state. This algorithm is nearly the same as the original Numerov method. We just use trigonometric fitting to extend the interval of periodicity of the traditional Numerov method to infinity. Then we use an embedded parameter to correct the error. This paper is arranged as follows: in the next section, we will give the derivation for the new method. In Section 3 we will analyze the stability and the local truncation error. In Section 4 we will test our new method by the Woods–Saxon potential. In the final section we will give a conclusion.
h2 (1 − 2c1 − c2 ) ω2 y(x) + y (x) h4 (1 − 12c1 ) ω4 y(x) − y (4) (x) + O(h6 ) = 0. 12 (5) From (5), we obtain the coefficients as
−
1 5 , c2 = . (6) 12 6 After substituting (4) and (6) into (2) and changing the notations such as y(x) → yk and y (x) → fk yk , we obtain the following iterative formula 5h2 αs + yk+2 = fk+1 yk+1 6 −1 2 h2 h − 1 − fk yk 1 − fk+2 , (7) 12 12 c1 =
where
5 2 2 1 2 2 αs = h ω + 2 1 + h ω cos(hω). 6 12
(8)
The local truncation error is LET(h) = −
h6 6 ω y(x) + y (6) (x) . 240
(9)
3. Stability and local truncation error 2. Derivation We consider the following structure which are the same for the Numerov method y(x + h) + y(x − h) − αs y(x) = h2 c1 y (x + h) + y (x − h) + c2 y (x) .
(2)
According to the stability theory by Lambert and Watson [6], the stability of (2) can be examined through the following test equation
y (x) = −ω y(x). 2
(3)
We apply (3) and its solutions to (2) and obtain the following equation αs = h2 ω2 c2 + 2 1 + h2 ω2 c1 cos(hω). (4) We use (4) to eliminate αs in (2) and convert the resultant expression to a Taylor series such as
It is easy to see that the traditional Numerov method is just the special case of (7) with αs = 2. We can show the interval of periodicity of (7) to be infinite by making the following substitutions: fk , fk+1 , fk+2 → −ω2 , and yk+2 → λ2 , yk+1 → λ and yk → 1, and hω → H to (7) and (8). The characteristic equation can be obtained as H2 2 λ − cos(H )λ + 1 = 0. 1+ (10) 12 For comparison, we let αs = 2 for (8) and obtain the characteristic equation 5H 2 H2 2 H2 − 2− λ+ 1+ λ 1+ (11) 12 6 12 for the traditional Numerov method. The characteristic roots are λ± = e±iH for (10), which is independent of H . So it is a P-stable and the interval of periodicity is infinite.
Z. Wang / Computer Physics Communications 167 (2005) 1–6
Although for the traditional Numerov method, the characteristic roots of (11) can √ be written as λ± = e±iθ(H ) , but it is only for H 6. So the corresponding interval of periodicity is 6. From (9), we can see that the local truncation error is ω dependent, so that it should be examined with more terms such as 1 6 6 h ω y(x) + y (6) (x) LTE(h) = − 240 11 8 8 h ω y(x) − y (8) (x) + 60480 13 − h10 ω10 y(x) + y (10) (x) 3628800 1 h12 ω12 y(x) − y (12) (x) + 23950080 17 h14 ω14 y(x) + y (14) (x) − 52306974720 + ···. (12) After rearranging, the above expression can be written as LTE(h) =
∞ (−1)n (n − 2)(2n + 3)h2n
3(2n)! 2n × ω y(x) − (−1)n y (2n) (x) . n=3
(13)
Because ω2n y(x) − (−1)n y (2n) (x) ≡ 0 for (3), the total local truncation error (3) is zero for the test equation. Therefore, the trigonometrically-fitted P-stable Numerov method not only extends the interval of periodicity to infinity but also provides an exact formula for the test-like problem. Can we use this property to improve the numerical accuracy for Schrödinger equation with periodic solution? We find the answer is “yes”. If we explain ω as the wave vector, which appears as a parameter in the characteristic root y(x) = e√iωx of the test equation (3), then we can let ω = |v(x) − Es |. Such that ω is position sensitive. If this assumption is correct, we can correct the error for every step by ω. We need to test our idea by the numerical test.
4. Numerical illustrations We will use the Woods–Saxon potential [15] v(x) = c0 z 1 − a(1 − z) (14)
3
Fig. 1. The Woods–Saxon potential defined by (14) is plotted in the range of x ∈ (0, 15).
with z = {exp[a(x − b)] + 1}−1 , c0 = −50, a = 5.0/3 and b = 7.0 to test our new method and compare it with the traditional Numerov method. The plot of the Woods–Saxon potential is given in Fig. 1. The test domain is 0 x 15, which is the same as in the literature [18]. The first problem is to find the eigenvalues in the bound state, for which −50 < Es < 0 and the √ boundary conditions are: y(0) = 0 and y(15) = e− Es Nk h/2 ≈ 0. The Schrödinger equation (1) is integrated forwards and backwards from the two boundary nodes, and an approximate eigenvalue can be found by the fact that the wave function and its derivative from the two sides at the matching point M are equal. − The connecting condition can be written as yM y˜M y˜M yM = 0, where yM and yM are obtained from are from the backward. the forward and y˜M and y˜M The matching√point is chosen to be x = 7.5. y0 = 0 and y1 = sin[ (50 + Es )h] are the initial wave functions for the forward integral; yNk = 0 and yNk −1 = √ − −E s (Nk −1)h/2 are for the backward, where N is e k the total number of steps, h is the step size. Like the other authors [15], the second problem is to find the real eigenvalues in the range 5 Es 1000 for which the phase shift δ(Es ) is just to π/2 in the resonant state with the same domain. The √ boundary conditions are: y(0) = 0 and y(15) = cos[ Es x]. The matching point√is also chosen to be x = 7.5. y0 = 0 and y1 = sin[ ES h] are the initial wave √ functions = cos( Es Nk h) and for the forward integral; y N k √ yNk −1 = cos( ES (Nk − 1)h) for the backwards. − We use the new method to calculate (yM y˜M m y˜M yM )/Es as a function of energy Es for the bound state with h = 1/32 and m = 6, and for the real res-
4
Z. Wang / Computer Physics Communications 167 (2005) 1–6
(a)
(b)
− y˜ y )/E 6 to energy E ∈ (−50, −5) with h = 1/32 and m = 6; (b) is the plot of (y y˜ − y˜ y ) to Fig. 2. (a) is the plot of (yM y˜M s M M M M M M s energy Es ∈ (5, 1005) with h = 1/64 and δ(Es ) = π/2.
Table 1 Error comparison between the new P-stable Numerov and the traditional Numerov method for the seven eigenvalues in the bound states [15] Ei
h
Es1
t (s)
Es2
t (s)
Es2 /Es1
−49.457788728
1/16 1/32 1/64 1/16 1/32 1/64 1/16 1/32 1/64 1/16 1/32 1/64 1/16 1/32 1/64 1/16 1/32 1/64 1/16 1/32 1/64
15 1 0 2.4 × 103 152 10 3.1 × 104 1.9 × 103 120 1.6 × 105 1.0 × 104 634 5.5 × 105 3.4 × 104 2.1 × 103 1.4 × 106 8.8 × 104 5.5 × 103 2.9 × 106 1.8 × 105 1.1 × 104
0.28 0.55 1.11 0.19 0.39 0.73 0.17 0.33 0.66 0.17 0.33 0.66 0.16 0.33 0.64 0.14 0.33 0.66 0.14 0.28 0.56
2 0 0 163 10 1 1.0 × 103 64 4 3.1 × 103 192 13 6.5 × 103 403 26 1.1 × 104 663 41 1.4 × 104 879 55
0.22 0.45 0.92 0.11 0.30 0.61 0.13 0.27 0.53 0.13 0.27 0.53 0.14 0.27 0.53 0.13 0.27 0.53 0.13 0.22 0.45
1 × 10−1 6 × 10−2 5 × 10−1 7 × 10−2 7 × 10−2 1 × 10−1 3 × 10−2 3 × 10−2 3 × 10−2 2 × 10−2 2 × 10−2 2 × 10−2 1 × 10−2 1 × 10−2 1 × 10−2 8 × 10−3 8 × 10−3 8 × 10−3 5 × 10−3 5 × 10−3 5 × 10−3
−46.290753954
−41.232607772
−34.672313205
−26.873448915
−18.094688282
−8.676081670
The true eigenvalues with nine decimal places are shown in the first column, the step-size with 1/16, 1/32 and 1/64 are in the second. The absolute deviations of the computed eigenvalues from the exact values, Es1 for Numerov’s method and Es2 for P-stable Numerov’s method, in 109 units and the consumed CPU time (in seconds) are shown in each sub-column. In the last column, the ratio of the absolute deviations of Es2 /Es1 .
onant states with h = 1/64, m = 0 and δ(Es ) = π/2, which are shown in Fig. 2 (a) and (b), respectively. From the figures, we find that there exists total nineteen eigenstates, thirteen for the bound states and six for resonant states, which are in agreement with those in [18]. In Table 1, we compare the absolute error, |Es | = Es − Ei for seven bound states with the step size of 1/16, 1/32, and 1/64, and the consuming CPU time.
In Table 2, we compare the four resonant states with the size of 1/32, 1/64, and 1/128. Where the true solution Ei to the Woods–Saxon potential problem were obtained by using analytic method which is correct to nine decimal places for the bound states and six decimal places for the resonant states [15]. In order to give an impression how great the improvement has been made to the resonant state, in Fig. 3 we give the plot of the log10 |Es × 106 |
Z. Wang / Computer Physics Communications 167 (2005) 1–6
5
Table 2 Error comparison between the new P-stable Numerov and the traditional Numerov method for the four eigenvalues in the resonant states [15] Ei 53.588872
163.215341
341.495874
989.701916
h 1/32 1/64 1/128 1/32 1/64 1/128 1/32 1/64 1/128 1/32 1/64 1/128
Es1 1.4 × 104 878 55 4.8 × 105 2.9 × 104 1.8 × 103 7.5 × 106 4.4 × 105 2.7 × 104 6.4 × 108 2.9 × 107 1.6 × 106
t (s)
Es2
t (s)
Es2 /Es1
0.23 0.47 0.95 0.19 0.48 0.97 0.19 0.38 0.75 0.19 0.38 0.75
3 1 0 8 1 0 15 2 1 49 6 1
0.14 0.38 0.75 0.16 0.31 0.55 0.27 0.31 0.63 0.41 0.53 0.94
2 × 10−4 6 × 10−4 3 × 10−3 2 × 10−5 5 × 10−5 3 × 10−4 2 × 10−6 5 × 10−6 4 × 10−5 8 × 10−8 2 × 10−7 1 × 10−6
The true eigenvalues, Ei , with six decimal places are shown in the first column, the step-size with 1/32, 1/64 and 1/128 are in the second. The absolute deviations of the computed eigenvalues from the exact values in 106 units and the consumed CPU time (in seconds) are shown in each sub-column. As the same as in Table 1, Es 1 is for Numerov’s method and Es 2 for P-stable Numerov’s method. In the last column, the ratio of the absolute deviations of |Es 2 /Es 1 |.
Fig. 3. Es = Es − Ei as a function of h is plotted for the new method, where Ei = 989.701916. In the figure, 10−6 and 10−7 are denoted for the plots of Es × 10−6 and Es × 10−7 obtained by the traditional Numerov method.
as a function of h for the last eigenvalue, Ei = 989.701916, obtained by the new method, and the plots of log10 |Es | and log10 |Es × 10−1 | by the Numerov method. From the tables and figures, it is easy to see the accuracy has been greatly improved. It is interesting to note (1) Es2 /Es1 decreases with the value of the eigenvalue. It seems that the larger Es is the higher the accuracy is. This result can be easily understood following the discussion of (13). In the resonant state, Es v(x) can be satisfied. Furthermore, if Es is large enough to satisfy
Fig. 4. log10 |Es2 /Es1 | as a function of Es − v(x) and the fitting function are plotted, where v(x) ≈ −50.
√ √ ω = |v(x) − Es | ≈ Es , then Shrödinger equation describes a plane-wave and (1) becomes a test-like equation as (3). In this case the P-stable Numerov method is an exact iterative formula. So it is not surprising to see from Table 2 that Es2 /Es1 can approach 10−8 for h = 1/32 for the last eigenvalue. (2) For the bound state, Es2 /Es1 for a given Es is essentially independent of step size. In Fig. 4, we plot log10 |Es 2 /Es 1 | as a function of Es − v(x). From the figure, we can see Es 2 /Es 1 ≈ −6 2 0.026 × 10−0.011Es +6.62×10 Es . It tells us that it is not necessary to use a method with a high algebraic order to obtain a high accurate oscillatory solution.
6
Z. Wang / Computer Physics Communications 167 (2005) 1–6
All the computations are carried out by the famous algorithm system Mathematica4.2 on an IBM PC-AT with AMD Athlon XP 2600+, 512M memory.
the progress of this work. This work is supported by National Natural Science Foundation of China (Grant no. 60371033).
5. Conclusion
References
In this paper we modify the traditional Numerov method by using trigonometric fitting to extend the interval of period to infinity. We also take advantage of the parameter to improve the accuracy. Although our new method has the same structure and algebra order as the traditional Numerov method, it is more accurate than the original. For the resonant state, it is very powerful and it can surpass the traditional method by 6 orders. The idea behind this method is easy to generalize to the other method. We think that the new method can replace the Numerov method for solving Schrödinger equation.
Acknowledgements Authors would like to thank Department of Physics, Shanghai University, for its continuous support during
[1] J.M. Blatt, J. Comput. Phys. 1 (1967) 382. [2] A.C. Allison, J. Comput. Phys. 6 (1970) 378. [3] G. Herzberg, Spectra of Diatomic Molecules, Van Nostrand, Toronto, 1950. [4] W. Cooley, Math. Comput. 15 (1961) 363. [5] J.R. Cash, A.D. Raptis, Comput. Phys. Comm. 33 (1984) 299. [6] J.D. Lambert, I.A. Watson, J. Inst. Math. Appl. 18 (1976) 189. [7] A.D. Raptis, J.R. Cash, Comput. Phys. Comm. 36 (1985) 113. [8] P.S. Williams, T.E. Simos, Int. J. Mod. Phys. C 11 (2000) 785. [9] L.Gr. Ixaru, M. Rizea, Comput. Phys. Comm. 19 (1980) 23. [10] A.D. Raptis, A.C. Allison, Comput. Phys. Comm. 14 (1978) 1. [11] A.D. Raptis, Computing 28 (1982) 373. [12] A.D. Raptis, J.R. Cash, Comput. Phys. Comm. 44 (1987) 95. [13] T.E. Simos, J. Comput. Phys. 148 (1999) 305. [14] G. Avdelas, T.E. Simos, Comput. Math. Appl. 31 (1996) 85. [15] A.C. Allison, A.D. Raptis, T.E. Simos, J. Comput. Phys. 97 (1991) 240. [16] Z. Wang, Y. Dai, N. Math. J. Chin. Univ. (Suppl.) 12 (2003) 146. [17] Z. Wang, Y. Dai, Int. J. Mod. Phys. C 14 (2003) 1087. [18] Z. Wang, Y. Ge, Y. Dai, D. Zhao, Comput. Phys. Comm. 160 (2004) 23.