J. Nig. Soc. Phys. Sci. 2 (2020) 12–25

Journal of the
Nigerian Society

of Physical
Sciences

Original Research

Block Third Derivative Trigonometrically-Fitted Methods for
Stiff and Periodic Problems

R. I. Abdulganiya,∗, O. A. Akinfenwab, O. A. Yusuffb, O. E. Enobaborc, S. A. Okunugab

aDistance Learning institute, University of Lagos, Nigeria
bDepartment of Mathematics, University of Lagos, Nigeria

cDepartment of Mathematics, Yaba College of Technology, Lagos, Nigeria

Abstract

This article constructed and implemented a family of a third derivative trigonometric fitted method of order k + 3 whose coefficients are functions
of frequency and step size for the integration of systems of first-order stiff and periodic Initial Value Problems. The Block Third Derivative
Trigonometric Fitted methods (BTDTFMs) are constructed via multistep collocation technique and applied in block form as simultaneous numer-
ical integrators which make them self-starting. The basic properties of the BTDTFMs are analyzed and presented. The accuracy and efficiency of
the methods based on number of function evaluation are established through some standard numerical examples which are either stiff or periodic
in nature.

Keywords: Convergence, Frequency, Stiff, Trigonometrically-Fitted.

Article History :
Received: 02 January 2020
Received in revised form: 16 February 2020
Accepted for publication: 17 February 2020
Published: 28 February 2020

c©2020 Journal of the Nigerian Society of Physical Sciences. All rights reserved.
Communicated by: T. Latunde

1. Introduction

A-stable methods are of great importance in solving stiff
problems but put stark limitations on the choice of suitable
methods for stiff problems. Dahlquist [1] established that the
most accurate A-stable linear multistep method has order 2.
The search for higher-order A-stable multistep method is car-
ried out in the following three major directions (Biala et al. [2]):

1. Use of higher derivatives of the solution;

2. Incorporating additional stages and hybrid points;

∗Corresponding author tel. no: +2348035195870
Email addresses: profabdulcalculus@gmail.com (R. I. Abdulganiy),

akinolu35@yahoo.com (O. A. Akinfenwa),
olaoluwa.yusuff@gmail.com (O. A. Yusuff),
enobaborosaretin@yahoo.com (O. E. Enobabor),
waleokunuga@gmail.com (S. A. Okunuga)

3. Incorporating super future points.

Many numerical methods have been developed along these
three directions which include but not limited to Gear [3], En-
right [4], Cash [5], Wu [6], Hojjati [7], Jator [8], Sahi et al. [9],
Ngwane and Jator [10], Akinfenwa et al. [11,12,13], Mehdizadeh
et al [14] and Abdulganiy et al. [15].

Method based on Simpson’s was explored by Akinfenwa et
al. [16] while methods bases on exponential fitting were consid-
ered in Okunuga [17], Vaquero and Vigo-Aguiar [18], Abhuli-
men et al. [19,20], Ehigie et al. [21] and Adesanya et al. [22]
to numerically approximate stiff problems. Although the expo-
nential fitting methods for solving stiff problems were easy to
implement, tedious mathematical analyses are involved in ob-
taining the A-stability.

Periodic problems, on the other hand, have received much
attention in the past few decades. The methods for solving peri-

12



Abdulganiy et al. / J. Nig. Soc. Phys. Sci. 2 (2020) 12–25 13

odic problems according to Yakubu et al., [23] can be classified
into two: methods with constant coefficients and methods hav-
ing coefficients depending on the frequency of the problems.
If a good estimate of the frequency or some suitable substitute
is known in advance, a linear multistep method can be used
for the integration of periodic problems (Vanden Berghe et al.,
[24] and Van de Vyver [25]). The idea of using basis function
which integrates a set of linearly independent function exactly
other than polynomial can be traced to the work of Gautchi [26]
and Lyche [27]. Many of such extension has been discussed in
Coleman and Duxbury [28], Ixaru et al. [29], Vanden Berghe
and his collaborators ([24], [30,31]), Simos [32,33], Tsitouras
and Simos [34], Nguyen et al. [35], Senu et al. [36], Jator
and his collaborators [37-41], Jator [42] and Abdulganiy and
his collaborators [43-45].

In what follows, a Continuous third Derivative Method with
Trigonometric Coefficients (CTDMTC) is derived via a mul-
tistep collocation technique for which the approximate inter-
polating function is a linear combination of polynomial and
trigonometric terms. The CTDMTC is used to generate a family
of Third Derivative Trigonometrically-Fitted Method (TDTFM)
and some other discrete methods as by-products that are com-
bined together into a Block Third Derivative Trigonometrically-
Fitted Method (BTDTFM). The BTDTFM is applied as a simul-
taneous numerical integrator for the integration of the system of
first-order Ordinary Differential Equations (ODEs)

y
′

(t) = f (t, y (t)) , y (t0) = y0, t ∈ [t0, tN ] (1)

where f satisfies the Lipschitz condition with a periodic so-
lution whose frequency is known in advance. Stiff or/and peri-
odic differential equation of the form in equation (1) frequently
arises in many area of Science and Engineering, a list of such
are provided in Jator and Agyingi [46], Varden Berghe et al.
[24] and Ndukum et al. [47]. The methods in this paper is dif-
ferent from the methods in Jator [48,49] and Akinfenwa [50]
in that they are used to solve system of first-order ODEs while
the methods in Jator [48,49] are direct numerical integrators of
general second-order ODEs and the methods in Akinfenwa [50]
are hybrid in nature respectively.

Some of the advantages of continuous methods include but
not restricted to providing defect control (Enright [51]) and
ability to produce complementary methods which are combined
and applied in block-by-block fashion (Onumanyi et al. [52],
Akinfenwa et al. [11,13]). The foremost block methods are
credited to Milne [53], Rosser [54], Shampine and Watts [55]
and Chu and Hamilton [56] and are implemented in predictor-
corrector modes that reduce the stability properties of the meth-
ods. The block method in this paper is implemented without
the use of predictors. According to Jator and Agyingi [46], for
a block method, it is needless to make a function evaluation
at the preliminary part of the original block since at all blocks
(with exception of the first block) the first function evaluation
is previously known from the previous block.

The rest of this paper is arranged as follows: the theoretical
and basic elements of the BTDTFM is discussed in section 2.
The basic properties of the methods are discussed in section 3.

In section 4, the estimation of the computational frequency is
discussed while numerical experiments are presented in section
5. Finally, section 6 concludes the paper.

2. Theoretical and Basic Elements of the BTDTFM

The k−step third derivative method has the form

yn+k−yn+k−1 = h
k∑

j=0

β j (u) fn+ j+h
2δk (u) gn+k+h

3γk (u) ln+k(2)

where u = ωh ,ω is the frequency, β j,δk,γk which depends
on frequency and step size are parameters to be determined
uniquely. Also, yn+k is the numerical approximation to exact
solution y (xn+k) and fn+ j = f

(
xn+ j , yn+ j

)
,

gn+k =
d f (x, y (x))

d x
|
x = xn+k
y = yn+k

ln+k =
d2 f (x, y (x))

d x2
|
x = xn+k
y = yn+k

In order to obtain (2), the exact solution y (x) is approxi-
mated on the interval [xn , xn+k] by the interpolating function
I (x) given by

I (x) =
k+1∑
j=0

a j x
j
+ ak+2sin(ωx) + ak+3cos(ωx) (3)

Imposing conditions in equations (4) and equation (4) below on
equation (3)

I
(
xn+ j

)
= yn+ j (4)


I
′
(
xn+ j

)
= fn+ j j = 0(1)k

I
′′
(
xn+ j

)
= gn+ j j = k

I
′′′

(
xn+ j

)
= ln+ j j = k

(5)

lead to a system of (k + 4) equations which is solved using a
Computer Algebra System (CAS) such as Maple. The deriva-
tion of the methods manually becomes more difficult for k ≥ 2.
Hence, the explorations of CAS become necessary. The contin-
uous form is developed by substituting the values of a j, j =
0 (1) (k + 3) into equation equation (3). After some algebraic
manipulations, the continuous form is obtained as the expres-
sion equation (6).

I (x) = yn+k−1+h
k∑

j=0

β j (u, x) fn+ j+h
2δk(u, x)gn+k+h

3γk(u, x)ln+k(6)

The continuous method in equation (6) is evaluated at x = xn+k
to generate the principal method given by

13



Abdulganiy et al. / J. Nig. Soc. Phys. Sci. 2 (2020) 12–25 14

yn+k = yn+k−1 + h
k∑

j=0

β j (cos (u) , sin(u)) fn+ j+h
2δk (cos (u) , sin (u) ) gn+k + h

3γk(cos (u) , sin(u))ln+k (7)

while the (k − 1) secondary methods are obtained as

yn+i = yn+k−1 + h
k∑

j=0

β j(cos (u) , sin(u)) fn+ j+h
2δk(cos (u) , sin(u))gn+k + h

3γk(cos (u) , sin(u))ln+k (8)

by evaluating equation (6) at x = xn+i , i = 1, 2, . . . , k − 1.
Following the steps discussed above, the secondary and principal methods for k = 2 and k = 3, their coefficients and the corre-
sponding power series conversion up to O

(
u10

)
are given as follows For k = 2,

yn = yn+1 + hβ0,1 (cos(u) , sin(u) ) fn + hβ1,1 (cos(u) , sin(u) ) fn+1

+h2δ2,1 (cos(u) , sin(u) )gn+1 + h3γ2,1 (cos(u) , sin(u) )ln+1
(9)

yn+2 = yn+1 + hβ0 (cos(u) , sin(u) ) fn + hβ1 (cos(u) , sin(u) ) fn+1 + hβ2 (cos(u) , sin(u) ) fn+2

+h2δ2 (cos(u) , sin(u) ) gn+2 + h3γ2 (cos(u) , sin(u) ) ln+2
(10)



β0,1 =
1
6

(
−14sin (u) u3 + 5u4 − 12cos (u) u2 − 6cos (2u) u2 + 12sin (2u) u + 6u2 − 24cos (u) + 12cos (2u) + 12

)
/
((
−2u − 2cos (2u) u − sin (2u) u2 − 4sin (u) + 2sin (2u) + 4cos (u) u + 4sin (u) u2 − 2u3

)
u
)

β1,1 =
1
6

((
−7u3 + 18u

)
sin(2u) +

(
−21u2 + 6

)
cos (2u) − 4u4 + 12cos (u) u2 − 3u2 − 12sin (u) u − 6

)
/
(((

1
2 u

2 − 1
)

sin (2u) + u3 − 2sin (u) u2 − 2cos (u) u + cos (2u) u + u + 2sin(u)
)

u
)

β2,1 =
1
6

((
8u3 − 12u

)
sin(2u) + u4 + 10sin(u) u3 − 12cos (u) u2 + 24cos (2u) u2 − 24cos (u) + 24

)
/
((
−4sin (u) u2 + sin (2u) u2 + 2u3 − 4cos (u) u + 2cos (2u) u + 4sin (u) − 2sin(2u) + 2u

)
u
)

δ2,1 =
1
6

(
8sin (u) u3 + 5sin(2u) u3 − 12cos (u) u2 − 12cos (2u) u2 − 12sin (u) u + 6sin (2u) u − 48cos (u) + 12cos (2u) + 36

)
/
((
−2u − 2cos (2u) u − sin (2u) u2 − 4sin (u) + 2sin (2u) + 4cos (u) u + 4sin (u) u2 − 2u3

)
u
)

γ2,1 =
((
−5u2 + 18

)
cos (2u) − 8cos (u) u2 + u2 − 16sin (u) u + 20sin (2u) u − 48cos (u) + 30

)
/
(
12cos (2u) u2 + 6

(
2u3 − 4sin (u) u2 + sin (2u) u2 − 4cos (u) u + 2u + 4sin (u) − 2sin (2u)

)
u
)

(11)



β0,1 = −
49

160 −
13

1680 u
2 −

277
1728000 u

4 −
14911

13970880000 u
6 + 3775394358914560000 u

8

β1,1 = −
13
10 +

67
2100 u

2 + 289378000 u
4 + 7703873180000 u

6 −
334687

1362160800000 u
8

β2,1 =
97

160 −
29

1200 u
2 −

7309
12096000 u

4 −
108337

13970880000 u
6 + 346729721794572800000 u

8

δ2,1 = −
33
80 +

23
1400 u

2 + 179403200 u
4 + 155712328480000 u

6 −
263267

3632428800000 u
8

γ2,1 =
23
240 −

1
2100 u

2 −
373

6048000 u
4 −

15901
6985440000 u

6 −
11203

222393600000 u
8

(12)



β0 =
1
6

((
−2u3 + 24u

)
sin(u) − u4 − 12cos (u) u2 + 24cos (u) − 24

)
/
((
−4sin (u) u2 + sin (2u) u2 + 2u3 − 4cos (u) + 2cos (2u) u + 4sin (u) − 2sin (2u) + 2u

)
u
)

β1 =
1
6

((
u3 − 6u

)
sin(2u) +

(
3u2 − 6

)
cos (2u) + 4 u4 + 12cos (u) u2 − 3u2 − 12sin (u) u + 6

)
/
(
u
((

1
2 u

2 − 1
)

sin (2u) + u3 − 2sin (u) u2 − 2cos (u) u + cos (2u) u + u + 2sin (u)
))

β2 =
1
6

(
22sin(u) u3 + sin(2u) u3 − 5u4 + 36cos (u) u2 − 6cos (2u) u2 − 24sin (u) u − 18u2 + 24cos (u) − 12cos (2u) − 12

)
/
((
−2u − 2cos (2u) u − sin (2u) u2 − 4sin (u) + 2sin (2u) + 4cos (u) u + 4sin (u) u2 − 2u3

)
u
)

δ2 =
1
6

(
−8sin (u) u3 + sin (2u) u3 − 12cos (u) u2 − 12sin (u) u + 6sin (2u) u + 12u2 − 48cos (u) + 12cos (2u) + 36

)
/
((
−2u − 2cos (2u) u − sin (2u) u2 − 4sin (u) + 2sin (2u) + 4cos (u) u + 4sin (u) u2 − 2u3

)
u
)

γ2 =
1
6 ((−u2 +6)cos(2u) +8cos(u) u2 +5u2−32sin(u) u+4sin(2u) u−48cos(u) +42)

((−4sin(u) u2 +sin(2u) u2 +2u3−4cos(u)u +2cos(2u) u+4sin(u) )−2sin(2u) +2u)u

(13)

14



Abdulganiy et al. / J. Nig. Soc. Phys. Sci. 2 (2020) 12–25 15

β0 = −
1

160 −
1

1680 u
2 −

53
1728000 u

4 −
15119

13970880000 u
6 −

115693
4358914560000 u

8

β1 =
3
10 +

1
700 u

2 + 41378000 u
4 + 3847873180000 u

6 + 1658171362160800000 u
8

β2 =
113
160 −

1
1200 u

2 −
941

12096000 u
4 −

46433
13970880000 u

6 −
2074607

21794572800000 u
8

δ2 =
17
80 +

1
4200 u

2 + 19403200 u
4 + 52192328480000 u

6 + 831191210809600000 u
8

γ2 =
7

240 +
1

2100 u
2 + 436048000 u

4 −
269

6985440000 u
6 −

84803
10897286400000 u

8

(14)

For k = 3,

yn=yn+2+hβ0,1 (ucos (u) , sin(u)) fn+hβ1,1(cos (u) , sin(u)) fn+1+hβ2,1(cos (u) , sin(u)) fn+2

+hβ3,1(cos (u) , sin(u)) fn+3+h2δ3,1(cos (u) , sin(u))gn+3+h3γ3,1(cos (u) , sin(u))ln+3
(15)

yn+1=yn+2+hβ0,2 (cos (u) , sin(u)) fn+hβ1,2 (ucos (u) , sin(u)) fn+1+hβ2,2(cos (u) , sin(u)) fn+2

+hβ3,2(cos (u) , sin(u)) fn+3+h2δ3,2(cos (u) , sin(u))gn+3+h3γ3,2(cos (u) , sin(u))ln+3
(16)

yn+3=yn+2+hβ0 (cos (u) , sin(u)) fn+hβ1 (cos (u) , sin(u)) fn+1+hβ2 (ucos (u) , sin(u)) fn+2

+hβ3 (cos (u) , sin(u)) fn+3+h2δ3(cos (u) , sin(u))gn+3+h3γ3(cos (u) , sin(u))ln+3
(17)



β0,1 =
1
3 (−16sin(u)u

3 + 17sin(2u)u3 + 6u4 − 30cos(u)u2 + 6cos(3u)u2 + 48cos(2u)u2 + 90sin(u)u
−18sin (3u) u − 54sin (2u) u + 12u2 + 21cos (u) − 21cos (3u) + 24cos (2u) − 24)/ρ

β1,1 =
1
3 (−81sin (u) u

3 − 17sin (3u) u3 + 24u4 − 117cos (u) u2 − 75cos (3u) u2 + 30sin (u) u + 126sin (3u) u
−60sin (2u) u + 48u2 − 78cos (u) + 78cos (3u) − 84cos (2u) + 84)/ρ

β2,1 =
1
3 (16sin (3u) u

3 + 81sin (2u) u3 + 6u4 − 54cos (u) u2 + 78cos (3u) u2 + 144cos (2u) u2 + 150sin (u) u
−102sin (3u) u − 120sin (2u) u + 12u2 + 57cos (u) + 57cos (3u) + 24cos (2u) − 24)/ρ

β3,1 =
1
3 (−11sin (u) u

3 − 11sin (3u) u3 − 44sin (2u) u3 + 6u4 + 21cos (u) u2 − 45cos (3u) u2 − 48cos (2u) u2

+12sin (u) u + 36sin (3u) u + 12sin (2u) u + 36cos(2u) + 12u2 + 57cos (u) + 57cos (3u) + 24cos (2u) − 24)/ρ

δ3,1 = (
(
−24u3 − 48

)
cos(2u) +

(
−6u2 + 33

)
cos (3u) − 6cos (u) u2 − 71sin (u) u + 28sin (2u) u + 29sin (3u) u

−33cos (u) + 48)/ρ

γ3,1 = (29 sin (3 u) u + 28 u sin (2 u) − 6 u2 cos (u) − 24 cos (2 u) u2 − 6 cos (3 u) u2 − 71 sin (u) u − 33 cos (u) − 48 cos (2 u)
+33 cos (3 u) + 48)/(

(
−6 u3 + 21 u

)
sin (3 u) + 18 (−cos (3 u) u +

(
3/2 u2 − 173

)
sin (2 u) + u3 − 3 u2 sin (u)

−5 u cos (u) + 4 cos (2 u) u + 2 u + 47 sin(u)6 )u)

(18)



β0,1 = −
121
405 −

289
42525 u

2 −
4397

53581500 u
4 + 518879123773265000 u

6 + 1540154587405481216140000 u
8

β1,1 = −
23
15 +

101
3150 u

2 + 18733969000 u
4 −

12401
833490000 u

6 −
48315783

300356456400000 u
8

β2,1 =
1
3 −

23
315 u

2 −
619

396900 u
4 + 5353916839000 u

6 + 7851770930035645640000 u
8

β3,1 = −
203
405 +

4061
85050 u

2 + 125353107163000 u
4 + 1200029247546530000 u

6 −
11234074463

8109624322800000 u
8

δ3,1 =
10
27 −

83
2835 u

2 −
3079

3572100 u
4 −

93587
8251551000 u

6 + 1115045320793908520000 u
8

γ3,1 = −
4
25 +

2
675 u

2 + 6112976750 u
4 + 548236876292500 u

6 + 45298619225267342300000 u
8

(19)

15



Abdulganiy et al. / J. Nig. Soc. Phys. Sci. 2 (2020) 12–25 16

β0,2 =
1
24 (44sin (u) u

3 + 17sin (2u) u3 − 6u4 + 60cos (u) u2 + 6cos (3u) u2 + 102cos (2u) u2 + 132sin (u) u
−210sin (2u) u − 18u2 + 384cos (u) − 168cos (2u) − 216)/ρ

β1,2 =
1
24 (−351sin (u) u

3 − 17sin (3u) u3 + 78u4 − 486cos (u) u2 − 54cos (3u) u2 − 216cos (2u) u2 − 18sin (u)
+66sin (3u) u + 486sin (2u) u + 108u2 − 1248cos (u) + 600cos (2u) + 684)/ρ

β2,2 =
1
24 (

(
−351u3 + 1278u

)
sin(2u) +

(
44u3 − 102u

)
sin (3u) − 78u4 + 432cos (u) u − 1134cos (2u) u2

+108cos (3u) u2 − 126u2 − 810sin (u) u + 384cos (u) + 246cos (2u) − 648)/ρ

β3,2 =
1
24 (−125sin (u) u

3 + 13sin (3u) u3 − 152sin (2u) u3 − 6u4 + 138cos (u) u2 + 18cos (3u) u2

−444cos (2u) u2 + 204sin (u) u + 186sin (3u) u + 1248sin (2u) u − 168cos(2u) − 1080)/ρ

δ3,2 =
1
4 (13sin (u) u

3 − sin (3u) u3 + 13sin (2u) u3 − 24cos (u) u2 + 24cos (2u) u2 − 87sin (u) u − 3sin (3u) u
+48sin (2u) u − 288cos (u) + 72cos (2u) + 216/ρ

γ3,2 =
1

24 (78cos (u) u
2 − 6cos (3u) u2 + 78cos (2u) u2 + 353sin (u) u + 13sin (3u) u − 340sin (2u) u − 6u2

+960cos (u) − 312cos (2u) − 648)/ρ

(20)



β0,2 =
1

90 +
1

800 u
2 + 1571984500 u

4 + 1053089293388480000 u
6 + 1456666991201425825600000 u

8

β1,2 = −
61

160 −
331

67200 u
2 −

14369
42336000 u

4 −
6318043

391184640000 u
6 −

611857117
1067934067200000 u

8

β2,2 = −1 +
19

3360 u
2 + 6131058400 u

4 + 63165719559232000 u
6 + 672527513429840000 u

8

β3,2 =
533

1440 −
19

9600 u
2 −

40501
127008000 u

4 −
23157647

1173553920000 u
6 −

8248324979
961140660480000 u

8

δ3,2 = −
11
48 −

1
2240 u

2 + 5834233600 u
4 + 42094139118464000 u

6 + 169075937320380220160000 u
8

γ3,2 =
11
240 +

47
33600 u

2 + 70321168000 u
4 + 491195592320000 u

6 −
1081123

19776556800000 u
8

(21)



β0 =
1
24 (−20sin (u) u

3 + sin (2u) u3 − 6u4 − 132cos (u) u2 + 6cos (2u) u2 + 324sin (u) u − 18sin (2u) u − 18u2

+384cos(u) − 24cos (2u) − 360)/ρ

β1 =
1
24 (81sin (u) u

3 − sin (2u) u3 + 30u4 + 522cos (u) u2 − 6cos (2u) u2 − 1218sin (u) u + 18sin (3u) u
+6sin (2u) u + 60u2 − 1272cos (u) + 24cos (3u) − 24cos (2u) + 1272)/ρ

β2 =
1
24 (20sin (u) u

3 − 81sin (2u) u3 − 114u4 − 432cos (u) u2 + 84cos (3u) u2 − 306cos (2u) u2 + 714sin (u) u
−186sin (3u) u + 642sin (2u) u − 66u2 − 192cos (u) − 192cos (3u) + 840cos (2u) − 456)/ρ

β3 =
1
24 (371sin (u) u

3 + 29sin (3u) u3 + 136sin(2u)u3 − 54u4 + 762cos (u) u2 + 66cos (3u) u2 − 276cos (2u) u2

−948sin (u) u + 186sin (2u) − 264u2 + 1080cos (u) + 168cos (3u) − 792cos (2u) − 456)/ρ

δ3 =
1
4 (−19sin (u) u

3 − sin (3u) u3 + 5sin (2u) u3 − 24cos (u) u2 − 31sin (u) u − 11sin (3u) + 32sin (2u) + 24u2

−264cos (u) − 24cos (3u) + 120cos (2u) + 168)/ρ

γ3 =
1

144 (116sin(2u) +
(
6u2 − 48

)
cos (3u) + 114cos(u)u2 − 30cos (2u) u2 + 54u2 − 433sin (u) u − 29sin (3u) u

−912cos (u) + 264cos (2u) + 696)/ρ

(22)



β0 =
1

810 +
229

1360800 u
2 + 35326790750 u

4 + 60449877921488960000 u
6 + 114774361732438497291200000 u

8

β1 = −
7

480 −
97

201600 u
2 −

6403
127008000 u

4 −
3738641

1173553920000 u
6 −

1490068837
9611406604800000 u

8

β2 =
1
3 −

1
1440 u

2 + 1513175200 u
4 + 2864595867769000 u

6 + 1711626160071291280000 u
8

β3 =
8813
12960 +

5483
5443200 u

2 −
35383

3429216000 u
4 −

77924501
31685955840000 u

6 −
42892337857

259507978329600000 u
8

δ3 = −
83
432 −

209
181440 u

2 −
1571

114307200 u
4 + 8447031056198528000 u

6 + 7008125718650265944320000 u
8

γ3 =
17

720 +
167

302400 u
2 + 3383190512000 u

4 + 8740511760330880000 u
6 + 12076795714417109907200000 u

8

(23)

16



Abdulganiy et al. / J. Nig. Soc. Phys. Sci. 2 (2020) 12–25 17

where

ρ = u(18u2sin(u) + 2sin(3u)u2 − 9u2sin (2u)

−6u3 + 30ucos (u) + 6cos (3u) u − 24cos (2u) u
−47sin (u) − 7sin (3u) + 34sin(2u) − 12u) (24)

It is interesting to note that as u → 0 in the power series
expansion of the BTDTFMs, classical third derivative methods
based on polynomial basis are recovered.

3. Basic Properties of BTDTFM

The basic properties of BTDTFM which includes the Local
Truncation Error (LTE), Order, Error Constant, Zero Stability,
Convergence and Region of Stability are discussed.

3.1. Local Truncation Error (LTE) and Order of BTDTFM

Theorem 1. The BTDTFM has a local truncation error of the
form LT E = Ck+4hk+4(ω2yk+2(xn) − y(k+4)(xn)) + O(hk+5).

Proof 1
Since the BTDTFM is made up of generalized linear multistep
methods with trigonometric coefficients, we associate the BT-
DTFM with a linear operator Lω

[
y (xn) ; h

]
for the principal

method and Lω
[
y (xn) ; h

]
for the secondary methods defined

respectively by



Lω
[
y (xn) ; h

]
= y (xn + kh)

−
(
yn+k−1 + h

∑k
j=0 β j (u)y

′

(xn + kh)
+h2δky

′′

(xn + kh) + h3γky
′′′

(xn + kh)
Lω

[
y (xn) ; h

]
= y (xn + ih)

−
(
yn+k−1 + h

∑k
j=0 β j (u)y

′

(xn + kh)
+h2δk(u)y

′′

(xn + kh) + h3γk(u)y
′′′

(xn + kh)

(25)

Assuming that y(xn) is sufficiently differentiable, expanding
with Taylor series expansions of y(xn +kh), y(xn +ih), y′(xn +kh),
y′′(xn + kh) and y′′′(xn + kh) about the point xn.
Substituting the coefficients β j (u) ,δk (u) ,γk (u) ,β j (u) δk(u) and
γk(u) into equations (7) and (8) respectively, after simplifica-
tion, we obtain the order, the Local Truncation Error as pre-
sented in Table 1.

�

3.2. Consistency of the BTDTFMs

Since each BTDTFM is of order p ≥ 5 > 1, (for each k
)as shown in Table 1, then it is consistent (Lambert [57] and
Fatunla [58]).

3.3. Stability of the BTDTFM

The combined methods (Principal and Secondary) are assemble
in the block form given by

A1Yµ+1 = A0Yµ+hB0 Fµ+hB1 Fµ+1+h
2Gµ+1C1+h

3 D1 Lµ+1(26)

Where

Yµ+1= (yn+1, yn+2, yn+3)
T , Yµ = (yn−2, yn−1, yn)

T ,

Fµ+1 = ( fn+1, fn+2, fn+3)
T ,

Fµ = ( fn−2, fn−1, fn)
T , Gµ+1 = (gn+1, gn+2, gn+3)

T , Lµ+1 = (ln+1, ln+2, ln+3)
T .

A0, A1, B0, B1, C1 and D1 are k×k matrices defined in canonical
form as follows:

For k= 2, we have

A1=
(
−1 0
−1 1

)
, A0=

(
0 −1
0 0

)
, B0=

(
0 β0,0
0 β0

)
,

B1=
(
β1,0 β2,0
β1 β2

)
,

C1=
(

0 δ2,0
0 δ2

)
and D1=

(
0 γ2,0
0 γ2

)
For k = 3, we have

A1 =

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

 , A0 =


0 0 0
0 0 −1
0 0 0

 , B0 =


0 0 β0,1
0 0 β0,0
0 0 β0

 ,
B1 =


β1,1 β2,1 β3,1
β1,0 β2,0 β3,0
β1 β2 β3

 ,
C1 =


0 0 δ3,1
0 0 δ3,0
0 0 δ3

 and D1 =


0 0 γ3,1
0 0 γ3,0
0 0 γ3


3.4. Zero Stability

Definition 2. A block method is zero stable if the roots of the
first characteristic polynomial have modulus less than or equal
to one and those of modulus one do not have multiplicity greater
than 2. i.e. ρ(R) = det(RA(1) − A(0)) = 0 satisfies |Ri| ≤ 1 and
for those roots with |Ri| = 1, the multiplicity does not exceed 2
Fatunla [58].

Proposition 3. The BTDTFM is zero stable.

Proof 2
From the normalized first characteristic polynomial of BTDTFM,
we have in canonical form that ρk (R) = det

[∑1
i=0 A1−iR

i
]
.

so that ρk (R) = 0 =⇒ −Rk−1 (1 + R) = 0. Consequently, the
roots R j, j = 1, 2, . . . k of ρk(R) satisfy

∣∣∣R j∣∣∣ = 1, the roots are
simple. Hence for each k = 2 and k = 3, the BTDTFM is Zero
stable �

17



Abdulganiy et al. / J. Nig. Soc. Phys. Sci. 2 (2020) 12–25 18

Table 1: Local Truncation Error of BTDTFM

k LTE Order ( p) Error Constant (C p+1)

2


(
−

1
1800 D

(6) (y) (x) − 11800 ω
2 D(4) (y)(x)

)
h6(

−
1

200 D
(6) (y) (x) − 1200 ω

2 D(4) (y)(x)
)

h6


[

5
5

]
−

1
1800 −

1
200

3


(
−

11
50400 D

(7) (y) (x) − 1150400 ω
2 D(5) (y)(x)

)
h7( 29

6300 D
(7) (y) (x) + 296300 ω

2 D(5) (y)(x)
)

h7(
−

59
50400 D

(7) (y) (x) − 5950400 ω
2 D(5) (y)(x)

)
h7

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


6
6
6

 − 1150400 296300 − 5950400
3.5. Convergence of BTDTFM

Theorem 4. Let Y be an approximation of the solution vector
Y for the system obtained from BTDTFM given by equations
(6) and (7) respectively, If en = |y (xn) − yn|, where the exact
solution is several times differentiable on [x0, xN ] and if ‖E‖ =∥∥∥Y − Y∥∥∥, then for sufficiently small h,BTDTFM is a (k + 3)th-
order convergent method. That is ‖E‖ = O(hk+3).

Proof 3
The proof can be readily obtained, similarly to the one given in
Abdulganiy et al. (2017). �

3.6. Linear Stability and Region of Absolute Stability of BT-
DTFM

According to Ndukum et al. [47], linear stability analysis of
trigonometrically fitted methods is more complicated than those
of the corresponding polynomial-based methods for which λandh
occur only in the combination z = λh . For trigonometrically
fitted methods, three parameters are considered since the step
length occurs in both z = λh and u = ωh . To analyze the linear
stability of BTDTFM, the block method (25) is applied to the
test equations y

′

= λy , y
′′

= λ2y and y
′′′

= λ3y . After simple
algebraic simplification and letting z = λh , we obtain

Yw+1 = M(z, u)Yw

where

M (z, u) =
[
A1 − B1z − C1z

2
− D1z

3
]−1

[A0 + B0z)] (27)

The rational function M (z, u) is called the amplification matrix
which determines the stability of the method.

Definition (Coleman and Ixaru, [59]: A region of stability is
a region in the zu−plane throughout which |ρ (z, u)| ≤ 1, where
ρ (z, u) is the spectral radius of M (z, u).

Since the stability matrix depends on two parameters z and u,
we plot the stability regions in the (z, u)−plane for both k = 2
and k = 3 respectively in Figures 1 and 2.

Definition (Ndukum etal., [47]): A BTDTFM with variable
coefficients (A0 (u) , A1 (u) , B0 (u) ,

B1(u), C1(u), D1(u)) with the stability function ρ (z, u) is said to
be A-stable at u = u0 if |ρ (z, u0)| < 1 for all z ∈ C−.

We plot the region |ρ (z, u)| < 1 (for each k) in the complex
plane via boundary locus method and look for the interval of

Figure 1: The z − u Plot for BTDTFM k = 2

Figure 2: The z − u Plot for BTDTFM k = 3

u for which z ∈ C−. For k = 2, we notice that the values
of u ∈ [2π,∞) are satisfactory and for k = 3, the values of
u ∈ [2.96π, 2.70π] are satisfactory. Figures 3 and 4 show the
plot of |ρ (z, u)| < 1 at u0 = 2π for k = 2 and u0 = 2.96π for
k = 3 respectively.

It is also worthy to note that the stability function of the BT-
DTFM for each k satisfies |ρ (z, u)| < 1 for all z ∈ C−, for some
u and limz→∞ ρ (z, u) = 0 is a property analogous to L−stability
for the classical methods which is essential for a method to per-
form well on highly stiff problems (Ndukum etal., [47]). Hence
the BTDTFMs are L−stable.

18



Abdulganiy et al. / J. Nig. Soc. Phys. Sci. 2 (2020) 12–25 19

Figure 3: RAS for BTDTFM k = 2

Figure 4: RAS for BTDTFM k = 3

4. Estimation of Computational Frequency

In fitted methods for solving periodic problems, the choice
of estimating computational frequency remains a hard nut to
crack. A rigorous theory for the exact computation of the fre-
quency is yet to be developed. However, some attempts have
been made by a number of researchers. For exponentially fit-
ted methods, the procedure for estimating frequency is given in
Ixaru et al. ([29], [60,61]), Vanden Berghe et al. ([30], [62,63])
and Van de Vyve.r [25]. The estimation of the computational
frequency of trigonometrically-fitted methods can be found in
Ramos and Vigo-Aguiar [64], Vigo-Aguiar and Ramos [65],
Ngwane and Jator [42], Jator [48,49] and Ndukum et al., [47].
It is interesting to note that all the aforementioned methods fo-
cus on the optimization of the local truncation error (LTE).

In the present study, the approach described in Vanden Berghe
et al. [30] for the exponentially fitted methods is adapted to the
BTDTFM with Trigonometric coefficients.

It is observed that the LTE of the principal and secondary
methods of the BTDTFM consist of a product of three factors.
The LTE of the principal method of k=2, for instance, consists
of the product of a general factor, a numerical factor in rational
form

(
−

1
200

)
, and a factor which involves two derivatives of the

solution vector. For easy reference, we introduce the follow-

ing functional D
[
y (x) ; ω

]
=

(
y(6) (x) + ω2y(4) (x)

)
. Assume ω

exists such that D is identically vanishes on the given interval,
then the principal method corresponding to that ω will be exact.
The reason according to Vanden Berghe et al. [30] is equiva-
lent to solving the differential equation y(6) (x) + ω2y(4) (x) = 0
and, indeed, sin (ωx) and cos (ωx) are solutions. In general,
no constant ω can be found such that D is identically vanished
but it makes logic to discourse the problem of obtaining ω for
which D are kept close to zero as possible for x in the given
interval with the result that upon the Taylor series expansion of
ω2 = −

y(6) (xn )
y(4) (xn )

about xn given by

D
[
y (x) ; ω

]
= − (x − xn)

(
y(7) (xn) −ω2y(5) (xn)

)
+O

(
(x − xn)

2
)

.

This implies that the bound of D behaves as h and conse-
quently, the bound of LTE for k = 2 in Table 1 behaves as h7

and hence increase the order by 1.

Remark 1. We remark that other steps and the technical de-
tails for implementations are as described in section 3 of Var-
den Berghe et al. [30].

5. Numerical Experiments

A number of numerical examples are provided in this sec-
tion to illustrate the accuracy and computational efficiency of
BTDTFM. We have calculated the absolute error at the end
point for stiff problem as |y (t) − y| and the maximum abso-
lute error of the approximate solution as Err = max |y (t) − y|
for periodic problems and the efficiency is plotted as Number
of Function Evaluation (NFE) against max |y (t) − y| .

5.1. Stiff Problems

Example 5.1 We consider as our first example a non-linear
system of first order differential equations in the range 0 ≤ t ≤
10

y
′

1 = µy1 + y
2
1 , y1 (0) = −

1
µ + 2

y
′

2 = −y2 , y2 (0) = 1

whose solution in closed form is given by y1 = −
e−2t
µ+2 , y2 =

e−2t , where µ = 10000 .

The new BTDTFM is used to solve this problem with ω = 1
and compared with the Second Derivative Multistep Method
(SDMM) in Hojjati et al. [7], Simpson’s 38 Method

(
S 3

8

)
in

Akinfenwa et al. [16], Multiderivative Hybrid Implicit Runge-
Kutta Method (MHIRK) in Akinfenwa et al. [13], L stable
Method (LM) in Mehdizabeh and Molayi [14], Third deriva-
tive block hybrid method (TDBHM) in Akinfenwa [50] and the

19



Abdulganiy et al. / J. Nig. Soc. Phys. Sci. 2 (2020) 12–25 20

Figure 5: Numerical solution of Example 5.1 for h = 0.1

results are displayed in Table 2 for the accuracy while Figure 5
show the numerical solution of BTDTFM against the time t for
h = 0.1 and ω = 1 .

It is evident from Table 2 that the newly constructed method
is more accurate than some of the previously known method
in the literature. In particular, BTDTFM k = 2 is far more
accurate with large step size compare to SSDM with smaller
step size and compete favorably with S 3

8
with same step size,

notwithstanding S 3
8

is of order 8.

Example 5.2

Consider the non-linear stiff system given by

y
′

1 (t) = −1002y1 + 1000y
2
2 , y1 (0) = 1

y
′

2 (t) = y1 − y2 (1 + y2) , y2 (0) = 1

Whose exact solution is given as y1 (t) = y22 , y2 (t) = e
−t

We solved this problem with ω = 1 . The results of BT-
DTFM k = 2 for this problem in comparison with Exponen-
tially fitted Gauss (EF-Gauss), Gauss-2s ,methods in Vaquero
and Vigo-Aguiar [18] and Hybrid Backward differentiation for-
mula (HBDF) of Jator and Agyingi [46] for h = 0.1 and h =
0.01 are displayed in Table 3. Table 4 presents the results of BT-
DTFM k = 2 for different values of h as compared to Extended
Continuous Block Backward Differentiation Formula (ECBBDF)
of Akinfenwa and Jator [12] and Continuous Block Backward
Differentiation Formula (CBBDF) of Akinfenwa et al. [11] for
0 ≤ t ≤ 10 while Table 5 shows the results of comparison with
the LM method in Mehdizabeh and Molayi [14] and the method
in Wu [6] for h = 0.05 in the interval 0 ≤ t ≤ 50 respec-
tively. The graphical representation of the numerical results of
BTDTFM at h = 0.05 is displayed in Figure 6.

It is seen from Tables 3-5 that for various values of h con-
sidered, BTDTFM does better in terms of accuracy than the
methods in [18], [46], [12], [11], [14] and [6] respectively.

Figure 6: Numerical solution of Example 5.2 for h = 0.05

Example 5.3

Consider a classical four-dimensional problem given by
y′1 (t)
y′2 (t)
y′3 (t)
y′4 (t)

 =

−104y1 (t) + 100y2 (t) − 10y3 (t) + y4 (t)
−1000y2 (t) + 10y3 (t) − 10y4 (t)

−y3 (t) + 10y4 (t)
−0.1y4 (t)

 ,

y1 (0)
y2 (0)
y3 (0)
y4 (0)

 =
1
1
1
1


The eigenvalue of the Jacobian matrix λ1 = −0.1, λ2 =

−1.0, λ3 = −1000 andλ4 = −10000 . The exact solution of
Example 5.3 is given as

y1 (t) = −
89990090

8999010009
e−0.1t +

818090
89901009

e−t

+
9989911

899010090
e−1000t +

89071119179
89990100090

e−10000t

y2 (t) =
9100

89991
e−0.1t −

910
8991

e−t +
9989911
9989001

e−1000t,

y3 (t) =
100

9
e−0.1t −

91
9

e−t, y4 (t) = e
−1.0t

This problem is considered to show the performance of BT-
DTFM on a 4 × 4 stiff problems.

Table 6 shows the numerical results of BTDTFM with ω =
1 in comparison with the Second Derivative Extended Back-
ward Differentiation Formula (SDEBDF) of Ehigie et al. [21]
and the method NMTD in Adesanya et al. [22] for h = 0.05
and 0.1 at t = 20 while Table 7 displays the results BTDTFM as
it compared with AB7 method in Abhulimen and Otunta [19],
CEGE method of Abhulimen and Omeike [20] for h = 0.05
and 0.1 at t = 1 . The graphical representation of the numerical
results of BTDTFM at h = 0.05 is shown in Figure 7.

The results in Tables 6 and 7 showed the superiority of BT-
DTFM in terms of accuracy over the methods it compared with

20



Abdulganiy et al. / J. Nig. Soc. Phys. Sci. 2 (2020) 12–25 21

Table 2: Comparison of results for Example 5.1

t h=0.0001
SSDM

h=0.1 S 3
8

h=0.1
MHIRK

h = 0.0001
LM

h = 0.1
T DBH M

h = 0.1
BT DT F Mk =
2

Err (y1)
Err (y2)

3 2.48E-11
2.47E-06

2.03E-19
1.44E-14

3.06E-15
3.08E-10

1.78E-20
2.08E-12

1.34E-17
2.98E-13

6.54E-19
6.56E-14

Err (y1)
Err (y2)

5 3.45E-14
2.30E-08

1.20E-20
3.21E-15

9.35E-17
6.94E-11

2.49E-19
4.66E-13

2.94E-18
1.79E-13

2.00E-20
1.48E-14

Err (y1)
Err (y2)

10 3.46E-18
3.15E-11

1.11E-20
4.38E-17

8.49E-21
9.36E-13

5.74E-20
6.35E-12

3.82E-20
2.84E-14

1.81E-24
2.00E-16

Table 3: Comparison of errors at t = 5 for Example 5.2

h Gauss-2s EF-Gauss HBDF BTDTFM k = 2
0.1 5.12E − 06 8.36E − 07 5.49E − 11 1.82E − 11
0.01 6.46E − 12 4.24E − 12 6.07E − 15 2.01E − 16

Table 4: Comparison of methods at t = 10 for Example 5.2

h y C BBDF4 C BBDF5 EC BBDF4 EC BBDF5 BTDTFM k = 2
0.02 y1

y2
4.88E-16
5.39E − 12

8.37E-18
9.16E − 14

2.48E-19
3.75E − 16

1.33E-20
1.35E − 16

5.76E-19
6.34E − 15

0.01 y1
y2

3.13E-17
3.45E − 13

3.39E-21
1.23E − 17

2.68E-19
2.93E − 15

2.87E-22
2.93E − 15

1.82E-20
2.00E − 16

Table 5: Comparison of result of at t = 50 for Example 5.2

h y LMM WU Method BTDTFM k = 2
0.05 y1

y2
4.13E − 25
1.29E − 22

2.3E − 21
1.4E − 18

4.89E − 51
1.27E − 29

Table 6: Results comparison at t = 20 for Example 5.3

h Method y1 − Error y2 − Error y3 − Error y4 − Error
0.05 S DE BDF

N MT D
BT DT F M
k = 2

5.31E-12
−

3.71 × E − 17

7.27E-11
−

3.75E − 16

5.90E-09
−

4.13E − 17

1.34E-09
−

3.41E − 15

0.1 DE BDF
N MT D
BTDTFM
k = 2

2.25E-10
8.96E-07
1.18E − 15

2.29E-09
2.06E-08
1.19E − 10

2.50E-07
8.96E-10
1.31E − 12

2.06E-08
8.06E-11
1.19E − 13

in the reviewed literature.

5.2. Problems with Periodic Solutions

Example 5.4

We consider the following system of couple differential equa-

tions which is well known as the two body problem

y
′′

1 (t) = −
y1
r3
, y1 (0) = 1, y

′

1 (0) = 0

y
′′

2 (t) = −
y2
r3
, y2 (0) = 0, y

′

2 (0) = 0

21



Abdulganiy et al. / J. Nig. Soc. Phys. Sci. 2 (2020) 12–25 22

Table 7: Results comparison at t = 1 for Example 5.3

h Method y1 − Error y2 − Error y3 − Error y4 − Error
0.05 AB7

CEGE
N MT D

BT DT F M k = 2

3.2E-02
3.5E-05
7.5E-11
2.2E − 12

3.2E-02
3.8E-04
8.6E-10
2.5E − 12

3.3E-01
3.5E-07
8.3E-08
2.5E − 09

3.7E-05
3.7E-08
2.9E-09
1.3E − 11

0.1 AB7
CEGE
N MT D
BTDTFM k = 2

2.5E-02
2.9E-05
1.4E-08
6.7E − 11

2.1E-01
2.7E-04
1.3E-08
7.5E − 10

2.4E-03
2.6E-06
2.1E-09
7.5E − 08

2.7E-05
2.6E-08
2.7E-11
4.0E − 14

Figure 7: Numerical solution of Example 5.3 for h = 0.05

where r =
√

y21 + y
2
2 and whose analytic solution is given by

y1 (t) = cos (t) , y2 (t) = sin (t)

The problem was considered in Senu et al. [36] and Abdul-
ganiy et al. [44] in the interval 0 ≤ t ≤ 10 with ω = 1.
The BTDTFM is compared with Block Hybrid Trigonometry
Method (BHTM) of Abdulganiy et al. [44] and a New Diago-
nally Implicit Runge-Kutta Nystrom Method for Periodic IVPs
(DIRKNNew) of Senu et al. [36] and the numerical results are
displayed in Table 8 below. While the graphical illustration of
BTDTFM is shown in Figure 8, Figure 9 shows the efficiency
of BTDTFM. It is evident from Table 8 that the newly devel-
oped method BTDTFM uses fewer number of function evalua-
tion and consequently produces good results and gives a better
approximation when compared with BHTM and DIRKNNew.
Also, Figure 9 hows that the method in this paper is efficient.

Example 5.5: Nonlinear Strehmel-Weiner problem

We consider the nonlinear second order IVP which was also
solved by Nguyen et al. [35] and Jator [37] in the interval 0 ≤
t ≤ 10 respectively

y”1 (t) = (y1 (t) − y2 (t))
3

+ 6368y1

Figure 8: Numerical solution of Example 5.4

Figure 9: Efficiency Curve of Example 5.4

× (t) − 6384y2 (t) + 42 cos (10t) ,
y1 (0) = 0.5, y

′

1 (0) = 0
y”1 (t) = − (y1 (t) − y2 (t))

3
+ 12768y1 (t) − 12784y2 (t)

+42 cos (10t) , y2 (0) = 0.5, y
′

2 (0) = 0

with analytic solution in closed form given by y1 (t) = y2 (t) =
cos (4t) − cos(10t)2 .

Numerical results of the maximum global errors of BTDTFM

22



Abdulganiy et al. / J. Nig. Soc. Phys. Sci. 2 (2020) 12–25 23

Table 8: Comparison of Numerical Results for Example 5.4

N
BTDTFM k = 2 BHTM DIRKNNew

Error NFE Error NFE Error NFE
100 2.84E − 29 201 5.13E − 27 1600 7.49 E − 4 5000
200 1.92E − 28 401 1.30E − 29 3200 5.62E − 7 10000
400 1.18E − 27 801 4.0E − 30 6400 1.00E − 9 25000
800 2.47E − 27 1601 7.43E − 28 12800 1.78E − 7 60000

Table 9: Comparison of Maximum Errors and Number of Function Evaluation

BTDTFM k = 2 SVBM TIRKN 3
N F E Err N F E Err N F E Err
751 2.25E-7 801 2.6E-7 907 2.5E-4
1126 2.95E-8 1201 1.6E-8 1288 6.6E-6
1501 7.01E-9 1601 2.8E-9 1682 7.0E-6

Table 10: Comparison of Number of steps used in the computation for ω = 10

Steps BTDTFM k = 2
Err

Vigo-Aguiar and Ramos
Err

NEWS5(4)-2
Err

898 5.5 2.5 4.9
1344 6.6 3.2 6.5
2123 7.8 3.4 7.8
2990 8.7 4.6 8.9
4690 9.9 5.3 9.9
7215 11.0 6.1 11.0

k = 2 were compared with Symmetric Boundary Value Method
(SBVM) of Jator [37] and Trigonometric Implicit Runge-Kutta
Methods (TIRKM) of Nguyen [35] are presented in Table 9
while Figure 10 is the graphical representation of the numerical
results of BTDTFM. Also, Figure 11 shows that the method in
this paper is more efficient

In Table 9 and Figure 11, we show that a fifth order BT-
DTFM has almost the same maximum errors with SVBM of
order 6 but uses fewer number of function evaluation and con-
sequently a more accurate and more efficient integrator for the
nonlinear Strehmel-Weiner problem.

Example 5.6: A Nonlinear Oscillator Lastly, we consider
the initial-value problem

ẍ = −100x + sin (x) , x (0) = 0, ẋ (0) = 1, t ∈ [0, 20π]

The analytical solution is not known but for comparison
purposes we use that x (20π) = 0.000392823991 given by Tsi-
touras and Simos [34]. The comparison of different number
of steps of BTDTFM k = 2 with NEW5(4)-2 of Tsitouras and
Simos [34] and method in Vigo-Aguiar and Ramos [65] with
ω = 10. The number of steps are obtained as contained in the
references from literature. The error is calculated as Err =

Figure 10: Numerical solution of Example 5.5

−log
(
|x

(
t f − x f

)
|
)

where x f is the approximation value ob-
tained at the final point t f .

Table 10 revealed that BTDTFM compete favourably with
methods in Tsitouras and Simos [34] and method in Vigo-Aguiar
and Ramos [65].

23



Abdulganiy et al. / J. Nig. Soc. Phys. Sci. 2 (2020) 12–25 24

Figure 11: Efficiency Curve of Example 5.5

Figure 12: Numerical solution of Example 5.6

6. Conclusion

A family of order k+3 Block Third Derivative Trigonometrically-
Fitted Methods are considered and implemented in this paper
for Stiff and Periodic Problems. The stability properties of the
method was analyzed and was found to be L-stable which in
fact is a requirement for an algorithm to solving stiff problems.
Numerical results on the representative examples show the ac-
curacy of the BTDTFM on linear and nonlinear stiff problems
and its efficiency on nonlinear periodic problems compared to
some accurate and efficient methods respectively in the litera-
ture. The present study is limited to trigonometric basis. Our
future research will consider functionally fitted fitted form as
suggested in Nguyen et al. [35].

Acknowledgments

We thank the referees for the positive enlightening com-
ments and suggestions, which have greatly helped us in making
improvements to this paper.

References

[1] G. G. Dahlquist, A special stability problem for linear multistep methods,
BIT 3 (1963) 27.

[2] T. A. Biala S. N. Jator, R. B. Adeniyi, P. L. Ndukum, “Block Hybrid
Simpson’s Method with Two Offgrid Points for Stiff Systems”, Interna-
tional Journal of Nonlinear Science 20 (2015) 3.

[3] C. W. Gear, Hybrid methods for initial value problems in ordinary differ-
ential equations. SIAM J. Numer. Anal.2 (1965) 69.

[4] W. H. Enright, “Second Derivative Multistep Methods for Stiff ordinary
differential equations”, SIAM J. Numer. Anal. 11 (1974) 321.

[5] J. R. Cash, On the Integration of Stiff Systems of ODEs Using Extended
Backward Differentiation Formulae. Numerische Mathematik 34 (1980)
235-246.

[6] X. U. Wu. A sixth-order A-stable explicit one-step method for stiff sys-
tems, Comput math. Applic.35 (1998) 59.

[7] G. Hojjati, M. Rahimi, S. M. Hosseini. New second derivative multistep
methods for stiff systems. Appl. Math. Model 30 (2006) 466.

[8] S. N. Jator, “Solving second order initial value problems by a hybrid mul-
tistep method without predictors”, Applied Mathematics and Computa-
tion 277 (2010) 4036.

[9] R. K Sahi, S. N. Jator, N. A. Khan, A Simpson’s-Type Second Deriva-
tive Method for Stiff System. International Journal of Pure and Applied
Mathematics 81 619.

[10] F. F. Ngwane, S. N. Jator, “Block hybrid-second derivative method for
stiff systems”, Intern. J. Pure Appl. Math. 4 (2012) 543.

[11] O. A. Akinfenwa, S. N. Jator, , N. Yao, Continuous block backward differ-
entiation formula for stiff ordinary differential equations. Computer and
Mathematics with Applications 65(2013)996.

[12] O. A. Akinfenwa, S. N. Jator, Extended continuous block backward
differentiation formula for stiff systems. Fasciculi Mathematici.(2015).
https://doi.org/10.1515/fascmath-2015-0010.

[13] O. A. Akinfenwa, S. A. Okunuga, B. I. Akinnukawe, U. O. Rufai, R.
I. Abdulganiy, Multiderivative hybrid implicit Runge-Kutta method for
solving Stiff system of a first order Differential equation. Far East Journal
of Mathematical Sciences 106(2018)543.

[14] M. Mehdizadeh, M. Molayi, A new class of L-stable hybrid one-step
methods for the numerical solution of ordinary differential equations.
Journal of Computer Science and Applied Mathematics 1(2015)39.

[15] R. I. Abdulganiy, O. A. Akinfenwa, S. A. Okunuga, A family of L0 stable
Third derivative Block Methods for solving systems of First Order Initial
Value Problems. J. of NAMP 36(2016)47.

[16] O. A. Akinfenwa, R. I. Abdulganiy S. A. Okunuga, V. Irechukwu, Simp-
son’s 38 Type Block Method For Stiff Systems of Ordinary Differential
Equations. Journal of the Nigerian Mathematical Society 36(2017)503.

[17] S. A. Okunuga, A fourth order composite two step method for stiff prob-
lems. International Journal of Computer Mathematics 72(1999)39.

[18] J. M. Vaquero and J. Vigo-Aguiar, Exponential fitted Runge-Kutta meth-
ods of collocation type based on Gauss, Radau, and Labatto tradi-
tional methods. Proceedings of the International Conference on Computa-
tional and Mathematical Methods in Science and Engineering (CMMSE
’07)289.

[19] C. E. Abhulimen, O. F. Otunta, A family of two step Exponentially fit-
ted Multiderivative methods for the numerical integration of stiff IVPs on
ODEs. International Journal of Numerical Mathematics 13(2007)1.

[20] C. E. Abhulimen, G. E. Omeike, A Six-Order Exponentially Fitted
Scheme for The Numerical Solution of Systems of Ordinary Differ-
ential Equations. Journal of Applied Mathematics and Bioinformatics
1(2011)175.

[21] J. O. Ehigie, S. A. Okunuga, A. B. Sofoluwe, A class of exponentially-
Fitted second derivative extended backward differentiation formula for
solving stiff problems. Fasciculi Mathematici 51(2013)71.

[22] A. O. Adesanya, R. O. Onsachi, , M. R. Odekunle, New algorithm for
first order stiff initial value problems. Fasciculi Mathematici 58(2017)2.

[23] D. G. Yakubu, M. Aminu, P. Tumba, M. Abdulhameed, An efficient fam-
ily of second derivative Runge-Kutta collocation methods for oscillatory
systems. Journal of the Mathematical Society 37(2018)111.

[24] G. Vanden Berghe, H. De Meyer, M. Van Daele, T. V. Hecke,
Exponentially-Fitted explicit Runge-Kutta methods. Computer Physics
Communications 123(1999)7.

[25] H. Van de Vyver, Frequency evaluation for exponentially fitted Runge-

24



Abdulganiy et al. / J. Nig. Soc. Phys. Sci. 2 (2020) 12–25 25

Kutta methods. Journal of Computational and Applied Mathematics
184(2005)442.

[26] W. Gautschi, Numerical Integration of Ordinary Differential Equa-
tions Based on Trigonometric Polynomials, Numerische Mathematik
3(1961)381.

[27] T. Lyche, Chebyshevian Multistep Methods for Ordinary Differential
Equations, Numerische Math.19(1972)65.

[28] J. P. Coleman, S. C. Duxbury, Mixed collocation methods for y′′ = f (x, y)
. Journal of computational and Applied Mathematics 126(2000)47.

[29] L. Gr. Ixaru , G. Vanden Berghe, H. De Meyer, “Frequency evaluation in
exponentially-fitted algorithms for ODEs”, Journal of Computational and
Applied Mathematics 140 (2002) 423.

[30] G. Vanden Berghe, Ixaru L. Gr., M. Van Daele, Optimal implicit exponen-
tially fitted Runge-Kutta methods. Computer Physics Communications
140 (2001) 346.

[31] G. Vanden Berghe & M. Van Daele, “Exponentially-fitted Numerov meth-
ods”, J. Comp. Appl. Math.200(2007) 140.

[32] T. E. Simos, An Exponentially-Fitted Runge-Kutta Method for the Nu-
merical Integration of Initial Value Problems with Periodic or Oscillating
Solutions. Computer Physics Communications 115 (1998) 1.

[33] T. E. Simos, Exponentially-Fitted Runge-Kutta-Nyström Method for the
Numerical Solution of Initial Value-Problems with Oscillating Solutions.
Applied Mathematics Letters 15 (2002) 217.

[34] Ch. Tsitouras, T. E. Simos, Optimized Runge-Kutta pairs for problem
with oscillatory solutions, Journal of computational and Applied Mathe-
matics 147 (2002) 397.

[35] H. S. Nguyen, R. B. Sidje, N. H. Cong, Analysis of trigonometric implicit
Runge-Kutta methods. Journal of computational and Applied Mathemat-
ics 198 (2007) 187.

[36] N. Senu, M. Suleimon, F. Ismail, M. Othman, A New Diagonally Implicit
Runge-Kutta-Nystrom Method for Periodic IVPs. WSEAS Transactions
on Mathematics 9(2010)679.

[37] S. N. Jator, Trigonometric symmetric boundary value method for oscillat-
ing solutions including the sine-Gordon and Poisson equations. Applied
& Interdisciplinary Mathematics 3(2016)1.

[38] F. F. Ngwane, S. N. Jator, Solving Oscillatory Problems Using a Block
Hybrid Trigonometrically Fitted Method with Two Off-Step Points. Texas
State University. San Marcos, Electronic Journal of Differential Equation
20(2013)119.

[39] F. F. Ngwane, S. N. Jator, Block hybrid method using trigonometric ba-
sis for initial problems with oscillating solutions. Numerical Algorithm
63(2013)713.

[40] S. N. Jator, S. Swindell, R. D. French, Trigonometrically Fitted Block
Numerov Type Method for y

′′
= f

(
x, y, y

′
)

. Numer Algor 62(2013)13.
[41] F. F. Ngwane, S. N. Jator, Trigonometrically-Fitted Second Derivative

Method for Oscillatory Problems. Springer Plus 3(2014)304.
[42] F. F. Ngwane, S. N. Jator, A Family of Trigonometrically Fitted Enright

Second Derivative Methods for Stiff and Oscillatory Initial Value prob-
lems. Journal of Applied Mathematics (2015)1.

[43] R. I. Abdulganiy, O. A. Akinfenwa, S. A. Okunuga, Maximal Order Block
Trigonometrically Fitted Scheme for the Numerical Treatment of Sec-
ond Order Initial Value Problem with Oscillating Solutions. International
Journal of Mathematical Analysis and Optimization(2017)168.

[44] R. I. Abdulganiy, O. A. Akinfenwa, S. A. Okunuga, G. O. Oladimeji, A
robust block hybrid trigonometry method for the numerical integration of
oscillatory second order nonlinear initial value problem. AMSE journals-

AMSE IIETA publication 2017 series 54(2017)497.
[45] R. I. Abdulganiy, O.A. Akinfenwa, S. A. Okunuga, Construction of L

stable second derivative trigonometrically fitted block backward differ-
entiation formula for the solution of oscillatory initial value problems,
African Journal of Science, Technology, Innovation and Development
10(2018)411.

[46] S. N. Jator and E. Agyingi, Block Hybrid k-Step Backward Differentiation
Formulas for Large Stiff Systems. International Journal of Computational
Mathematics. (2014).https://doi.org/ 10.1155/2014/162103.

[47] P. L. Ndukum, Biala, T. A., S. N. Jator, R. B. Adeniyi, On a family of
trigonometrically fitted extended backward differentiation formulas for
stiff and oscillatory initial value problems. Numer Algor, https://doi.org/
10. 1007/s11075-016-0148-1.

[48] S. N. Jator Implicit third derivative Runge-Kutta-Nyström method with
trigonometric coefficients. Numerical Algorithms 70(2015)133.

[49] S. N. Jator, Block third derivative method based on trigonometric
polynomials for periodic initial-value problems. Afrika Matematika 27
(2016)365.

[50] O. Akinfenwa, Third derivative hybrid block integrator for solution of
stiff systems of initial value problems. Afrika Matematika 28(2017)629.

[51] W. H. Enright, Continuous numerical methods for ODEs with defect con-
trol. Journal of Computational and Applied Mathematics 125(2000)159.

[52] P. Onumanyi, U.W. Sirisena, S. N. Jator, “Continuous finite difference
approximations for solving differential equations”, International Journal
of Computer Mathematics 72 (1999) 15.

[53] W. E. Milne. Numerical Solution of Differential Equations. John Wiley
and Sons, New York, NY, USA, 1953.

[54] J. B. Rosser. A Runge-Kutta for all seasons, SIAM Review 9 (1967) 417.
[55] L. F. Shampine, H. A. Watts. Block implicit one-step methods. Math.

Comp. 23 (1969) 731.
[56] M. T. Chu, H. Hamilton, Parallel solution of ODE’s by multi-block meth-

ods. SIAM Journal on Scientific and Statistical Computing 8 (1987) 342.
[57] J. D. Lambert. Computational methods in ordinary differential equations,

John Wiley, New York, 1973.
[58] S. O. Fatunla, Numerical methods for initial value problems in ordinary

differential equation. United Kingdom: Academic Press Inc, (1988).
[59] J. P. Coleman, L. Gr. Ixaru, P-Stability and Exponential Fitting Methods

for y
′′

= f (x, y) , IMAJ. Num. An. 16 (1996) 179.
[60] L. Gr. Ixaru , M. Rizea, H. De Meyer & G. Vanden Berghe, “Weights of

the exponential fitting multistep algorithms for ODEs”, Journal of Com-
putational and Applied Mathematics 132 (2001) 83.

[61] L. Gr. Ixaru, G. Vanden Berghe, H. De Meyer, “Exponentially-fitted Vari-
able two step BDF algorithms for first order ODEs”, Computer Physics
Communication 150 (2003) 116.

[62] G. Vanden Berghe, H. De Meyer, M. Van Daele & T. Van Hecke, “Ex-
ponentially fitted Runge-Kutta methods”, Journal of Computational and
Applied Mathematics 125 (2000) 107.

[63] G. Vanden Berghe, L. Gr. Ixaru & H. De Meyer, “Frequency Determina-
tion and Step-length Control for Exponentially-Fitted Runge-Kutta Meth-
ods”, Journal of Computational and Applied Mathematics 132 (2001) 95.

[64] H. Ramos, “J. Vigo-Aguiar, On the frequency choice in trigonometrically
fitted methods”, Applied Mathematics Letters 23 (2010) 1378.

[65] J. Vigo-Aguiar & H. Ramos, “On the choice of the frequency in trigono-
metrically fitted methods for periodic problems”, Journal of computa-
tional and Applied Mathematics 277 (2015) 94.

25