Research Article | Open Access
A New Trigonometrically Fitted Two-Derivative Runge-Kutta Method for the Numerical Solution of the Schrödinger Equation and Related Problems
A new trigonometrically fitted fifth-order two-derivative Runge-Kutta method with variable nodes is developed for the numerical solution of the radial Schrödinger equation and related oscillatory problems. Linear stability and phase properties of the new method are examined. Numerical results are reported to show the robustness and competence of the new method compared with some highly efficient methods in the recent literature.
We are interested in initial value problems related to systems of first-order ordinary differential equations (ODEs) of the form whose solutions have a pronounced oscillatory character. Such problems often occur in many applied fields such as celestial mechanics, physical chemistry, quantum mechanics, and electronics. The investigation of the numerical solution of (1) has been the subject of extensive research activity during the last few decades [1–16]. Such special optimized methods fall into two classes. The first class consists of the methods with constant coefficients. These methods can be applied to any problem with periodic solution. Methods of the second class have coefficients depending on the frequency of problem. When a good estimate of the frequency is known in advance one can use procedures such as exponential/trigonometric fitting or phase fitting. For example, Chen et al.  considered phase-fitted and amplification-fitted RK methods whose updates are also phase fitted and amplification fitted. Their methods are shown to be more efficient than the codes in the literature for some typical test problems and for the Lotka-Volterra system and a two-gene regulatory network in biology as well. You et al.  investigated the trigonometrically fitted Scheifele methods for oscillatory problems. A good theoretical foundation of the exponential fitting techniques for multistep methods was presented by Gautschi  and Lyche .
Recently, Simos et al.  introduced an exponentially fitted explicit Runge-Kutta method which integrates exactly the model problem. Van de Vyver  constructed a new fourth-order explicit Runge-Kutta method based on Simos’ approach  for the numerical integration of the Schrödinger equation. And another exponentially fitted Runge-Kutta method with four stages was constructed by Vanden Berghe et al.  which exactly integrates differential initial value problems whose solutions are linear combinations of functions of the form. More recently, Chan and Tsai  presented a family of two-derivative Runge-Kutta (TDRK) methods. An advantage of the two-derivative Runge-Kutta methods over classical Runge-Kutta methods is that they can reach higher order with fewer function evaluations.
Inspired by the approaches of Simos  and of Chan and Tsai , we will construct in this paper a new two-derivative Runge-Kutta method for the numerical integration of the Schrödinger equation and related oscillatory problems. The remainder of this paper is organized as follows: in Section 2, we introduce the scheme of two-derivative Runge-Kutta methods and present (up to fifth) order conditions for TDRK methods. In Section 3, we construct a new trigonometrically fitted TDRK method and analyze its linear stability and phase properties. In Section 4, some numerical examples are given to show the effectiveness and competence of our new method compared to the selected methods in the recent literature. Section 5 is devoted to conclusive remarks.
2. Two-Derivative Runge-Kutta Methods
We start with a special form of explicit two-derivative Runge-Kutta (TDRK) methods proposed by Chan and Tsai : where. The coefficients of the scheme can be expressed by the Butcher tableau or simplify by. This explicit TDRK method involves only one evaluation of the functionandevaluations of the functionper step.
According to Chan and Tsai , the conditions for order up to five are listed as follows.(i)Order 2 requires (ii)Order 3 requires, in addition, (iii)Order 4 requires, in addition, (iv)Order 5 requires, in addition, In practice, the following simplifying assumption is useful: Choosingand solving (4)–(8) yield a three-stage TDRK method of Chan and Tsai  which is given by the following tableau:
3. A Fifth-Order TDRK Method with Frequency-Dependent Coefficients
In order to adapt TDRK methods to the problem (1) whose solutions share an oscillatory feature withbeing an accurate estimate of the principal frequency, we allow the coefficients to depend on. In particular, we consider in this paper the three-stage explicit TDRK method given by the following Butcher tableau:
3.1. Construction of the New Method
If we require the TDRK method (10) to integrate exactly; that is, then we have where Solving (4)–(6) for, the simplifying conditions (8) and (12), we obtain the coefficients expressed in terms of in which Assuming thatis a constant and letting we obtain. It is easy to check that The limit of the above equation asverifies the second equation of (7). Thus we obtain a new TDRK method of order five given by (14) with. We denote this method by TFTDRKV5.
For small values ofthe above formula is subject to heavy cancelations and in that case the following series expansions should be used:
3.2. Stability and Phase Analysis of the New Method
In this section, we investigate the linear stability of the new method. We consider the following test equation: When applied to (20) the TDRK method (2) produces the difference equation where withthe identity matrix and.
Definition 1. For the TDRK method (2) with stability function, the region in the-plane is called the imaginary stability region of the method. And any closed curve defined byis a stability boundary of the method.
The imaginary stability region of the TFTDRKV5 is plotted in Figure 1.
Definition 2 (see ). For the TDRK method (2) with stability function, the quantities are called the phase-lag (dispersion) and amplification factor error (dissipation), respectively. If the method is said to have phase-lag order and dissipation order , respectively.
Denoting the ratio, we obtain the following expressions for the phase-lag and the amplification of the method TFTDRKV5: Thus, the method TFTDRKV5 has a phase-lag of order six and a dissipation of order seven.
4. Numerical Results
In this section we carry out six numerical results to illustrate the performance of the new method. The criterion used in the numerical comparison is the decimal logarithm of the maximum global error (LOG10(ERROR)) versus the computational effort measured by the CPU seconds (CPU SECONDS) required by each problem. The methods used for comparison are listed in the following.(i)EFRK4: the fourth-order exponentially fitted Runge-Kutta method given by Vanden Berghe et al. in .(ii)RK4V: the fourth-order optimized Runge-Kutta method given by Van de Vyver in .(iii)RK5S: the fifth-order trigonometrically fitted Runge-Kutta method derived by Anastassi and Simos in .(iv)RK4S: the fourth-order exponentially fitted Runge-Kutta method derived by Simos in .(v)EFRK5: the fifth-order trigonometrically fitted Runge-Kutta method derived by Sakas and Simos in .(vi)TFTDRKV5: the fifth-order TFTDRK method with one evaluation of functionand three evaluations of functionper step derived in Section 3.1 of this paper.
Problem 1. Consider the numerical integration of the Schrödinger equation with the well-known Woods-Saxon potential where, andThe domain of numerical integration is. It is appropriate to chooseas follows (see ): In this experiment we consider the resonance problem (). The numerical results are compared with the analytical solution of the Woods-Saxon potential, rounded to six decimal places. In Figures 2, 3, 4, and 5, we plot the error ofversus the computational effort measured by CPU seconds required by each method for, respectively.
Problem 2. We consider the inhomogeneous equation whose exact solution is In this test we choose. The numerical results given in Figure 6 are computed with the stepsizes, for all the methods on the interval.
Problem 3. We consider the linear problem studied in  whose exact solution is In this test we take. The numerical results given in Figure 7 are computed with the stepsizes, for all the methods on the interval.
Problem 4. We consider the following “almost periodic” orbit problem studied by Franco and Palacios :
whereandThe analytic solution to the problem is given by
In this test we take. The numerical results given in Figure 8 are computed with the stepsizes, for all the methods on the interval.
Problem 5. We consider the two-dimensional problem studied in  where The analytical solution is In this test we choose. The numerical results given in Figure 9 are computed with the stepsizes, for all the methods on the interval.
Problem 6. We consider the Fermi-Pasta-Ulam problem studied in . This problem consists of a chain of springs, where nonlinear springs alternate with stiff harmonic springs. The problem can be described by a Hamiltonian system with Hamiltonian function
where. This is equivalent to the nonlinear oscillatory problem of the form
Since the analytic solution is not available, we study the numerical total energyon the interval. In this test we choosewithandand we take the initial values
the remaining initial values being zero.
The numerical results given in Figure 10 are computed with the stepsizes , for all the methods.
In Figures 2–10, we can see that the new method TFTDTRK5V is the most efficient among the six methods we select for comparison.
A new trigonometrically fitted two-derivative Runge-Kutta method of algebraic order five is derived in this paper. We also analyze the linear stability and phase properties of the new method. Numerical results are reported to show the robustness and competence of the new method compared with some highly efficient methods in the recent literature. To explain the superiority of the new method, we observe that, compared with RK methods of a specific number of stages, TDRK methods have the possibility of gaining a higher algebraic order than traditional RK methods with the same number of stages. We conclude that, for problems with oscillatory solutions, trigonometrically fitted TDRK methods are more accurate and more efficient than the adapted (exponentially or trigonometrically fitted) RK methods.
The authors are grateful to the editors and the anonymous referees for their constructive comments and valuable suggestions which have helped to improve the paper. This research was partially supported by NSFC (no. 11101357 and no. 11171155), the foundation of Shandong Outstanding Young Scientists Award Project (no. BS2010SF031), the foundation of Scientific Research Project of Shandong Universities (no. J13LI03), NSF of Shandong Province, China (no. ZR2011AL006), the Fundamental Research Fund for the Central Universities (no. Y0201100265), and the foundation of Scientific Research Project of Weifang (no. 20121103).
- L. G. Ixaru and M. Rizea, “A numerov-like scheme for the numerical solution of the Schrödinger equation in the deep continuum spectrum of energies,” Computer Physics Communications, vol. 19, no. 1, pp. 23–27, 1980.
- T. E. Simos, “A family of fifth algebraic order trigonometrically fitted Runge-Kutta methods for the numerical solution of the Schrödinger equation,” Computational Materials Science, vol. 34, no. 4, pp. 342–354, 2005.
- H. van de Vyver, “An embedded phase-fitted modified Runge-Kutta method for the numerical integration of the radial Schrödinger equation,” Physics Letters A, vol. 352, no. 4-5, pp. 278–285, 2006.
- T. E. Simos, “An embedded Runge-Kutta method with phase-lag of order infinity for the numerical solution of the Schrödinger equation,” International Journal of Modern Physics C, vol. 11, no. 6, pp. 1115–1133, 2000.
- T. E. Simos and J. V. Aguiar, “A modified phase-fitted Runge-Kutta method for the numerical solution of the Schrödinger equation,” Journal of Mathematical Chemistry, vol. 30, no. 1, pp. 121–131, 2001.
- Z. Kalogiratou, Th. Monovasilis, and T. E. Simos, “Computation of the eigenvalues of the Schrödinger equation by exponentially-fitted Runge-Kutta-Nyström methods,” Computer Physics Communications, vol. 180, no. 2, pp. 167–176, 2009.
- A. A. Kosti, Z. A. Anastassi, and T. E. Simos, “Construction of an optimized explicit Runge-Kutta-Nyström method for the numerical solution of oscillatory initial value problems,” Computers & Mathematics with Applications, vol. 61, no. 11, pp. 3381–3390, 2011.
- A. A. Kosti, Z. A. Anastassi, and T. E. Simos, “An optimized explicit Runge-Kutta-Nyström method for the numerical solution of orbital and related periodical initial value problems,” Computer Physics Communications, vol. 183, no. 3, pp. 470–479, 2012.
- J. M. Franco, “Runge-Kutta methods adapted to the numerical integration of oscillatory problems,” Applied Numerical Mathematics, vol. 50, no. 3-4, pp. 427–443, 2004.
- H. Van de Vyver, “An embedded exponentially fitted Runge-Kutta-Nyström method for the numerical solution of orbital problems,” New Astronomy, vol. 11, no. 8, pp. 577–587, 2006.
- E. Hairer, S. P. Nørsett, and G. Wanner, Solving Ordinary Differential Equations. I. Nonstiff Problems, vol. 8 of Springer Series in Computational Mathematics, Springer, Berlin, Germany, 2nd edition, 1993.
- J. M. Farto, A. B. González, and P. Martín, “An algorithm for the systematic construction of solutions to perturbed problems,” Computer Physics Communications, vol. 111, no. 1–3, pp. 110–132, 1998.
- H. Ramos and J. Vigo-Aguiar, “On the frequency choice in trigonometrically fitted methods,” Applied Mathematics Letters, vol. 23, no. 11, pp. 1378–1381, 2010.
- G. Vanden Berghe, L. Gr. Ixaru, and H. De Meyer, “Frequency determination and step-length control for exponentially-fitted Runge-Kutta methods,” Journal of Computational and Applied Mathematics, vol. 132, no. 1, pp. 95–105, 2001.
- J. Vigo-Aguiar and H. Ramos, “Dissipative Chebyshev exponential-fitted methods for numerical solution of second-order differential equations,” Journal of Computational and Applied Mathematics, vol. 158, no. 1, pp. 187–211, 2003.
- G. Vanden Berghe and M. Van Daele, “Exponentially-fitted Obrechkoff methods for second-order differential equations,” Applied Numerical Mathematics, vol. 59, no. 3-4, pp. 815–829, 2009.
- Z. Chen, X. You, X. Shu, and M. Zhang, “A new family of phase-fitted and amplification-fitted Runge-Kutta type methods for oscillators,” Journal of Applied Mathematics, vol. 2012, Article ID 236281, 27 pages, 2012.
- X. You, Y. Zhang, and J. Zhao, “Trigonometrically-fitted Scheifele two-step methods for perturbed oscillators,” Computer Physics Communications, vol. 182, no. 7, pp. 1481–1490, 2011.
- W. Gautschi, “Numerical integration of ordinary differential equations based on trigonometric polynomials,” Numerische Mathematik, vol. 3, pp. 381–397, 1961.
- T. Lyche, “Chebyshevian multistep methods for ordinary differential equations,” Numerische Mathematik, vol. 19, pp. 65–75, 1972.
- T. E. Simos, “An exponentially-fitted Runge-Kutta method for the numerical integration of initial-value problems with periodic or oscillating solutions,” Computer Physics Communications, vol. 115, no. 1, pp. 1–8, 1998.
- H. Van de Vyver, “Stability and phase-lag analysis of explicit Runge-Kutta methods with variable coefficients for oscillatory problems,” Computer Physics Communications, vol. 173, no. 3, pp. 115–130, 2005.
- G. Vanden Berghe, H. De Meyer, M. Van Daele, and T. Van Hecke, “Exponentially fitted Runge-Kutta methods,” Journal of Computational and Applied Mathematics, vol. 125, no. 1-2, pp. 107–115, 2000.
- R. P. K. Chan and A. Y. J. Tsai, “On explicit two-derivative Runge-Kutta methods,” Numerical Algorithms, vol. 53, no. 2-3, pp. 171–194, 2010.
- P. J. van der Houwen and B. P. Sommeijer, “Explicit Runge-Kutta (-Nyström) methods with reduced phase errors for computing oscillating solutions,” SIAM Journal on Numerical Analysis, vol. 24, no. 3, pp. 595–617, 1987.
- Z. A. Anastassi and T. E. Simos, “Trigonometrically fitted Runge-Kutta methods for the numerical solution of the Schrödinger equation,” Journal of Mathematical Chemistry, vol. 37, no. 3, pp. 281–293, 2005.
- D. P. Sakas and T. E. Simos, “A fifth algebraic order trigonometrically-fitted modified Runge-Kutta Zonneveld method for the numerical solution of orbital problems,” Mathematical and Computer Modelling, vol. 42, no. 7-8, pp. 903–920, 2005.
- J. M. Franco and M. Palacios, “High-order P-stable multistep methods,” Journal of Computational and Applied Mathematics, vol. 30, no. 1, pp. 1–10, 1990.
- J. M. Franco, “A class of explicit two-step hybrid methods for second-order IVPs,” Journal of Computational and Applied Mathematics, vol. 187, no. 1, pp. 41–57, 2006.
- L. Galgani, A. Giorgilli, A. Martinoli, and S. Vanzini, “On the problem of energy equipartition for large systems of the Fermi-Pasta-Ulam type: analytical and numerical estimates,” Physica D, vol. 59, no. 4, pp. 334–348, 1992.
Copyright © 2013 Yanwei Zhang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.