Journal of Applied Mathematics

Volume 2013 (2013), Article ID 937858, 9 pages

http://dx.doi.org/10.1155/2013/937858

Research Article

A New Trigonometrically Fitted Two-Derivative Runge-Kutta Method for the Numerical Solution of the Schrödinger Equation and Related Problems

Yanwei Zhang,¹ Haitao Che,² Yonglei Fang,¹ and Xiong You³

¹Department of Mathematics and Information Science, Zaozhuang University, Zaozhuang 277160, China

²School of Mathematics and Information Science, Weifang University, Weifang, Shandong 261061, China

³Department of Applied Mathematics, Nanjing Agricultural University, Nanjing 210095, China

Received 4 July 2013; Revised 28 September 2013; Accepted 4 October 2013

Academic Editor: Jesus Vigo-Aguiar

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.

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.

Figure 1: Stability region of method TFTDRKV5.

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.

Figure 2: Efficiency curves for Problem 1:.

Figure 3: Efficiency curves for Problem 1:.

Figure 4: Efficiency curves for Problem 1:.

Figure 5: Efficiency curves for Problem 1:.

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.

Figure 6: Efficiency curves for Problem 2.

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.

Figure 7: Efficiency curves for Problem 3.

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.

Figure 8: Efficiency curves for Problem 4.

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.

Figure 9: Efficiency curves for Problem 5.

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.

Figure 10: Efficiency curves for Problem 6.

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).

References

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.