J. Nig. Soc. Phys. Sci. 5 (2023) 1017

Journal of the
Nigerian Society

of Physical
Sciences

An Accuracy-preserving Block Hybrid Algorithm for the
Integration of Second-order Physical Systems with Oscillatory

Solutions

Joshua Sundaya,∗, Joel N. Ndama, Lydia J. Kwarib

aDepartment of Mathematics, University of Jos, Jos 930003, Nigeria
bDepartment of Mathematics, Federal College of Education, Pankshin 933105, Nigeria

Abstract

It is a known fact that in most cases, to integrate an oscillatory problem, higher order A-stable methods are often needed. This is because such
problems are characterized by stiffness, chaos and damping, thus making them tedious to solve. However, in this research, an accuracy-preserving
relatively lower order Block Hybrid Algorithm (BHA) is proposed for solution of second-order physical systems with oscillatory solutions.
The sixth order algorithm was derived using interpolation and collocation of power series within a single step interval [tn, tn+1]. In order to
circumvent the Dahlquist-barrier and also obtain an accuracy-preserving algorithm, four off-step points were incorporated within the single step
interval. A number of special cases of oscillatory problems were solved using the proposed method and the results obtained clearly showed that
it outperformed other existing methods we compared our results with even though the BHA is of lower order relative to such methods. Some of
the second-order physical systems considered were the Kepler, Bessel and damped problems. Some important properties of the BHA were also
analyzed and the results of the analysis showed that it is consistent, zero-stable and convergent.

DOI:10.46481/jnsps.2023.1017

Keywords: Accuracy-preserving, Algorithm, Block hybrid method, Oscillation, Physical systems, Second-order

Article History :
Received: 28 August 2022
Received in revised form: 05 December 2022
Accepted for publication: 11 December 2022
Published: 14 January 2023

c© 2023 The Author(s). Published by the Nigerian Society of Physical Sciences under the terms of the Creative Commons Attribution 4.0 International license

(https://creativecommons.org/licenses/by/4.0). Further distribution of this work must maintain attribution to the author(s) and the published article’s title, journal citation, and DOI.

Communicated by: T. Latunde

1. Introduction

Second-order physical systems with oscillatory solutions find
applications in diverse areas of human endeavors like engineer-
ing and sciences. Such systems are applied in vibration of
mass-spring systems, astrophysics, control theory, mechanics,
circuit theory, biology among others [1, 2]. In this research,
an accuracy-preserving BHA shall be derived for the direct in-

∗Corresponding author tel. no: +234 7034884482
Email address: sundayjo@unijos.edu.ng (Joshua Sunday)

tegration of second order oscillatory physical systems of the
form

y′′(t) = f
(
t, y(t), y′(t)

)
(1)

subject to the initial conditions

y(t0) = y0, y
′(t0) = y

′

0 (2)

on the interval t ∈ [t0, tN ], where f : < × <m → <m, N is
an integer and m the dimension of equation (1). Equation (1)
is assumed to satisfy the uniqueness theorem stated in Theorem
1.1.

1



Sunday et al. / J. Nig. Soc. Phys. Sci. 5 (2023) 1017 2

Theorem 1.1 [3]
Let,

y(n)(t) = f
(
t, y(t), y′(t), ..., y(n−1)(t)

)
, y(k)(t0) = yk (3)

where k = 0, 1, ..., (n − 1),y(0) = y, y(t) and f are scalars.
Let < be the region defined by the inequalities t0 ≤ t ≤ t0 +
a, |sk − yk| ≤ bk, k = 0, 1, ..., n − 1, where yk ≥ 0 for k > 0.
Suppose the function f (t, s0, s1, ..., sn−1) in (3) is non-negative,
continuous and non-decreasing in t, and continuous and non-
decreasing in sk for each k = 0, 1, ..., n−1 in the region <. If in
addition f (t, y0, y1, ..., yn−1) , 0 in < for t > t0, then the initial
value problem (3) has at most one solution in <”.

The derivation will be carried out via a continuous scheme
based on linear multistep method by incorporating four off-step
points. The implementation of the algorithm will be effected in
a block-by-block mode, thus making it self-starting (i.e. with-
out the need for predictors).

Equations of the form (1) can be solved by first transform-
ing them into their equivalent system of first order differen-
tial equations and then employing an appropriate method [4-
6]. However, one of the setbacks of such approach is that some
of the vital properties and characteristics of the higher order
differential equations are lost in the course of the conversion.
Besides, coding such methods are often cumbersome since in
most cases subroutines have to be incorporated to provide the
starting values.

Some methods have also been proposed in literature for the
direct solution of special second order differential equations.
Such equations are termed ‘special’ because they do not depend
on y′, in order words they are of the form y′′ = f (x, y). These
methods often require fewer function evaluations and less mem-
ory space [7-9].

Over the years, several authors have proposed different meth-
ods for the direct solution of oscillatory problems of the form
(1). Ref[2] proposed an eleventh order block hybrid method
for the direct solution of system of second order differential
equations including Hamiltonian systems. The method was
derived from continuous scheme via hybrid method approach
with several off-step points. The method was implemented in a
block manner. Ref[10] formulated a continuous explicit hybrid
method for the solution of second order differential equations.
The authors interpolated the basis function at both grid and off-
grid points while the differential systems were collocated at se-
lected points. The authors further derived starting values of the
same order with the methods by adopting Taylor series expan-
sion to circumvent the inherent disadvantage of starting values
of lower order. Ref[11] developed an order eight implicit block
method for the solution of second order differential equations.
The authors adopted the Hermite polynomial as basis function
to construct the method that comprises first and second deriva-
tives. The basic properties of the method were analysed and the
method was implemented on some linear and nonlinear second
order differential equations. Other researchers that also devel-
oped direct methods for solving problems of the form (1) are
[12-23].

2. Derivation and Implementation of the BHA

2.1. Derivation of the BHA

The accuracy-preserving BHA shall be derived by seeking
a continuous approximate solution Y (t) to the second order os-
cillatory problems of the form (1) on the interval [tn, tn+1]. This
is expressed compactly in a vector form as,

Y (t) =
[
1 t t2 ... t7

]


σ0
σ1
σ2
.
.
.
σ7


(4)

where σ0,σ1,σ2, ...,σ7are uniquely determined parameters in
<n. Equation (4) is interpolated at the points (tn+τ, yn+τ) , τ =
r, s and collocated at the points (tn+κ, fn+κ) , κ = 0, p, q, r, s, 1,
where τ and κare the interpolation and collocation points re-
spectively. The off-step points p, q, r, s are taken as p = 1/5,
q = 2/5, r = 3/5 and s = 4/5. The parameters σ0,σ1,σ2, ...,σ7are
determined by imposing the conditions

Y (tn+τ) = yn+τ, τ = r, s (5)

Y′′ (tn+κ) = fn+κ, κ = 0, p, q, r, s, 1 (6)

to obtain the system of equations



1 tn+r t2n+r t
3
n+r t

4
n+r t

5
n+r t

6
n+r t

7
n+r

1 tn+s t2n+s t
3
n+s t

4
n+s t

5
n+s t

6
n+s t

7
n+s

0 0 2 6tn 12t2n 20t
3
n 30t

4
n 42t

5
n

0 0 2 6tn+p 12t2n+p 20t
3
n+p 30t

4
n+p 42t

5
n+p

0 0 2 6tn+q 12t2n+q 20t
3
n+q 30t

4
n+q 42t

5
n+q

0 0 2 6tn+r 12t2n+r 20t
3
n+r 30t

4
n+r 42t

5
n+r

0 0 2 6tn+s 12t2n+s 20t
3
n+s 30t

4
n+s 42t

5
n+s

0 0 2 6tn+1 12t2n+1 20t
3
n+1 30t

4
n+1 42t

5
n+1





σ0
σ1
σ2
.
.
.
σ7



=



yn+r
yn+s
fn
fn+p
fn+q
fn+r
fn+s
fn+1


(7)

The system (7) is solved using SWP 5.5 for the parameters σ′j s
and then substituted into equation (4) to get the continuous form

2



Sunday et al. / J. Nig. Soc. Phys. Sci. 5 (2023) 1017 3

Y (t) =
[
1 t t2 ... t7

]


1 tn+r t2n+r t
3
n+r t

4
n+r t

5
n+r t

6
n+r t

7
n+r

1 tn+s t2n+s t
3
n+s t

4
n+s t

5
n+s t

6
n+s t

7
n+s

0 0 2 6tn 12t2n 20t
3
n 30t

4
n 42t

5
n

0 0 2 6tn+p 12t2n+p 20t
3
n+p 30t

4
n+p 42t

5
n+p

0 0 2 6tn+q 12t2n+q 20t
3
n+q 30t

4
n+q 42t

5
n+q

0 0 2 6tn+r 12t2n+r 20t
3
n+r 30t

4
n+r 42t

5
n+r

0 0 2 6tn+s 12t2n+s 20t
3
n+s 30t

4
n+s 42t

5
n+s

0 0 2 6tn+1 12t2n+1 20t
3
n+1 30t

4
n+1 42t

5
n+1



−1 

yn+r
yn+s
fn
fn+p
fn+q
fn+r
fn+s
fn+1


(8)

On simplification and evaluation of the continuous scheme (8), we obtain

Y (t) =
∑
τ=r,s

ατ(x)yn+τ + h
2

 1∑
j=0

β j(t) fn+ j + βk(t) fn+k

 , k = p, q, r, s (9)
where

αr (x) = −5x + 4
αs(x) = −3 + 5x
β0(x) = −

625
1008 x

7 + 12548 x
6 −

425
96 x

5 + 12532 x
4 −

137
72 x

3 + 12 x
2 −

397
6300 x +

1
375

βp(x) =
3125
1008 x

7 −
875
72 x

6 + 177596 x
5 −

1925
144 x

4 + 256 x
3 −

14059
50400 x +

127
3000

βq(x) = −
3125
504 x

7 + 162572 x
6 −

1475
48 x

5 + 2675144 x
4 −

25
6 x

3 −
157
1260 x +

29
375

βr (x) =
3125
504 x

7 −
125
6 x

6 + 122548 x
5 −

325
24 x

4 + 259 x
3 −

5813
25200 x +

161
1500

βs(x) = −
3125
1008 x

7 + 1375144 x
6 −

1025
96 x

5 + 1525288 x
4 −

25
24 x

3 −
1

1575 x +
4

375
β1(x) =

625
1008 x

7 −
125
72 x

6 + 17596 x
5 −

125
144 x

4 + 16 x
3 −

107
50400 x −

1
3000


(10)

and x is expressed as

x =
t − tn

h
(11)

Solving (9) for the independent solution at the grid points gives the method

Yn+ j =
1∑

i=0

( jh)i

i !
y(i)n + h

2

 1∑
j=0

η j(x) fn+ j + ηk(x) fn+k

 , k = p, q, r, s (12)
where

η0(x) = −
625
144 x

6 + 1258 x
5 −

2125
96 x

4 + 1258 x
3 −

137
24 x

2 + x
ηp(x) =

3125
144 x

6 −
875
12 x

5 + 887596 x
4 −

1925
36 x

3 + 252 x
2

ηq(x) = −
3125

72 x
6 + 162512 x

5 −
7375
48 x

4 + 267536 x
3 −

25
2 x

2

ηr (x) =
3125
72 x

6 − 125x5 + 612548 x
4 −

325
6 x

3 + 253 x
2

ηs(x) = −
3125
144 x

6 + 137524 x
5 −

5125
96 x

4 + 152572 x
3 −

25
8 x

2

η1(x) =
625
144 x

6 −
125
12 x

5 + 87596 x
4 −

125
36 x

3 + 12 x
2


(13)

On the evaluation of (12) at x = p, q, r, s, 1, we obtain the discrete accuracy-preserving BHA whose coefficients are presented in
Table 1.

Therefore, the accuracy-preserving BHA is given explicitly as,

yn+ 15 = yn +
1
5 hy

′

n + h
2
(

1231
126000 fn +

863
50400 fn+ 15 −

761
63000 fn+ 25 +

941
12600 fn+ 35 −

341
126000 fn+ 45 +

107
252000 fn+1

)
yn+ 25 = yn +

2
5 hy

′

n + h
2
(

71
3150 fn +

544
7875 fn+ 15 −

37
1575 fn+ 25 +

136
7875 fn+ 35 −

101
15750 fn+ 45 +

8
7875 fn+1

)
yn+ 35 = yn +

3
5 hy

′

n + h
2
(

123
3500 fn +

3501
28000 fn+ 15 −

9
3500 fn+ 25 +

87
2800 fn+ 35 −

9
875 fn+ 45 +

9
5600 fn+1

)
yn+ 45 = yn +

4
5 hy

′

n + h
2
(

376
7875 fn +

1424
7875 fn+ 15 +

176
7875 fn+ 25 +

608
7875 fn+ 35 −

16
1575 fn+ 45 +

16
7875 fn+1

)
yn+1 = yn + hy

′

n + h
2
(

61
1008 fn +

475
2016 fn+ 15 +

25
504 fn+ 25 +

125
1008 fn+ 35 +

25
1008 fn+ 45 +

11
2016 fn+1

)


(14)

y
′

n+ 15
= y

′

n + h
(

19
288 fn +

1427
7200 fn+ 15 −

133
1200 fn+ 25 +

241
3600 fn+ 35 −

173
7200 fn+ 45 +

3
800 fn+1

)
y
′

n+ 25
= y

′

n + h
(

14
225 fn +

43
150 fn+ 15 +

7
255 fn+ 25 +

7
255 fn+ 35 −

1
75 fn+ 45 +

1
450 fn+1

)
y
′

n+ 35
= y

′

n + h
(

51
800 fn +

219
800 fn+ 15 +

57
400 fn+ 25 +

57
400 fn+ 35 −

21
800 fn+ 45 +

3
800 fn+1

)
y
′

n+ 45
= y

′

n + h
(

14
225 fn +

64
255 fn+ 15 +

8
75 fn+ 25 +

64
255 fn+ 35 +

14
255 fn+ 45 + (0) fn+1

)
y
′

n+1 = y
′

n + h
(

19
288 fn +

25
96 fn+ 15 +

25
144 fn+ 25 +

25
144 fn+ 35 +

25
96 fn+ 45 +

19
288 fn+1

)


(15)

Equations (14) and (15) together form the proposed BHA.
3



Sunday et al. / J. Nig. Soc. Phys. Sci. 5 (2023) 1017 4

Table 1. Coefficients of the BHA

yn y′n fn fn+p fn+q fn+r fn+s fn+1
yn+p 1

1
5

1231
126000

863
50400

−761
63000

941
12600

−341
126000

107
252000

yn+q 1
2
5

71
3150

544
7875

−37
1575

136
7875

−101
15750

8
7875

yn+r 1
3
5

123
3500

3501
28000

−9
3500

87
2800

−9
875

9
5600

yn+s 1
4
5

376
7875

1424
7875

176
7875

608
7875

−16
1575

16
7875

yn+1 1 1
61

1008
475

2016
25
504

125
1008

25
1008

11
2016

y′n+p 0 1
19
288

1427
7200

−133
1200

241
3600

−173
7200

3
800

y′n+q 0 1
14
225

43
150

7
255

7
255

−1
75

1
450

y′n+r 0 1
51
800

219
800

57
400

57
400

−21
800

3
800

y′n+s 0 1
14
225

64
255

8
75

64
255

14
55 0

y′n+1 0 1
19
288

25
96

25
144

25
144

25
96

19
288

2.2. Implementation of the BHA

Unlike the conventional predictor-corrector methods that require starting values, the proposed BHA is self-starting. It is imple-
mented in block mode over a non- overlapping one-step interval [tn, tn+1], n = 0, 1, 2, ..., N−1. Firstly, we set the data inputs namely
initial conditions, step size and differential equation to be solved. Thus at n = 0, the block method is applied simultaneously on the
subinterval [t0, t1] to get

[
yκi, y

′

κi

]T
, i = 1, 2, ....5. Applying the values of the known previous block, we get the values of the next

subinterval [t1, t2] and the process goes on in this manner until the final subinterval [tN−1, tN ] is obtained. The code was written in
MATLAB 2021a.

3. Properties of the BHA

The properties of the proposed BHA shall be analysed in this section.

3.1. Order and Error Constant of the BHA

Definition 3.1
The linear difference operators associated with the BHA in equations (14) and (15) are defined as

`p {y(tn); h} = y (tn + ph) −
(
αr (t)y (tn + rh) + αs(t)y (tn + sh) + h2

∑1
j=0

(
β j(t) fn+ j + βk(t) fn+k

))
, k = p, q, r, s

`q {y(tn); h} = y (tn + qh) −
(
αr (t)y (tn + rh) + αs(t)y (tn + sh) + h2

∑1
j=0

(
β j(t) fn+ j + βk(t) fn+k

))
, k = p, q, r, s

`r {y(tn); h} = y (tn + rh) −
(
αr (t)y (tn + rh) + αs(t)y (tn + sh) + h2

∑1
j=0

(
β j(t) fn+ j + βk(t) fn+k

))
, k = p, q, r, s

`s {y(tn); h} = y (tn + sh) −
(
αr (t)y (tn + rh) + αs(t)y (tn + sh) + h2

∑1
j=0

(
β j(t) fn+ j + βk(t) fn+k

))
, k = p, q, r, s

`1 {y(tn); h} = y (tn + h) −
(
αr (t)y (tn + rh) + αs(t)y (tn + sh) + h2

∑1
j=0

(
β j(t) fn+ j + βk(t) fn+k

))
, k = p, q, r, s


(16)

`
′

p {y(tn); h} = hy
′ (tn + ph) −

(
α
′

r (t)y (tn + rh) + α
′

s(t)y (tn + sh) + h
2 ∑1

j=0

(
β
′

j(t) fn+ j + β
′

k(t) fn+k
))
, k = p, q, r, s

`
′

q {y(tn); h} = hy
′ (tn + qh) −

(
α
′

r (t)y (tn + rh) + α
′

s(t)y (tn + sh) + h
2 ∑1

j=0

(
β
′

j(t) fn+ j + β
′

k(t) fn+k
))
, k = p, q, r, s

`
′

r {y(tn); h} = hy
′ (tn + rh) −

(
α
′

r (t)y (tn + rh) + α
′

s(t)y (tn + sh) + h
2 ∑1

j=0

(
β
′

j(t) fn+ j + β
′

k(t) fn+k
))
, k = p, q, r, s

`
′

s {y(tn); h} = hy
′ (tn + sh) −

(
α
′

r (t)y (tn + rh) + α
′

s(t)y (tn + sh) + h
2 ∑1

j=0

(
β
′

j(t) fn+ j + β
′

k(t) fn+k
))
, k = p, q, r, s

`
′

1 {y(tn); h} = hy
′ (tn + h) −

(
α
′

r (t)y (tn + rh) + α
′

s(t)y (tn + sh) + h
2 ∑1

j=0

(
β
′

j(t) fn+ j + β
′

k(t) fn+k
))
, k = p, q, r, s


(17)

Suppose y(t)is differentiable, the terms ` and `′ are expanded as Taylor series about the point tn to get

`k {y(tn); h} = ckhp+2y(p+2) (tn) + O
(
hp+3

)
, k = p, q, r, s

`
′

k {y(tn); h} = c
′

kh
p+2y( p+2) (tn) + O

(
hp+3

)
, k = p, q, r, s

 (18)
where p in equation (18) denotes the order of the BHA and ` and `′ the error constants.

4



Sunday et al. / J. Nig. Soc. Phys. Sci. 5 (2023) 1017 5

Table 2. Order and error constant of the BHA

y(tn+ j) Order Error Constant
j = 1/5 6 −2.1058 × 10−8

j = 2/5 6 −5.1471 × 10−8

j = 3/5 6 −8.0571 × 10−8

j = 4/5 6 −1.0836 × 10−7

j=1 6 −1.4550 × 10−7

Therefore, the proposed BHA is of order 6.

Definition 3.2 [24]
A linear multistep method for a second order problem is said to be of order pif it satisfies c0=c1=c2=. . . =cp=cp+1=0, cp+2 , 0,

where,

c0 =
∑k

j=0 α j

c1 =
∑k

j=0

(
jα j −β j

)
.
.
.

cp =
∑k

j=0

[
1
p! j

pα j −
1

(p−1)! j
p−1β j

]
, p = 2, 3, ..., q + 1


(19)

The parameter cp+2 , 0 is referred to as the error constant of the method with the local truncation error defined as tn+k =
cp+2hp+2y(p+2)(tn) + O

(
hp+3

)
”.

The BHA is expanded in Taylor series about the point tn to obtain



∑
∞
j=0

( 15 h)
j

j! y
j
n − yn −

1
5 hy

′

n −
1231

126000 h
2y
′′

n −
∑
∞
j=0

h j+2
j! y

j+2
n

{
863

50400

(
1
5

) j
−

761
63000

(
2
5

) j
+ 94112600

(
3
5

) j
−

341
126000

(
4
5

) j
+ 107252000 (1)

j
}

∑
∞
j=0

( 25 h)
j

j! y
j
n − yn −

2
5 hy

′

n −
71

3150 h
2y
′′

n −
∑
∞
j=0

h j+2
j! y

j+2
n

{
544

7875

(
1
5

) j
−

37
1575

(
2
5

) j
+ 1367875

(
3
5

) j
−

101
15750

(
4
5

) j
+ 87875 (1)

j
}

∑
∞
j=0

( 35 h)
j

j! y
j
n − yn −

3
5 hy

′

n −
123
3500 h

2y
′′

n −
∑
∞
j=0

h j+2
j! y

j+2
n

{
3501

28000

(
1
5

) j
−

9
3500

(
2
5

) j
+ 872800

(
3
5

) j
−

9
875

(
4
5

) j
+ 95600 (1)

j
}

∑
∞
j=0

( 45 h)
j

j! y
j
n − yn −

4
5 hy

′

n −
376
7875 h

2y
′′

n −
∑
∞
j=0

h j+2
j! y

j+2
n

{
1424
7875

(
1
5

) j
+ 1767875

(
2
5

) j
+ 6087875

(
3
5

) j
−

16
1575

(
4
5

) j
+ 167875 (1)

j
}

∑
∞
j=0

(h) j

j! y
j
n − yn − hy

′

n −
61

1008 h
2y
′′

n −
∑
∞
j=0

h j+2
j! y

j+2
n

{
475

2016

(
1
5

) j
+ 25504

(
2
5

) j
+ 1251008

(
3
5

) j
+ 251008

(
4
5

) j
+ 112016 (1)

j
}


=



0

0

0

0

0


(20)

Table 2 displays the order and error constant of the proposed BHA.

3.2. Consistency of the BHA

Definition 3.3 [24]
“If a method has orderp ≥ 1, then it is consistent”.
Thus, the proposed BHA is consistent since it is of order 6.

3.3. Zero-Stability of the BHA

Definition 3.4 [25]
“If no root of the characteristic polynomial has a modulus greater than one and every root with modulus one is simple, then

such a method is called zero-stable”.
In order words, the BHA is zero-stable if its first characteristic polynomial ρ(z)satisfies |zd| ≤ 1, d = 1, 2, 3, ..., n. Therefore,

5



Sunday et al. / J. Nig. Soc. Phys. Sci. 5 (2023) 1017 6

Figure 1. Region of absolute stability of the BHA

ρ(z) =

∣∣∣∣∣∣∣∣∣∣∣∣∣∣


z 0 0 0 0
0 z 0 0 0
0 0 z 0 0
0 0 0 z 0
0 0 0 0 z

 −


0 0 0 0 1
0 0 0 0 1
0 0 0 0 1
0 0 0 0 1
0 0 0 0 1



∣∣∣∣∣∣∣∣∣∣∣∣∣∣
= z4(z − 1) = 0

(21)

On the evaluation of equation (21), we obtain z = 0, 0, 0, 0, 1.
Therefore, the proposed BHA is zero-stable.

3.4. Convergence of the BHA

Theorem 3.1 [25]
“For a method to be convergent, it is necessary and suffi-

cient that it be consistent and zero-stable.”
The BHA is convergent since it is consistent and zero-stable.

3.5. Region of Absolute Stability of the BHA

Definition 3.5 [26]
“The region of absolute stability of a method is the region

in the complex z plane, where z = λh.”
In such region, the approximate solution y′′ = −λ2y satisfies

yi → 0 as i → 0 for any initial condition.
Definition 3.6 [26]
“If the region of absolute stability of method contains the

whole of the left half-plane, that is Re (hλ) < 0, then the method
is said to be A-stable.”

Using the boundary locus method, the stability polynomial
for the proposed BHA is

h(w) = −h10
((

8.13 × 10−10
)

w5 +
(
1.00 × 10−8

)
w4

)
− h8

((
5.02 × 10−8

)
w5 +

(
5.04 × 10−6

)
w4

)
− h6

((
1.32 × 10−6

)
w5 +

(
7.18 × 10−4

)
w4

)
− h4

((
3.71 × 10−2

)
w4 −

(
2.00 × 10−4

)
w5

)
− h2

((
2.00 × 10−2

)
w5 +

(
6.27 × 10−1

)
w4

)
+ w5 − 2w4 (22)

For more details on boundary locus method, see [27-28]. The
region of absolute stability of the BHA is shown in Figure 1.

Note that the absolute stability of the proposed BHA in Fig-
ure 1 is the region that lies inside the closed contour. From the
figure, it is clear that the BHA is A-stable.

4. Numerical Examples

The newly derived BHA shall be applied in integrating some
second-order physical systems with oscillatory solutions. The
following abbreviations shall be used in the Tables 3-5 and Fig-
ures 3-10.

t: Point of evaluation
h: Step-size
NS: Number of steps taken
Err: Maximum global error
ErrH: Maximum global error for the Hamiltonian
ode 45: MATLAB inbuilt solver
VSSHM: Variable step-size hybrid method of order seven

by [13]
BHM: Order eleven block hybrid method by [2]
IBM: Order eight implicit block method by [11]
BHA: Newly proposed sixth order block hybrid algorithm
Let y(tn) and yn be the exact and approximate solutions re-

spectively, then the maximum global error is computed as Err =
max |y(tn) − yn| and the maximum global error for the Hamilto-
nian as ErrH = max |Hn − H0|.

Example 1. The system below describes the perturbed Ke-
pler problem

y
′′

1 (t) = −
y1 (t)

(y21 (t)+y22 (t))
3/2 −

(e2 +2e)y1 (t)
(y21 (t)+y22 (t))

5/2 , y1(0) = 1, y
′

1(0) = 0

y
′′

2 (t) = −
y2 (t)

(y21 (t)+y22 (t))
3/2 −

(e2 +2e)y2 (t)
(y21 (t)+y22 (t))

5/2 , y2(0) = 0, y
′

2(0) = 1 + e


(23)

whose exact solution is given by

y1(t) = cos
(
et + t

)
y2(t) = sin

(
et + t

) } (24)
where e is the eccentricity taken as e = 10−3. The problem (23)
has an angular momentum

L = y1(t)y
′

2(t) − y2(t)y
′

1(t) (25)

as a first integral and also has the Hamiltonian

H =
1
2

(
y
′ 2
1 (t) + y

′ 2
2 (t)

)
−

1(
y21(t) + y

2
2(t)

)1/2 −
(
e2 + 2e

)
3
(
y21(t) + y

2
2(t)

)3/2
(26)

The author in [2] developed an eleventh order block hybrid
method to compute the maximum global errors in the solu-
tion as well as the maximum global error in the Hamiltonian
of Example 1 at different step-sizes. The newly derived BHA
was also used to solve the same problem. The results obtained
clearly showed that even though the BHA is of the sixth or-
der, it performed better than the eleventh order block hybrid
method proposed by [2], see Table 3. Figure 2 shows the plots
for mean anomaly against eccentric anomaly at different eccen-
tricities (e=0.01, 0.1, 0.2, 0.4, 0.6 and 0.9) using the newly de-
rived BHA. On the other hand, Figure 3 shows the exact plots
for mean anomaly against eccentric anomaly at the same ec-
centricities. The plots generated using the BHA (see Figure 2)

6



Sunday et al. / J. Nig. Soc. Phys. Sci. 5 (2023) 1017 7

Figure 2. Mean anomaly against eccentric anomaly (for Example 1) at different
values of eccentricity using the proposed BHA

Figure 3. Exact plots for mean anomaly against eccentric anomaly (for Example
1) at different values of eccentricity

Figure 4. Showing the trajectories in the phase plane flow (for Example 1) using
the proposed BHA at N=160

converges to those of the exact plots in Figure 3 at the respec-
tive eccentricities. It is important to state that the Figures 2
and 3 show the relationship between mean anomaly (that is, the
angle between line drawn from the sun to perihelion and to a
point moving in the orbit) and eccentric anomaly (that is, the
angle that describes the position of a body that is moving along
an elliptic orbit). Furthermore, Figures 4, 5 and 6 respectively
depict the trajectories in the phase plane flow using proposed
BHA, BHA developed by [2] and the exact flow N = 160. For
more details on Kepler equations, see [29-35].

Example 2. The nonlinear system below describes the aero-
dynamic damping of a pendulum

θ′′ +
g
L

sin θ +
ca
mL

(
θ′
)2 sgn (θ′) = 0 (27)

and the linear system below describes the viscous damping of a
pendulum

θ′′ +
cd
m

(
θ′
)2

+
g
L
θ = 0 (28)

To determine the quadratic damping in equations (27) and (28)
respectively, the proposed BHA is employed in simulating the

Figure 5. Showing the trajectories in the phase plane flow (for Example 1) using
BHM at N=160

Figure 6. Showing the exact trajectories in the phase plane flow (for Example
1) at N=160

two equations and the simulation results are compared with that
of author [1]. The author transformed (27) and (28) to their
equivalent first order systems and then applied MATLAB solver
(ode 45). Let g = 981cms−2, ca = cd = 14, L = 5cm, m = 10g,
θ(0) = 1.57and θ′(0) = 0. The initial value θ is in radians. Also,
sgn in equation (27) represents signum function; ca and cd are
the coefficients of aerodynamic damping and viscous damping
respectively. Figure 7 shows a simple pendulum with two forces
acting on the mass; that is, the weight mg and the tension T . The
net force is found to be F = mg sin θ.

From simulation results presented in Figure 8, it is obvious
that quadratic damping causes a rapid initial amplitude reduc-
tion but is less effective as time increases and amplitude de-
creases. On the other hand, linear damping is less effective ini-
tially but is more effective as time increases. The plots confirm
that the proposed BHA is computationally reliable as they agree
with those of ode 45.

Example 3. Consider the nonlinear oscillatory system

y′′(t) +
[

13 − 12
−12 13

]
y(t) =

∂U
∂y

, y(0) =
[
−1
1

]
,

y′(0) =
[
−5
5

]
(29)

where U(y) = y1(t)y2(t) (y1(t) + y2(t))
3and the exact solution

7



Sunday et al. / J. Nig. Soc. Phys. Sci. 5 (2023) 1017 8

Table 3. Juxtaposition of maximum global error for Example 1 on [0, 100]
h BHA-Err BHM-Err BHA-ErrH BHM-ErrH

1/2 4.39 × 10−15 1.74 × 10−13 2.62 × 10−17 5.27 × 10−15

1/4 3.46 × 10−14 8.62 × 10−13 5.89 × 10−17 1.23 × 10−14

1/8 5.09 × 10−14 3.49 × 10−12 9.22 × 10−17 2.11 × 10−14

1/16 8.13 × 10−14 6.94 × 10−12 2.16 × 10−16 4.19 × 10−14

1/32 1.67 × 10−13 1.92 × 10−12 6.65 × 10−16 5.06 × 10−14

Figure 7. A simple pendulum

Figure 8. Solution curves for Example 2 depicting quadratic and linear damping

given by

y(t) =
[
−sin 5t − cos 5t
sin 5t + cos 5t

]
(30)

The system (29) has the Hamiltonian

H (y(t) + y′(t)) = 12 y
′T y′ + 12 y

T ,[
13 − 12
−12 13

]
y(t) + U(y)

 (31)
The eleventh order block hybrid method formulated by [2] was
employed in computing the maximum global errors in the solu-
tion as well as the maximum global error in the Hamiltonian for
Example 3 at different step sizes. The proposed BHA was also
used to solve the same example. The results obtained clearly
showed that even though the BHA is of the sixth order, it out-

Figure 9. Solution curves for Example 4 showing the Xand Y displacements

performed the block hybrid method of order eleven derived by
[2], see Table 4.

Example 4. The following nonlinear system describes the
trajectory motion of an object in terms of X and Y displace-
ments

mx′′(t) + cd x
′(t)

((
x′(t)

)2
+

(
y′(t)

)2)(1/2)
= 0 (32)

my′′(t) + cd y
′(t)

((
x′(t)

)2
+

(
y′(t)

)2)(1/2)
= −w (33)

The system was solved using the following parameters of aero-
dynamic drag cd = {0.1, 0.3, 0.5} N/m. Let x(0) = 0, y(0) =
0,x′(0) = 100ms−1, y′(0) = 10ms−1,m = 10kg,w = mg and
g = 9.81ms−2.

The proposed BHA was directly employed in simulating the
system and the results compared with that of the author in [36]
who applied ode 45 solver. It is obvious from the plots in Fig-
ure 9 that as cd is increases, the maximum value of the displace-
ment in the ydirection decreases. The value of the displacement
in the x direction, for which the displacement in the y direction
is maximum, also decreases. The curves generated via the pro-
posed BHA converge to those generated using the ode45 solver,
implying that the BHA is computationally robust and accurate.

Example 5. Consider the Bessel’s nonlinear oscillatory
problem

t2y′′(t) + ty′(t) +
(
t2 − 0.25

)
y(t) = 0, y(1) =

√
2
π

sin(1),

y′(1) =
2 cos(1) − sin(1)

√
2π

(34)

whose exact solution is given by

y(t) =

√
2
πt

sin(t) (35)

The newly derived sixth order BHA was applied on Example 5
in the interval [1, 8]. The results obtained were compared with

8



Sunday et al. / J. Nig. Soc. Phys. Sci. 5 (2023) 1017 9

Table 4. Juxtaposition of maximum global error for Example 3 on [0, 100]
h BHA-Err BHM-Err BHA-ErrH BHM-ErrH

1/2 1.98 × 100 4.58 × 100 1.04 × 103 2.95 × 103

1/4 8.41 × 10−10 7.54 × 10−8 5.24 × 10−10 3.13 × 10−8

1/8 3.55 × 10−12 2.20 × 10−11 3.33 × 10−13 5.61 × 10−12

1/16 2.23 × 10−14 4.95 × 10−14 6.19 × 10−13 1.62 × 10−12

1/32 5.70 × 10−15 9.15 × 10−14 4.89 × 10−14 3.78 × 10−12

Table 5. Juxtaposition of maximum global error for Example 5 on [1, 8]
NS BHA-Err IBM-Err VSSHM-Err
67 1.1984 × 10−17 7.7716 × 10−16 6.5286 × 10−11

82 6.1189 × 10−18 1.8874 × 10−15 1.3679 × 10−11

112 2.8346 × 10−18 1.1601 × 10−15 1.1897 × 10−12

Figure 10. Accuracy curves for Example 5

those of order eight implicit block method (IBM) by [11] and
order seven variable step-size hybrid method (VSSHM) devel-
oped by [13] using the same number of steps. The results of
the new sixth order BHA performed better than the other two
methods, see Table 5. Figure 10 also showed that the BHA is
more accurate than the IBM and VSSHM.

For more details on oscillatory problems, refer to [37-45].

5. Conclusion

An attempt has been made in this research to derive an
accuracy-preserving BHA for the solution of second-order phys-
ical systems with oscillatory solutions. Problems considered
included Kepler, Bessel and damped systems. The results gen-
erated showed that the BHA is accuracy-preserving and com-
putationally reliable. Simulation results obtained also showed
that the new method was in agreement with the MATLAB in-
built solver (ode 45). Properties such as order, error constant,
zero-stability, stability, consistence and convergence of the new
BHA were also validated. Finally, the proposed BHA is recom-
mended for the integration of physical systems with oscillatory
solutions. Further research will focus on the solution of physi-
cal models for higher order partial differential equations.

Acknowledgements

We would like to appreciate the University of Jos for pro-
viding the facilities for conducting this research. We also appre-
ciate the Reviewers and Editors for all the valuable comments

and suggestions, which helped us to improve the quality of the
manuscript.

References

[1] H. Soraya, “Pendulum with aerodynamic and viscous damping”, Journal
of Applied Information, Communication and Technology 3 (2016) 43.

[2] S. N. Jator & K. L. King, “Integrating oscillatory general second-order
initial value problems using a block hybrid method of order 11”, Mathe-
matical Problems in Engineering (2018) 3750274.

[3] D. V. V. Wend, “Uniqueness of solution of ordinary differential equa-
tions”, The American Mathematical Monthly 74 (1967) 948.

[4] L. Brugnano and D. Trigiante, Solving ODEs by Multistep Initial and
Boundary Value Methods, Gordon Breach, Amsterdam, (1998).

[5] M. I. Modebei, O. O. Olaiya & I. P. Ngwongwo, “Computational study
of some three-step hybrid integrators for solution of third order ordinary
differential equations”, Journal of Nigerian Society of Physical Sciences
3 (2021) 459.

[6] F. O. Obarhua & O. J. Adegboro, “Order four continuous numerical
method for solving general second order ordinary differential equations”,
Journal of Nigerian Society of Physical Sciences 3 (2021) 42.

[7] E. Hairer, “A one-step method of order 10 for y′′ = f (x, y)”, IMA Journal
of Numerical Analysis 2 (1982) 83.

[8] R. D’Ambrosio, M. Ferro & B. Paternoster, “Two-step hybrid collocation
methods fory′′ = f (x, y)”, Applied Mathematics Letters 22 (2009) 1076.

[9] S. D. Yakubu, Y. A. Yahaya & K. O. Lawal,”3-point block hybrid linear
multistep methods for solution of special second order ordinary differen-
tial equations”, Journal of the Nigerian Mathematical Society 40 (2021)
149.

[10] F. O. Obarhua & S. J. Kayode, “Continuous explicit hybrid method for
solving second rder ordinary differential equations”, Pure and Applied
Mathematics Journal 9 (2020) 26.

[11] R. Allogmany & F. Ismail, “Direct solution of using three point
block method of order eight with applications”, Journal of King Saud
University-Science 33 (2021) 101337.

[12] B. Y. Guo & J. P. Yan, “Legendre-Gauss collocation method for initial
value problems of second order ordinary differential equations”, Applied
Numerical Mathematics 59 (2009) 1386.

[13] S. N. Jator, “Solving second order initial value problems by hybrid
method without predictor”, Appl. Math. Comput. 217 (2010) 4036.

[14] H. Ramos, Z. Kalogiratou, T. Monovasilis & T. E. Simos, “A trigonomet-
rically fitted optimized two-step hybrid block method for solving initial-
value problems of the form y′′ = f (x, y, y′) with oscillatory solutions”,
Proceedings of the International Conference on Numerical Analysis and
Applied Mathematics 2014.

[15] F. F. Ngwane & S. N. Jator, “Block hybrid method using trigonometric
basis for initial value problems with oscillating solutions”, Numerical Al-
gorithms 63 (2013) 713.

[16] B. Wang, A. Iserles & X. Wu, “Arbitrary-Order Trigonometric Fourier
Collocation Methods for Multi-Frequency Oscillatory Systems”, Foun-
dations of Computational Mathematics 16 (2016) 158.

[17] S. N. Jator & H. B. Oladejo, “Block Nystrom method for singular differ-
ential equations of the Lane-Emden type and problems with highly os-
cillatory solutions”, International Journal of Applied and Computational
Mathematics 3 (2017) 1385.

[18] J. Sunday, “Optimized two-step second derivative methods for the so-
lutions of stiff ystems”, Journal of Physics Communications 6 (2022)
055016.

9



Sunday et al. / J. Nig. Soc. Phys. Sci. 5 (2023) 1017 10

[19] H. Ramos, S. Mehta & J. Vigo-Aguiar, “A unified approach for the devel-
opment of k-step block Falkner-type methods for solving general second
order initial value problems in ODEs”, Journal of Computational and Ap-
plied Mathematics 318 (2017) 550.

[20] O. Adeyeye & O. Zurmi, “Maximal order block method for the solution of
second order ordinary differential equations”, IAENG International Jour-
nal of Applied Mathematics 46 (2016) 4.

[21] J. Sunday, G. M. Kumleng, N. M. Kamoh, J. A. Kwanamu, Y. Skwame
& O. Sarjiyus, “Implicit four-point hybrid block integrator for the simu-
lations of stiff models”, Journal of the Nigerian Society of Physical Sci-
ences 4 (2022) 287.

[22] V. J. Shaalini & S. E. Fadugba, “A new multistep method for solving delay
differential equations using Lagrange interpolation”, Journal of Nigerian
Society of Physical Sciences 3 (2021) 159.

[23] M. Kida, S. Adamu, O. O. Aduroja & T. P. Pantuvo, “Numerical solution
of stiff and oscillatory problems using third derivative trigonometrically
fitted block method”, Journal of Nigerian Society of Physical Sciences 4
(2022) 34.

[24] S. O. Fatunla, “Numerical integrators for stiff and highly oscillatory dif-
ferential equations”, Mathematics of Computation 34 (1980) 373.

[25] J. D. Lambert, Numerical methods for ordinary differential systems: The
initial value problem, John Wiley & Sons LTD, United Kingdom (1991).

[26] J. D. Lambert, Computational methods in ordinary differential equations,
John Wiley & Sons, Inc., New York (1973).

[27] L. J. Kwari, J. Sunday, J. N. Ndam & A. A. James, “On the numerical
approximation and simulation of damped and undamped Duffing oscilla-
tors”, Science Forum (Journal of Pure and Applied Sciences) 21 (2021)
503.

[28] M. R. Odekunle, A. O. Adesanya & J. Sunday, “4-Point Block Method
for the Direct Integration of First-Order Ordinary Differential Equations”,
International Journal of Engineering Research and Applications 5 (2012)
1182.

[29] R. H. Ibrahim & A. H. Saleh, “Re-evaluation solution methods for Ke-
pler’s equation of an elliptical orbit”, Iraqi Journal of Science 60 (2019)
2269.

[30] R. I. Abdulganiy, H. Ramos, A. O. Akinfenwa & S. A. Okunuga, “A
functionally-fitted block numerov method for solving second-order initial
value problems with oscillatory solutions”, Mediterr. J. Math 18 (2021)
259.

[31] S. O. Fatunla, “Numerical integrators for stiff and highly oscillatory dif-
ferential equations”, Math Comput. 34 (1980) 373.

[32] S. Faydaoglu & T. Ozis, “Periodic solutions for certain non-smooth oscil-
lators with high nonlinearities” Appl. Comput. Math. 20 (2021) 366.

[33] F. Bakhtiari-Nejad & M. Nazari, “Nonlinear vibration analysis of
isotropic cantilever plate with visco-elastic laminate”, Nonlinear Dynam-
ics 56 (2009) 325.

[34] J. Sunday, A. Shokri & D. Marian, “Variable Step Hybrid Block Method
for the Approximation of Kepler Problem”, Fractal and Fractional 6
(2022) 343.

[35] J. Guckenheimer & P. Holmes, Nonlinear oscillations, dynamical systems
and bifurcations of vector fields, Springer-Verlag, (1983).

[36] P. Harihara & N. D. Childs, “Solving problems in dynamics and vibra-
tions using MATLAB”, Department of Mechanical Engineering, Texas A
and M University College Station (2020) 1.

[37] B. S. H. Kashkari & M. I. Syam, “Optimization of one-step block method
with three hybrid points for solving first-order ordinary differential equa-
tions”, Results in Physics 12 (2019) 592.

[38] J. Sunday, A. Shokri, R. O. Akinola, K. V. Joshua & K. Nonlaopon, “A
convergence-preserving non-standard finite difference scheme for the so-
lutions of singular Lane-Emden equations”, Results in Physics 42 (2022)
106031.

[39] S. S. Shyam, B. Omar, A. Hijaz & C. Yu-Ming, “Second order differential
equation: oscillation theorems and applications”, Mathematical Problems
in Engineering (2020) 1.

[40] L. Tongxing, V. R. Yuriy & T. Shuhong, “Oscillation of second-order
nonlinear differential equations with damping”, Mathematica Slovaca 64
(2014) 1227.

[41] J. Sunday, A. Shokri, J. A. Kwanamu & K. Nonlaopon, “Numerical Inte-
gration of Stiff Differential Systems using Non-Fixed Step-Size Strategy”,
Symmetry 14 (2022) 1575.

[42] F. J. Agocs, W. J. Handley, A. N. Lasenby & M. P. Hobson, “Efficient
method for solving highly oscillatory ordinary differential equations with
applications to physical systems”, Physical Review Research 2 (2020) 1.

[43] P. F. Rowat & A. I. Selverston, “Modeling the gastric mill central pattern
generator of the lobster with a relaxation-oscillator network”, Journal of
Neurophysiology 70 (1993) 1030.

[44] S. Nourazar & A. Mirzabeigy, “Approximate solution for nonlinear Duff-
ing oscillator with damping effect using the modified differential trans-
form method”, Scientia Iranica B 20 (2013) 364.

[45] J. Sunday, C. Chigozie, E. O. Omole & J. B. Gwong, “A pair of three-step
hybrid block methods for the solutions of linear and nonlinear first-order
systems”, European Journal of Mathematics and Statistics 3 (2022) 14.

10