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

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 IVPs of the type 𝑦^′′=𝑓𝑥,𝑦,𝑦𝑥₀=𝑦₀,𝑦^′𝑥₀=𝑦^′₀, (1.1) 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 [1–3]. 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 IVPs 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 Stö‌rmer-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 [8–10]. This is determined by the application of the symmetric multistep method to a test equation given by 𝑦^′′(𝑡)=−𝜆²𝑦(𝑡). If 𝜆²ℎ²∈(0,𝜒²₀) 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𝑦^′′=−𝜆²𝑦. (2.1) Based on the theory developed in [7, 18], when a symmetric multistep method, given by ^𝑘_𝑗=0𝑎_𝑗𝑦_𝑛+𝑖=ℎ^2𝑘_𝑗=0𝛽_𝑗𝑓_𝑛+𝑗 (2.2) is applied to the scalar test equation (2.1), a difference equation is obtained of the form: ^𝑘_𝑖=0𝑎_𝑖+𝜒²𝛽_𝑖𝑦_𝑛+𝑖=0, (2.3) where 𝜒=𝜆ℎ, ℎ is the step length, and 𝑦_𝑛 is the computed approximation to 𝑦(𝑛ℎ),𝑛=0,1,2,…. The general solution of the above difference equation can be written as𝑦_𝑛=^𝑘_𝑗=0𝐴_𝑗𝜉^𝑛_𝑗, (2.4) where 𝜉_𝑗, 𝑗=1(1)𝑘 are the distinct roots of the polynomial:𝑃𝜉;𝜒²=𝜌(𝜉)+𝜒²𝜎(𝜉), (2.5) where 𝜌 and 𝜎 are polynomials given by 𝜌(𝜉)=^𝑘_𝑖=0𝑎_𝑖𝜉^𝑖;𝜎(𝜉)=^𝑘_𝑖=0𝛽_𝑖𝜉^𝑖. (2.6) 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 𝜉_𝑖,𝑖=1(1)𝑘.

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 (0,𝜒²₀), if, for all 𝜒²∈(0,𝜒²₀),𝜒²=𝜆²ℎ²,𝜉_𝑖,𝑖=1(1)𝑘 are complex and satisfy: 𝜉₁=𝜉₂=1,𝜉_𝑗≤1𝑗=3(1)𝑘. (2.7)

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

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:^𝑘_𝑗=−𝑘𝛼_𝑗𝑦_𝑛+𝑗=^𝑚_𝑖=1ℎ^2𝑖2𝛽_𝑖0𝑦_𝑛^(2𝑖)+𝛽_𝑖1𝑦^(2𝑖)_𝑛+𝑠+𝑦^(2𝑖)_𝑛−𝑠. (2.8) With the help of the associated formula,𝑦_𝑛±𝑠=^𝑘_𝑗=−𝑘𝑎_𝑗𝑦_𝑛±𝑗+ℎ^2𝑘_𝑗=−𝑘𝑏_𝑗𝑦^′′_𝑛∓𝑗+ℎ^4𝑘_𝑗=−𝑘𝑏_𝑗𝑦⁽⁴⁾_𝑛±𝑗. (2.9) 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: 𝑦_𝑛+1−2𝑦_𝑛+𝑦_𝑛−1=^𝑚_𝑖=1ℎ^2𝑖2𝛽_𝑖0𝑦_𝑛^(2𝑖)+𝛽_𝑖1𝑦^(2𝑖)_𝑛+𝑠+𝑦^(2𝑖)_𝑛−𝑠, (2.10) where 𝑠 is a constant; the basic idea behind our approach is to approximate 𝑦_𝑛±𝑠 by the expression involving the quantities 𝑦_𝑛,𝑦_𝑛+1, and 𝑦_𝑛−1. We introduce 𝑃-stable methods of 𝑂(ℎ⁴), 𝑂(ℎ⁶), and 𝑂(ℎ⁸) with minimal phase-lag errors using (2.10) with the associated formula (2.9). This gives 𝑦_𝑛±𝑠=𝑎₀𝑦_𝑛+𝑎₁𝑦_𝑛±1+𝑎₋₁𝑦_𝑛∓1+ℎ²𝑏₀𝑦_𝑛^′′+𝑏₁𝑦^′′_𝑛±1+𝑏₋₁𝑦^′′_𝑛∓1. (2.11) Applying the method given by (2.10) and (2.11) with 𝑖=1 to the scalar test problem (2.1), we obtain the following difference equation: 𝐴𝜒𝑦_𝑛+1−2𝐵𝜒𝑦_𝑛+𝐴𝜒𝑦_𝑛−1=0,𝜒=𝜆ℎ, (2.12) where ℎ is the step length, and 𝐴(𝜒) and 𝐵(𝜒) are polynomials in 𝜒.

The characteristic equation associated with (2.12) is obtained:𝜉²−2𝑅𝜒𝜉+1=0,where𝑅𝜒=𝐴𝜒𝐵𝜒. (2.13) The roots of (2.13) are complex and of modulus one, see [7, 19], if𝑅𝜒<1. (2.14)

Under this condition the roots of (2.13) can be written as 𝜉_1,2=𝑒^{±𝑖𝜃(𝜒)} 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 cos(𝜒²), [20–22]. 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 |𝑅(𝜒²)|<1 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 |𝜉₁|≤1 and |𝜉₂|≤1 for all values of 𝜒, 𝜒=𝜆ℎ.

Definition 2.4 (see [6]). The numerical method (2.4) has an interval of periodicity (0,𝜒²) if, for all 𝜒²∈(0,𝜒²), 𝜉₁ and 𝜉₂ satisfy 𝜉₁=𝑒^{𝑖𝜃(𝜒)},𝜉₂=𝑒^{−𝑖𝜃(𝜒)}, (2.15) 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: 𝑡=𝜒−𝜃𝜒=𝜒−cos⁻¹𝑅𝜒². (2.16) If 𝑡=𝑂(𝜒^𝑞+1) as 𝜒→0, the order of phase-lag is 𝑞.

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

Theorem 2.5 (see [19]). A method which has the characteristic equation (2.13) has an interval of periodicity (0,𝜒²₀) if, for all 𝜒²∈(0,𝜒²₀), 𝑅𝜒²<1. (2.17)

Definition 2.6. A region of stability is a region of the plane, throughout which |𝑅(𝜒)|<1. Any closed curve defined by |𝑅(𝜒)|=1 is a stability boundary.

When 𝑅_𝑛𝑚(𝜈²;𝜃)<1, the roots of (2.13) are distinct and lie on the unit circle.

When 𝑅_𝑛𝑚(𝜈²;𝜃)>1, 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 (0,𝜒²₀), one gets cos𝜃𝜒=𝑅𝜒²,where𝜒²∈0,𝜒²₀. (2.18)

Remark 2.8 (see [27]). If the phase-lag order is 𝑞=2𝑠, then 𝑇=𝑐𝜒^2𝑠+1+𝑂𝜒^2𝑠+3⟹cos𝜒−𝑐𝜒=cos𝜒−cos𝜒−𝑡=𝑐𝜒^2𝑠+2+𝑂𝜒^2𝑠+3. (2.19)

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

Definition 2.10. A family of phase-fitted methods with the stability function 𝑅(𝜉,𝜒), where 𝜒=𝜆ℎ, is 𝑃-stable if, for each value of 𝜆, the quantity |𝑅(𝜒²)|<1 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 𝐴(𝜒)>0,𝑅𝜒<1⟹𝐴𝜒−𝐵𝜒>0,𝐴𝜒+𝐵𝜒>0. (2.21)

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𝑦_𝑛+1−2𝑦_𝑛+𝑦_𝑛−1=ℎ²2𝛽₁₀𝑦_𝑛^′′+𝛽₁₁𝑦^′′_𝑛+𝑠+𝑦^′′_𝑛−𝑠. (2.22)

Therefore, we introduce 𝑃-stable methods of higher order with minimal phase-lag errors using the formula:𝑦_𝑛±𝑠=𝑎₀𝑦_𝑛+𝑎₁𝑦_𝑛±1+𝑎₋₁𝑦_𝑛∓1+ℎ²𝑏₀𝑦_𝑛^′′+𝑏₁𝑦^′′_𝑛±1+𝑏₋₁𝑦^′′_𝑛∓1, (2.23) with the coefficients: 𝛽₁₀=12−𝛽₁₁,𝑎₁=12𝑠+1−𝑎₀,𝑎₋₁=121−𝑎₀−𝑠, (2.24) where 𝛽₁₁,𝑎₀,𝑏₀,𝑏₁, and 𝑏₋₁ are free parameters.

Also, the local truncation error of the method is given by1𝐸=12−𝛽₁₁𝑎₁+𝑎₋₁+2𝑏₀+𝑏₁+𝑏₋₁ℎ⁴𝑦_𝑛⁽⁴⁾𝑥_𝑛+𝑂ℎ⁶,LTE=112−𝛽₁₁𝑎₁+𝑎₋₁+2𝑏₀+𝑏₁+𝑏₋₁. (2.25) We denote this family of methods 𝑀₂(𝑎₀,𝑏₀,𝑏₁,𝑏₋₁). For an 𝑀₂(𝑎₀,𝑏₀,𝑏₁,𝑏₋₁) applied to the test equation (2.1), setting 𝜒=𝜆ℎ, we obtain (2.13), and we get 𝐴𝜒=1+𝜒²𝛽₁₁𝑎₁+𝑎₋₁−𝜒⁴𝛽₁₁𝑏₁+𝑏₋₁,𝐵𝜒=1−𝜒²𝛽₁₀+𝛽₁₁𝑎₀+𝜒⁴𝛽₁₁𝑏₀. (2.26) Then we get 𝐴𝜒+𝐵𝜒=2+𝜒²𝛽₁₁𝑎₀+𝑎₁−𝑎₋₁−𝛽₁₀+𝜒⁴𝛽₁₁𝑏₀−𝑏₁−𝑏₋₁,𝐴𝜒−𝐵𝜒=𝜒²𝛽₁₁𝑎₀+𝑎₁+𝑎₋₁+𝛽₁₀+𝜒⁴𝛽₁₁𝑏₀+𝑏₁+𝑏₋₁. (2.27) As given in (2.14), the method represented by (2.22), (2.23) with (2.24) will be 𝑃-stable provided 𝐴(𝜒)±𝐵(𝜒)>0 for all 𝜒∈(0,∞); consequently, it is easy to prove that 𝐴𝜒±𝐵𝜒>0if𝛽₁₁𝑏₀+𝑏₁+𝑏₂𝛽≤0,≤0,₁₁𝑏₀−𝑏₁−𝑏₂𝛽≥0,₁₁1−𝑎₀≥14. (2.28) Following the phase-lag, which is denoted by 𝑃(𝜒), as a leading coefficient in the expression of |(𝜃(𝜒)−𝜒)/𝜒|, then𝑃𝜒=1𝜒²𝐴𝜒cos𝜒−𝐵𝜒=𝜒²−𝛽₁₁𝑏₀+𝑏₁+𝑏₋₁−12𝛽₁₁𝑎₁+𝑎₋₁+124+𝜒⁴12𝛽₁₁𝑏₁+𝑏₋₁+1𝛽24¹¹𝑎₁+𝑎₋₁+16!+𝜒⁶−1𝛽24¹¹𝑏₁+𝑏₋₁−1𝛽6!¹¹𝑎₁+𝑎₋₁+18!+𝜒⁸1𝛽6!¹¹𝑏₁+𝑏₋₁+1𝛽8!¹¹𝑎₁+𝑎₋₁+1.10! (2.29) The method has phase-lag error of 𝑂ℎ⁴if𝛽₁₁𝑏₀+𝑏₁+𝑏₋₁𝛽<0,₁₁1−𝑎₀>14, (2.30) 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). 𝑏₀=𝑏₁+𝑏₋₁,𝑏₀>1/24,𝛽₁₁<0,𝑎₀>1
Let 𝑎₀=11/5,𝑏₀=1/4, and 𝛽₁₁=−5/12

Case 2 (2.1). 𝑏₀=0,𝛽₁₁≠0,𝑎₀=1−(1/12𝛽₁₁)+2(𝑏₁+𝑏₋₁), and 𝑏₀=0,and𝛽₁₁(𝑏₁+𝑏₋₁)=−1/12

Case 3 (2.2). 𝑏₀=−(𝑏₁+𝑏₋₁),𝛽₁₁≠0,𝑎₀=1−1/12𝛽₁₁
Let 𝑎₀=11/5, and 𝑏₀=1/4

Case 4 (2.3). 𝑏₀=0,𝑎₀=1−(1−24𝛽₁₁(𝑏₁+𝑏₋₁))/12𝛽₁₁ for 𝛽₁₁(𝑏₁+𝑏₋₁)≤−1/12

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:𝛽₂₁=1𝑠²,𝛽₁₁=6𝑠²−112𝑠²,𝑎₁=12𝑠+1−𝑎₀,𝑎₋₁=121−𝑎₀−𝑠,𝑏₁=13𝑎12⁰−2𝑏₀+(𝑏𝑠−1)(𝑠+1)(𝑠+3)₋₁=13𝑎12⁰−2𝑏₀−(𝑠−1)(𝑠+1)(𝑠−3) (2.31) where 𝛽₁₁, 𝑠,𝑎₀, and 𝑏₀ are free parameters.

Also, the local truncation error of the method is given by 1𝐸=−1360𝛽12¹¹𝑎₁+𝑎₋₁+12𝑏₁+𝑏₋₁ℎ⁶𝑦_𝑛⁽⁶⁾𝑥_𝑛+𝑂ℎ⁸. (2.32) For an 𝑀₄(𝑎₀,𝑏₀,𝑏₁,𝑏₋₁), setting 𝜒=𝜆ℎ, (2.13) is obtained with the coefficient given by, and we get𝐴𝜒±𝐵𝜒≥0⟹𝑃-stablemethod⟹1−𝑎₀−𝑠²≥0,4𝑏₀−𝑎₀+1−𝑠²𝑎≥0,₀−1+𝑠²𝑃𝜒=1≤0,𝜒²𝐴𝜒cos𝜒−𝐵𝜒=𝜒²−𝛽₁₁𝑏₀+𝑏₁+𝑏₋₁−12𝛽₁₁𝑎₁+𝑎₋₁+124+𝜒⁴12𝛽₁₁𝑏₁+𝑏₋₁+1𝛽24¹¹𝑎₁+𝑎₋₁−1.6! (2.33)

The coefficient of 𝜒²⇒(1/2)𝛽₁₁𝑠²=1/24. the coefficient of 𝜒⁴⇒𝛽₁₁(5𝑎₀−5+12𝑏₀+6𝑠²)=1/30,

We can obtain a symmetric 𝑃-stable second derivative, two-step method shown by (2.22) if 𝑎₀<1−3𝑠² and we have the following.(1)𝑏₀=0, with phase-lag error of order: 𝑂ℎ⁴if𝛽₁₁=112𝑠²,𝑂ℎ⁶if𝛽₁₁≥1270𝑠². (2.34)(2)𝑏₀≠0,𝑏₀≥(𝑎₀−1+𝑠²)/4, with phase-lag error of order: 𝑂ℎ⁴if𝛽₁₁=112𝑠²,𝑂ℎ⁶if𝛽₁₁4𝑏₀−3𝑠²≤1.90 (2.35)

2.2.3. Sixth-Order Scheme

The two-step 𝑃-stable methods of order six with minimal phase-lag errors using the formula𝑦_𝑛+1−2𝑦_𝑛+𝑦_𝑛−1=ℎ²2𝛽₁₀𝑦_𝑛^′′+𝛽₁₁𝑦^′′_𝑛+𝑠+𝑦^′′_𝑛−𝑠+ℎ⁴2𝛽₂₀𝑦_𝑛⁽⁴⁾+𝛽₂₁𝑦⁽⁴⁾_𝑛+𝑠+𝑦⁽⁴⁾_𝑛−𝑠+ℎ⁶2𝛽₃₀𝑦_𝑛⁽⁶⁾+𝛽₃₁𝑦⁽⁶⁾_𝑛+𝑠+𝑦⁽⁶⁾_𝑛−𝑠, (2.36) associated with the formula 𝑦_𝑛±𝑠=𝑎₀𝑦_𝑛+𝑎₁𝑦_𝑛±1+𝑎₋₁𝑦_𝑛∓1+ℎ²𝑏₀𝑦_𝑛^′′+𝑏₁𝑦^′′_𝑛±1+𝑏₋₁𝑦^′′_𝑛∓1. (2.37) Similarly, we can obtain the following results: 𝛽₁₀=12−𝛽₁₁,𝛽₂₁=𝛽₃₁=𝛽₁₁,𝛽₂₀=124−𝛽₁₁12𝑠²+1,𝑎₀=125𝑏₀−32,𝑎45₁=125𝑏₀−77−45Rootof2025𝑧²+4860𝑏₀−3465𝑧+2916𝑏²₀+301−4158𝑏₀,𝑏₋₁=11𝑏30⁰−29Rootof2025𝑧²+4860𝑏₀−3465𝑧+2916𝑏²₀+301−4158𝑏₀−14,𝑎405₋₁=Rootof2025𝑧²+4860𝑏₀−3465𝑧+2916𝑏²₀+301−4158𝑏₀,𝑏₁1=−6𝑏₀−28−2819Rootof2025𝑧²+4860𝑏₀−3465𝑧+2916𝑏²₀+301−4158𝑏₀. (2.38) The 𝑃-stability conditions will be 𝐴𝜒=1+𝜒²𝛽₁₁1−125𝑏₀−𝜒⁴+𝜒⁶15𝛽₂₁𝑏₀+𝛽₃₁1−125𝑏₀−15𝜒⁸𝛽₃₁𝑏₀,𝐵𝜒=1−𝜒²𝛽₁₀+125𝛽₁₁𝑏₀+𝜒⁴𝛽₂₀−125𝛽₂₁𝑏₀+𝛽₁₁𝑏₀+𝜒⁶𝛽₂₁𝑏₀. (2.39)

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:𝑦_𝑛+2−2𝑦_𝑛+1+2𝑦_𝑛−2𝑦_𝑛−1+𝑦_𝑛+2=ℎ²2𝛽₁₀𝑦_𝑛^′′+𝛽₁₁𝑦^′′_𝑛+𝑠+𝑦^′′_𝑛−𝑠+ℎ⁴2𝛽₂₀𝑦_𝑛⁽⁴⁾+𝛽₂₁𝑦⁽⁴⁾_𝑛+𝑠+𝑦⁽⁴⁾_𝑛−𝑠. (3.1) Associated with𝑦_𝑛±𝑠=𝑎₀𝑦_𝑛+𝑎₁𝑦_𝑛±1+𝑎₋₁𝑦_𝑛∓1+𝑎₂𝑦_𝑛±2+𝑎₋₂𝑦_𝑛∓2+ℎ²𝑏₀𝑦_𝑛^′′+𝑏₁𝑦^′′_𝑛±1+𝑏₋₁𝑦^′′_𝑛∓1+𝑏₂𝑦^′′_𝑛±2+𝑏₋₂𝑦^′′_𝑛∓2, (3.2) the characteristic equation will be 𝐴𝜒𝜉⁴+𝐵𝜒𝜉³+𝐶𝜒𝜉²+𝐵𝜒𝜉+𝐴𝜒=0, (3.3) where 𝐴𝜒=1+𝑎₂+𝑎₋₂−𝜒²𝑏₂+𝑏₋₂𝜒²𝛽₁₁−𝜒²𝛽₂₁,𝐵𝜒=−2+𝑎₁+𝑎₋₁−𝜒²𝑏₁+𝑏₋₁𝜒²𝛽₁₁−𝜒²𝛽₂₁,𝐶𝜒=21+𝜒²𝛽₁₀−𝜒⁴𝛽₂₀+𝜒²𝑎₀−𝑏₀𝜒²𝛽₁₁−𝜒²𝛽₂₁. (3.4) The 𝑃-stability conditions are𝑃₁𝜒𝜒𝜒𝜒𝑃=2𝐴−2𝐵+𝐶≥0,₂𝜒=12𝐴𝜒−2𝐶𝜒≥0,𝑃²₂𝜒−4𝑃₁𝜒𝑃₃𝜒𝑃≥0,₃𝜒=2𝐴𝜒+2𝐵𝜒+𝐶𝜒≥0. (3.5)

Following the phase-lag, denoted by 𝑃(𝜒), is leading coefficient in the expression of 𝑃𝜒=2𝐴𝜒cos2𝜒+2𝐵𝜒cos𝜒+𝐶𝜒𝜒𝜒𝜒𝜒4𝐴sin2𝜒+2𝐵sin=𝑂⁴. (3.6) The local truncation error of the method is given byLET=31180−𝛽₁₁1𝑎12¹+𝑎₋₁+43𝑎₂+𝑎₋₂+𝑏₁+𝑏₋₁+13𝑏₂+𝑏₋₂−2𝛽₂₀−𝛽₂₁𝑎₁+𝑎₋₁+4𝑎₂+𝑎₋₂+𝑏₀+𝑏₁+𝑏₋₁+𝑏₂+𝑏₋₂ℎ⁶𝑦_𝑛⁽⁶⁾𝑥_𝑛+𝑂ℎ⁸. (3.7)

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 𝑞=4, the coefficients of (3.1) and (3.2) will be𝛽₂₁=7−1122𝑠²𝛽₁₁−𝛽₂₀,𝛽₁₀=1−𝛽₁₁,𝑎₀=1,𝑏_𝑖𝑎=𝑖,𝑖=1,2,₁1=−6𝑠³−𝑠²+4𝑠−48,𝑎₋₁=16𝑠³−𝑠²,𝑎−4𝑠−48₂=1𝑠12³+𝑠²−𝑠−48,𝑎₋₂1=−𝑠12³−2𝑠²−𝑠−48. (3.8) Also, the local truncation error of the method is given byLET=2⁷−46!−6!2𝛽₁₁𝑎₁+𝑎₋₁+16𝑎₂+16𝑎₋₂−𝑏4!¹+𝑏₋₁+4𝑏₂+4𝑏₋₂.2! (3.9) The coefficients of the characteristic equation will be 𝐴𝜒𝑠=1+²3𝜒²𝛽₁₁−𝜒²𝛽₂₁,𝐵𝜒=−2−4𝑠3𝜒²𝛽₁₁−𝜒²𝛽₂₁,𝐶𝜒=21+𝜒²−𝜒⁴𝛽₂₀−𝜒²𝛽₂₁. (3.10) The 𝑃-stability conditions are𝑃₁𝜒=8+𝜒²2𝑠²3𝛽₁₁+83𝑠𝛽₁₁+2+𝜒⁴−2𝑠²3𝛽₂₁−83𝑠𝛽₂₁−2𝛽₂₀−2𝛽₂₁2=8+3𝜒²𝛽₁₁𝑠(𝑠+4)+3+𝜒⁴−2𝛽₂₀−23𝛽₂₁3+4𝑠+𝑠²𝑃≥0,₂𝜒=8+−4+143𝑠²+4𝛽₁₁−1+𝑠²𝜒²𝑃≥0,₃𝜒=13𝑠²−7+61+𝛽₁₁𝑠²𝜒²𝑃≥0,₄𝜒=64+−1281+𝛽₁₁+4483𝑠²+64𝛽₁₁𝑠²𝜒²−649𝛽₁₁−13+6𝑎₀−6𝛽₁₁+19𝑠²448𝛽₁₁−1+𝑎₀−128𝛽₁₁𝑠²31+𝛽₁₁+16𝛽²₁₁𝑠⁴𝜒⁴≥0. (3.11) We can obtain a symmetric 𝑃-stable second derivative, four-step methods with minimal phase-lag errors by suitable choices of the free parameters, see [8–10].

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: 𝐴𝜒𝜉⁴+𝐵𝜒𝜉³+𝐶𝜒𝜉²+𝐵𝜒𝜉+𝐴𝜒=0, (3.12) with the following coefficients: 𝛽₁₀=1−𝛽₁₁,𝛽₂₀1=−180𝑠²31−105𝑠²+75𝑠⁴𝛽₁₁,𝛽₂₁=1180𝑠²31−15𝑠⁴𝛽₁₁,𝑎₁=𝑠870528−102𝑠²−29𝑠³+9𝑠⁴,𝑎₋₁=𝑠870−528+102𝑠²−29𝑠³−9𝑠⁴,𝑎₂=𝑠870−93+102𝑠²+29𝑠³−9𝑠⁴,𝑎₋₂=𝑠87093−102𝑠²+29𝑠³+9𝑠⁴,𝑏₁=𝑠290024−31𝑠²+7𝑠⁴,𝑏₀=𝑠²105−𝑠²,𝐴𝜒=1+𝑎₂+𝑎₋₂−𝜒²𝑏₂+𝑏₋₂𝜒²𝛽₁₁−𝜒²𝛽₂₁,𝐵𝜒=−2+𝑎₁+𝑎₋₁−𝜒²𝑏₁+𝑏₋₁𝜒²𝛽₁₁−𝜒²𝛽₂₁,𝐶𝜒=21+𝜒²𝛽₁₀−𝜒⁴𝛽₂₀+𝜒²𝑎₀−𝑏₀𝜒²𝛽₁₁−𝜒²𝛽₂₁. (3.13)(3.14)

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

Let 𝛽₁₀=1−𝛽₁₁, 𝛽₂₁=(7/12)−(1/2𝑠²)𝛽₁₁−𝛽₂₀, let 𝑎₁=(2/3)(16+𝑠−𝑎₀)−(𝑠²/6)(𝑠+1), let 𝑎₂=(1/12)(𝑠(𝑠²+2𝑠−1)+2𝑎₀−62), let 𝑎₋₁=(1/6)(𝑠(𝑠²−𝑠−4)−(4𝑎₀+56)),let 𝑎₋₂=(1/12)(𝑠(1+2𝑠−𝑠²)+2𝑎₀+58), and let 𝑏_𝑖=𝑖forall𝑖 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.