Trigonometrically-fitted second derivative method for oscillatory problems

Abstract A continuous Trigonometrically-fitted Second Derivative Method (CTSDM) whose coefficients depend on the frequency and stepsize is constructed using trigonometric basis functions. A discrete Trigonometrically-fitted second derivative method (TSDM) is recovered from the CTSDM as a by-product and applied to solve initial value problems (IVPs) with oscillating solutions. We discuss the stability properties of the TSDM and present numerical experiments to demonstrate the efficiency of the method. AMS Subject Classification 65L05; 65L06


Introduction
In this paper, we consider the subclass of first order differential equation with periodic or oscillating solutions where f : × m → m , y, y 0 ∈ m . Oscillatory IVPs frequently arise in areas such as classical mechanics, celestial mechanics, quantum mechanics, and biological sciences. Several numerical methods based on the use of polynomial basis functions have been developed for solving this class of important problems (see Lambert , Hairer et al. in (Hairer and Wanner 1996), Hairer 1982, andSommeijer 1993). Other methods based on exponential fitting techniques which take advantage of the special properties of the solution that may be known in advance have been proposed (see Simos 1998, Vanden Berghe et al. 2001a, Vanden Berghe et al. 2009, Vigo-Aguiar et al. 2003, Franco 2002, Fang et al. 2009, Nguyen et al. 2007, Ozawa 2005, Jator et al. 2012, and Ngwane et al. 2012b). In the spirit of 2005, the motivation governing the exponentially-fitted methods is inherent to the fact that if the frequency or a reasonable estimate of it is known in advance, these methods will be more advantageous than the polynomial based methods.
*Correspondence: fifonge@yahoo.com 1 Department of Mathematics, USC Salkehatchie, Allendale, SC 29810, USA Full list of author information is available at the end of the article The aim of this paper is to construct a TSDM. This construction is done by initially developing a CTSDM which then provides a discrete method that is applied as a TSDM which takes the frequency of the solution as a priori knowledge. In particular, CTSDM consists of a sum of continuous functions while TSDM is a by-product of CTSDM. The coefficients of the TSDM are functions of the frequency and the stepsize, hence the solutions provided by the proposed method are highly accurate if (1) has periodic solutions with known frequencies. We adopt the approach given in Jator et al. in (Ngwane and Jator 2012a;Jator et al. 2012), where the TSDM is used to obtain the approximation y n+1 to the exact solution y(x n+1 ) on the interval [x n , x n+1 ], h = x n+1 − x n , n = 0, . . . , N − 1, on a partition [a, b], where a, b ∈ R, h is the constant stepsize, n is a grid index and N > 0 is the number of steps. We note that second derivative methods with polynomial basis functions were proposed to overcome the Dahlquist 1956 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 Enright 1974).
This paper is organized as follows. In Section "Development of method", we obtain a trigonometric basis representation U(x) for the exact solution y(x) which is used to generate a TSDM for solving (1). The analysis and implementation of the TSDM are discussed in Section "Error analysis and stability". Numerical examples are given in http://www.springerplus.com/content/3/1/304 Section "Numerical examples" to show the accuracy and efficiency of the TSDM. Finally, we give some concluding remarks in Section "Conclusion".

Development of method
In this section, our objective is to construct a CTSDM which produces a discrete method as a by-product. The method has the form where u = wh, β j (u), γ j (u), j = 0, 1, are coefficients that depend on the stepsize and frequency. We note that y n+j is the numerical approximation to the analytical solution y(x n+j ), and with j = 0, 1. In order to obtain equation (2) we proceed by seeking to approximate the exact solution y(x) on the interval [x n , x n + h] by the interpolating function U(x) of the form where a 0 , a 1 , a 2 , a 3 and a 4 are coefficients that must be uniquely determined. We then impose that the interpolating function in (3) coincides with the analytical solution at the point x n to obtain the equation We also demand that the function (3) satisfies the differential equation (1) at the points x n+j , j = 0, 1 to obtain the following set of three equations: Equations (4) and (5) lead to a system of five equations which is solved by Cramer's rule to obtain a j , j = 0, 1, 2, 3, 4. Our continuous CTSDM is constructed by substituting the values of a j into equation (3). After some algebraic manipulation, the CTSDM is expressed in the form where w is the frequency, β 0 (w, x), β 1 (w, x), γ 0 (w, x), and γ 1 (w, x), are continuous coefficients. The continuous method (6) is used to generate the method of the form (2).
Thus, evaluating (6) at x = x n+1 and letting u = wh, we obtain the coefficients of (2) as follows: Error analysis and stability

Local truncation error
We note that when u → 0 the coefficients given by (7) are vulnerable to heavy cancellations and hence the following Taylor series expansion must be used (see Simos 1998).
In fact, for practical computations when u is small, it is better to use the series expansion (8) (see Calvo et al. 2009).
Thus the Local Truncation Error (LTE) of (2) subject to (8) is obtained as Remark 1. The method (2) specified by (8) is a fourthorder method and reduces to a one-step conventional second derivative method as u → 0 (see Lambert 1973, p. 201).

Stability
Proposition 1. The TSDM (2) applied to the test equations y = λy and y = λ 2 y yields Proof. We begin by applying (2) to the test equations y = λy and y = λ 2 y which are expressed as f (x, y) = λy and g(x, y) = λ 2 y respectively; letting q = hλ and u = wh, we obtain a linear equation which is used to solve for y n+1 with (11) as a consequence. http://www.springerplus.com/content/3/1/304 Remark 2. The rational function M(q; u) is called the stability function which determines the stability of the method.

Definition 1. A region of stability is a region in the q
The TSDM method (2) specified by (7) is given by (12) Definition 2. The method (12) is zero-stable provided the roots of the first characteristic polynomial have modulus less than or equal to unity and those of modulus unity are simple (see Lambert 1991). (12) is consistent if it has order p > 1 (see (Fatunla 1991)). (12) is consistent as it has order p > 1 and zero-stable, hence it is convergent since zerostability + consistency = convergence.

Remark 3. The TSDM
Remark 4. In the q − u plane the TSDM (12) is stable for q ≤ 0, and u ∈[−2π, 2π], since from above |M(q; u)| ≤ 1, q ≤ 0. Figure 1 is a plot of the stability region and Figure 2 shows the zeros and poles of M(q; u). We note from   Nguyen et al. 2007).

Definition 4. The TSDM with the stability function (11) is said to be
Remark 6. We observe from definition 1, remarks 4, 5, and Figure 2, that TSDM is A-stable. In particular, |M(q; iy)| = 1 ∀y ∈ R, and by the maximum principle,  Hairer et al. 1996, p.43, 53). Moreover, the real part of the zeros of |M(q; u)| must be negative, while the real part of the poles of |M(q; u)| must be positive.

Implementation
In the spirit of Ngwane et al. in (Ngwane and Jator 2012a;2012b), the TSDM (12) is implemented to solve (1) without requiring starting values and predictors. For instance, if we let n = 0 in (12), then y 1 is obtained on the subinterval [x 0 , x 1 ], as y 0 is known from the IVP. Similarly, if n = 1, then y 2 is obtained on the sub-interval [x 1 , x 2 ], as y 1 is known from the previous computation, and so on; until we reach the final sub-interval [x N−1 , x N ]. We note

Numerical examples
In this section, we give numerical examples to illustrate the accuracy (small errors) and efficiency (fewer number of function evaluations (NFEs)) of the TSDM. We find the approximate solution on the partition π N , where π N : a = x 0 < x 1 < x 2 < ... < x n < x n+1 < . . . < x N = b, and we give the errors at the endpoints calculated as Error=y N − y(x N ). We note that the method requires only two function evaluations per step and in general requires (2N + 2) NFEs on the entire interval. All computations were carried out using a written code in Matlab. Ozawa 2005. y 1 = − y 1 r 3 , y 2 = − y 2 r 3 , r = y 2 1 + y 2 2 , y 1 (0) = 1 − e, y 1 (0) = 0, y 2 (0) = 0, y 2 (0)

Example 1. Consider the given two-body problem which was solved by
where e, 0 ≤ e < 1 is an eccentricity. The exact solution of this problem is where k is the solution of the Kepler's equation k = x + e sin(k). We choose ω = 1.

The analytic solution is given by
We choose ω = 1.01 and for more on frequency choice see Ramos et al. 2010.
We compare the end-point global errors for TSDM with the fourth order methods in Ixaru et al. 2004. We see from Table 2 that the results produced by TSDM are better than Simos' method used in (Ixaru and Berghe 2004), as TSDM produces better error magnitude while using less number of steps and fewer number of function evaluations. TSDM is very competitive to the method used by Ixaru et al. 2004.
Example 3. We consider the following inhomogeneous IVP by Simos 1998.
where the analytic solution is given by The exponentially-fitted method in Simos 1998 is fourth order and hence comparable to our method, TSDM. We  Example 4. Linear Kramarz problem We consider the following second-order IVP, (see Nguyen et al. 2007[p. 204]) Exact : y(t) = (2 cos(t), − cos(t)) T .
We use this example to show the efficiency of TSDM on linear systems. Nguyen et al. 2007 used the "trigonometric implicit Runge-Kutta", TIRK3, method to solve the above linear Kramarz problem. Clearly, TSDM performs better as seen in Table 4.
Example 5. We consider the IVP (see Vigo-Aguiar et al. 2003) y + K 2 y = K 2 x, y(0) = 10 −5 , y (0) where K = 314.16, and we choose ω = 314.16. The analytic solution is given by Exact : y(x) = x + 10 −5 (cos(Kx) − cot(K) sin(Kx)).    This problem demonstrates the performance of TSDM on a well-known oscillatory problem. We compare the results from TSDM with the Dissipative Chebyshev exponential-fitted methods, CHEBY24 and CHEBY1 used in Vigo-Aguiar et al. 2003. We see that TSDM uses fewer number of function evaluations with better accuracy than CHEBY24 that is designed to use fewer number of steps. Integrating in the interval [0, 1] with a stepsize equal to the total length of the interval, we obtain an error of order 10 −21 . Hence TSDM is a more efficient integrator. We note that compared with the methods CHEBY24 and CHEBY1 which use stepsizes considerably larger than those used in multistep methods, TSDM is very competitive and superior to both CHEBY24 and CHEBY1.
We choose β = −3 and β = −1000 in order to illustrate the phenomenon of stiffness. Given the initial conditions y 1 (0) = 2 and y 2 (0) = 3, the exact solution is β-independent and is given by This example is chosen to demonstrate the performance of TSDM on stiff problems. We compute the solutions to        Example (6) with β = −3, −1000. We obtain better absolute errors than Nguyen et al. (2007). This efficiency is achieved using fewer number of steps and less number of function evaluations than Nguyen et al. (2007). For example when β = −3, our method generates a solution with error magnitude 10 −6 involving just 6 steps and 28 function evaluations, whereas (Nguyen et al. 2007) attains the same error magnitude using 10 steps and 47 function evaluations. When β = −1000, TSDM generates solutions with comparable error magnitude. We see that TSDM is competitive and better than the method in Nguyen et al. (2007) which is of order six and is thus expected to do better.

An implementation in predictor-corrector mode
In this section, we also implement our CTSDM in a predictor-corrector mode. The predictor is given by where and the corrector is given by equations (6) and (7). We note that when u → 0 we use the following Taylor series expansion (see Simos 1998) ⎧ ⎪ ⎨ ⎪ ⎩ α 0 = 1 − u 2 6 + u 4 120 − u 6 5040 + u 8 362880 − u 10 39916800 + . . .
As we expected, the predictor-corrector (PreCor) mode runs faster than the TSDM but is less accurate compared to the TSDM. We illustrate this by applying the predictorcorrector to Example 2 and Example 3. We plot the efficiency curves showing the accuracy versus the CPU computation time, and the accuracy versus the NFEs.

Estimating the frequency
A preliminary testing indicates that a good estimate of the frequency can be obtained by demanding that LTE = 0, and solving for the frequency. That is, solve for ω given that − h 5 720 w 2 y (3) (x n ) + y (5) (x n ) = 0, where y (j) , j = 2, . . . , 5 denote derivatives. We used this procedure to calculate ω for the problem given in example (5) and obtained ω ≈ ± 314.16, which approximately gives the known frequency ω = 314.16. Hence, this procedure is interesting and will be seriously considered in our future research.
We note that estimating the frequency and the choice of the frequency in trigonometrically-fitted methods is challenging and has grown in interest. Existing references on how to estimate the frequency and on the choice of the frequency include G. Vanden Berghe et al. 2001b, andRamos et al. 2010.

Conclusion
We have proposed a TSDM for solving oscillatory IVPs. The TSDM is A-stable and hence, an excellent candidate for solving stiff IVPs. This method has the advantages of being self-starting, having good accuracy with order 4, and requiring only two functions evaluation at each integration step. We have presented representative numerical examples that are linear, non-linear, stiff and highly oscillatory. These examples show that the TSDM is more accurate and efficient than those in Nguyen et al. 2007, Simos 1998, Ixaru et al. 2004, and Ozawa 2005. Details of the numerical results are displayed in Tables 1, 2 , 3, 4, 5, 6, 7, 8, 9 and 10 and the efficiency curves are presented in Figures 3,4,5,6,7,8,9,10,11 and 12. Our future research will incorporate a technique for accurately estimating the frequency as suggested in subsection "Estimating the frequency" as well as implementing the method in a variable step mode. http://www.springerplus.com/content/3/1/304