J. Nig. Soc. Phys. Sci. 4 (2022) 797

Journal of the
Nigerian Society

of Physical
Sciences

A One-Step Block Hybrid Integrator for Solving Fifth Order
Korteweg-de Vries Equations

Olumide O. Olaiyaa, Mark I. Modebeia,b,∗, Saheed A. Belloc

aMathematics Programme, National Mathematical Centre, Abuja, Nigeria
bDepartment of Mathematics, University of Abuja, Nigeria

cNational Water Resources Institute, Kaduna, Nigeria

Abstract

Fifth-order Korteweg-de Vries (KdV) equations, arise in modeling waves phenomena such as the propagation of shallow water waves over a flat
surface, gravity-capillary waves and sound waves in plasmas. In this work, a one-step block hybrid linear multistep method was derived using the
collocation technique, to solve fifth-order KdV models via the Method of Line (MoL). The consistency, stability and convergence of the method
were established. The efficiency of the method can be seen from comparison of the exact solutions of problems and other methods cited from
literature.

DOI:10.46481/jnsps.2022.797

Keywords: Korteweg-de Vries (KdV) equations, Fifth-order PDE, Linear multistep, Block Method, Convergence

Article History :
Received: 03 May 2022
Received in revised form: 12 July 2022
Accepted for publication: 25 July 2022
Published: 19 August 2022

c© 2022 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

In this work, the numerical solution of fifth order PDE of
the form;{

a(x)ytt + b(x)yxxxxx = G(x, t, y, y2, yx, yt, yxx, yxxx),
y(x, b) = g1(x), y(a, t) = g3(t)

(1)

defined on the domain Ω = {(x, t) : a < x < b, c < t < d}, and
the boundary ∂Ω. where a, b, c, d are real numbers.

Formulation of models into PDEs of the type in equation
(1) arise in sciences, for example, fifth-order Korteweg-de Vries
(KdV) type equations, models different wave phenomena which

∗Corresponding author tel. no: +2348031576262
Email address: gmarc.ify@gmail.com (Mark I. Modebei)

includes the propagation of shallow water waves over a flat sur-
face, gravity-capillary waves and sound wave propagation in
plasmas see [1, 2, 3]. Prominent examples includes the ”good”
Boussinesq equation, the Kaup-Kupershmidt equation and an
extended KdV5 equation, [4]. These equations model phe-
nomenon like laser optics, nonlinear dispersive waves in a wide
range of applications, water waves, and plasma: generally, these
are PDEs with higher order spatial derivatives, [4, 5, 6]. Some
of the notable general (linear and nonlinear) fifth order KdV
equations include though not limited to: the Kawahara equa-
tion [4]; the Lax equation [7]; the Caudrey-Dodd-Gibbon (C-D-
G) equation [8]; the fifth-order KdV equation [9]; the Sawada-
Kotera (SK) equation [10]. They take the general form

Ayt + By
2yx + Cyxyxx + Dyyxxx + Eyxxxxx = 0 (2)

1



Olaiya et al. / J. Nig. Soc. Phys. Sci. 4 (2022) 797 2

with appropriate initial-boundary conditions, where A, B, C, D, E
are constants. the form (2) is referred to as the generalized fifth-
order Korteweg-de Vries (GFKdV), [11, 12].

These GFKdV equations models diverse important physi-
cal phenomena. This equation do not just show the movement
of long waves in shallow water under gravity and in a one-
dimensional nonlinear lattice, but at the same time, it is a sig-
nificant mathematical model for magneto-sound propagation in
plasmas [13] and a chain of coupled nonlinear oscillators [14].

Obtaining the exact solution of the GFKdV equation ap-
pears to be unreachable in most cases especially for the nonlin-
ear case, save the special case of solitary waves in [15]. Some
methods which are semi numerical and complete numerical in
literature that were derived for solving fifth order KdV types in-
clude, Modified Variational Iteration Technique [11], Modified
Variational Iteration Algorithm-II (MVIA-II) [12], He’s Varia-
tional Iteration Method, [16], homotopy analysis method, [17],
Group analysis method [18], among others.

In this work, we derive a 1-Step Block Hybrid Integrator
(1SBHI), assembled in block form using collocation technique,
for the solution of fifth-order Korteweg-de Vries (KdV) type
equations using the MOLs approach, see [19, 20, 21]. The
linear multistep method (LMM) derived using the collocation
technique is conveniently used to approximate, with high ac-
curacy, continuous Initial-Boundary Value Problems (IBVPs)
of both the ordinary and partial differential equations [20]. It
has several advantages among which are the fact that they do
not require any staring value to kick-start their implementation,
they are self-starting [22, 23]. They also have the advantage
of producing smaller global errors (at the end of the range of
integration) than those produced by the step-by-step methods
due to the presence of accumulated errors at each step in the
step-by-step method, see [24, 25, 26, 27, 28].

The Method of Lines (MOLs) approach is a very impor-
tant tool used for solving Partial Differential Equations (PDEs),
where a PDE is transformed into a system of Ordinary Differ-
ential Equations (ODEs) whereby appropriate derivatives are
replaced by finite difference approximations. In this work, we
transform the PDEs into a system of ODEs such that only one of
the spatial derivatives is substituted with central difference ap-
proximations. The system that results is solved using 1SHBI.
In the case of this work, we discretize the space variable t with
meshing

∆t =
d − c

M
, tm = c + m∆t, m = 0, 1, . . . M

and then define the vector

y = [y1,1, y1,2, y2,1, . . . , yn,m−1, yn−1,m, . . . , yn−1,m−1]
T

and

G = [G1,1, G1,2, G2,1, . . . , Gn,m−1, Gn−1,m, . . . , Gn−1,m−1]
T

where ym ≈ y(x, tm). The central differences approximation for-
mulae for the first and second derivatives are

yt ≈
ym+1 − ym−1

(2∆t)

ytt ≈
ym+1 − 2ym + ym−1

(∆t)2

The solution y(x, t), of equation (1), where (x, t) is in the rectan-
gle [a, b]×[c, d], is obtained by solving the system of fifth-order
ODEs in the spatial variable x. Here the mesh of x is spaced is:

πN = a = x0 < x1 < .. . < xN−1 < xN = b,

with the step-size h given as

h =
b − a

N
, xn = a + nh, n = 0, 1, . . . , N

Then, the equation (1) has the semidiscretized form

dy5n,m
d x5

≈−
an,m
bn,m

+
1

bn,m

(
ym+1 − 2ym + ym−1

(∆t)2
+ Gn,m

)
(3)

where

Gn,m =
xn, tm, yn,m, y2n,m, dyn,md x , ym+1 − um−1(2∆t) , dy

2
n,m

d x2
,

dy3n,m
d x3


and bn,m , 0 for n = 0, 1, . . . , N and m = 0, 1, . . . , M. The
equivalent form for equation (3) is the following system of fifth
order ODE given by

y(5) = f (x, y, y′, y′′, y′′′, y(4)), a < x < b (4)

subject to the boundary conditions:

y(a) = y0, y’(a) = y
′
0, y”(a) = y

′′
0 ,

y(b) = y1n,m, y
′(b) = y2n,m

where f (x, y, y′, y′′, y′′′, y(4)) = AY + g and A is a k by k matrix
of constants with (k = (N − 1)(M − 1)) which results from the
semi-discretized system of equation (3), expressed in the form
the equation (4) and hence solved using 1SBHI. Here g is a
vector of constants. Note that j in y jn,m, for j = 1, 2, are only
index notation for different yn,m. In what follows, the derivation
of the approximate solution of equation (4) is discussed.

It is noteworthy to see that this technique is easy to de-
rive and implement for solving this class of problem defined
in equation (1).

2. Derivation of the Method

Suppose the approximate solution yn(x) of equation (4) have
the continuous form

y(x) ≈ u(x) =
10∑
j=0

a j x
j (5)

where its fifth derivative is given by

y(v)(x) ≈ u(v)(x) =
10∑
j=5

j( j−1)( j−2)( j−3)( j−4)a j x
j−5(6)

Evaluating equation (5) at the points x = xn+v j , j = 0(1)4 where
v0 = 0, v1 =

1
6 , v2 =

1
3 , v3 =

2
3 , v4 =

5
6 and equation (6) at the

2



Olaiya et al. / J. Nig. Soc. Phys. Sci. 4 (2022) 797 3

points x = xn+v j , j = 0(1)5 where v5 = 1, the following system
of equations are obtained.

u(xn) = yn; u
′(xn) = y

′
n; u

′′(xn) = y
′′
n ; u

′′′(xn) = y
′′′
n ;

u(iv)(xn) = y
(iv)
n ; u

(v)(xn) = fn; u
(v)(xn+v) = fn+v j, j = 1(1)5.

(7)

Solving the above system simultaneously with the use of CAS
in Mathematica 12.0, the coefficients a j, j = 0(1)10 are ob-
tained (though not included here due to their cumbersome terms).
They are in turn substituted into equation (5), to obtain the fol-
lowing continuous method:

yn+vk =
4∑

j=0

αkj h
iy( j)n + h

5
5∑

j=0

βkj fn+v j

∣∣∣∣∣
k=1(1)5

(8)

where α j and β j are shown in the Table 1.
The first derivative of equation (8) is given as

y′n+vk =
4∑

i=1

αih
i−1y(i)n + h

4
5∑

j=0

βi fn+vi

∣∣∣∣∣
k=1(1)5

(9)

here, α j and β j are shown in the Table 2.
The second derivative of equation (8) is given as

y′′n+vk =
4∑

i=2

αih
i−2y(i)n + h

3
5∑

j=0

βi fn+vi

∣∣∣∣∣
k=1(1)5

(10)

where α j and β j are shown in the Table 3
The third derivative of equation (8) is given as

y′′′n+vk =
4∑

i=3

αih
i−3y(i)n + h

2
5∑

j=0

βi fn+vi

∣∣∣∣∣
k=1(1)5

(11)

here α j and β j are shown in the Table 4.
The fourth derivative of equation (8) is given as

y(iv)n+vk = y
(iv)
n + h

5∑
j=0

βi fn+vi

∣∣∣∣∣
k=1(1)5

(12)

where the coefficients α j and β j are shown in the Table 5.

2.1. Characteristics of the method
The formula in equation (8) (and its derivatives in equations

(9)-(12)) is a continuous schemes associated with the linear dif-
ferential operator defined by

L[z(x + vkh); h]

=

4∑
i=0

αkj h
jz( j)(x) −

5∑
j=0

h5βiz
(5)(x + v jh)

∣∣∣∣∣
k=1(1)5

(13)

Expanding equation (13) in Taylor series, the constants C′i s
written as linear combination of the derivatives of z(x) up to
p + 1 derivative as

L[z(x); h]

= C0z(x)+C1hz
′(x)+C2h

2z′′(x)+· · ·+C ph
pz(p)(x)+O(h( p+1))

(14)

Figure 1. The region of stability

where p is the order of equation (13) and consequently the
formula (13). By definition, the LMM (8) is of order p if
C0 = C1 = C2 = · · · = C p+4 = 0, and C p+5 , 0 in which

L[z(x); h] = C p+5h
p+5z(p+5)(x) + O(h(p+6)) (15)

In this case, C p+5 is the principal error constant, see [29].

Definition 2.1. A linear difference operator L[z(x); h] of order
p is consistent if p > 1.

Thus, the order of the formula in equation (8) is p = 6 with the
following error constants

C p+5 = (−6.87857 × 10
−13,−2.36466 × 10−11,

− 5.40035 × 10−10,−1.47456 × 10−9,−3.4408 × 10−9)T

The order of the formulas in equations (9)-(12) is p = 6 and
their error constants are obtained similarly.

2.2. Zero-stability and convergence

A numerical method in equations (8)-(12) is zero-stable pro-
vided as h → 0, the solutions is bounded.

Following [20], the block method in (8)-(12) may be rewrit-
ten in a matrix form as

A0Yµ = A1Yµ−1 + h
5(Fµ + Fµ−1) (16)

where where

Yµ =
(
Y 0µ, Y

1
µ, . . . , Y

5
µ

)T
, Yµ−1 =

(
Y 0µ−1, Y

1
µ−1, . . . , Y

5
µ−1

)T

3



Olaiya et al. / J. Nig. Soc. Phys. Sci. 4 (2022) 797 4

Table 1. Coefficients of α j, j = 0(1)4, and β j, j = 0(1)5 in equation (8)

k α0 α1 α2 α3 α4 β0 β1 β2 β3 β4 β5
1 1 16

1
72

1
1296

1
31104 86197η1 58392η1 −30815η1 14415η1 −9112η1 1883η1

2 1 13
1
18

1
162

1
1944 7963η2 10528η2 −4365η2 1985η2 −1248η2 257η2

3 1 23
2
9

4
81

2
243 622η3 1392η3 −215η3 180η3 −112η3 23η3

4 1 56
25
72

125
1296

625
31104 6685η4 16600η4 −375η4 2375η4 −1368η4 275η4

5 1 1 12
1
6

1
24 401η5 1056η5 105η5 195η5 −96η5 19η5

η01 =
1

112870195200 ; η
0
2 =

1
440899200 ; η

0
3 =

2
3444525 ; η

0
4 =

625
4514807808 ; η

0
5 =

1
201600

Table 2. Coefficients of α j, j = 1(1)4, and β j, j = 0(1)5 in equation (9)

k α1 α2 α3 α4 β0 β1 β2 β3 β4 β5
1 1 16

1
18

1
162 100795η

1
1 82832

1
1 −42175

1
1 19525

1
1 −12320

1
1 2543

1
1

2 1 13
1
18

1
162 2209η

1
2 3518η

1
2 −1300η

1
2 595η

1
2 −374η

1
2 77η

1
2

3 1 23
2
9

4
81 1316η

1
3 3328η

1
3 −140η

1
3 425η

1
3 −256η

1
3 52η

1
3

4 1 56
25
72

125
1296 7027η

1
4 19040η

1
4 2225η

1
4 3325η

1
4 −1712η

1
4 335η

1
4

5 1 1 12
1
6 35η

1
5 98η

1
5 25η

1
5 25η

1
5 −10η

1
5 2η

1
5

η11 =
1

4702924800 ; η
1
2 =

1
9185400 ; η

1
3 =

2
1148175 ; η

1
4 =

125
188116992 ; η

1
5 =

1
4200 ;

Table 3. Coefficients of α j, j = 2, 3, 4, and β j, j = 0(1)5 in equation (10)

k α2 α3 α4 β0 β1 β2 β3 β4 β5
1 1 16

1
72 40501η

2
1 42484η

2
1 −20455η

2
1 9335η

2
1 −5876η

2
1 0

2 1 13
1

18 1646η
2
2 3304η

2
2 −985η

2
2 470η

2
2 −296η

2
2 61η

2
2

3 1 13
1
9 467η

2
3 1328η

2
3 190η

2
3 205η

2
3 −112η

2
3 22η

2
3

4 1 56
25
72 497η

2
4 1444η

2
4 485η

2
4 395η

2
4 −164η

2
4 31η

2
4

5 1 1 12 222η
2
5 648η

2
5 315η

2
5 270η

2
5 −72η

2
5 17η

2
5

η21 =
1

87091200 ; η
2
2 =

1
680400 ; η

2
3 =

1
42525 ; η

2
4 =

125
3483648 ; η

2
5 =

1
8400 ;

Table 4. Coefficients of α j, j = 3, 4, and β j, j = 0(1)5 in equation (11)

k α3 α4 β0 β1 β2 β3 β4 β5
1 1 16 5561η

3
1 37504η

3
1 −16325η

3
1 7295η

3
1 −4576η

3
1 941η

3
1

2 1 13 1843η
3
2 5032η

3
2 −830η

3
2 515η

3
2 −328η

3
2 68η

3
2

3 1 23 254η
3
3 776η

3
3 395η

3
3 235η

3
3 −104η

3
3 19η

3
3

4 1 56 269η
3
4 800η

3
4 575η

3
4 475η

3
4 −128η

3
4 25η

3
4

5 1 1 239η34 696η
3
4 600η

3
4 555η

3
4 −24η

3
4 34η

3
4

η31 =
1

3628800 ; η
3
2 =

1
113400 ; η

3
3 =

2
14175 ; η

3
4 =

25
145152 ; η

3
5 =

1
4200 ;

Table 5. Coefficients of α4 and β j, j = 0(1)5 in equation (12)

.

k α4 β0 β1 β2 β3 β4 β5
1 1 4991η41 12824η

4
1 −4365η

4
1 1885η

4
1 −1176η

4
1 241η

4
1

2 1 287η42 1248η
4
2 245η

4
2 45η

4
2 −32η

4
2 7η

4
2

3 1 43η43 112η
4
3 180η

4
3 155η

4
3 −48η

4
3 8η

4
3

4 1 43η44 120η
4
4 175η

4
4 225η

4
4 8η

4
4 5η

4
4

5 1 13η45 32η
4
5 55η

4
5 55η

4
5 32η

4
5 13η

4
5

η41 =
1

86400 ; η
4
2 =

1
5400 ; η

4
3 =

1
675 ; η

4
4 =

5
3456 ; η

4
5 =

1
200 ;

Y 0µ = (yn+ 16 , yn+ 13 , yn+ 23 , yn+ 56 , yn+1)
T

...

Y 5µ = (y
(v)
n+ 16

, y(v)
n+ 13

, y(v)
n+ 23

, y(v)
n+ 56

, y(v)n+1)
T

Y 0µ−1 = (yn− 56 , yn− 23 , yn− 13 , yn− 16 , yn)
T

...

Y 5µ−1 = (y
(v)
n− 56

, y(v)
n− 23

, y(v)
n− 13

, y(v)
n− 16

, y(v)n )
T

Fµ = ( fn+ 16 , fn+ 13 , fn+ 23 , fn+ 56 , fn+1)
T

Fµ−1 = ( fn− 56 , fn− 23 , fn−13 , fn− 16 , fn)
T

Hence, as h → 0 the method in equation (16) becomes

A0Yµ = A1Yµ−1 (17)

where A0 is the identity matrix of order 25, A0 = I25, and A1 is

4



Olaiya et al. / J. Nig. Soc. Phys. Sci. 4 (2022) 797 5

Table 6. Absolute Error for Problem 1
x = 1.0 x = 1.5

t 1SHBI mVIA-I 1SHBI mVIA-I
0.01 0.00000 0.00000 2.14572 × 10−18 0.00000
0.02 1.11257 × 10−18 8.88178 × 10−16 2.22141 × 10−18 8.88178 × 10−16

0.03 8.47894 × 10−17 1.19904 × 10−14 4.71125 × 10−17 1.95399 × 10−14

0.04 5.41381 × 10−17 8.79296 × 10−14 3.12412 × 10−17 1.43884 × 10−13

0.05 9.12414 × 10−16 4.19229 × 10−13 1.49787 × 10−17 6.90114 × 10−13

Comparing absolute errors of 1SHBI with the mVIA-I in [11] for Example 1.

Table 7. Error for Example 2
x yex yapprox Error
0.12566 -0.07426755 -0.07426785 2.9720E-7
0.25132 0.05130957 0.05130487 4.7001E-6
0.31415 0.11392287 0.11391146 1.1407E-5
0.43982 0.23757397 0.23753066 4.3304E-5
0.50265 0.29812884 0.29805541 7.3430E-5
0.62831 0.41551891 0.41534181 1.7710E-4
0.75398 0.52644076 0.52607803 3.6273E-4
0.81681 0.57893577 0.57843927 4.9650E-4
0.94247 0.67698465 0.67611561 8.6903E-4
1.00530 0.72216309 0.72104522 1.1178E-3

0 ≤ x ≤ 2π, t = 1, h = 0.0628319.

Table 8. Maximum Error for Problem 2 at t = 1
Method N L1 L2 L∞

DGFEM 10 0.11E-1 0.12E-1 0.17
20 0.39E-2 0.43E-2 0.61E-2
40 0.12E-3 0.14E-3 0.19E-3
80 0.36E-5 0.40E-5 0.57E-5

1SHBI 10 1.52E-6 3.23E-6 1.17E-7
20 2.15E-8 3.24E-6 3.21E-8
40 5.18E-8 2.79E-9 5.27E-8
80 2.32E-10 4.71E-10 8.11E-10

DGFEM is Discontinuous Galerkin Finite Element Method
developed in [32]. The p-norm in Table 8 is given as

Lp = ||yi − y(xi)||p =

 inf∑
i=1

|yi − y(xi)|
p


1
p

L∞ = ||yi − y(xi)||∞ = max
1≤i≤N

{|yi − y(xi)|}

a 25 × 25 matrix given by
A11

A22
A33

A44
A55


with the A11, . . ., A55 being 5 × 5 matrices respectively, given

yexact

ynumerical

Error

Figure 2. The shaded regions shows the surface plots for the numerical solution,
the analytic solution and the residue (error) for Example 1.

by

Aii =


0 1 0 0 0
0 1 0 0 0
0 1 0 0 0
0 1 0 0 0
0 1 0 0 0

 , i = 1, . . . , 5
5



Olaiya et al. / J. Nig. Soc. Phys. Sci. 4 (2022) 797 6

Table 9. Numerical Results with t = 1 for Example 3
x Exact Approximate Error

0.1 1.999999999998E-10 1.999999999998E-10 0.000000000000000
0.2 1.999999999992E-10 1.999999999991E-10 1.292469707114105E-25
0.3 1.999999999982E-10 1.999999999981E-10 3.618915179919496E-25
0.4 1.999999999968E-10 1.999999999967E-10 6.720842476993354E-25
0.5 1.999999999950E-10 1.999999999949E-10 1.214921524687259E-24
0.6 1.999999999928E-10 1.999999999927E-10 1.757758801675183E-24
0.7 1.999999999902E-10 1.999999999901E-10 2.429843049374518E-24
0.8 1.999999999872E-10 1.999999999871E-10 2.998529720504725E-24
0.9 1.999999999838E-10 1.999999999837E-10 3.127776691216136E-24
1.0 1.999999999800E-10 1.999999999799E-10 2.403993655232236E-24

Table 10. Maximum Error for Example 3 at t = 1
Method GA with RK MQ with RK IMQ with RK 1SHBI

L2 1.9774E-21 5.9692E-12 1.9774E-21 5.96847E-24
L∞ 1.2705E-21 3.2896E-12 1.2705E-21 3.12778E-24

The matrices Aii have the characteristic polynomial |A j j−λI5| =
0, for j = 1, . . . , 5, that is, λ4(λ − 1) = 0. The characteristic
polynomial has the root λr = 0, for r = 1, . . . , 4 and λ5 = 1.
The method is zero-stable, since the roots of the characteristic
polynomial have one unity while others are zero (see [29]). For
convergence, the following theorem apply.

Theorem 2.1. [30]. A linear multistep method is convergent
provided it is consistent (having order p ≥ 1) and zero-stable.

The derived method has p = 6, and equally zero-stable and
hence it is convergent.

2.3. Linear stability analysis
The linear stability for any given h > 0 is concern with

the behaviour of the underlined problem not just the numerical
method. Here, the stability properties for the intending numer-
ical method is analyzed by considering the linear test equation
for µ > 0 of the form

y(v) = −µ5y (18)

the test in equation (18) is appropriate since it has a bounded
solution as µ → 0. Here, µ is the frequency. set δ = µh and then
applying to the test equation in (18), and after some manipula-
tions, the following equation results

Zµ = M(δ)Zµ−1, δ = µh (19)

where M(δ) =
A0 + δB0

(A1 −δB1)
is the required stability matrix.

The amplification matrix in equation (19) results on the domi-
nant eigenvalues given as

M(δdominant)

=
6(−2)4/531/5

(
−1 − 4δ1/6 + 5δ1/3 − 5δ2/3 + 4δ5/6 −δ

)1/5
(
11 − 56δ1/6 + 285δ1/3 + 285δ2/3 − 56δ5/6 + 11δ

)1/5
(20)

having the stability region in the Figure 1.
These roots (containing the real and imaginary parts) are

plotted and as shown in Figure 1. Here we study the bounded-
ness of the solution under consideration through the eigenval-
ues of the stability matrix M(δ) for which these eigenvalues δ is
such that |δ| < 1 for the derived method to be stable. If δ is real,
then the absolute stability region is reduced to a real interval
which is called an interval of stability, see [31]. Figure 1 shows
the stability region for the proposed method, with the stability
interval (0, 8.7893).

2.4. Implementation

The system formed by the semi-discretization given by equa-
tion (4) is solved using the unified block scheme formed by
equations (8)-(12) with 5N equations and 5N unknowns solve
simultaneously, whose solution provides a set of approximate
values of (1) using codes written in Mathematica 12.0, enhanced
by the feature NSolve[] for linear problems while nonlinear
problems were solved by Newton’s method enhanced by the
feature FindRoot[]. Following [20], the following algorithm
summarizes the computational procedures.
We begin by noting that the solution of the problem in equation
(1) is sought in the subintervals DN = a = x0 < x1 < .. . <
xN = b and DM = c = t0 < t1 < .. . < tM = d, where h =

b−a
N ,

k = d−cM are constant step-size.
Step 1: Use the block method setting n = 0, m = 0, to obtain
V1 on the rectangle [y0,m, y1,m](n,m)∈[a,b]×[c,d], similarly, for, n = 1
so that V2 is obtained on the rectangle [y1,m, y2,m](n,m)∈[a,b]×[c,d],
and on the rectangle
[y3,m, y4,m](n,m)∈[a,b]×[c,d] . . . [yn=N−1,m, yn=N,m](n,m)∈[a,b]×[c,d] for m =
0, 1, 2, . . . , M we thus obtain V3, V4, . . . , VD.
Step 2: Solve the unified block given by the system
V1

⋃
V2

⋃
, · · · , VD obtained in step 1. Step 3: The solution of

equation (1) is approximated by the solutions in step 2 as
yn,m = [y(x1, t); y(x2, t) . . . y(xn, t)]T for t > 0, n = 1, 2, . . . N.

6



Olaiya et al. / J. Nig. Soc. Phys. Sci. 4 (2022) 797 7

yexact

ynumerical

Error

Figure 3. Surface plots showing the numerical solution, the analytic solution
and the residue for Example 2.

3. Numerical Examples

In this section, we test the efficiency of the derived method
and thus compare only with the exact solution of the given prob-
lem extracted from literature.

Example 1.. Consider the fifth order KdV equation,

yt + yyx − yyxxx + yxxxxx = 0
y(x, 0) = ex,

}
(21)

with exact solution y(x, t) = ex−t.

After discretization of the ”t” variable, we obtain,

ym+1 − ym−1
(2∆t)

+ ym
dym
d x
− ym

d3ym
d x3

+
d5ym
d x5

= gm,

x ∈ [0, 1] m = 1, . . . M − 1 (22)

Figure 4. Surface plots for the numerical solution for Example 4

where ∆t, tm, m = 0, 1, . . . , M, y, ym(x) ≈ y(x, tm), and gm = 0,
expressed as

y(v) = f (x, y, y′, y′′′) = Ay + g (23)

A is a k by k matrix.
To show the efficiency of 1SHB1, the above Table 6 shows

the absolute error for different values of t the values x = 1.0 and
1.5 respectively. It can be observe that the 1SHBI outperformed
mVIA-I for the example considered. The Figure 2 shows a good
agreement of the numerical and exact solutions.

Example 2.. Consider a time dependent biharmonic fifth order
equation discussed in [32].

yt + yx + yxxxxx = 0, 0 ≤ x ≤ 2π, t ≥ 0
y(x, 0) = sin (x),
y(0, t) = y(2π, t)

 (24)
The analytic solution for this problem is y(x, t) = sin(x − 2t).

Upon discretization of the time variable, we obtain,

ym+1 − ym−1
(2∆t)

+
dym
d x

+
d5ym
d x5

= gm,

0 ≤ x ≤ 1, m = 1, . . . M − 1 (25)

7



Olaiya et al. / J. Nig. Soc. Phys. Sci. 4 (2022) 797 8

where ∆t, tm, m = 0, 1, . . . , M, y, ym(x) ≈ y(x, tm), g and gm are
expressed in the form

y(v) = f (x, y, y′, y′′, y′′′, y(iv)) = Ay + g (26)

A is a k × k matrix (k = (N − 1)(M − 1)) and gm = 0.
Table 7 compares the errors obtained for Example 2. Figure 3
shows the representation of the numerical solution, analytical
solution and residue (errors), for Example 2.
Table 8, for different subinterval N, shows the maximum errors
L∞ and the Lp with p = 1, 2 − norm obtained, in compari-
son with the 1SHBI and the method in [32]. This shows the
superiority of the 1SHBI.

Example 3.. Consider the nonlinear fifth order KdV equation
called Lax’s cases of generalized KdV equation discussed [33].

yt + 45y2yx + 15yxyxx + 15yyxxx + yxxxxx = 0,
0 ≤ x ≤ 1, t ≥ 0

y(x, 0) = 2k2 sec h2(k(x − x0))

 (27)
The exact solution for this problem is
y(x, t) = 2k2 sec h2(k(x − 16k4t − x0))

Upon discretization of the time variable, we obtain,

ym+1 − ym−1
(2∆t)

+ 45ymym
dym
d x

+ 15
dym
d x

d2ym
d x2

+ 15ym
d3ym
d x3

+
d5ym
d x5

= gm, 0 ≤ x ≤ 1, m = 1, . . . M − 1

(28)

where ∆t, tm, m = 0, 1, . . . , M, y, ym(x) ≈ y(x, tm), g and gm are
expressed in the form

y(v) = f (x, y, y′, y′′, y′′′, y(iv)) = Ay + g (29)

A is as expressed in Example 1 and gm = 0.
In [33], the authors presented a numerical solution of equa-

tion (27) using a meshless method of lines. This method uses
Meshless, Multiquadric (MQ), Inverse multiquadric (IMQ) and
Gaussian (GA) as Radial basis function (RBF) for spatial deriva-
tives and Runge-Kutta method as a time integrator.
The constant step-size used for Table 9 and Table 10 is h = 0.1.

Figure 4 shows the graphical representation with surface
plot for the numerical solution (shaded region), analytical so-
lution (shaded region) and residue or error (shaded region), for
Example 3.

The method presented in [33] for the numerical solution of
Example 4 uses a meshless method of lines. This method uses
Meshless, Multiquadric (MQ), Inverse multiquadric (IMQ) and
Gaussian (GA) as Radial basis function (RBF) for spatial deriva-
tives and Runge-Kutta method as a time integrator.

4. Conclusion

In this work, a one step hybrid block integrator (ISHBI) was
derived consisting of continuous linear multistep methods. This

was used to solve certain fifth- order PDEs subject to appro-
priate initial-boundary conditions. To show the robustness of
method derived, fifth-order KdV PDEs were solved were trans-
formed into a system of fifth-order ODEs using the method of
lines. From the numerical experiments performed, It’s evident
that the 1SHBI performed well in terms of accuracy as com-
pared to exact solutions of the problems considered in litera-
ture.

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] M. Kazeminia, S. Soleimani-Amiri & S.A. Zahedi, “Exact and numerical
solutions for nonlinear higher order modified KdV equations by using
variational iteration method”, Adv. Stud. Theor. Phys. 4 (2010) 437.

[2] G. Amit, J. Singh & D. Kumar, “A reliable algorithm for KdV equations
arising in warm plasma”, Nonlin. Eng. 5 (1) (2016) 7.

[3] D. Kumar, J. Singh & S. Kumar, “A fractional model of Navier-Stokes
equation arising in unsteady flow of a viscous fluid“, J. Assoc. Arab Univ.
Basic Appl. Sci. 17 (2015) 14.

[4] G. Amit, J. Singh & D. Kumar, “Numerical simulation of fifth order KdV
equations occurring in magneto-acoustic waves”, Ain. Shams. Eng. J.
(2017), http://dx.doi.org/10.1016/j.asej.2017.03.004.

[5] T. Kakutani. & H. Ono, “Weak non-linear hydromagnetic waves in a cold
collision free plasma”, J. Phys. Soc. Jpn. 26 (130) (1969) 5.

[6] P. Saucez, A. Vande Wouwer, W.E. Schiesser & P. Zegeling, “Method of
lines study of nonlinear dispersive waves”, Journal of Computational and
Applied Mathematics 168 (2004) 413.

[7] P. D. Lax, “Integrals of nonlinear equations of evolution and solitary
waves”, Commun. Pure Appl. Math. 21 (1968) 467.

[8] P. Caudrey, R. Dodd, & J. Gibbon, “A new hierarchy of Korteweg-de
Vries equations”, Proc. R. Soc. Lond. A Math. Phys. Sci. 351 (1976) 407.

[9] S. Handibag & B. Karande, “Existence the solutions of some fifth-order
KdV equation by Laplace decomposition method”, Am. J. Comput. Math.
3 (2013) 80, https://doi.org/10.4236/ajcm.2013.31013.

[10] K. Sawada & T. Kotera, “A method for finding n-soliton solutions of the
KdV equation and KdV-like equation”, Prog. Theor. Phys. 51 (5) (1355).

[11] H. Ahmad, T. A. Khan, P. S. Stanimirovic & I. Ahmad, “Modified
Variational Iteration Technique for the Numerical Solution of Fifth Or-
der KdV-type Equations”, J. Appl. Comput. Mech. 6 (2020) 1220,
https://doi.org/10.22055/JACM.2020.33305.2197

[12] H. Ahmad, T. A. Khan, & Shao-Wen Yao, “An efficient approach for
the numerical, solution of fifth-order KdV equations”, De Gruyter Open
Mathematics 18 (2020) 738.

[13] Kawahara, T., “Oscillatory solitary waves in dispersive media”, Journal
of the Physical Society of Japan 33 (1972) 260.

[14] K. A. Gorshkov, L. A. Ostrovsky, V. V. Papko & A. S. Pikovsky, “On the
existence of stationary multisolitons”, Physics Letters 74 (1979) 177.

[15] Y. Yamamoto, & E. I. Takizawa, “On a solution on non-linear time-
evolution equation of fifth order”, Journal of the Physical Society of Japan
50 (1981) 1421.

[16] M. O. Miansari, M. E. Miansari, A. Barari & D. D. Ganji, “Application of
He’s Variational Iteration Method to nonlinear Helmholtz and fifth-order
KdV equations”, J. Appl. Math. Stat. Inform. 5 (2009) 5.

[17] S. Abbasbandy & F. S. Zakaria, “Soliton solutions for the fifth-order
KdV equation with the homotopy analysis method”, Nonlinear Dyn. 2008
(2008) 83.

[18] G. Wang, A. H. Kara, K. Fakhar, J. V. Guzman & A. Biswas, “Group
analysis, exact solutions and conservation laws of a generalized fifth order
KdV equation”, Chaos, Solitons Fractals 86 (2016) 8.

[19] T. A. Biala & S. N. Jator, “Block Unification Algorithm for 2D and 3D
Elliptic PDES”, Journal of the Nig. Math. Soc. 36 (2017) 319.

8



Olaiya et al. / J. Nig. Soc. Phys. Sci. 4 (2022) 797 9

[20] M. I. Modebei, R. B. Adeniyi & S. N. Jator, “Numerical Approximations
of Fourth-order PDEs Using Block Unification Method”, Journal of the
Nigerian Mathematical Society 39 (2020) 47.

[21] M. I. Modebei & R. B. Adeniyi, “A six-step Block Unification Integra-
tor for Numerical Solution of Fourth Order Boundary Value Problems”,
General Letters in Mathematics 5 (2018) 71.

[22] O. O. Olaiya, R. A. Azeez & M. I. Modebei, “Efficient Hybrid Block
Method For The Numerical Solution Of Second-order Partial Differential
Problems via the Method of Lines”, Journal of the Nigerian Society of
Physical Sciences 3 (2021) 26.

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

[24] M. I. Modebei, “Optimized hybrid block integrator for Numerical solu-
tion of general third order Ordinary differential equations”, Journal of the
Nigerian MAthematical Society 41 (1) (2022) 49.

[25] M. I. Modebei, O. O. Olaiya & A. C. Onyekonwu, “A 3-step fourth deriva-
tives method for numerical integration of third order ordinary differen-
tial equations”, International Journal of Mathematical Analysis and Opti-
mization: Theory and Applications 7 (2021) 32.

[26] 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.
[27] V. O. Atabo & S. O. Adee, “A New Special 15-Step Block Method for

Solving General Fourth Order Ordinary Differential Equations”, Journal
of the Nigerian Society of Physical Sciences 3 (2021) 308.

[28] M. Kida, S. Adamu, O. O. Aduroja & T. P. Pantuvo, “Numerical Solution
of Stiff and Oscillatory Problems using Third Derivative Trigonometri-
cally Fitted Block Method”, Journal of the Nigerian Society of Physical
Sciences 4 (2022) 34.

[29] J. D. Lambert, Computational Methods in Ordinary Differential Equa-
tions, John Wiley, New York (1973).

[30] P. Henrici, Discrete Variable Methods for ODEs, John Wiley, New York,
USA (1962).

[31] M. A. Rufai & H. Ramos, “Numerical solution of second-order singu-
lar problems arising in astrophysics by combining a pair of one-step
hybrid block Nyström methods”, Astrophys Space Sci. 365 (2020) 96,
https://doi.org/10.1007/s10509-020-03811-8

[32] Y. Cheng & C Shu, “A Discontinuous Galerkin Finite Element Method for
Time Dependent Partial Differential Equations with Higher Order Deriva-
tives”, Mathematics of Computation 77 (2008) 699.

[33] S. T. Mohyud-Din, E. Negahdary & M. Usman, “A Meshless Numeri-
cal Solution of the Family of Generalized Fifth-order Korteweg-de Vries
Equations”, Mathematics Faculty Publications Paper 10 (2012), htt p :
//ecommons.udayton.edu/mth f acpub/10.

9