`Journal of Applied MathematicsVolume 2011, Article ID 407151, 15 pageshttp://dx.doi.org/10.1155/2011/407151`
Research Article

## -Stable Higher Derivative Methods with Minimal Phase-Lag for Solving Second Order Differential Equations

Department of Mathematics, Faculty of Science, King Abdul-Aziz University, Jeddah, Saudi Arabia

Received 2 July 2011; Accepted 3 September 2011

Copyright © 2011 Fatheah A. Hendi. 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

Some new higher algebraic order symmetric various-step methods are introduced. For these methods a direct formula for the computation of the phase-lag is given. Basing on this formula, calculation of free parameters is performed to minimize the phase-lag. An explicit symmetric multistep method is presented. This method is of higher algebraic order and is fitted both exponentially and trigonometrically. Such methods are needed in various branches of natural science, particularly in physics, since a lot of physical phenomena exhibit a pronounced oscillatory behavior. Many exponentially-fitted symmetric multistepmethods for the second-order differential equation are already developed. The stability properties of several existing methods are analyzed, and a new -stable method is proposed, to establish the existence of methods to which our definition applies and to demonstrate its relevance to stiff oscillatory problems. The work is mainly concerned with two-stepmethods but extensions tomethods of larger step-number are also considered. To have an idea about its accuracy, we examine their phase properties. The efficiency of the proposed method is demonstrated by its application to well-known periodic orbital problems. The new methods showed better stability properties than the previous ones.

#### 1. Introduction

In the recent years, there is much activity concerned with the numerical solution of the of the type where is a function in which the first derivative does not appear explicitly. In addition the function also satisfies the Lipschitz condition of the first order with respect to . There are two main groups of methods for the numerical solution of the form given in (1.1) with periodic of an oscillatory solutions, see [13]. The first consists of methods in which the coefficients that determine the numerical schemes depend on the particular problem to be solved. The second consists of methods with constant coefficients. These are of great interest in many applications such as orbital mechanics, theoretical physics, and electronics.

Numerous second-category numerical techniques have been obtained for the solution of (1.1). Methods of this category must be -stable, and this applies especially on the problem which has highly oscillatory solutions. The -stability property was first introduced by Awoyemi [1]. On the other hand, sixth-order -stable methods are obtained in Cash et al. [4]. An important contribution for these methods is shown in Hairer and Wanner [5] where the lower-order -stable methods are developed.

Indeed, standard numerical methods may require huge number of time steps to track the oscillations, such method should be chosen very carefully, and the best choice is strictly application dependent.

The aim of this paper is to develop an efficient free parameter class of -stable method with minimal phase-lag. The derivation of this method allows free parameters to lead to an efficient implementation method. The numerical tests show that these new classes of methods are more efficient than the other well-known -stable methods. This is due to the value of the free parameter, stability properties and the order of phase-lag are depending on one or more off points than the other well-known -stable methods. Consequently we will describe the basic theory of stability, (phase-lag of symmetric multistep methods) and develop higher-orders -stable methods.

A very well-known family of multistep methods for the solution of (1.1) is those family known as Strmer-Cowell methods. These methods have been used many times for long-term integrations of planetary orbits (see Quinlan and Tremaine [6] and references therein), and present a problem, called orbital instability when the number of steps exceeds 2. To solve the problem of orbital instability, the works of Lambert and Watson [7] have introduced the symmetric multistep methods and showed that the symmetric methods have nonvanishing interval of periodicity which is the interval of guaranteed periodic solution, see [810]. This is determined by the application of the symmetric multistep method to a test equation given by . If where is the step length of the integration, then this interval is called interval of periodicity). In addition the authors in [6] have produced higher-order symmetric multistep methods based on the work given by Lambert and Watson [7]. It was shown that the linear symmetric multistep methods developed in [6, 7, 11] are much simpler than the hybrid Runge-Kutta methods. For the above reasons of simplicity and accuracy in long-time integration of periodic initial value problems (and especially orbital problems), we give attention to this family of methods.

In the present paper we construct an exponentially fitted explicit symmetric multistep method, based on the symmetric higher algebraic-order method developed by Quinlan and Tremaine in [6] and Ixaru et al. [12]. Stability analysis is presented in Section 2. In Section 3 we present the construction of the new method.

Moreover, In this paper definitions of the periodicity interval and -stability, which are designed for linear multistep methods with constant coefficients given by [6, 13, 14] are modified. The modification is performed to provide a basis for linear stability analysis of exponential fitting methods for the special class of ordinary differential equations of second order in which the first derivative does not appear explicitly. The stability properties of several existing methods are analyzed, and a new P-stable method is proposed, to establish the existence of methods to which our definition applies and to demonstrate its relevance to stiff oscillatory problems. The work is mainly concerned with multistep methods, but extensions to methods of larger step number are also considered.

#### 2. Formulation

##### 2.1. Stability and Periodicity

For problems with oscillatory solutions, linear stability analysis is based on the test equation of the form where is a real constant. This test equation is previously introduced by Lambert and Watson [7] as well as the interval of its periodicity, in order to investigate the periodic stability properties of numerical method for solving the initial value problem given in (1.1). Stability means that the numerical solutions remain bounded as we move further away from the starting point, see Coleman [15], Simos [16], and Simos and Williams [17].

In order to investigate the periodic stability properties of numerical methods for solving the initial-value problem (1.1), in [7] they introduce the scalar test equation Based on the theory developed in [7, 18], when a symmetric multistep method, given by is applied to the scalar test equation (2.1), a difference equation is obtained of the form: where , is the step length, and is the computed approximation to . The general solution of the above difference equation can be written as where , are the distinct roots of the polynomial: where and are polynomials given by We note here that the roots of the polynomial (2.5) are perturbations of the roots of (2.6). We denote and the perturbations of the principal roots of .

Based on [7] when a symmetric multistep method is applied to the scalar test equation (2.1), a difference equation (2.3) is obtained. The characteristic equation associated with (2.3) is given by (2.5). The roots of the characteristic polynomial (2.5) are denoted as .

According to Lambert and Watson [7] the following definitions are that given by:

We have the following definitions.

Definition 2.1 (see [7]). The numerical method (2.3) has an interval of periodicity , if, for all are complex and satisfy:

Definition 2.2 (see [7]). The method given by (2.3) is -stable if its interval of periodicity is .

As a modification of this method this paper aimed to develop a family of symmetric -stable two step, four step and six step methods involving higher-order derivatives with minimal phase-lag error in the form: With the help of the associated formula, This latter formula depends on an offpoint . The precise choice enables getting a -stable formula with large interval of periodicity.

##### 2.2. Construction of Two-Step -Stable Higher-Order Derivative with Phase Fitted Schemes

Considering the two-step, -stable formula that depends on one of point involving higher-order derivative in the general form: where is a constant; the basic idea behind our approach is to approximate by the expression involving the quantities , and . We introduce -stable methods of , , and with minimal phase-lag errors using (2.10) with the associated formula (2.9). This gives Applying the method given by (2.10) and (2.11) with to the scalar test problem (2.1), we obtain the following difference equation: where is the step length, and and are polynomials in .

The characteristic equation associated with (2.12) is obtained: The roots of (2.13) are complex and of modulus one, see [7, 19], if

Under this condition the roots of (2.13) can be written as where . The numerical solution of (2.13) is bounded if both roots are unequal and their magnitude less than one or equal to one.

The periodicity condition requires those roots to lie on the unit circle, that is, is then a rational approximation or , [2022]. A method is said to be -stable if the interval of periodicity is infinite.

For a given method (i.e., a given ), one has to find a restriction which must be placed on the step length to ensure that the condition is satisfied.

The following definition is taken from the work by [20, 23]:

Definition 2.3 (see [19, 24]). The method (2.4) is unconditionally stable and for all values of , .

Definition 2.4 (see [6]). The numerical method (2.4) has an interval of periodicity if, for all , and satisfy where is a real function of . For any method corresponding to the characteristic equation (2.13), the phase-lag is defined by Hairer and Wanner in [5] as the leading term in the following expression: If as , the order of phase-lag is .

Due to Definition 2.2, the method (2.10) is -stable if its interval of periodicity is .

Theorem 2.5 (see [19]). A method which has the characteristic equation (2.13) has an interval of periodicity if, for all ,

Definition 2.6. A region of stability is a region of the plane, throughout which . Any closed curve defined by is a stability boundary.

When , the roots of (2.13) are distinct and lie on the unit circle.

When , the method is unstable since the corresponding difference equation has an unbounded solution, see [25, 26].

Theorem 2.7 (see [19]). For a method which has an interval of periodicity , one gets

Remark 2.8 (see [27]). If the phase-lag order is , then

Definition 2.9 (see [15]). A symmetric two-step method phase fitted has phase-lag of order infinity: So, we see that the methods represented by (2.10) and (2.11) are phase fitted if .

Definition 2.10. A family of phase-fitted methods with the stability function , where , is -stable if, for each value of , the quantity is satisfied for all values of and for all values of except possibly a set of exceptional values of determined by the chosen value of .

By the last definition we have assumed that ,

###### 2.2.1. Two-Step -Stable Involve Second Derivative with Minimal Phase-Lag Errors

Consider the symmetric family of second derivative, two-step methods of (2.10) and (2.11) as

Therefore, we introduce -stable methods of higher order with minimal phase-lag errors using the formula: with the coefficients: where , and are free parameters.

Also, the local truncation error of the method is given by We denote this family of methods . For an applied to the test equation (2.1), setting , we obtain (2.13), and we get Then we get As given in (2.14), the method represented by (2.22), (2.23) with (2.24) will be -stable provided for all ; consequently, it is easy to prove that Following the phase-lag, which is denoted by , as a leading coefficient in the expression of , then The method has phase-lag error of satisfying the -stability conditions. we can obtain a symmetric -stable second derivative, two-step method given by (2.22) with phase-lag error of order , as the following cases.

Case 1 (1.1).
Let , and

Case 2 (2.1). , and

Case 3 (2.2).
Let , and

Case 4 (2.3). for

###### 2.2.2. Fourth-Order Scheme

As in Section 2.2.1 we introduced -stable methods of higher order with minimal phase-lag errors using the formulae (2.22) and (2.23) with the coefficients: where , , and are free parameters.

Also, the local truncation error of the method is given by For an , setting , (2.13) is obtained with the coefficient given by, and we get

The coefficient of . the coefficient of ,

We can obtain a symmetric -stable second derivative, two-step method shown by (2.22) if and we have the following.(1), with phase-lag error of order: (2), with phase-lag error of order:

###### 2.2.3. Sixth-Order Scheme

The two-step -stable methods of order six with minimal phase-lag errors using the formula associated with the formula Similarly, we can obtain the following results: The -stability conditions will be

The method has phase-lag error of .

#### 3. Construction of Four-Step -Stable Higher-Order Derivative with Minimal Phase-Lag Errors

Consider the symmetric four-step methods in the form: Associated with the characteristic equation will be where The -stability conditions are

Following the phase-lag, denoted by , is leading coefficient in the expression of The local truncation error of the method is given by

##### 3.1. Second Derivative of Second-Order Scheme

We can establish some choices of the parameters to obtain symmetric four step -stable methods involve second derivative with minimal phase-lag errors tabulated in Table 1.

Table 1
##### 3.2. Second Derivative of Fourth-Order Scheme

Let , the coefficients of (3.1) and (3.2) will be Also, the local truncation error of the method is given by The coefficients of the characteristic equation will be The -stability conditions are We can obtain a symmetric -stable second derivative, four-step methods with minimal phase-lag errors by suitable choices of the free parameters, see [810].

##### 3.3. Fourth Derivative of Sixth-Order Scheme

Consider the symmetric four-step methods involve fourth derivative in the form (3.1) with (3.2) to obtain the characteristic equation: with the following coefficients:

Now we have some cases for -stable, symmetric four-step methods involve fourth derivative one is as follows.

Let , , let , let , let ,let , and let with , and as free parameters.

#### 4. Conclusions

In this paper a higher algebraic order exponentially fitted free-parameters method is developed. We have given explicitly the way for the construction of the method. Stability analysis of the new method is also presented. The numerical results, so far obtained in this paper, show the efficiency of the newly derived integrator of order five. We also observed that, for an exponentially fitted problems, our integrator do not use small step lengths, as may be required by many multistep methods before good accuracy is obtained. We exploit the freedom in the selection of the free parameters of one family with the purpose of obtaining specific class of the highest possible phase-lag order, which are also characterized by minimized principal truncation error coefficients. Finally, the new integrator derived in this paper is capable of handling stiff problems for which exponential fitting is applicable.

#### References

1. D. O. Awoyemi, “A P-stable linear multistep method for solving general third order ordinary differential equations,” International Journal of Computer Mathematics, vol. 80, no. 8, pp. 987–993, 2003.
2. C. E. Abhulimen and S. A. Okunuga, “Exponentially-fitted second derivative multistep methods for stiff Initial value problems in ordinary differential equations,” JESA: Journal of Engineering Science And Applications, vol. 5, no. 1-2, pp. 36–49, 2008.
3. C. E. Abhulimen and F. O. Otunta, “A family of two step exponentially fitted multiderivative methods for the numerical integration of stiff IVPs in ODEs,” IJNM: International Journal of Numerical Mathematics, vol. 2, no. 1, pp. 1–21, 2007.
4. J. R. Cash, A. D. Raptis, and T. E. Simos, “A sixth-order exponentially fitted method for the numerical solution of the radial schrodinger equation,” Journal of Computational Physics, vol. 91, no. 2, pp. 413–423, 1990.
5. E. Hairer and G. Wanner, Solving Ordinary Differential Equations. II, vol. 14 of Springer Series in Computational Mathematics, Springer, Berlin, Germany, 1991.
6. G. D. Quinlan and S. Tremaine, “Symmetric multistep methods for the numerical integration of planetary orbits,” The Astronomical Journal, vol. 100, pp. 1694–1700, 1990.
7. 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.
8. A. D. Raptis and T. E. Simos, “A four-step phase-fitted method for the numerical integration of second order initial value problems,” BIT Numerical Mathematics, vol. 31, no. 1, pp. 160–168, 1991.
9. T. E. Simos, “A family of four-step trigonometrically-fitted methods and its application to the Schrödinger equation,” Journal of Mathematical Chemistry, vol. 44, no. 2, pp. 447–466, 2008.
10. 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.
11. L. Gr. Ixaru and G. Vanden Berghe, Exponential Fitting, vol. 568 of Mathematics and Its Applications, Kluwer Academic Publishers, Dordrecht, The Netherlands, 2004, With 1 CD-ROM (Windows, Macintosh and UNIX.
12. L. Gr. Ixaru, G. Vanden Berghe, and H. De Meyer, “Frequency evaluation in exponential fitting multistep algorithms for ODEs,” vol. 140, no. 1-2, pp. 423–434.
13. C. E. Abhulimen and F. O. Otunta, “A new class of exponential fitting numerical integration for initial value problems in ordinary differential equations,” Nigerian Mathematical society, vol. 28, pp. 1–14, 2009.
14. Ch. Tsitouras and T. E. Simos, “High algebraic, high phase-lag order embedded Numerov-type methods for oscillatory problems,” Applied Mathematics and Computation, vol. 131, no. 1, pp. 201–211, 2002.
15. J. P. Coleman, “Numerical methods for ${y}^{″}=f\left(x,y\right)$,” in Proceedings of the 1st International Colloquium on Numerical Analysis, D. Bainov and V. Civachev, Eds., pp. 27–38, VSP, Utrecht, The Netherlands.
16. T. E. Simos, “Atomic structure coputations,” in Atomic Structure Computations in Chemical Modelling: Applications and Theory, A. Hinchliffe, Ed., pp. 38–142, UMIST, The Royal Society of Chemistry, 2000.
17. T. E. Simos and P. S. Williams, “On finite difference methods for the solution of the Schrödinger equation,” Journal of Computer Chemistry, vol. 23, pp. 513–554, 1999.
18. J. D. Lambert, Numerical Methods for Ordinary Differential Systems, John Wiley & Sons, Chichester, UK, 1991.
19. C. E. Abhulimen, On exponentially fitted multiderivative numerical methods for stiff initial value problems in ordinary differential equations, Ph.D. thesis, Department of Mathematics, University of Benin, Benin City, Nigeria, 2006.
20. J. P. Coleman and L. Gr. Ixaru, “P-stability and exponential-fitting methods for ${y}^{″}=f\left(x,y\right)$,” IMA Journal of Numerical Analysis, vol. 16, no. 2, pp. 179–199, 1996.
21. K. Eriksson, C. Johnson, and A. Logg, “Explicit time-stepping for stiff ODEs,” SIAM Journal on Scientific Computing, vol. 25, no. 4, pp. 1142–1157, 2003.
22. T. E. Simos and J. Vigo-Aguiar, “Symmetric eighth algebraic order methods with minimal phase-lag for the numerical solution of the Schrödinger equation,” Journal of Mathematical Chemistry, vol. 31, no. 2, pp. 135–144, 2002.
23. J. Vigo-Aguiar, J. Martín-Vaquero, and R. Criado, “On the stability of exponential fitting BDF algorithms,” Journal of Computational and Applied Mathematics, vol. 175, no. 1, pp. 183–194, 2005.
24. T. E. Simos and P. S. Williams, “A P-stable hybrid exponentially-fitted method for the numerical integration of the Schrödinger equation,” Computer Physics Communications, vol. 131, pp. 109–119, 2000.
25. C. E. Abhulimen and F. O. Otunta, “A Sixth order multiderivative multistep methods for stiff system of differential equations,” IJNM: International Journal of Numerical Mathematics, vol. 2, no. 1, pp. 248–268, 2006.
26. J. C. Butcher, Numerical Methods for Ordinary Differential Equations, John Wiley & Sons, Chichester, UK, 2003.
27. G. Psihoyios and T. E. Simos, “Exponentially and trigonometrically fitted explicit advanced step-point (EAS) methods for initial value problems with oscillating solutions,” International Journal of Modern Physics C, vol. 14, no. 2, pp. 175–184, 2003.