Abstract
We present a block hybrid functionally fitted Runge–Kutta–Nyström method (BHFNM) which is dependent on the stepsize and a fixed frequency. Since the method is implemented in a block-by-block fashion, the method does not require starting values and predictors inherent to other predictor-corrector methods. Upon deriving our method, stability is illustrated, and it is used to numerically solve the general second-order initial value problems as well as hyperbolic partial differential equations. In doing so, we demonstrate the method’s relative accuracy and efficiency.
1. Introduction
We consider the numerical solution of the second-order IVPwhere is a smooth function and is the dimension of the system and its application to hyperbolic partial differential equations. Frequently, whether analyzing or numerically solving systems like (1), the system is reduced to an equivalent system of first-order IVP of dimension to leverage methods developed for first-order systems (cf. Lambert [1, 2], Hairer et al. [3], Hairer [4], and Brugnano et al. [5, 6] in the numerical setting).
There are a wide range of methods that avoid order-reduction provided they have the special form . Importantly, these methods require less storage space and fewer function evaluations (cf. Hairer [4], Hairer et al. [7], Simos [8], Lambert et al., and [9], Twizell et al. [10]). Fewer methods have been developed which both avoid the order-reduction and solve the general second-order IVP (1) (see Vigo-Aguiar et al. [11], Awoyemi [12], Chawla and Sharma [13], Mahmoud and Osman [14], Franco [15], and Jator [16, 17]). Still other methods have been developed for classes of second-order IVP which use qualitative features of solutions (see Coleman and Duxbury [18], Coleman and Ixaru [19], Simos [20], Vanden et al. [21], Vigo-Aguiar et al. [11], Fang et al. [22], Nguyen et al. [23], Ramos and Vigo-Aguiar [24], Franco and Gomez [25], and Ozawa [26]).
Our objective is to present a block hybrid functionally fitted Runge–Kutta–Nyström method (BHFNM) that is implemented in a block-by-block fashion, which does not suffer from the disadvantages of requiring starting values and predictors inherent to predictor-corrector 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 [16, 29]. However, the BHFNM 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. It is also the case that the proposed method performs better than those given in Ngwane and Jator [30, 31]. In this paper, we propose a BHFNM which is of order and its application is extended to solving oscillatory systems and hyperbolic partial differential equations. It is in the same spirit as those presented by Ngwane and Jator [30, 31].
The organization of this article is as follows. In Section 2, we derive the BHFNM for solving (1). The analysis and implementation of the BHFNM are discussed in Section 3. Numerical examples are given in Section 4 to show the accuracy and efficiency of the BHFNM. Finally, the conclusion of the paper is given in Section 5.
2. Development of the BHFNM
To simplify the presentation we introduce the following notation: let denote the stepsize in the method, for any real number , and be a set of positive real numbers that are to be chosen, which will be used to specify intermediate points used in the method described below. Further, let .
In order to numerically integrate (1), we define the BHFNM which is applied on the partition , in which the step is given by combining the main methods: for each ,andwhere , , , , for and are coefficients to be determined and . We note that is the numerical approximation to the analytical solution evaluated at .
The coefficients of the method are chosen so that the method integrates the IVP (1) exactly whenever the solutions are members of the linear space , and they will depend on the frequency as well as the step-length .
In order to derive the main methods and additional methods, we initially seek a continuous local approximation on the interval , expressed in vector form aswhere is a vector coefficients to be uniquely determined and , , , and are basis functions. In order to determine the coefficients given by , we impose the following conditions on (5):and we obtain a system of equations which is expressed as
We note that is the numerical approximation to the analytical solution , is the numerical approximation to , and is an approximation to , for .
System (7) is solved with the aid of Mathematica to obtain the coefficients given by , which are substituted into (5) to give the continuous scheme given bywhich after simplification takes the formwhere , , and for are continuous coefficients. The first derivative of (9) is given by
Remark 1. We note that in the derivation of the BHFNM, the basis functions , , and are chosen because they are simple to analyze. Nevertheless, other possible bases are possible (see Nguyen et al. [23] and Hoang et al. [23]). We note that the zeros of the Chebyshev’s polynomial of the first kind or the zeros of the Legendre’s polynomial can also be chosen. However, we elect not to use them in this case because upon testing them the results were found to be the same as for the chosen points.
2.1. Specification of the Method
The continuous form (8) and its first derivative, which are equivalent to the forms (9) and (10), are used to generate four discrete methods and four additional methods. The discrete and additional methods are then applied as a BHFNM for solving (1). We choose , and evaluating (9) at , , gives the four discrete methods , which takes the form of the main methods (3). Evaluating (10) at gives the additional methods , with , which takes the form of the additional methods (4). The coefficients and the corresponding Taylor series equivalence of , where are, respectively, given in terms of as follows:
Remark 2. We note that the Taylor series expansions in (11) through (18) must be used when because the corresponding trigonometric coefficients given in these equations are vulnerable to heavy cancelations (see [8]).
3. Properties of the Method
3.1. Local Truncation Error
We define the local truncation errors (LTEs) of (3) and (4) specified by the coefficients in (11) through (18) asAssuming that is sufficiently differentiable, we can expand the terms in as a Taylor series about the point to obtain the expressions for the LTEs and for as
Remark 3. The local truncation error constants of the BHFNM formulated from and and specified by the coefficients in (11) through (18) are given, respectively, by , , , , , , and orders ; hence the BHFNM has order at least 5.
3.2. Block Form
BHFNM is formulated from the four discrete hybrid formulas stated in (3) which are provided by the continuous one-step hybrid trigonometrically fitted method with three off-grid points given by (9) and its first derivative (10). We define the following vectors:where . The methods in (3) specified by coefficients (11)-(18) are combined to give the BHFNM, which is expressed aswhere , , , and are square matrices of dimension eight whose elements characterize the method and are given by the coefficients of (3).
3.3. Stability
The linear-stability of the BHFNM 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 (23) to the test equation yieldswhere the matrix is the amplification matrix which determines the stability of the method.
Definition 4. A region of stability is a region in the -plane, throughout which the spectral radius and any closed curve given by define the stability boundary of the method (see [22]). We note that the plot for the stability region of the BHFNM method is given in Figure 1.
Remark 5. It is observed that in the -plane the BHFNM is stable for and (see Figure 1).
3.4. Implementation
The methods given by (11)-(18) are combined to form the block method BHFNM (23), which is used to solve (1) without requiring starting values and predictors. BHFNM is implemented in a block-by-block fashion using a Mathematica 10.0 code, enhanced by the features and (see Keiper and Gear [32]) for linear and nonlinear problems, respectively. In particular, the BHFNM is applied over , by choosing , and using (23) to simultaneously obtain the values of
4. Numerical Examples
The comparison of the following methods are given in this section:(i)Block Nyström method of order 5 based on polynomial basis and given in [33].(ii)Block Hybrid Trigonometrically Fitted Algorithm of order 5 given in [30].(iii)Block hybrid trigonometrically fitted Runge–Kutta–Nyström method of order 3 given in [31].(iv)Block hybrid functionally fitted Runge–Kutta–Nyström method (BHFNM) of order 5 given in this paper.
Example 6. We consider the following general second-order IVP given in [34] and the analytical solution is
This example is solved using the Block Nyström method based on polynomial basis given in [33], Block Hybrid Trigonometrically Fitted Algorithm of order 5 given in [30], Block hybrid trigonometrically fitted Runge–Kutta–Nyström method of order 3 given in [31], and the BHFNM given in this paper. The errors obtained which are compared for different step sizes and displayed in Table 1 show that BHFNM is more accurate.
Example 7. We consider the oscillatory system that is solved in [34]with and .
The exact solution of this system is
This example is solved using the Block Nyström method based on polynomial basis given in [33], Block Hybrid Trigonometrically Fitted Algorithm of order 5 given in [30], Block hybrid trigonometrically fitted Runge–Kutta–Nyström method of order 3 given in [31], and the BHFNM given in this paper. The errors obtained are compared for different step sizes and displayed in Table 2 show that BHFNM is more accurate.
4.1. Application to a PDE
and are calculated as follows: and , , where , .
Example 8. We consider the telegraph equation which was also solved in [33]The analytical solution is given byand the initial conditions are defined to match with the analytical solution.
This PDE is solved by first discretization the spatial variable via the finite difference method to obtainwhere , , , , , , , and , which can be written in the formsubject to the initial conditions , , , , , where , and S is a , matrix arising from the semidiscretized system.
This problem is solved using the BHFNM and the numerical results are displayed in Table 3 and Figure 2. The table gives and errors, while a comparison of the exact and numerical solutions is given in the figure plotted in 3 dimensions.
(a) Exact solution
(b) Numerical solution
Example 9. We consider the dissipative nonlinear wave equation given in [35]with initial conditionsThe analytical solution is given byand the initial conditions are defined to match with the analytical solution.
This PDE is solved by first discretization the spatial variable via the finite difference method to obtainwhere , , , , , , , and , which can be written in the formsubject to the initial conditions , , , , where , and S is a , , matrix arising from the semidiscretized system.
This problem is solved using the BHFNM and the numerical results are displayed in Table 4 and Figure 3. The table gives and errors, while a comparison of the exact and numerical solutions is given in the figure plotted in 3 dimensions.
(a) Exact solution
(b) Numerical solution
Example 10. We consider the dissipative nonlinear wave equation given in [35]with initial conditionsThe analytical solution is given byand the initial conditions are defined to match with the analytical solution.
This PDE is solved by first discretization the spatial variable via the finite difference method to obtainwhere , , , , , , , and , which can be written in the formsubject to the initial conditions , , , , where , and S is a ,, matrix arising from the semidiscretized system.
This problem is solved using the BHFNM and the numerical results are displayed in Table 5 and Figure 4. The table gives and errors, while a comparison of the exact and numerical solutions are given in the figure plotted in 3 dimensions.
(a) Exact solution
(b) Numerical solution
4.2. Estimating the Frequency
Some techniques for estimating the frequency are given in [24, 36–39]. In the spirit of [36], where the technique for estimating the frequency was to determine the roots of the polynomial associated with the LTE, we rewrite the principal terms given in (20) and (21) as follows:where , , are constants, , the positive represents the derivative, and is a linear differential operator and therefore invertible, which facilitates the simplification of . We estimate by setting any member of and solving for to obtain two possible solutions of as follows:This estimate has been tested on Example 6 on , for , where
We observe that when was used, an appropriate estimate for the parameter was provided at each step, in which case the computational time was increased and results were slightly less accurate than those given by the exact frequency, . We then noticed that, for this example, the results produced by setting were similar to those given by . Therefore, for this example, the appropriate choice for the frequency is . We do not claim that this approach will work in all cases; nevertheless, we will rigorously conduct research on this topic in our future work. Details of the numerical results are given in Table 6.
5. Conclusion
We have presented a BHFNM with coefficients that depend on a fixed frequency and the stepsize and that is implemented in a block-by-block fashion, which does not suffer from the disadvantages of requiring starting values and predictors inherent to predictor-corrector methods. It can be used to solve general second-order IVP, in particular, those that arise from the semidiscretization of hyperbolic PDE. Numerical studies have been presented to establish accuracy of the method. Since the exact value of the parameter is not always available, we have attempted to estimate an appropriate value of the parameter , in the spirit of [36], where the technique for estimating the frequency was to determine the roots of the polynomial associated with the LTE. We do not claim that this approach will work for all cases; hence we are determined to conduct rigorous research on this topic in our future work. We also note that the methods presented were not designed to incorporate a variable stepsize implementation, since we did not consider additional methods which are required to provide the error estimates as we proceed on the interval of interest. Nevertheless, our future research will include a variable stepsize implementation of the methods.
Data Availability
All references including data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.