Abstract
In some cases, high-order methods are known to provide greater accuracy with larger step-sizes than lower order methods. Hence, in this paper, we present a Block Hybrid Method (BHM) of order 11 for directly solving systems of general second-order initial value problems (IVPs), including Hamiltonian systems and partial differential equations (PDEs), which arise in multiple areas of science and engineering. The BHM is formulated from a continuous scheme based on a hybrid method of a linear multistep type with several off-grid points and then implemented in a block-by-block manner. The properties of the BHM are discussed and the performance of the method is demonstrated on some numerical examples. In particular, the superiority of the BHM over the Generalized Adams Method (GAM) of order 11 is established numerically.
1. Introduction
General second-order differential equations frequently arise in several areas of science and engineering, such as celestial mechanics, quantum mechanics, control theory, circuit theory, astrophysics, and biological sciences. Several of these differential equations have oscillatory solutions, for instance, the initial value problem (IVP)where , is an integer, , is the dimension of the system, is a real nonsingular diagonalizable constant matrix having large eigenvalues, and has been studied (see [1–9]). Special cases of (1) have also been extensively discussed in the literature ([10–14], Hairer [15]).
The method proposed in this paper can solve (1), as well as the general second-order IVPwhere , is an integer, and is the dimension of the system. Equation (2) is conventionally solved by converting it into an equivalent first-order system of double dimension and then solved using standard methods that are available in the literature for solving systems of first-order IVPs (see Lambert [16], Hairer et al. [17], and Brugnano et al. [18]). In general, these methods are implemented in a step-by-step fashion in which, on the partition , an approximation is obtained at only after an approximation at has been computed, where for some constant step-size and integer ,
Several methods have been proposed for directly solving the special second-order IVP, where the function does not depend on . It has been shown that such methods have the advantages of requiring less storage space and fewer number of function evaluations (see Hairer [15], Hairer et al. [19], Lambert et al. [20], and Twizell et al. [21]). Nevertheless, these methods are implemented as predictor-corrector methods and require starting values as well. In which case, the cost of implementation is increased, especially, for higher order methods.
In this paper, we propose a high-order BHM of order 11 which provides greater accuracy even when larger step-sizes are used. We show that the BHM of order 11 is applied to directly solve Hamiltonian systems as well as the general form (2). The BHM is formulated from a continuous scheme based on a hybrid method of a linear multistep type with several off-grid points and then implemented in a block-by-block manner. In this way, the method is self-starting and implemented without predictors (see Jator and Oladejo [22], Jator et al. [9], Jator [7], and Ngwane and Jator [6, 23, 24]). It is shown that the method can also be used to directly solve non-Hamiltonian systems with embedded first derivatives as well as partial differential equations via the method of lines. In particular, the superiority of the BHM over the Generalized Adams Method (GAM) of order 11 is established numerically.
The article is organized as follows. In Section 2, we derive our hybrid method and specify the coefficients as well. We discuss the properties and implementation of the method in Section 3. Numerical examples are given in Section 4 and concluding remarks are given in Section 5.
2. Derivation
We propose a BHM for the IVP (2) which can be solved by advancing the from to by fixing two interpolation points and a set of distinct collocation points . We then choose the coefficients of the method, such that the method integrates the IVP exactly, whenever the solutions are members of the linear space . Thus, we initially seek a continuous local approximation on the interval , which is expressed in vector form aswhere are coefficients in to be uniquely determined. To determine these coefficients, we interpolate at and and then collocate at . In particular, we determine these coefficients by imposing the following conditions on (4):to obtain the following system of equations represented in matrix formThe system is solved with the aid of Mathematica to obtain the coefficients, , which are substituted into (4) to obtain the continuous formThe continuous scheme (7) is simplified and evaluated at , , to give the set of formulas (8), whose coefficients are specified in Tables 1 and 2. In the same vein, the first derivative of the continuous scheme (7) is simplified and evaluated at , , to give the set of formulas (9), whose coefficients are also specified in Tables 1 and 2.
In order to numerically integrate (1), we advance the BHM from to on the partition , whereby the step is given by combining the following methods:where , , , , , are coefficients given in Tables 1 and 2. We note that , , and .
3. Properties of the Method
Local Truncation Error. We define the local truncation errors (LTEs) of (8) and (9), specified by the coefficients in Tables 1 and 2 asAssuming that is sufficiently differentiable, we can expand the terms in and as Taylor series about the point to obtain the expressions for the LTEs as follows: , ,
where is the order and and are the error constants.
Remark 1. The local truncation error constants of the BHM formulated from (8) and (9) and specified by the coefficients in Tables 1 and 2 as well as choosing are displayed in Table 3.
Stability. BHM is formulated from (8) and (9) and defined in vector form as follows:where . The methods in (8) and (9) specified by the coefficients from Tables 1 and 2 are combined to give the BHM, which is expressed aswhere , and are square matrices of dimension twenty whose elements characterize the method and are given by the coefficients of (8) and (9).
The linear stability of the BHM is discussed by applying the method to the test equation , where is the frequency and is the damping (see [2]). Letting and , it is easily shown as in [2] that the application of (13) to the test equation yieldswhere the matrix is the amplification matrix which determines the stability of the method.
Definition 2. The region of stability of the BHM is region throughout which the spectral radius (see [2]).
Remark 3. It is observed that in the -plane, the BHM is stable for as demonstrated in example 4.6 and confirmed in Figure 1(b).
(a)
(b)
Implementation. The methods specified by (8) and (9) are combined to form the BHM (37), which is then implemented in a block form on the nonoverlapping interval . This approach has an advantage over the usual predictor-corrector methods, since it is self-starting. Initially, for , the block method is simultaneously applied on the subinterval to obtain . Using the known values from the previous block, we obtain values on the next subinterval , and the process is repeated until we reach the final interval . In the spirit of [22], our code was written in Mathematica and the steps are given in Algorithm 1.
| ||||||||||||||||
4. Numerical Examples
We solve a variety of oscillatory general second-order IVPs and Hamiltonian systems to illustrate the accuracy of our hybrid method. This is done via the maximum global error, where the maximum global error of the approximate solution is calculated as , where is the exact solution and is the approximate solution on the partition . We have also included the maximum global errors for the Hamiltonian () as well as figures showing the global error of the Hamiltonian () for each problem.
4.1. Application to Hamiltonian Systems
The following methods have been compared in this subsection.(i)Block Hybrid Method (BHM) of order 11 given in this paper,(ii)Generalized Adams method (GAM) of order 11 in [18].
Example 4. We consider the following oscillatory nonlinear system that was studied in [10]:with , and initial conditionssuch that the analytic solution of this system is .
The Hamiltonian is
This example is solved using both our BHM and GAM and the errors in the solution and Hamiltonian obtained are compared for different step-sizes. The results displayed in Table 4 show that our method is more accurate. Details of the results are also displayed in Figure 1.
Example 5. We consider the Fermi-Pasta-Ulam problem that has been studied in [10, 25]wherewith initial conditionsThe remaining initial conditions are equal to zero, and the Hamiltonian is
This example is solved using both our BHM and GAM and the errors in the Hamiltonian obtained are compared for different step-sizes. The results displayed in Table 5 show that our method is more accurate. Details of the results are also displayed in Figure 2.
(a) BHM for
(b) Adams for
Example 6. We consider the pendulum oscillator given bywith initial conditionsand Hamiltonian
This example is solved using both our BHM and GAM and the errors in the Hamiltonian obtained are compared for different step-sizes. The results displayed in Table 6 show that our method is more accurate. Details of the results are also displayed in Figures 3 and 4.
(a) BHM for
(b) Adams for
(a) BHM for
(b) Adams for
Example 7. We consider the perturbed Kepler’s problem given bywhere the exact solution of this problem isand the Hamiltonian isThe system also has the angular momentum as a first integral. We let .
This example is solved using both our BHM and GAM and the errors in the solution, Hamiltonian, and momentum obtained are compared for different step-sizes. The results displayed in Table 7 show that our method is more accurate. Details of the results are also displayed in Figure 5.
(a) BHM for
(b) Adams for
(c) Exact Flow
4.2. Oscillatory Problems with Present
The following methods have been compared in this subsection.(i)Block Hybrid Method (BHM) of order 11 given in this paper,(ii)Generalized Adams method (GAM) of order 11 in [18].
Example 8. We consider the oscillatory system that is solved in [2]: with and .
The exact solution of this system is .
This example is solved using both our BHM and GAM and the errors in the solution and the first derivative obtained are compared for different step-sizes. The results displayed in Tables 8 and 9 show that our method is more accurate. We note that this example was chosen to show that the BHM performs well on an oscillatory general second-order system of IVPs.
Example 9. We consider the following general second-order IVP given in [2]: and the analytical solution is
This example is solved using both our BHM and GAM and the errors obtained are compared for different step-sizes. The results displayed in Table 10 show that our method is more accurate.
4.3. Efficiency Curves for Four Methods
In this subsection, we have solved examples 4.1–4.6 and generated efficiency curves for four methods. In particular, we have compared three methods given in the literature to the BHM given in this paper by plotting versus as well as versus . The plots given in Figure 7 show that the BHM is the most efficient in terms of the size of , since smaller values of give greater accuracy. However, the GAM performs better in terms of CPU time. Nevertheless, the BHM remains competitive pertaining to the CPU time. We note that the high accuracy of the BHM is due to the fact that the error constants for all the members of the BHM are smaller than those of the GAM. For instance, the error constant for the main method for the BHM is , while the error constant for the GAM is . In this subsection, the following methods have been compared and the details of these comparisons are given in Figure 7:(i)Block Hybrid Method (BHM) of order 11 given in this paper,(ii)Generalized Adams method (GAM) of order 11 given in [18],(iii)Block Continuous Hybrid Integrator of order 9 given in [26],(iv)the 10-step block Falkner-type method of order 11 given in [27].
4.4. Application to a PDE
Example 10. We consider the Telegraph equation which was also solved in [22]:The analytical solution is given byand the initial conditions are defined to match 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 L is an matrix arising from the semidiscretized system and is a vector of constants (Figure 8).
This example is chosen to show that the BHM performs well on PDEs such as the Telegraph equation and this is confirmed by the results displayed in Figure 6.
(a) BHM, , EH
(b) BHM, , EM
(c) GAM, , EH
(d) GAM, , EM
(a) Example 4
(b) Example 4
(c) Example 5
(d) Example 5
(e) Example 6
(f) Example 6
(g) Example 7
(h) Example 7
(i) Example 8
(j) Example 8
(k) Example 9
(l) Example 9
(a)
(b)
4.5. Application to a Stiff Problem
Example 11. We consider the stiff second-order IVP (see [28]), where is an arbitrary parameter.
This problem was chosen to justify that the stability of the BHM is achieved when . In Table 11, we give the absolute errors at selected values of , which indicate that, choosing , the method is stable since, for this value of , . However, for , ; hence the method becomes unstable.
5. Conclusion
We have proposed a BHM of order 11 for solving systems of general second-order IVPs and Hamiltonian systems. The BHM is formulated from a continuous scheme based on a hybrid method of a linear multistep type with several off-grid points and implemented in a block-by-block manner. In this way, starting values and predictors which are generally considered drawbacks in the implementation of predictor-corrector methods are avoided. The properties of the BHM are discussed and the performance of the method is demonstrated on several numerical examples. In particular, the superiority in terms of accuracy of the BHM over comparable methods given in the literature is established numerically. Our future research will be focused on developing methods for oscillatory systems arising from the semidiscretization of hyperbolic partial differential equations.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.