Abstract
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.
1. Introduction
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. [17] 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. [18] 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 [19] and Lyche [20].
Recently, Simos et al. [21] introduced an exponentially fitted explicit Runge-Kutta method which integrates exactly the model problem. Van de Vyver [22] constructed a new fourth-order explicit Runge-Kutta method based on Simos’ approach [21] 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. [23] which exactly integrates differential initial value problems whose solutions are linear combinations of functions of the form. More recently, Chan and Tsai [24] 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 [21] and of Chan and Tsai [24], 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 [24]: 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 [24], 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 [24] 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:
It can be seen that when, this new method reduces to the fifth-order TDRK method (9) of Chan and Tsai [24]. The local truncation error of the above method is given by
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 [25]). 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 [23].(ii)RK4V: the fourth-order optimized Runge-Kutta method given by Van de Vyver in [22].(iii)RK5S: the fifth-order trigonometrically fitted Runge-Kutta method derived by Anastassi and Simos in [26].(iv)RK4S: the fourth-order exponentially fitted Runge-Kutta method derived by Simos in [21].(v)EFRK5: the fifth-order trigonometrically fitted Runge-Kutta method derived by Sakas and Simos in [27].(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 [10]): 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 [9] 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 [28]:
or equivalently
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 [29] 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 [30]. 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
where
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.
5. Conclusions
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.
Acknowledgments
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).