Mathematical Problems in Engineering
Volume 2018 (2018), Article ID 3750274, 15 pages
https://doi.org/10.1155/2018/3750274

Research Article

Integrating Oscillatory General Second-Order Initial Value Problems Using a Block Hybrid Method of Order 11

Samuel N. Jator and Kindyl L. King

Department of Mathematics and Statistics, Austin Peay State University, Clarksville, TN 37044, USA

Correspondence should be addressed to Samuel N. Jator; jators@apsu.edu

Received 5 March 2018; Revised 30 April 2018; Accepted 14 May 2018; Published 11 June 2018

Academic Editor: Zhangxin Chen

Copyright © 2018 Samuel N. Jator and Kindyl L. King. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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 IVP

where ,   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 as

where 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 form

The system is solved with the aid of Mathematica to obtain the coefficients, , which are substituted into (4) to obtain the continuous form

The 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.

Table 1: Coefficients for (8) and (9).

0
0

0

1

1
1/5

1

1

2
2/5

1

1

3
3/5

1

1

4
4/5

1

1

5
5/5

0

1

Table 2: Coefficients for (8) and (9).

6
6/5

1

1

7
7/5

1

1

8
8/5

1

1

9
9/5

1

1

10
10/5

1

1

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 as

Assuming 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.

Table 3: Coefficients for (8) and (9) including order and error constant and .

0
0
11

1
1/5

2
2/5

3
3/5

4
4/5

5
1

6
6/5

7
7/5

8
8/5

9
9/5

10
2

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 as

where , 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 yields

where 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)

Figure 1: (a) The stability region for the BHM plotted in the -plane and (b) the stability interval plotted in the (-plane when () for the BHM.

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.

procedure ENTER PARTITIONS(, variables)

    For ,  

    Generate block from (13)                         ▹ Initially on

    Solve                         ▹ Block-by-block with respect to the

variables(solutions) on .

    Obtain . ▹ No need for starting

values and predictors

end Procedure

Algorithm 1: BHM algorithm.

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 conditions

such 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.

Table 4: Comparison of global error for example 4.1 on the interval .

BHM-Error
GAM-Error
BHM-HE
GAM-HE

1/2

1/4

1/8

1/16

1/32

Example 5. We consider the Fermi-Pasta-Ulam problem that has been studied in [10, 25]

where

with initial conditions

The 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.

Table 5: Comparison of global error for example 4.2 on the interval .

h
BHM-HE
GAM-HE

1/8

1/16

1/32

1/64

1/128

(a) BHM for

(b) Adams for

Figure 2: Results for 4.1. The logarithm of the global error of the Hamiltonian, on the interval .

Example 6. We consider the pendulum oscillator given by

with initial conditions

and 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.

Table 6: Comparison of global error for example 4.3 on the interval .

h
BHM-HE
GAM-HE

1/2

1/4

1/8

1/16

1/32

(a) BHM for

(b) Adams for

Figure 3: Results for example 4.2. The logarithm of the global error of the Hamiltonian, on the interval .

(a) BHM for

(b) Adams for

Figure 4: Results for example 4.3. The logarithm of the global error of the Hamiltonian, on the interval .

Example 7. We consider the perturbed Kepler’s problem given by

where the exact solution of this problem is

and the Hamiltonian is

The 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.

Table 7: Comparison of global error for example 4.4 on the interval .

h
BHM-Error
GAM-Error
BHM-HE
GAM-HE
BHM-EL
GAM-EL

1/2

1/4

1/8

1/16

1/32

(a) BHM for

(b) Adams for

(c) Exact Flow

Figure 5: Results for example 4.3.

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.

Table 8: Comparison of global error for example 4.5 on the interval .

h
BHM-Error
GAM-Error
BHM-Error
GAM-Error

1

1/2

1/4

1/8

1/16

Table 9: Comparison of global error of the first derivatives for example 4.5 on the interval .

h
BHM-Error
GAM-Error
BHM-Error
GAM-Error

1

1/2

1/4

1/8

1/16

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.

Table 10: Comparison errors for example 4.6 on .

h
BHM-Errors
Adams-Errors

2

1

1/2

1/4

1/8

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 by

and 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 obtain

where ,  ,  ,  ,  , ,  ,  , and , which can be written in the form

subject 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

Figure 6: Results for example 4.4: the logarithm of the global error of the Hamiltonian and the logarithm of the global error of Momentum on the interval .

(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

Figure 7: Efficiency curves for examples 4.1-4.6. The labels (a), (c), (e), (g), (i), and (k) represent the versus , where . The labels (b), (d), (f), (h), (j), and (l) represent the versus .

(a)

(b)

Figure 8: Errors for various values of N and M for example 4.7.

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.

Table 11: Demonstrating stability using absolute errors for example 4.7. on the interval .

BHM-Err
BHM-Err

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.

References

B. Wang, K. Liu, and X. Wu, “A Filon-type asymptotic approach to solving highly oscillatory second-order initial value problems,” Journal of Computational Physics, vol. 243, pp. 210–223, 2013.

K. Liu and X. Wu, “Multidimensional {ARKN} methods for general oscillatory second-order initial value problems,” Computer Physics Communications, vol. 185, no. 7, pp. 1999–2007, 2014.

M. M. Chawla and S. R. Sharma, “Families of three-stage third order Runge-Kutta-Nystrom methods for ,” The Journal of the Australian Mathematical Society. Series B. Applied Mathematics, vol. 26, no. 3, pp. 375–386, 1985.

H. Ramos, Z. Kalogiratou, T. Monovasilis, and T. E. Simos, “A trigonometrically fitted optimized two-step hybrid block method for solving initial-value problems of the form with oscillatory solutions,” in Proceedings of the International Conference on Numerical Analysis and Applied Mathematics 2014, ICNAAM 2014, 2014.

B.-y. Guo and J.-p. Yan, “Legendre-Gauss collocation method for initial value problems of second order ordinary differential equations,” Applied Numerical Mathematics, vol. 59, no. 6, pp. 1386–1408, 2009.

F. F. Ngwane and S. N. Jator, “Block hybrid method using trigonometric basis for initial value problems with oscillating solutions,” Numerical Algorithms, vol. 63, no. 4, pp. 713–725, 2013.

S. N. Jator, “Implicit third derivative Runge-Kutta-Nyström method with trigonometric coefficients,” Numerical Algorithms, vol. 70, no. 1, pp. 133–150, 2015.

S. N. Jator, “A continuous two-step method of order 8 with a block extension for ,” Applied Mathematics and Computation, vol. 219, no. 3, pp. 781–791, 2012.

S. N. Jator, S. Swindell, and R. French, “Trigonometrically fitted block Numerov type method for ,” Numerical Algorithms, vol. 62, no. 1, pp. 13–26, 2013.

B. Wang, A. Iserles, and X. Wu, “Arbitrary-Order Trigonometric Fourier Collocation Methods for Multi-Frequency Oscillatory Systems,” Foundations of Computational Mathematics, vol. 16, pp. 158–181, 2016.

H. S. Nguyen, R. B. Sidje, and N. H. Cong, “Analysis of trigonometric implicit Runge-Kutta methods,” Journal of Computational and Applied Mathematics, vol. 198, no. 1, pp. 187–207, 2007.

Y. Fang, Y. Song, and X. Wu, “A robust trigonometrically fitted embedded pair for perturbed oscillators,” Journal of Computational and Applied Mathematics, vol. 225, no. 2, pp. 347–355, 2009.

T. E. Simos, “Dissipative trigonometrically-fitted methods for second order IVPs with oscillating solution,” International Journal of Modern Physics C, vol. 13, no. 10, pp. 1333–1345, 2002.

M. Van Daele and G. Vanden Berghe, “Geometric numerical integration by means of exponentially-fitted methods,” Applied Numerical Mathematics, vol. 57, no. 4, pp. 415–435, 2007.

E. Hairer, “A one-step method of order 10 for ,” IMA Journal of Numerical Analysis (IMAJNA), vol. 2, no. 1, pp. 83–94, 1982.

J. D. Lambert, Numerical Methods for Ordinary Differential Systems. The Initial Value Problem, John Wiley & Sons, New York, NY, USA, 1993.

E. Hairer and G. Wanner, Solving Ordinary Differential Equations II, Springer, Berlin, Germany, 2nd edition, 1996.

L. Brugnano and D. Trigiante, Solving ODEs By Multistep Initial and Boundary Value Methods, Gordon Breach, Amsterdam, 1998.

E. Hairer, S. P. Nørsett, and G. Wanner, Solving Ordinary Differential Equations I: Nonstiff Problems, vol. 8 of Springer Series in Computational Mathematics, Springer, Berlin, Germany, 1993.

J. D. Lambert and A. Watson, “A sysmetric multistep method for periodic initial value problem,” Mathematical Problems in Engineering, vol. 18, pp. 189–202, 1976.

E. H. Twizell and A. Q. M. Khaliq, “Multiderivative methods for periodic initial value problems,” SIAM Journal on Numerical Analysis, vol. 21, no. 1, pp. 111–122, 1984.

S. N. Jator and H. B. Oladejo, “Block Nyström method for singular differential equations of the Lane-Emden type and problems with highly oscillatory solutions,” International Journal of Applied and Computational Mathematics, vol. 3, no. suppl. 1, pp. S1385–S1402, 2017.

S. N. Jator, “On a class of hybrid methods for ,” International Journal of Pure and Applied Mathematics, vol. 59, no. 4, pp. 381–395, 2010.

S. N. Jator, “Solving second order initial value problems by a hybrid multistep method without predictors,” Applied Mathematics and Computation, vol. 217, no. 8, pp. 4036–4046, 2010.

E. Hairer and C. Lubich, “Long-time energy conservation of numerical methods for oscillatory differential equations,” SIAM Journal on Numerical Analysis, vol. 38, no. 2, pp. 414–441, 2000.

O. A. Akinfenwa, “Ninth order block piecewise continuous hybrid integrators for solving second order ordinary differential equations,” International Journal of Differential Equations and Applications, vol. 12, pp. 49–67, 2013.

H. Ramos, S. Mehta, and J. Vigo-Aguiar, “A unified approach for the development of -step block Falkner-type methods for solving general second-order initial-value problems in {ODE}s,” Journal of Computational and Applied Mathematics, vol. 318, pp. 550–564, 2017.

R. D'Ambrosio, M. Ferro, and B. Paternoster, “Two-step hybrid collocation methods for ,” Applied Mathematics Letters, vol. 22, no. 7, pp. 1076–1080, 2009.