Functionally Fitted Block Method for Solving the General Oscillatory Second-Order Initial Value Problems and Hyperbolic Partial Differential Equations

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.


Introduction
We consider the numerical solution of the second-order IVP   =  (, ,   ) , where  : R × R  × R  → R  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 2 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).
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 5 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.

Development of the BHFNM
To simplify the presentation we introduce the following notation: let ℎ > 0 denote the stepsize in the method,   =  0 + ℎ 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  0 =  ∪ {0}.
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 ⟨1, ,  2 ,  3 ,  4 , sin(), cos()⟩, 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 [  ,  +1 ], expressed in vector form as where  = [ℓ 0 , ℓ 1 , . . ., ℓ 6 ] ⊤ is a vector coefficients to be uniquely determined and Ω  () =   ,  = 0, 1, 2, 3, 4, Ω 5 () = sin , and Ω 6 () = cos  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  ∈  0 .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 by which after simplification takes the form where  0 (),  0 (), and   () for  ∈  0 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 Ω  () =   ,  = 0, 1, 2, 3, 4, Ω 5 () = sin , and Ω 6 () = cos  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  = {, , V, 1} = {1/4, 1/2, 3/4, 1}, 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  → 0 because the corresponding trigonometric coefficients given in these equations are vulnerable to heavy cancelations (see [8]).

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  = 0, . . ., .The methods in (3) specified by coefficients ( 11)-( 18) are combined to give the BHFNM, which is expressed as where  0 ,  1 ,  0 , and  1 are square matrices of dimension eight whose elements characterize the method and are given by the coefficients of (3).

Stability.
The linear-stability of the BHFNM is discussed by applying the method to the test equation   = − 2 , where  is a real constant (see [18]).Letting Υ = ℎ, it is easily shown as in [19] that the application of ( 23) to the test equation yields where the matrix (Υ 2 ; ) 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 ((Υ 2 ; )) ≤ 1 and any closed curve given by ((Υ 2 ; )) = 1 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.

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].
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.
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.Example 8. We consider the telegraph equation which was also solved in [33] The analytical solution is given by and 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 obtain   2 sin(  ) sin())(sin(  ) cos()), which can be written in the form subject to the initial conditions U(0) = sin(  ),  = 1, . . ., , U  (0) = 0,  = 1, . . ., , where f(, U, U  ) = SU + G, 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  2 errors, while a comparison of the exact and numerical solutions is given in the figure plotted in 3 dimensions.
Example 10.We consider the dissipative nonlinear wave equation given in [35] with initial conditions  (, 0) = sin () ,   (, 0) = −Υ sin () . (40) The analytical solution is given by and the initial conditions are defined to match with the analytical solution.
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,  = 1.We then noticed that, for this example, the results produced by setting  =  0 were similar to those given by  = 1.Therefore, for this example, the appropriate choice for the frequency is  =  0 .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.

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.

( 21 ) 3 .
Remark The local truncation error constants of the BHFNM formulated from (2) and (3) and specified by the coefficients in (