About this Journal Submit a Manuscript Table of Contents
Discrete Dynamics in Nature and Society
Volume 2012 (2012), Article ID 324989, 20 pages
http://dx.doi.org/10.1155/2012/324989
Research Article

New 4(3) Pairs Diagonally Implicit Runge-Kutta-Nyström Method for Periodic IVPs

Department of Mathematics and Institute for Mathematical Research, Universiti Putra Malaysia, 43400 Serdang, Selangor, Malaysia

Received 14 February 2012; Revised 24 June 2012; Accepted 8 July 2012

Academic Editor: Taher S. Hassan

Copyright © 2012 Norazak Senu 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

New 4(3) pairs Diagonally Implicit Runge-Kutta-Nyström (DIRKN) methods with reduced phase-lag are developed for the integration of initial value problems for second-order ordinary differential equations possessing oscillating solutions. Two DIRKN pairs which are three- and four-stage with high order of dispersion embedded with the third-order formula for the estimation of the local truncation error. These new methods are more efficient when compared with current methods of similar type and with the L-stable Runge-Kutta pair derived by Butcher and Chen (2000) for the numerical integration of second-order differential equations with periodic solutions.

1. Introduction

In many scientific areas of engineering and applied sciences such as celestial mechanics, quantum mechanics, elastodynamics, theoretical physics and chemistry, and electronics, oscillatory problems of second-order ordinary differential equations (ODEs) can be found. An oscillatory problems of second-order ODEs have the following form:

An -stage Runge-Kutta-Nyström (RKN) method for the numerical integration of the IVP is given by where

The RKN parameters , and are assumed to be real and is the number of stages of the method. Introduce the -dimensional vectors , and and matrix , where , and , respectively. The latter contains the class of Diagonally Implicit RKN (DIRKN) methods for which all the entries in the diagonal of are equal. An embedded pair of DIRKN methods is based on the method of order and the other DIRKN method of order and can be expressed in Butcher notation by the table of coefficients (see Table 1).

tab1
Table 1: m-stage DIRKN pair.

Several authors in their papers have developed numerical methods for this class of problems, for example, van der Houwen and Sommeijer [1, 2], Sideridis and Simos [3], García et al. [4], and Senu et al. [5]. Next, Van de Vyver [6] and Senu et al. [7] obtained explicit RKN method with minimal phase lag. Franco [8] have developed explicit hybrid method for periodic IVPs. For implicit RKN methods, see for example, van der Houwen and Sommeijer [2], Sharp et al. [9]. Imoni et al. [10] and Al-Khasawneh et al. [11] have developed general purpose of DIRKN methods with variable stepsize which is not related to dispersion property. Another classes of numerical methods for solving (1.1) are exponentially fitted or phase fitted in which the period or frequency is known in advance (see e.g., [1221]).

Most of the numerical methods developed for solving (1.1) are in constant stepsize (see [13, 6, 22, 23]). In this paper, the development of efficient DIRKN methods with reduced phase-lag in variable stepsize is studied. The strategies introduced in Dormand et al. [24] and Simos [18] for finding the optimized pair is used and also new implementation code is discussed in this paper.

In this paper, dispersion relations are developed and used together with algebraic conditions to be solved explicitly. In the following section, the construction of the new 4(3) pairs of DIRKN method is described. Its coefficients are given using the Butcher tableau notation. Finally, numerical tests on second-order differential equation problems possessing an oscillatory solutions are performed.

2. Analysis of Phase-Lag and Stability

In this section, we will discuss the analysis of phase-lag for RKN method. The first analysis of phase-lag was carried out by Bursa and Nigro [25], then followed by Gladwell and Thomas [26] for the linear multistep method, and for explicit and implicit Runge-Kutta(-Nyström) methods by van der Houwen and Sommeijer [1, 2].

The phase-lag analysis of the method (1.2)-(1.3) is investigated using the homogeneous test equation

By applying the general method (1.2)-(1.3) to the test equation (2.1) yields where . Here is the stability matrix of the RKN method and its characteristic polynomial is the stability polynomial of the RKN method. Solving the difference system (2.2), the computed solution is given by The exact solution of (2.1) is given by Equations (2.4) and (2.5) led to the following definition.

Definition 2.1 (phase-lag). Apply the RKN method (1.2)-(1.3) to (1.1). Next, we define the phase-lag . If , then the RKN method is said to have phase-lag order . Additionally, the quantity is called amplification error. If , then the RKN method is said to have dissipation order .

Let us denote From Definition 2.1, it follows that Let us denote and in the following form: where is diagonal element in the Butcher tableau.

Based on the functions and defined as (2.8), a few properties of the functions and are summarized in the following theorem which is introduced by Van der Houwen and Sommeijer [2]. The development of dispersion relations is according to the following theorem.

Theorem 2.2. (1) The functions and are consistent, dispersive, and dissipative of orders , and , respectively,
(2) An RKN method of algebraic order , dispersion of order , and dissipation order possess functions R and S that are consistent, dispersive, and dissipative of orders , and .
(3) If , then the order of consistency and dispersion of and is equal.

Proof (see Van der Houwen and Sommeijer [2]). Based on the above theorem, the dispersion relations are developed. For the dispersion relation of order six () in terms of and is and for the dispersion relations up to order eight () for are given by
The following quantity is used to determine the dissipation constant of the formula.(i) for (ii) for

Notice that the fourth-order method is already dispersive order four and dissipative order five due to consistency of the method. Furthermore, dispersive order is even and dissipative order is odd.

We next discuss the stability properties of method for solving (1.1) by considering the scalar test problem (2.1) applied to the method (1.2)-(1.3), from which the expression given in (2.2) is obtained. Eliminating and in (2.2), we obtain a difference equation of the form The characteristic equation associated with (2.15) is given as in (2.3). Since our concerned here is with the analysis of high-order dispersive RKN method, we therefore drop the necessary condition of periodicity interval that is, . Hence, the method derived will be with empty interval of periodicity. We now consider the interval of absolute stability of RKN method. We therefore need the characteristic equation (2.3) to have roots with modulus less than one so that approximate solution will converge to zero as tends to infinity. For convenience, we note the following definition adopted for method (2.2).

Definition 2.3. An interval is called the interval of absolute stability of the method (2.2) if for all , .

3. Construction of the Method

In the following, we will derive a three-stage fourth-order and a four-stage fourth-order accurate DIRKN method with dispersive order six and eight, respectively, by taking into account the dispersion relation in Section 2. The RKN parameters must satisfy the following algebraic conditions to find fourth-order accuracy as given in Hairer and Wanner [27]: For most methods, the need to satisfy

The following strategies are used for developing our new efficient pairs.(1) The high-order DIRKN formula with high order of dispersion. Our aim is to find the ratio (phase-lag order/algebraic order) as high as possible and the dissipation constant is “small”. (2) The following quantities as in [24] should be as small as possible:(a), (b), where are called error coefficients for and , respectively.(3) The strategy to control the error is based on the phase-lag procedure first introduced by Simos [18] and also see [19, 20]. A local error estimation at the point is determined by the expressions where and are solutions using the third-order formula. These local error estimations can be used to control the step size by the standard formula [28, 29] where 0.9 is a safety factor, represents the local error estimation at each step and Tol is the accuracy required which is the maximum allowable local error. If Est Tol, then the step is accepted, and we applied the accepted procedure of performing local extrapolation (or higher-order mode) meaning that the more accurate approximation will be used to advance the integration. If Est Tol, then the step is rejected and the will be updated using formula (3.8).

3.1. The Three-Stage DIRKN Formula
3.1.1. The Fourth-Order Formula

In this section, we derive the fourth-order formula with dispersive order six and dissipative order five. Sharp et al. [9] stated that fourth-order method with dispersive order eight does not exist. Therefore, the method of algebraic order four ( with dispersive order six ( and dissipative order five ( is now considered. From algebraic conditions (3.1)–(3.4) and (3.7), it formed eleven equations with thirteen unknowns to be solved. We let and be a free parameter. Therefore, the following solution of one-parameter family is obtained:

From the above solution, we are going to derive a method with dispersion of order-six. The six order dispersion relation (2.10) needs to be satisfied and this can be written in terms of RKN parameters which correspond to the above family of solution yields the following equation: and solving for yields The first two values will give us a nonempty stability interval while the others will produce the methods with empty stability interval. Taking the first two values of gives us two fourth-order diagonally implicit RKN methods with dispersive order six. For it will give the method with PLTE for and , respectively. The dissipation constant and the stability interval are and , respectively.

3.1.2. The Third-Order Formula

Based on the values of and in Section 3.1.1, we now derive the three-stage third-order embedded formula. Solving equations (3.1)–(3.3) simultaneously yields a solution for and in terms of , and is the same as the fourth-order formula in Section 3.1.1. With this solution, and are functions in terms of . Next, we plot the graph for and against . We only consider with and lying between , and , respectively. From numerical experiments, the optimal pair and giving and , respectively, and . We denote this pair as DIRKN4(3)6New method (see Table 2).

tab2
Table 2: The DIRKN4(3)6New method.
3.2. The Four-Stage DIRKN Formula
3.2.1. The Fourth-Order Formula

To derive four-stage fourth-order () with dispersive order eight (), (3.1)–(3.4) and (3.7) together with the dispersion relation of order six equation (2.11) are solved simultaneously and will yield the following solution:

By substituting the above solution to the dispersion relation of order eight (), (2.12) gives us expression in terms of and solving for will give us the values , , , , , , and . Numerical results show that choosing will give us smallest dissipation constant hence more accurate scheme compared to the others. The dissipation constant and the stability interval are and , respectively.

3.2.2. The Third-Order Formula

Based on the values of and in Section 3.2.1, here we solve third-order embedded formula for the values of and obtaining where and are free parameters. The and are functions in terms of and while and are in terms of . By setting , we have and in terms of .

Similarly, we plot the graph for and against . Consider with and lying between and , respectively, while with and lying between and , respectively. From numerical experiments, the optimal pair chosen are and . With these values will give and . The pair we denote by DIRKN4(3)8New method (see Table 3).

tab3
Table 3: The DIRKN4(3)8New method.

4. Problems Tested

In order to evaluate the effectiveness of the new embedded methods, we solved several problems which have oscillatory solutions. The code developed uses constant and variable step size mode and the results obtained are compared with the methods proposed in [10, 11, 21, 29, 30]. Table 4 and Figures 1, 2, 3, 4, 5, 6, 7, and 8 show the numerical results for all methods used. These codes have been denoted by the following:(i)DIRKN4(3)8New: a new 4(3) pair with phase-lag order 8 derived in this paper. (ii)DIRKN4(3)6New: a new 4(3) pair with phase-lag order 6 derived in this paper. (iii)RKND: a fourth-order explicit RKN derived by Dormand [29]. (iv)FPRKN: a phase-fitted fourth-order explicit RKN derived by Papadopoulos et al. [21]. (v)DIRKN4(3)Imoni: a 4(3) RKN pair of three-stage derived by Imoni et al. [10]. (vi)DIRKN4(3)Raed: a 4(3) RKN pair derived by Al-Khasawneh et al. [11]. (vii)DIRK4(3)Butcher: a 4(3) Runge-Kutta pair with six-stage, L-stable, and FSAL property derived by Butcher and Chen [30].

tab4
Table 4: Numerical results for Problem 1 to Problem 8 with .
324989.fig.001
Figure 1: Efficiency curves for Problem 1.
324989.fig.002
Figure 2: Efficiency curves for Problem 2.
324989.fig.003
Figure 3: Efficiency curves for Problem 3.
324989.fig.004
Figure 4: Efficiency curves for Problem 4.
324989.fig.005
Figure 5: Efficiency curves for Problem 5.
324989.fig.006
Figure 6: Efficiency curves for Problem 6.
324989.fig.007
Figure 7: Efficiency curves for Problem 7.
324989.fig.008
Figure 8: Efficiency curves for Problem 8.

Problem 1 (homogenous). One has Exact solution

Problem 2 (inhomogeneous problem studied by Allen and Wing [31]). One has Exact solution

Problem 3 (problem studied by van der Houwen and Sommeijer [1]). One has where .
Exact solution is . Numerical result is for the case .

Problem 4 (inhomogeneous system studied by Franco [8]). One has Exact solution

Problem 5 (The Duffing's equations as given in [32]). One has The exact solution computed by the Galerkin method and given by with and .

Problem 6 (Inhomogeneous problem studied by Lambert and Watson [33]). One has Exact solution is is chosen to be and parameters and are 20 and 0.1, respectively.

Problem 7 (an almost periodic orbit problem given in stiefel and bettis [34]). One has Exact solution , .

Problem 8 (linear Strehmel-Weiner problem studied in Cong [35]). One has Exact solution

5. Numerical Results

In this section, we evaluate the effectiveness of the new DIRKN pairs derived in the previous section when they are applied to the numerical solution of eight problems which are model and nonmodel problems for constant and variable step size.

For constant step size, Table 4 shows the absolute maximum error for fourth-order DIRKN4(3)6New, DIRKN4(3)8New, PFRKN, and RKND methods when solving Problems 18 with three different step values. From numerical results in Table 4, we observed that the new DIRKN4(3)8New method is more accurate compared with PFRKN and RKND method which is not related to the dispersion of the method. Also the new DIRKN4(3)6New method is more accurate compared with RKND method while comparable with PFRKN method for certain problem. Notice that all the methods are of the same algebraic order.

For variable step size, Figures 18 show the decimal logarithm of the maximum global error for the solution (MAXE) versus the function evaluations and also the decimal logarithm of the maximum global error for the solution (MAXE) versus the time taken. From Figures 18, we observed that DIRKN4(3)6New and DIRKN4(3)8New performed better compared to DIRKN4(3)Raed, DIRKN4(3)Imoni, and DIRK4(3)Butcher for integrating second-order differential equations possessing an oscillatory solution in terms of function evaluations. This is due to the fact that when using DIRK4(3)Butcher, the second-order system of ODEs needs to be transformed to a first-order system and hence the dimension is doubled. Furthermore, the DIRK4(3)Butcher method has five effective stages per step. This means that the function evaluations per step for DIRK4(3)Butcher are more than the function evaluations per step used in the DIRKN4(3)6New and DIRKN4(3)8New methods which have three and four function evaluations per step, respectively. Meanwhile the DIRKN4(3)Raed, DIRKN4(3)Imoni methods are less accurate compared with the DIRKN4(3)6New and DIRKN4(3)8New when phase lag and dissipation is considered. Therefore the new methods converges faster and consequently less steps are needed for a specified value of tolerance even though DIRKN4(3)Raed and DIRKN4(3)Imoni have the same number of stages per step with DIRKN4(3)8New and DIRKN4(3)6New method respectively.

6. Conclusion

In this paper, we have derived two 4(3) pairs, namely, DIRKN4(3)6New and DIRKN4(3)8New which have dispersive order six and eight, respectively, with “small” dissipation constant which is suitable for problems with oscillating solutions and moderate frequency. From the results shown in Table 4 and Figures 18, we conclude that the new methods are more efficient for integrating second-order equations possessing an oscillatory solution when compared to the RKND derived in [29], PFRKN derived in [21], DIRK4(3)Butcher as in [30] and also with the others the same type of pairs for example, DIRKN4(3)Imoni pair derived in [10] and DIRKN4(3)Raed derived in [11]. The DIRKN4(3)6New method is the most efficient method since it has three function evaluations per step.

Acknowledgments

This work is partially supported by IPTA Fundamental Research Grant, Universiti Putra Malaysia (Project no. 05-10-07-385FR) and UPM Research University Grant Scheme (RUGS) (Project no. 05-01-10-0900RU).

References

  1. 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. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  2. P. J. van der Houwen and B. P. Sommeijer, “Diagonally implicit Runge-Kutta-Nyström methods for oscillatory problems,” SIAM Journal on Numerical Analysis, vol. 26, no. 2, pp. 414–429, 1989. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  3. A. B. Sideridis and T. E. Simos, “A low-order embedded Runge-Kutta method for periodic initial value problems,” Journal of Computational and Applied Mathematics, vol. 44, no. 2, pp. 235–244, 1992. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  4. A. García, P. Martín, and A. B. González, “New methods for oscillatory problems based on classical codes,” Applied Numerical Mathematics, vol. 42, no. 1–3, pp. 141–157, 2002. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  5. N. Senu, M. Suleiman, and F. Ismail, “An embedded explicit Runge–Kutta–Nyström method for solving oscillatory problems,” Physica Scripta, vol. 80, no. 1, 2009. View at Publisher · View at Google Scholar
  6. H. Van de Vyver, “A symplectic Runge-Kutta-Nyström method with minimal phase-lag,” Physics Letters A, vol. 367, no. 1-2, pp. 16–24, 2007. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  7. N. Senu, M. Suleiman, F. Ismail, and M. Othman, “A zero-dissipative Runge-Kutta-Nyström method with minimal phase-lag,” Mathematical Problems in Engineering, vol. 2010, Article ID 591341, 15 pages, 2010. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  8. J. M. Franco, “A 5(3) pair of explicit ARKN methods for the numerical integration of perturbed oscillators,” Journal of Computational and Applied Mathematics, vol. 161, no. 2, pp. 283–293, 2003. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  9. P. W. Sharp, J. M. Fine, and K. Burrage, “Two-stage and three-stage diagonally implicit Runge-Kutta Nyström methods of orders three and four,” IMA Journal of Numerical Analysis, vol. 10, no. 4, pp. 489–504, 1990. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  10. S. O. Imoni, F. O. Otunta, and T. R. Ramamohan, “Embedded implicit Runge-Kutta Nyström method for solving second-order differential equations,” International Journal of Computer Mathematics, vol. 83, no. 11, pp. 777–784, 2006. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  11. R. A. Al-Khasawneh, F. Ismail, and M. Suleiman, “Embedded diagonally implicit Runge-Kutta-Nystrom 4(3) pair for solving special second-order IVPs,” Applied Mathematics and Computation, vol. 190, no. 2, pp. 1803–1814, 2007. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  12. Z. A. Anastassi and T. E. Simos, “An optimized Runge-Kutta method for the solution of orbital problems,” Journal of Computational and Applied Mathematics, vol. 175, no. 1, pp. 1–9, 2005. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  13. T. E. Simos, “Closed Newton-Cotes trigonometrically-fitted formulae of high order for long-time integration of orbital problems,” Applied Mathematics Letters, vol. 22, no. 10, pp. 1616–1621, 2009. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  14. S. Stavroyiannis and T. E. Simos, “Optimization as a function of the phase-lag order of nonlinear explicit two-step P-stable method for linear periodic IVPs,” Applied Numerical Mathematics, vol. 59, no. 10, pp. 2467–2474, 2009. View at Publisher · View at Google Scholar
  15. T. E. Simos, “Exponentially and trigonometrically fitted methods for the solution of the Schrödinger equation,” Acta Applicandae Mathematicae, vol. 110, no. 3, pp. 1331–1352, 2010. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  16. T. E. Simos, “New stable closed Newton-Cotes trigonometrically fitted formulae for long-time integration,” Abstract and Applied Analysis, vol. 2012, Article ID 182536, 15 pages, 2012. View at Publisher · View at Google Scholar
  17. T. E. Simos, “Optimizing a hybrid two-step method for the numerical solution of the Schrödinger equation and related problems with respect to phase-lag,” Journal of Applied Mathematics, vol. 2012, Article ID 420387, 17 pages, 2012. View at Publisher · View at Google Scholar
  18. T. E. Simos, “Embedded Runge-Kutta methods for periodic initial value problems,” Mathematics and Computers in Simulation, vol. 35, no. 5, pp. 387–395, 1993. View at Publisher · View at Google Scholar
  19. G. Avdelas and T. E. Simos, “Block Runge-Kutta methods for periodic initial-value problems,” Computers & Mathematics with Applications, vol. 31, no. 1, pp. 69–83, 1996. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  20. G. Avdelas and T. E. Simos, “Embedded methods for the numerical solution of the Schrödinger equation,” Computers & Mathematics with Applications, vol. 31, no. 2, pp. 85–102, 1996. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  21. D. F. Papadopoulos, Z. A. Anastassi, and T. E. Simos, “A phase-fitted Runge-Kutta-Nyström method for the numerical solution of initial value problems with oscillating solutions,” Computer Physics Communications, vol. 180, no. 10, pp. 1839–1846, 2009. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  22. 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. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  23. 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. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  24. J. R. Dormand, M. E. A. El-Mikkawy, and P. J. Prince, “Families of Runge-Kutta-Nyström formulae,” IMA Journal of Numerical Analysis, vol. 7, no. 2, pp. 235–250, 1987. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  25. L. Bursa and L. Nigro, “A one-step method for direct integration of structural dynamic equations,” International Journal for Numerical Methods in Engineering, vol. 15, pp. 685–699, 1980.
  26. I. Gladwell and R. M. Thomas, “Damping and phase analysis for some methods for solving second-order ordinary differential equations,” International Journal for Numerical Methods in Engineering, vol. 19, no. 4, pp. 495–503, 1983. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  27. E. Hairer and G. Wanner, “A theory for Nyström methods,” Numerische Mathematik, vol. 25, no. 4, pp. 383–400, 1975/76. View at Publisher · View at Google Scholar
  28. J. C. Butcher, Numerical Methods for Ordinary Differential Equations, John Wiley & Sons, Chichester, UK, 2nd edition, 2008. View at Publisher · View at Google Scholar
  29. J. R. Dormand, Numerical Methods for Differential Equations, CRC Press, Boca Raton, Fla, USA, 1996.
  30. J. C. Butcher and D. J. L. Chen, “A new type of singly-implicit Runge-Kutta method,” Applied Numerical Mathematics, vol. 34, no. 2-3, pp. 179–188, 2000. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  31. R. C. Allen, Jr. and G. M. Wing, “An invariant imbedding algorithm for the solution of inhomogeneous linear two-point boundary value problems,” Journal of Computational Physics, vol. 14, pp. 40–58, 1974. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  32. H. Van de Vyver, “A Runge-Kutta-Nyström pair for the numerical integration of perturbed oscillators,” Computer Physics Communications, vol. 167, no. 2, pp. 129–142, 2005.
  33. J. D. Lambert and I. A. Watson, “Symmetric multistep methods for periodic initial value problems,” Journal of the Institute of Mathematics and its Applications, vol. 18, no. 2, pp. 189–202, 1976. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  34. E. Stiefel and D. G. Bettis, “Stabilization of Cowell's method,” Numerische Mathematik, vol. 13, no. 2, pp. 154–175, 1969. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  35. N. H. Cong, “A-stable diagonally implicit Runge-Kutta-Nyström methods for parallel computers,” Numerical Algorithms, vol. 4, no. 3, pp. 263–281, 1993. View at Publisher · View at Google Scholar · View at Zentralblatt MATH