A Family of Trigonometrically Fitted Enright Second Derivative Methods for Stiff and Oscillatory Initial Value Problems

A family of Enright’s second derivative formulas with trigonometric basis functions is derived using multistep collocation method. The continuous schemes obtained are used to generate complementary methods. The stability properties of the methods are discussed.Themethodswhichcanbeappliedinpredictor-correctorformareimplementedinblockformassimultaneousnumericalintegratorsovernonoverlappingintervals.Numericalresultsobtainedusingtheproposedblockformrevealthatthenewmethodsareefficientandhighlycompetitivewithexistingmethodsintheliterature.


Introduction
Many real life processes in areas such as chemical kinetics, biological sciences, circuit theory, economics, and reactions in physical systems can be transformed into systems of ordinary differential equations (ODE) which are generally formulated as initial value problems (IVPs). Some classes of IVPs are stiff and/or highly oscillatory as described by the following model problem: where ( ) ∈ R and is × real matrix with at least one eigenvalue with a very negative real part and/or very large imaginary part, respectively (see Fatunla [1]). Many conventional methods cannot solve these types of problems effectively.
Stiff systems have been solved by several authors including Lambert [2,3], Gear [4,5], Hairer [6], and Hairer and Wanner [7]. Different methods including the Backward Differentiation Formula (BDF) have been used to solve stiff systems. Second derivative methods with polynomial basis functions were proposed to overcome the Dahlquist [8] barrier theorem whereby the conventional linear multistep method was modified by incorporating the second derivative term in the derivation process in order to increase the order of the method, while preserving good stability properties (see Gear [9], Gragg and Stetter [10], and Butcher [11]).
Obrechkoff [20] proposed a general multiderivative method for solving systems of ordinary differential equations. Special cases of Obrechkoff method have been developed by many others including Cash [21] and Enright [22]. The methods by Enright [22] have order = + 2 for a step method.
In this paper, we propose a numerical integration formula which more effectively copes with stiff and/or oscillatory 2 Journal of Applied Mathematics IVPs. We will construct a continuous form of the second derivative multistep method (CSDMM) using a multistep collocation technique such that Enright's second derivative methods (ESDM) will be recovered from the derived continuous methods. The aim of this paper is to derive a family of Enright's second derivative formulas with trigonometric basis functions using multistep collocation method. Many methods for solving IVPs are implemented in a step-by-step fashion in which, on the partition , an approximation is obtained at +1 only after an approximation at has been computed, where : = 0 < 1 < ⋅ ⋅ ⋅ < < +1 < ⋅ ⋅ ⋅ < = , +1 = + ℎ, = 1, . . . , , ℎ = ( − )/ , ℎ is the step size, is a positive integer, and is the grid index. We implement ESDM in block form.
In Section 2, we present a derivation of the family of Enright methods. Error analysis and stability are discussed in Section 3. The implementation of the ESDM and numerical examples to show the accuracy and efficiency of the ESDM are given in Section 4. Finally, we conclude in Section 5.

Derivation of the Family of Methods
We consider the first-order differential equation where is assumed to satisfy the conditions to guarantee the existence of a unique solution of the initial-value problem.

Block Specification and Implementation of the Methods.
We consider a general procedure for the block implementation of the methods in matrix form (see Fatunla [26]). First we define the following vectors: where + = ( + ℎ), + = ( + ℎ, ( + ℎ)), and + = ( ( , ( ))/ )| + + . The integration on the entire block will be compactly written as which forms a nonlinear equation because of the implicit nature, and hence we employ the Newton iteration for the evaluation of the approximate solutions. We use Newton's approach for the implementation of implicit schemes to get the following solution of the block: The × matrices 0 , 1 , 0 , 1 , and 1 are defined as follows.
Define the local truncation error of (23) as follows:   Table 1. (ii) From the local truncation error constant computation, it follows that the order ( ) of the method (23) is = + 2.

Stability
Definition 6. The block method (23) is zero-stable, provided the roots of the first characteristic polynomial have modulus less than or equal to one and those of modulus one are simple (see [2]).

Remark 9. The block method (23) is consistent as it has order
> 1 and zero-stable; hence it is convergent since zerostability + consistency = convergence.
Proof. We begin by applying (23) to the test equations = and = 2 which are expressed as ( , ) = and ( , ) = 2 , respectively; letting = ℎ and = ℎ, we obtain a linear equation which is used to solve for +1 with (23) as a consequence.
Remark 11. The rational function ( ; ) is called the stability function which determines the stability of the method.  Remark 15. Figures 1, 3, 5, and 7 are plots of the stability region of ( ; ) for case = 1, 2, . . . , 4, respectively. We note from these figures that the stability region of ( ; )

Numerical Illustration
In this section we consider some standard problems: stiff, oscillatory, linear, and nonlinear systems that appear in

Re Im
Poles and zeros for k = 3  of the ESDM once for all integration. Some of the methods of orders 4 and 6 in the literature have been compared to EM2 and EM4, respectively. We find the approximate solution on the partition , and we give the errors at the endpoints calculated as Error = − ( ). We denote the Max| − ( )| by Err, the number of steps by , and the number of function evaluations by NFEs. We will write an error of the form Err = × 10 − as (− ).
EM2 is fourth-order and hence comparable to the exponentially fitted method by Simos [25] which is also of fourthorder. PC1 and PC2 denote the predictor-corrector mode for = 1 and = 2, respectively. The efficiency curves in Figure 9 show the computational efficiency of the two methods (Simos and EM2) by considering the NFEs over integration steps for each method. Hence for this example, EM2 performs better than Simos. We see from Table 2 that ESDM is efficient for each case.
Example 17. We consider the nonlinear Duffing equation which was also solved by Ixaru and Vanden Berghe [30]:  We compare the errors produced by our EM2 with the fourth-order methods by Ixaru and Vanden Berghe [30].  We see from Table 3 that the results produced by ESDM on Example 17 are very good. In fact, EM2 produces results that are better than Simos' method used in [30], as it produces better error magnitude while using less number of steps and fewer number of function evaluations. This example once more shows us that the ESDM produces good results and in particular EM2 is very competitive to the method used by Ixaru and Vanden Berghe [30].
Example 18 (a nearly sinusoidal problem). We consider the following IVP on the range 0 ≤ ≤ 10 (see the study by Nguyen et al. [16]): We choose = −3 and = −1000 in order to illustrate the phenomenon of stiffness. Given the initial conditions 1 (0) = 2 and 2 (0) = 3, the exact solution is -independent and is given by We choose this example to demonstrate the performance of ESDM on stiff problems. We compute the solutions to Example 18 with = −3, −1000. We compare EM4 of order six to the method by Nguyen et al. [16] which is also of order six. For both = −3 and = −1000, EM4 clearly obtains better absolute errors compared to Nguyen et al. [16]. This efficiency is achieved using fewer number of steps and less number of function evaluations than Nguyen et al. [16].

A Predictor-Corrector Mode Implementation of ESDM.
We can implement our ESDM in a predictor-corrector (PC) mode. The predictor for = 1, 2, 3, 4 is given below while the corrector for each case is given by the main method (10). PC1 and PC2 denote the PC mode for = 1 and 2, respectively: = 2: = 3, 4: The coefficients , = 0, 1, . . . , − 1, and + −1 for = 3 are given below while those for = 4 are given in Appendix 3 of the supplementary material. We observe that the PC implementation performs poorly relative to the block implementation; see Table 2. (44)

4.2.
Estimating the Frequency. Though we are mainly interested in problems where is taken as the exact frequency of the analytical solution and is known in advance, it is important to note that the exact frequency may be unknown for some problems. A preliminary testing indicates that a good estimate of the frequency can be obtained by demanding that the LTE in Theorem 3 equals zero and solving for the frequency. That is, solve for given that +3 ℎ +3 ( 2 ( +1) ( ) + ( +3) ( )) = 0, where ( ) , = + 1, + 3, denotes derivatives. We used this procedure to estimate for the problem given in Example 16 and obtained ≈ ±9.999996, which approximately gives the known frequency = 10. Hence, this procedure is interesting and will be seriously considered in our future research.
If a problem has multiple frequencies, then is approximatively calculated so that it is an indicative frequency (see Nguyen et al. [16]). We note that estimating the frequency when it is unknown as well as finding the frequency for problems for which the frequency varies over time is very challenging. This challenge and the choice of the frequency in trigonometrically fitted methods have grown in interest. Existing references on how to estimate the frequency and the choice of the frequency include Vanden Berghe et al. [31] and Ramos and Vigo-Aguiar [32].

Conclusion
In this paper, we have proposed a family of Enright methods using trigonometric bases for solving stiff and oscillatory IVPs. The ESDM is zero-stable and produced good results on stiff IVPs. This method has the advantages of being selfstarting and having good accuracy properties. ESDM has order ( +2) similar to that in Enright [22]. We have presented representative numerical examples that are linear, nonlinear, stiff, and highly oscillatory. The need that the frequency be known in advance might be a shortcoming, yet these examples show that the ESDM is not only promising but more accurate and efficient than those in Nguyen et al. [16], Simos [25], Ixaru and Vanden Berghe [30], and Ozawa [15]. Details of the numerical results are displayed in Tables 2, 3, 4, and 5, and the efficiency curves are presented in Figures 9, 10, 11, 12, and 13. Our future research will incorporate a technique for accurately estimating the frequency as suggested in    Section 4.2 as well as implementing the method in a variable step mode.