Research Article  Open Access
F. F. Ngwane, S. N. Jator, "A Trigonometrically Fitted Block Method for Solving Oscillatory SecondOrder Initial Value Problems and Hamiltonian Systems", International Journal of Differential Equations, vol. 2017, Article ID 9293530, 14 pages, 2017. https://doi.org/10.1155/2017/9293530
A Trigonometrically Fitted Block Method for Solving Oscillatory SecondOrder Initial Value Problems and Hamiltonian Systems
Abstract
In this paper, we present a block hybrid trigonometrically fitted RungeKuttaNyström method (BHTRKNM), whose coefficients are functions of the frequency and the stepsize for directly solving general secondorder initial value problems (IVPs), including Hamiltonian systems such as the energy conserving equations and systems arising from the semidiscretization of partial differential equations (PDEs). Four discrete hybrid formulas used to formulate the BHTRKNM are provided by a continuous onestep hybrid trigonometrically fitted method with an offgrid point. We implement BHTRKNM in a blockbyblock fashion; in this way, the method does not suffer from the disadvantages of requiring starting values and predictors which are inherent in predictorcorrector methods. The stability property of the BHTRKNM is discussed and the performance of the method is demonstrated on some numerical examples to show accuracy and efficiency advantages.
1. Introduction
In what follows, we consider the numerical solution of the general secondorder IVPs of the formwhere , is an integer, and is the dimension of the system. Problems of form (1) frequently arise in several areas of science and engineering such as classical mechanics, celestial mechanics, quantum mechanics, control theory, circuit theory, astrophysics, and biological sciences. Equation (1) is traditionally solved by reducing it into a system of firstorder IVPs of double dimension and then solved using the various methods that are available for solving systems of firstorder IVPs (see Lambert [1, 2], Hairer and Wanner in [3], Hairer [4], and Brugnano and Trigiante [5, 6]).
Nevertheless, there are numerous methods for directly solving the special secondorder IVPs in which the first derivative does not appear explicitly and it has been shown that these methods have the advantages of requiring less storage space and fewer number of function evaluations (see Hairer [4], Hairer et al. [7], Simos [8], Lambert and Watson, and [9], Twizell and Khaliq [10]). Fewer methods have been proposed for directly solving secondorder IVPs in which the first derivative appears explicitly (see VigoAguiar and Ramos [11], Awoyemi [12], Chawla and Sharma [13], Mahmoud and Osman [14], Franco [15], and Jator [16, 17]). It is also the case that some of these IVPs possess solutions with special properties that may be known in advance and take advantage of when designing numerical methods. In this light, several methods have been presented in the literature which take advantage of the special properties of the solution that may be known in advance (see Coleman and Duxbury [18], Coleman and Ixaru [19], Simos [20], Vanden Berghe et al. [21], VigoAguiar and Ramos [11], Fang et al. [22], Nguyen et al. [23], Ramos and VigoAguiar [24], Franco and Gómez [25], and Ozawa [26]). However, most of these methods are restricted to solving special secondorder IVPs in a predictorcorrector mode.
Our objective is to present a BHTRKNM that is implemented in a blockbyblock fashion; in this way, the method does not suffer from the disadvantages of requiring starting values and predictors which are inherent in predictorcorrector methods (see Jator et al. [27], Jator [16], and Ngwane and Jator [28]). We note that multiderivative trigonometrically fitted block methods for have been proposed in Jator [29] and Jator [16]. However, the BHTRKNM proposed in this paper avoids the computation of higher order derivatives which have the potential to increase computational cost, especially, when applied to nonlinear systems. In this paper, we propose a BHTRKNM which is of order 3 and its application is extended to solving oscillatory systems, PDEs, and Hamiltonian systems including the energy conserving equation.
The organization of this article is as follows. In Section 2, we derive the BHTRKNM for solving (1). The analysis and implementation of the BHTRKNM are discussed in Section 3. Numerical examples are given in Section 4 to show the accuracy and efficiency of the BHTRKNM. Finally, the conclusion of the paper is given in Section 5.
2. Development of the BHTRKNM
In order to numerical integrate (1) we define the BHTRKNM as consisting of the following four discrete formulas:where , , , and are coefficients that depend on the steplength and frequency . In general, the frequency is chosen near the exact frequency of the true solution (see [30]). The coefficients of the method are chosen so that the method integrates the IVP (1) exactly where the solutions are members of the linear space .
The main method has the formwhere , , and , , and are to be determined coefficient functions of the frequency and stepsize. In order to derive the main method and additional methods we initially seek a continuous local approximation on the interval of the formwhere , , and , , are continuous coefficients. The first derivative of (4) is given byWe assume that is the numerical approximation to the analytical solution , is the numerical approximation to , and is an approximation to , .
The following theorem shows how the continuous method (4) is constructed. This is done by requiring that on the interval from to the exact solution is locally approximated by function (4) with (5) obtained as a consequence.
Theorem 1. Let , , , and be basis functions and let be a vector, where is the transpose. Define the matrix byand is obtained by replacing the th column of by the vector . Let the following conditions be satisfied:then the continuous representations (4) and (5) are equivalent to the following:
Proof. To prove this theorem, we use the approach given in Jator [17] with appropriate notational modification. We start by requiring that the method (4) be defined by the assumed basis functionswhere , , and are coefficients to be determined. Substituting (10) into (4) we get which is simplified to and expressed aswhere By imposing conditions (7) on (13), we obtain a system of five equations which can be expressed as where is a vector whose coefficients are determined via Cramer’s rule as where is obtained by replacing the th column of by . In order to obtain the continuous approximation, we use the elements of to rewrite (13) aswhose first derivative is given by
Remark 2. We note that, in the derivation of the BHTRKNM, the basis functions , , , and are chosen because they are simple to analyze. Nevertheless, other possible bases are possible (see Nguyen et al. [23]).
2.1. Specification of the Method
The continuous methods (8) and (9) which are equivalent to forms (4) and (5) are used to generate two discrete methods and two additional methods. The discrete and additional methods are then applied as a BHTRKNM for solving (1). We choose and evaluating (8) at and , respectively, gives the two discrete methods and which takes the form of the main method. Evaluating (9) at and , respectively, gives the additional methods and . The coefficients and their corresponding Taylor series equivalence of , , , and are, respectively, given as follows:
Remark 3. We note that the Taylor series expansions in (19) through (22) must be used when because the corresponding trigonometric coefficients given in these equations are vulnerable to heavy cancelations (see [8]).
2.2. Block Form
BHTRKNM is formulated from the four discrete hybrid formulas stated in (2) which are provided by the continuous onestep hybrid trigonometrically fitted method with one offgrid point given by (4) and its first derivative (5). We define the following vectors: where , . The methods in (2) specified by the coefficients (19)–(22) are combined to give the BHTRKNM, which is expressed aswhere , , , and are matrices of dimension four whose elements characterize the method and are given by the coefficients of (2).
3. Error Analysis and Stability
3.1. Local Truncation Error (LTE)
We define the local truncation error of (24) aswhere and is linear different operator.
Suppose that is sufficiently differentiable. Then, a Taylor series expansion of the terms in (25) about the point gives the following expression for local truncation error:where , , are constant coefficients (see [17]).
Definition 4. The block method (24) has algebraic order at least provided there exists a constant such that the local truncation error satisfies , where is the maximum norm.
Remark 5. (i) The local truncation error constants of of the block method (24) are given, respectively, by , where .
(ii) From the local truncation error constant computation, it follows that the method (24) has order at least three.
3.2. Stability
The linearstability of the BHTRKNM is discussed by applying the method to the test equation , where is a real constant (see [18]). Letting , it is easily shown as in [19] that the application of (24) to the test equation yieldswhere the matrix is the amplification matrix which determines the stability of the method. In the spirit of [22], the spectral radius of can be obtained from the characteristics equationwhere , , and are rational functions. We let in the following definition.
Definition 6. A region of stability is a region in the plane, throughout which and any closed curve given by defines the stability boundary of the method (see [22]). We note that the plot for the stability region of the BHTRKNM method is given in Figure 1.
Remark 7. It is observed that, in the plane, the BHTRKNM is stable for and (see Figure 1).
3.3. Implementation
The main method and the additional methods specified by (19)–(22) are combined to form the block method BHTRKNM (24), which is used to solve (1) without requiring starting values and predictors. BHTRKNM is implemented in a blockbyblock fashion using a Mathematica 10.0 code, enhanced by the feature for linear problems while nonlinear problems were solved by Newton’s method enhanced by the feature (see Keiper and Gear [33]). Mathematica can symbolically compute derivatives and so the entries of the Jacobian matrix which involve partial derivatives are automatically generated. In what follows, we summarize how BHTRKNM is applied.
Step 1. Choose , , and the number of blocks . Using (24), , and , the values of are simultaneously obtained over the subinterval , as and are known from the IVP (1).
Step 2. For and , the values of are simultaneously obtained over the subinterval , as and are known from the previous block.
Step 3. The process is continued for and to obtain the numerical solution to (1) on the subintervals .
In order to illustrate the efficiency of our method, we solved a variety of problems including oscillatory systems, PDEs such as the Telegraph equation, and Hamiltonian systems. The following methods are selected for comparison:(i)BHTRKNM given in this paper.(ii)ARKN: adapted RungeKuttaNyström method in [34] which has order five.(iii): difference scheme in [32].(iv)ESDIRK: explicit singly diagonally implicit RungeKutta method in [26].(v)FESDIRK: functionally fitted ESDIRK in [26].(vi)EFRK: exponentially fitted RungeKutta method (Method (b)) in Simos [8].(vii)N4: fourthorder standard RungeKuttaNyström method in [35].
4. Numerical Examples
In this section, numerical experiments are performed using a code in Mathematica 10.0 to illustrate the accuracy and efficiency of the method.
Example 1. We consider the following inhomogeneous IVP by Simos [8].where the analytical solution is given by
This example was solved using the order 3 BHTRKNM and the endpoint errors () obtained were compared to the order 4 exponentially fitted method given in Simos [8]. In Table 1 it is shown that BHTRKNM is more efficient than the method in Simos [8]. We also compare the computational efficiency of the two methods in Figure 2 by considering the FNEs (number of function evaluations) over integration steps for each method. This example illustrates that the BHTRKNM performs better.

Example 2. We consider the nonlinear Duffing equation which was also solved by Simos [8] and Ixaru and Vanden Berghe [31]:The analytical solution is given bywhere , , , , , , and . We choose .
We compare the endpoint global errors for our method with those of Simos [8] and Ixaru and Vanden Berghe [31]. We see from Table 2 that the results produced by our method are competitive to those given in Simos [8] and Ixaru and Vanden Berghe [31]. Hence our method is more accurate and efficient as demonstrated in Figure 3.
Example 3. We consider the following twobody problem which was solved by Ozawa [26] on :where is an eccentricity. The exact solution of this problem iswhere is the solution of Kepler’s equation . We choose .
We show in Table 3 that the results obtained using the BHTRKNM method are more accurate than the explicit singly diagonally implicit RungeKutta (ESDIRK) and the functionally fitted ESDIRK (FESDIRK) methods given in Ozawa [26]. In Figure 4, we also illustrate the efficiency advantage of the BHTRKNM method over those in Ozawa [26].
Example 4. We consider the stiff secondorder IVP (see [16] and references herein) ; where is an arbitrary parameter.
This problem was chosen to demonstrate the stability of the BHTRKNM (Figure 5). As mentioned in Remark 7, the method is stable when and .
(a)
(b)
4.1. Problems Where Appears Explicitly
Example 5. We consider the harmonic oscillator with frequency and small perturbation that was solved in Franco [15] and Guo and Yan [34]:where the analytical solution is given bywhere , , and . Guo and Yan [34] solved this problem using ARKN method. The results in Table 4 show that the BHTRKNM is competitive with the order 5 RungeKuttaNyström method.

4.2. Hyperbolic PDE
Example 6. We consider the given Telegraph equation (see Ding et al. [32]).
The exact solution is given by .
In order to solve this PDE using the BHTRKNM, we carry out the semidiscretization of the spatial variable using the secondorder finite difference method to obtain the following secondorder system in the second variable .where , , , , , , and , which can be written in the formsubject to the boundary conditions and where , and is , matrix arising from the semidiscretized system, and is a vector of constants.
The boundary conditions are chosen accordingly. This example was chosen to demonstrate that the BHTRKNM can be used to solve the Telegraph equation. In Table 5, the results produced by the BHTRKNM using and space step are compared to scheme (3.12) ( and ), time step , and space step , given in Ding et al. [32]. It is obvious from Table 5 that the BHTRKNM is more accurate than the method given in [32]. Moreover, the errors produced by BHTRKNM method using and space step are given in Figure 6.
