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

Journal of the
Nigerian Society

of Physical
Sciences

Approximate solution of space fractional order diffusion
equations by Gegenbauer collocation and compact finite

difference scheme

K. Issa∗, A. S. Olorunnisola, T. O. Aliu, A. D. Adeshola

Department of Mathematics and Statistics, Kwara State University, Malete, Kwara State, Nigeria.

Abstract

In this paper, approximation of space fractional order diffusion equation are considered using compact finite difference technique to discretize the
time derivative, which was then approximated via shifted Gegenbauer polynomials using zeros of (N − 1) degree shifted Gegenbauer polynomial
as collocation points. The important feature in this approach is that it reduces the problems to algebraic linear system of equations together with
the boundary conditions gives (N + 1) linear equations. Some theorems are given to establish the convergence and the stability of the proposed
method. To validate the efficiency and the accuracy of the method, obtained results are compared with the existing results in the literature. The
graphical representation are also displayed for various values of β− Gegenbauer polynomials. It can be observe in the tables of the results and
figures that the proposed method performs better than the existing one in the literature.

DOI:10.46481/jnsps.2023.1368

Keywords: Caputo fractional derivative; Stability; fractional diffusion equation; Gegenbauer polynomials; Compact finite difference method
(CFDM).

Article History :
Received: 18 February 2023
Received in revised form: 26 April 2023
Accepted for publication: 10 May 2023
Published: 22 May 2023

© 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: Oluwatimilehin Oluwadare

1. Introduction

In the last decade, significant interest has been given to the
study of fractional calculus due to its application in applied
Mathematics. Some of its application can be traced to modeling
of some physical phenomena such as communication, robotics,
transportation systems, finance, signal processing, genetic algo-
rithms and the damping visco-elasticity as discussed in [1–4].
In recent time, the theory and applications of fractional calcu-

∗Corresponding author tel. no: +2348036554437
Email address: issa.kazeem@kwasu.edu.ng (K. Issa )

lus have expanded enormously, its first appearance was in a let-
ter written by Gottfried Wilhelm Leibniz in 1695, as discussed
in [5]. One of the popular areas of fractional calculus is the
space fractional order diffusion equations (SFODEs). Several
researchers have studied SFODEs and applied different meth-
ods to find the approximate solution, [6, 7] employed finite dif-
ference method (FDM) to discretize the time derivative, then
used Gegenbauer polynomial as approximating polynomial,
[8] developed tau method approach to find an approximant to
SFODEs, [9] employed second kind shifted Chebyshev polyno-
mials for solving SFODEs, [10] proposed generalized shifted
Chebyshev polynomials for solving space fractional optimal

1



Issa et al. / J. Nig. Soc. Phys. Sci. 5 (2023) 1368 2

control problems, new aspects of fractional Biswas-Milovic
model with Mittag-Leffler law was reported in [11]. [12] em-
ployed Riesz derivative to solve high-order solvers. [13] solved
fractional Rayleigh-Stokes model in their paper titled ”Numer-
ical solution of fractional Rayleigh-stokes model arising in a
heated generalized second grade fluid”. [14] investigate non-
linear variable order fractional reaction-diffusion equation with
Mittag-Leffler kernel. [15] developed lines and splines meth-
ods for solving SFODEs. [16] employed spline method com-
bined with Richardson extrapolation to solve SFODEs, [17]
discussed shifted Legendre spectral method for fractional-order
multi-point boundary value problems, [18] analysis the initial-
boundary value problems for the linear time-fractional diffu-
sion equations with a uniformly elliptic spatial differential op-
erator of the second order and the Caputo type time-fractional
derivative acting in the fractional Sobolev spaces, [19] proposed
spatio-temporal homogenization function method for solving
1D fourth-order fractional diffusion-wave equations. [20] in-
vestigate Caputo-Fabrizio and Fractal fractional derivative op-
erator in the of studied HIV/AIDS fractional-order model is
studied, [21] studied the fractional-order COVID-19 epidemic
model using Laplace homotopy analysis method and [22] in-
vestigated the numerical solution of SFODEs using Chebyshev
collocation method of the fourth kind and compact finite differ-
ence scheme. The main aim of this paper is to extend the work
reported in [22] by introducing shifted Gegenbauer polynomi-
als as approximating polynomial UN (x, T ) for solving SFODEs
using shifted Gegenbauer collocation method to find the so-
lution, since its solution generalizes the results of some other
orthogonal polynomials such as Legendre, Jacobi polynomi-
als Pα,βn (x) of certain kinds with α = β and shifted Chebyshev
polynomials of certain kinds. In this paper, we want to em-
ploy compact finite difference method (CFDM) to discretize the
time derivative, then the discretized SFODEs will be collocated
at (N − 1) points. These equations together with the boundary
conditions generate (N + 1) system of linear equations.
This paper is organized as follows. In section 2, we define
Caputo fractional order derivative, review some notable or-
thogonal polynomials, state the problem under consideration
and some useful lemmas from the literature are recalled (see
2.1,2.2,2.3). Shifted Gegenbauer interms of differential opera-
tor is presented in section 3, we present the formulation of the
scheme in section 3.1. We analyze the convergence and the
stability of the proposed method in section 3.2, in section 4 nu-
merical examples are presented with computational results and
graphical representation to show the efficiency and the accu-
racy of the proposed method, concluding remarks are given in
section 5.

2. Preliminaries

2.1. Caputo fractional order derivative operator
The Caputo fractional derivative operator Dω of order ω is

given as

Dω f (t) =
1

Γ(n −ω)

∫ t
0

f (i)(t)
(t − x)ω+1−i

d x, i − 1 < ω < i, i ∈ N
(1)

where ω > 0 is the order of the derivative. The linearity prop-
erty of the Caputo fractional derivative is:

Dω(λ f (t) + µg(t)) = λDω f (t) + µDωg(t), (2)

where, λ & µ are constants. The following results are obtained:

Dωt j =


0, for j ∈ N0 , j < dωe

Γ( j+1)
Γ( j+1−ω) t

j−ω, for j ∈ N0 , j ≥ dωe
(3)

where dωe is the smallest integer greater than or equal to ω.

2.2. Some of orthogonal polynomials
Some of the notable orthogonal polynomials ψ(t) are de-

fined here:
ψ(t) is an orthogonal polynomial with respect to the weight
function ω(t) in an interval [a, b], if the inner product of ψ(t)
is zero.
That is,

〈ψm(t),ψn(t)〉 =
∫ b

a
ω(t)ψm(t)ψn(t)dt =


0, m 6= n

λn, m = n
(4)

Some of these orthogonal polynomials ψi(t) are:

2.2.1. Legendre polynomials
An orthogonal polynomial ψm(t) defined in an interval

[−1, 1] is said to be a Legendre polynomial if the weight func-
tion ω(t) = 1. The recurrence relation of a legendre polynomial
is given as:

Pk+1(u) =
2k + 1
k + 1

uPk(u) −
k

k + 1
Pk−1(u), k ≥ 1, (5)

with P0(u) = 1, P1(u) = u.

2.2.2. Chebyshev polynomials
The prominent one listed with their respective weight func-

tions ω(t) in the interval [−1, 1] are found in [23] as:

ψm =

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

Tm(t) = cos(mt), ω(t) =
1

√
1 − t2

,

Um =
sin (m + 1) t

sin (t)
, ω(t) =

√
1 − t2,

Vm =
cos

(
m + 12

)
t

cos
(

t
2

) , ω(t) = √ 1 + t
1 − t

,

Wm =
sin

(
m + 12

)
t

sin
(

t
2

) , ω(t) = √ 1 − t
1 + t

.

(6)

2



Issa et al. / J. Nig. Soc. Phys. Sci. 5 (2023) 1368 3

The inner product for the third and fourth kinds Chebyshev
polynomials in the interval [−1, 1] are defined as:

〈ψm,ψn〉 =

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

〈Vm, Vn〉 =
∫ 1
−1

√
1 + t
1 − t

Vm(t)Vn(t)dt =

0, m 6= nπ, m = n ,
〈Wm, Wn〉 =

∫ 1
−1

√
1 − t
1 + t

Wm(t)Wn(t)dt =

0, m 6= nπ, m = n .
(7)

while the fifth kind was reported in [14, 24].

2.2.3. Shifted Gegenbauer polynomials
Gegenbauer polynomials C(β)i (t) defined in the interval

[−1, 1] with respect to weight function ω(t) = (1 − t2)(β−
1
2 ) can

be determined using

C(β)i (t) =
i∑

j=0

(−1) jΓ(2β + 2i − j)Γ(β + 12 )

(i − j)! Γ(2β)Γ( j + 1)Γ(i − j + β + 12 )
ti− j. (8)

The recurrence relation is given as

C(β)i (t) =
1
i

[
2(i + β− 1)tC(β)i−1(t) − (i + 2β− 2)C

(β)
i−2(t)

]
, i ≥ 2,

(9)
where Cβ0 (t) = 1, C

β
1 (t) = 2βt. To use this polynomial in the

interval t ∈ [a, b] a shifted Gegenbauer polynomial is defined
by introducing the variable λ = 2t−(a+b)b−a . Hence, the shifted
Gegenbauer polynomial in t is obtained as:

C(β)∗i (t) =
1
i

[
2(i + β− 1)

(
2t − (a + b)

b − a

)
C(β)∗i−1 (t) (10)

− (i + 2β− 2)C(β)∗i−2 (t)
]
, i ≥ 2

where C(β)∗0 (t) = 1, C
(β)∗
1 (t) = 2β

(
2t−(a+b)

b−a

)
.

The analytic form of the shifted Gegenbauer polynomial
C(β)∗i (t) is given as:

C(β)∗i (t) =
i∑

j=0

(−1) jΓ(2β + 2i − j)Γ(β + 12 )

(i − j)! Γ( j + 1)Γ(i − j + β + 12 )Γ(2β)
ti− j. (11)

The orthogonality condition is

〈C(β)∗m (t), C
(β)∗
n (t)〉 =

∫ 1
0

(
t − t2

)(α− 12 ) C(β)∗m (t)C(β)∗n (t)dt (12)
=


0, for m 6= n

π21−4αΓ(n+2α)
n![Γ(α)]2 (n+α)

, for m = n

Let f (x) be a square integrable function in [0, 1], expressing it
in terms of the shifted Gegenbauer polynomials as follows:

f (x) =
∞∑

i=0

αiC
(β)∗
i (x), (13)

where αi is defined as:

αi =
i! [Γ(β)]2(i + β)
π21−2βΓ(i + 2β)

∫ 1
−1

(1 − x2)β−
1
2 f

(
x + 1

2

)
Cβi (x)d x, (14)

or,

αi =
n! [Γ(β)]2(i + β)
π21−4βΓ(i + 2β)

∫ 1
0

(x − x2)β−
1
2 f (x)C(β)∗i (x)d x, (15)

(N + 1) terms of C(β)∗n (x) are considered in the approximation.
Then Eq. (13) becomes

f (x) =
N∑

i=0

αiC
(β)∗
i (x). (16)

2.3. Statement of the problem
The problem under consideration and its initial and bound-

ary conditions are:

∂u(x, t)
∂t

= a(x)
∂ωu(x, t)
∂xω

+b(x, t) 0 < x < ς, 0 ≤ t ≤ υ, 1 < ω < 2
(17)

u(x, 0) = g(x), 0 < x < ς (18)

u(0, t) = µ0(t), 0 ≤ t ≤ υ
u(ς, t) = µ1(t), 0 ≤ t ≤ υ

(19)

Eqs. (18) and (19) are initial and boundary conditions respec-
tively.

When ω = 2, Eq. (17) becomes a classical diffusion equation,
that is:

∂u(x, t)
∂t

= a(x)
∂2u(x, t)
∂x2

+ b(x, t). (20)

Here, we try to reduce the fractional derivative to a system of
ordinary differential equations by utilizing the shifted Gegen-
bauer polynomials and then solve the resulting system of equa-
tion using compact finite difference technique.

We recall some lemmas to justify our finding
Let η denote an open bounded domain in R2 space and L2(η)
represents a Hilbert space with the inner product

〈p(x), q(x)〉 =
∫
η

p(x)q(x)d x, (21)

together with the Euclidean norm ‖u(x)‖= 〈u(x),υ(x)〉2, and the
Sobolev space as

H s(η) = u�L2(η),
dsu
d xs

�L2(η)

Suppose that b(x) > 0 for 0 ≤ x ≤ 1

Lemma 2.1. For any p, q ∈ H
β
2 (η), 〈a D

β
x p, q〉 = 〈a D

β
2
x p,x D

β
2
b q〉,

〈x D
β
b p, q〉 = 〈x D

β
2
b p,a D

β
2
x q〉, for 1 < β < 2,

〈a D
β
x p,x D

β
b p〉 = cos(βπ)‖a D

β
x p‖

2= cos(βπ)‖x D
β
b p‖

2, f orβ > 0
3



Issa et al. / J. Nig. Soc. Phys. Sci. 5 (2023) 1368 4

Lemma 2.2. For the functions f (x) and 〈a D
β
x f (x)〉�H

β(η),
∃ δτ > 0 3

‖ f (x) +
δ

2
b(xn)a D

β
x f (x)‖≤ ‖ f (x)‖, f or 1 < β < 2

see [22] for details.

3. Shifted Gegenbauer CFDM

Here, we employed the formula for fractional derivative
f (x) derived in [7]:

Assume ω > 0, now using Caputo linearity property, we have:

Dω fN (x) =
N∑

i=0

αi D
ω
(
C(β)∗i (x)

)
. (22)

By linearity property and Eq. (3) we obtain:

Dω
(
C(β)∗i (x)

)
= 0, i = 0, 1, . . . ,dωe− 1,ω > 0 (23)

Dω
(
C(β)∗i (x)

)
=

i∑
j=0

(−1) jΓ(2β + 2i − j)Γ(β + 12 )

(i − j)! Γ( j + 1)Γ(i − j + β + 12 )Γ(2β)

(24)

× Dωxi− j, i ≥ dωe

Applying Eq. (3) to Eq. (24), we obtain:

Dω
(
C(β)∗i (x)

)
=

i−dωe∑
j=0

(−1) jΓ(2β + 2i − j)Γ(β + 12 )

Γ( j + 1)Γ(i − k + β + 12 )Γ(2β)Γ(i + 1 − j −σ)

(25)

× xi− j−ω

Substituting (25) into (22) to obtain:

Dω(gN (x)) = (26)
N∑

i=dωe

i−dωe∑
j=0

αi
(−1) jΓ(2β + 2i − j)Γ(β + 12 )

Γ( j + 1)Γ(i − j + β + 12 )Γ(2β)Γ(i + 1 − j −σ)
xi− j−ω

Let

Ni, j =
(−1) jΓ(2β + 2i − j)Γ(β + 12 )

Γ( j + 1)Γ(i − j + β + 12 )Γ(2β)Γ(i + 1 − j −σ)
(27)

Simplifying Eq. (26) gives

Dw( fN (x)) =
N∑

i=dωe

i−dωe∑
j=0

αi Ni, j x
i− j−ω. (28)

See [7] for for more details

3.1. Formulation of the scheme

Solving Eq.(17) based on CFDM, using Gegenbauer col-
location, given two integers M, N > 0 and two mesh points
τm−1 = (m − 1)δτ f or m = 1, 2, ..., M + 1, δτ =

T
M .

Introducing Taylor’s expansion for the time-discretizing
that is:

∂u(xn, tm)
∂t

= δτu(xn, tm) +
δτ

2
∂2u(xn, tm)

∂t2
+ 0(δτ2) (29)

where δτu(xn, tm) =
umn −u

m−1
n

δτ
, substituting Eq.(29) into Eq.(3)

δτu(xn, tm)+
δτ

2
∂2u(xn, tm)

∂t2
+0(δτ2) = a(xn)

∂wu(xn, tm)
∂xw

+b(xn, tm)
(30)

Differentiating Eq.(17)

∂2u(x, t)
∂t2

= a(xn)δτ
∂wu(xn, tm)

∂xw
+ δτb(xn, tm) (31)

substituting Eq.(31) in Eq.(30) gives:

δτu(xn, tm) = a(xn)
∂wu(xn, tm)

∂xw
+ b(xn, tm) (32)

−
δτ

2

[
a(xn)δτ

∂wu(xn, tm)
∂xw

+ δτb(xn, tm)
]

+ · · ·

Simplifying Eq.(32) gives

umn − u
m−1
n =

δτ

2
a(xn)

∂wumn
∂xw

+
δτ

2
a(xn)

∂wum−1n
∂xw

+
δτ

2
(bmn + b

m−1
n )

(33)

+ Rn(x)(δτ)3

where,Rn(x) is the truncation term, let u(xn, tm) = U mn

U mn −
δτ

2
a(xn)

∂wU mn
∂xw

= U m−1n +
δτ

2
a(xn)

∂wU m−1n
∂xw

(34)

+
δτ

2
(bmn + b

m−1
n ) + R

n(x)(δτ)3,

to obtain full discrete form, it is required to approximate the Ca-
puto derivative in ∂

w U mn
∂xw . The approximant uN (x, t) is constructed

using the Gegenbauer collocation technique.

uN (x, t) =
N∑

k=0

uk(t)C
(β)∗
k (x) (35)

Using Eqs. (34),(35) and the roots of shifted Gegenbauer poly-
nomial of degree N − 1, we obtain the simplified form as

N∑
k=0

uk(t)C
(β)∗
k (xn) −

δτ

2
a(xn)

N∑
k=dwe

k−dwe∑
i=0

uk(t)Nk,i x
k−i−w

=

N∑
k=0

um−1k C
(β)∗
k (xn) +

δτ

2
a(xn)

N∑
k=dwe

k−dwe∑
i=0

uk(t)Nk,i x
k−i−w

+
δτ

2
(b(xn, tm) + b(xn, tm−1))

(36)

4



Issa et al. / J. Nig. Soc. Phys. Sci. 5 (2023) 1368 5

with the boundary conditions:

U(0, t) =
N∑

n=0

(−1)nΓ(n + 2β)
n! Γ(2β)

un(t) = µ0,

U(1, t) =
N∑

n=0

Γ(n + 2β)
n! Γ(2β)

un(t) = µ1

(37)

Eqs. (36) and (37) gives (N +1) system of linear equations, then
obtain uN, n = 0, 1, · · · , N. To get the initial condition U 0n of
Eq.(36), the initial condition of the problem U(x, 0) combining
Eq. (15) will be used to obtain U 0n .

3.2. Convergence and Stability Analysis

In this section, we consider convergence and stability of the
proposed method. Let η denote an open bounded domain in
R2 space and L2(η) represents a Hilbert space with the inner
product

〈p(x), q(x)〉 =
∫
η

p(x)q(x)d x, (38)

together with the Euclidean norm ‖u(x)‖= 〈u(x),υ(x)〉2, and the
Sobolev space as

H s(η) = u ∈ L2(η),
dsu
d xs
∈ L2(η)

Suppose that b(x) > 0 for 0 ≤ x ≤ 1 the following lemmas
are required to investigate the stability and convergence of the
method.

Lemma 3.1. Let U j ∈ H1(η), j = 1, 2, ..., J be the solution of
equation (34)and U 0 be the initial condition, then the following
inequality holds

‖U j‖≤ ‖U 0‖+ max
0≤τ≤N

δτ

2
(‖b ji‖+‖b

j−1
i ), (39)

where, U j = u(x, t j)

Proof: By Mathematical induction on j.
For j = 1 in Eq. (34), we obtain

U 1−
δτ

2
a(xi)a D

ω
x U

1
= U 0 +

δτ

2
a(xi) a D

ω
x U

0
+
δτ

2
(b1 +b0), (40)

multiplying Eq. (40) by U 1 and integrating on η, we get:

‖U 1‖2−
δτ

2
a(xi)〈a D

ω
x U

1, U 1〉 = 〈U 0, U 1〉 +
δτ

2
a(xi)〈a D

ω
x U

0, U 1〉

+
δτ

2
(〈b1, U 1〉 + 〈b0, U 1〉).

(41)

Since cos
(
β
2 π

)
< 0, using Lemma 2.1, we have

〈a D
β
xU

1, U 1〉 = 〈a D
β
2
x U

1,x D
β
2
b U

1
〉 = cos

(
β

2
π
)
‖a D

β
2
x U

1
‖

2< 0
(42)

since cos( β2 π) < 0, for 1 < β < 2.
From the left side of Eq. (41), we obtain

‖U 1‖2≤ ‖U 1‖2−
δτ

2
a(xi)〈a D

β
xU

1, U 1〉 (43)

Using Cauchy-Schwartz inequality, the right side of Eq. (41)
becomes

‖〈U 0, U 1〉 +
δτ

2
a(xi)〈a D

ω
x U

0, U 1〉‖

≤ ‖U 0 +
δτ

2
a(xi) a D

ω
x U

0
‖‖U 1|≤ ‖U 0‖‖U 1‖

(44)

Using Eqs. (42),(43), and (44), to obtain

‖U 1‖≤ ‖U 0‖+ max
0≤n≤N

δτ

2
(‖b ji‖+‖b

j−1
i ‖) (45)

Assume Eq. (39) is true for k = 1, 2, ..., j − 1, then

‖U k‖≤ ‖U 0‖+ max
0≤n≤N

δτ

2
(‖bki ‖+‖b

k−1
i ‖) (46)

Now, multiplying Eq. (34) by U j and integrating on η

‖U j‖2−
δτ

2
a(xi)〈a D

ω
x U

j, U j〉 = 〈U j−1, U j〉 +
δτ

2
a(xi)〈a D

ω
x U

j−1, U j〉

+
δτ

2
(〈b j, U j〉 + 〈b j−1, U j〉)

(47)

From the previous procedure, Eq. (47) becomes

‖U j‖≤ ‖U j−1‖+ max
0≤i≤N

δτ

2
(‖b ji‖+‖b

j−1
i ‖)

Theorem 3.2. The approximation method introduced in equa-
tion (34) is unconditionally stable.

Proof Let U ji , j = 1, 2, ..., J be the approximate solution
obtained in Eq.(34) with the initial condition U 0i = u(xi, 0), then
the error ε j = u(xi, t j) − U

j
i satisfies

ε j −
δτ

2
p(xi)

∂βε j

∂xβ
−
δτ

2
p(xi)

∂βε j−1

∂xβ
= ε j−1. (48)

According to lemma 3.1,

‖ε j‖≤ ‖ε0‖ for, j = 1, 2, . . . , M,

hence the proof.

4. Numerical Experiments

In this section, we validate the efficiency and the accuracy
of the proposed method (PM) on some selected examples from
the literature by comparing the maximum error �N with the ex-
isting results in the latest literature. Where

�N = max
0≤i≤N

|U(xi, T ) − uN (xi, T )| (49)

5



Issa et al. / J. Nig. Soc. Phys. Sci. 5 (2023) 1368 6

Example 4.1
Considering SFODE [7, 9, 22]

∂u(x, t)
∂t

= Γ

(
6
5

)
x

9
5
∂

9
5 u(x, t)

∂x
9
5

+ 3x2(2x − 1) exp(−t)

u(x, 0) = x3(x−1 − 1),
u(0, t) = u(1, t) = 0, t > 0.

(50)

with the closed form solution u(x, t) = x3(x−1 − 1) exp(−t).

Comparing Eq. (50) with (17), we have

a(x) = Γ
(

6
5

)
x

9
5 , b(x, t) = 3x2(2x − 1) exp(−t), ω =

9
5

For approximation of degree 3, that is N = 3, Eqs. (35) and
(36) gives

u3(x, t) =
3∑

k=0

umk (t)C
(β)∗
k (x) (51)

um0 C
(β)∗
0 (xn) + u

m
1 C

(β)∗
1 (xn) + u

m
3 C

(β)∗
3 (xn)

+
δτ

2
a(xn)

[
um2 N2,0 x

0.2
+ um3 (N3,0 x

1.2
+ N3,1 x

0.2)
]

= um−10 C
(β)∗
0 (xn) + u

m−1
1 C

(β)∗
1 (xn) + u

m−1
3 C

(β)∗
3 (xn)

+
δτ

2
a(xn)

[
um−12 N2,0 x

0.2
+ um−13 (N3,0 x

1.2
+ N3,1 x

0.2)
]

+
δτ

2
(b(xn, tm) + b(xn, tm−1)),

(52)

respectively.
At x = x0, Collocating at zeros of C

(β)∗
2 (x), we obtain

um0 C
(β)∗
0 (x0) + u

m
1 C

(β)∗
1 (x0) + u

m
3 C

(β)∗
3 (x0)

+
δτ

2
a(x0)

[
um2 N2,0 x

0.2
+ um3 (N3,0 x

1.2
+ N3,1 x

0.2)
]

= um−10 C
(β)∗
0 (x0) + u

m−1
1 C

(β)∗
1 (x0) + u

m−1
3 C

(β)∗
3 (x0)

+
δτ

2
a(x0)

[
um−12 N2,0 x

0.2
+ um−13 (N3,0 x

1.2
+ N3,1 x

0.2)
]

+
δτ

2
(b(x0, tm) + b(x0, tm−1))

(53)

we have C(β)∗0 (x0) = 1. Let G1 = C
(β)∗
1 (x0), G2 = C

(β)∗
3 (x0),

H1 =
δτ
2 a(x0)N2,0 x

0.2, H2 =
δτ
2 a(x0)(N3,0 x

1.2 + N3,1 x0.2)

the equation becomes

um0 + u
m
1 G1 + u

m
3 G2 + u

m
2 H1 + u

m
3 H2

= um−10 + u
m−1
1 G1 + u

m−1
3 G2 + u

m−1
2 H1 + u

m−1
3 H2

(54)

at x = x1, we get

um0 C
(β)∗
0 (x1) + u

m
1 C

(β)∗
1 (x1) + u

m
3 C

(β)∗
3 (x1)

+
δτ

2
a(x1)

[
um2 N2,0 x

0.2
+ um3 (N3,0 x

1.2
+ N3,1 x

0.2)
]

= um−10 C
(β)∗
0 (x1) + u

m−1
1 C

(β)∗
1 (x1) + u

m−1
3 C

(β)∗
3 (x1)

+
δτ

2
a(x1)

[
um−12 N2,0 x

0.2
+ um−13 (N3,0 x

1.2
+ N3,1 x

0.2)
]

+
δτ

2
(b(x1, tm) − b(x1, tm−1)).

(55)

Let G11 = C
(β)∗
1 (x0), G22 = C

(β)∗
3 (x0), H11 =

δτ
2 a(x0)N2,0 x

0.2,
H22 =

δτ
2 a(x0)(N3,0 x

1.2 + N3,1 x0.2)
Eq. (55) becomes

um0 + u
m
1 G11 + u

m
3 G22 + u

m
2 H11 + u

m
3 H22

= um−10 + u
m−1
1 G11 + u

m−1
3 G22 + u

m−1
2 H11 + u

m−1
3 H22.

(56)

From Eq. (37) we obtain

um0 (t)−
Γ(1 + 2β)

Γ(2β)
um1 (t) +

Γ(2 + 2β)
2Γ(2β)

um2 (t)−
Γ(3 + 2β)

6Γ(2β)
um3 (t) = µ0,

(57)
and

um0 (t) +
Γ(1 + 2β)

Γ(2β)
um1 (t) +

Γ(2 + 2β)
2Γ(2β)

um2 (t) +
Γ(3 + 2β)

6Γ(2β)
um3 (t) = µ1.

(58)
Eqs. (54),(56)-(58) can be written as

Aum = Bum−1 +
δτ

2

(
bmn − b

m−1
n

)
, (59)

or
um = A−1 Bum−1 +

δτ

2

(
bmn − b

m−1
n

)
, (60)

where

A =


1 G1 H1 G2 + H2
1 G11 H11 G22 + H22
1 −Γ(1+2β)

Γ(2β)
Γ(2+2β)
2Γ(2β) −

Γ(3+2β)
6Γ(2β)

1 Γ(1+2β)
Γ(2β)

Γ(2+2β)
2Γ(2β)

Γ(3+2β)
6Γ(2β)

 ,

B =


1 G1 H1 G2 + H2
1 G11 H22 G33 + H22
0 0 0 0
0 0 0 0

 , um =

um0
um1
um2
um3


Table 1 shows the computational errors at different values

of β while Figure 1a is the corresponding figure. Table 2 dis-
play the comparison of the absolute errors at β = 0.5 with the
existing results in the literature at β = 0.5.

Example 4.2

Consider the space fractional order diffusion problem [7, 9]

∂u(x, t)
∂t

=
Γ(2.2)

6
x2.8

∂
3
2 u(x, t)

∂x
3
2

− x3(1 + x) exp(−t)

u(x, 0) = x3

u(0, t) = 0, u(1, t) = exp(−t), t > 0.

The exact solution is u(x, t) = x3e−t.

Table 3 display the comparison of the maximum errors at dif-
ferent values of β and at β = 1 in [7] is equivalent to the results
obtained in [9]. Figure 2 is the graphical representation of the
exact solution and its corresponding approximate solutions at
different values of β while Figure 3 display the absolute errors
at various values of β.

6



Issa et al. / J. Nig. Soc. Phys. Sci. 5 (2023) 1368 7

Figure 1: Example 4.1, relationship between exact and approximate solution at
N = 3, & T = 1

Example 4.3
Consider the space fractional diffusion problem [7]

∂u(x, t)
∂t

= Γ

(
3
2

)
x

1
2
∂

3
2 u(x, t)

∂x
3
2

− 2x sin(t + 1) + (x2 + 1) cos(t + 1)

u(x, 0) = (x2 + 1) sin(1)
u(0, t) = sin(t + 1), u(1, t) = 2 sin(t + 1), t > 0.

The exact solution is u(x, t) = (x2 + 1) sin(t + 1)

Comparison of the maximum errors are display in Table 4. Fig-
ure 4 is the graphical representation of the absolute errors at
various values of β.

5. Discussion of results and conclusion

5.1. Discussion of Results
Table 1 shows the absolute errors of the proposed method

for example 4.1 with respect to the values of β at N = 3 and
T = 1. The absolute errors of the proposed method is compared
with the results in [22],[7] and [9] as shown in Table 2. The ap-
proximate solutions at β = 1 and β = 1.5 correspond to the

Table 1: Example 4.1, absolute errors for various values of β

x β = 0.5 β = 1 β = 1.5

0 0 2.2959 × 10−41 4.4409 × 10−16

0.1 1.0738 × 10−9 1.2324 × 10−9 1.1618 × 10−9

0.2 1.7375 × 10−9 1.9252 × 10−9 1.8596 × 10−9

0.3 2.0553 × 10−9 2.1777 × 10−9 2.1705 × 10−9

0.4 2.0916 × 10−9 2.0890 × 10−9 2.1718 × 10−9

0.5 1.91075 × 10−9 1.7583 × 10−9 1.9406 × 10−9

0.6 1.5770 × 10−9 1.2849 × 10−9 1.5542 × 10−9

0.7 1.1547 × 10−9 7.6791 × 10−10 1.0897 × 10−9

0.8 7.0821 × 10−10 3.0655 × 10−10 6.2437 × 10−10

0.9 3.0187 × 10−10 2.1667 × 10−17 2.3540 × 10−10

1 0 2.1667 × 10−17 8.8818 × 10−16

Table 2: Exapmle 4.1, comparison of absolute errors

x [9] [22] [7] PM

0 0 4.77 × 10−6 0 0

0.1 5.33 × 10−6 3.17 × 10−9 5.46 × 10−6 1.07 × 10−9

0.2 8.06 × 10−6 5.85 × 10−9 8.51 × 10−6 1.74 × 10−9

0.3 8.72 × 10−6 7.97 × 10−9 9.60 × 10−6 2.06 × 10−9

0.4 7.84 × 10−6 9.44 × 10−9 9.18 × 10−6 2.09 × 10−9

0.5 5.96 × 10−6 1.02 × 10−8 7.69 × 10−6 1.91 × 10−9

0.6 3.59 × 10−6 1.01 × 10−8 5.60 × 10−6 1.58 × 10−9

0.7 1.29 × 10−6 9.12 × 10−9 3.33 × 10−6 1.15 × 10−9

0.8 4.32 × 10−7 7.17 × 10−9 1.34 × 10−6 7.08 × 10−10

0.9 4.04 × 10−6 4.16 × 10−9 8.39 × 10−8 3.02 × 10−10

1 0 7.55 × 10−17 0 0

Table 3: Example 4.2, maximum error for various values of β

β [7] PM

β = 0.5 2.25 × 10−06 6.30 × 10−09

β = 1 8.83 × 10−06 1.32 × 10−08

β = 1.5 8.35 × 10−06 3.69 × 10−08

7



Issa et al. / J. Nig. Soc. Phys. Sci. 5 (2023) 1368 8

Table 4: Maximum error for Example 4.3 relative to values of β

β [7] PM

β = 0.5 1.56 × 10−05 8.23 × 10−09

β = 1 1.36 × 10−05 2.01 × 10−08

β = 1.5 11.05 × 10−05 1.06 × 10−09

Figure 2: Example 4.2, relationship between exact and approximate solutions
at N = 3, & T = 1

Figure 3: Absolute errors for Example 4.2 at N = 3, & T = 1

Figure 4: Absolute errors for Example 4.3 at N = 3, & T = 1

approximate solutions when using second kind shifted Cheby-
shev polynomial U∗j (x) and Jacobi ( j, 1, 1, x) (that’s P

1,1
j (x) as

an approximation polynomials. Tables 2-4 show the compari-
son of the errors at various values of β relative to the existing
results in the literature. Figures 1-2 are the comparison of the
solution for various values of β, while figures 3 and 4 represent
the absolute errors of the proposed method at various values of
β.

5.2. Conclusion
This paper studied space fractional order diffusion equation

by proposing a new technique for finding the approximate so-
lution using compact finite difference method to discretize the
time derivative, then use shifted Gegenbauer polynomials as ap-
proximating polynomial. The computational results are plotted
and detailed in the tables to justify the accuracy of the propose
method. It can be observe in the tables of the results and figures
that the proposed method performs better than the existing one
in the literature.

References

[1] J. T. Machado, V. Kiryakova & F. Mainardi, “Recent history of fractional
calculus”, Commun. Nonlinear Sci. Numer. Simul., 16 (2011) 1140.

[2] A. E. Abouelrega, S. Yao & H. Ahmad, “Analysis of a functionally graded
thermopiezoelectric finite rod excited by a moving heat source”, Results
in Physics 19 (2020) 103389.

[3] A. Atangana & J. Gómez-Aguilar. “Fractional derivatives with no-index
law property: application to chaos and statistics”, Chaos Solitons Fractals
114 (2018) 516.

[4] A. A. Hamoud & K. P. Ghadle, “Some new existence, uniqueness and
convergence results for fractional volterra-fredholm integro-differential
equations”, J. Appl. Comput. Mech. 5 (2019) 58.

[5] M. Dehghan & A. Saadatmandi, “Chebyshev finite difference method for
Fredholm integro-differential equation”, International Journal of Com-
puter Mathematics 85 (2008) 123.

[6] V. O. Bohaienko, “A fast finite-difference algorithm for solving space-
fractional filtration equation with a generalised Caputo derivative”, Com-
putational and Applied Mathematics 38 (2019) 1.

[7] K. Issa, B. M. Yisa & J. Biazar, “Numerical solution of space frac-
tional diffusion equation using shifted Gegenbauer polynomials”, Com-
putational Methods for Differential Equations 10 (2022) 431.

[8] A. Saadatmandi & M. Dehghan, “A tau approach for solution of the space
fractional diffusion equation”, Computers & Mathematics with Applica-
tions 62 (2011) 1135.

8



Issa et al. / J. Nig. Soc. Phys. Sci. 5 (2023) 1368 9

[9] N. H. Sweilam, A. M. Nagy & A. A. El-Sayed, “Second kind shifted
Chebyshev polynomials for solving space fractional order diffusion equa-
tion”, Chaos, Solitons & Fractals, 73 (2015) 141.

[10] H. Hassani, J. T. Machado & E. Naraghirad, “Generalized shifted Cheby-
shev polynomials for fractional optimal control problems”, Communica-
tions in Nonlinear Science and Numerical Simulation 75 (2019) 50.

[11] J. Singh, D. Kumar & D. Baleanu, “New aspects of fractional Biswas-
Milovic model with Mittag-Leffler law”, Mathematical Modelling of Nat-
ural Phenomena 14 (2019) 303.

[12] K. M. Owolabi & A. Atangana, “High-order solvers for space-fractional
differential equations with Riesz derivative”, Discrete & Continuous Dy-
namical Systems-S 12 (2019) 567.

[13] O. Nikan, A. Golbabai, J. A. Machado, & T. Nikazad, “Numerical solu-
tion of the fractional Rayleigh-Stokes model arising in a heated general-
ized second-grade fluid”, Engineering with Computers, 37 (2021) 1751.

[14] P. Pandey & J. F. Gömez-Aguilar, “On solution of a class of nonlinear
variable order fractional reaction-diffusion equation with Mittag-Leffler
kernel”, Numer. Methods Partial Differential Eq. 37 (2021) 998.

[15] Y. Salehi, M. T. Darvishi & W. E. Schiesser, “Numerical solution of space
fractional diffusion equation by the method of lines and splines”, Applied
Mathematics and Computation 336 (2018) 465.

[16] Z. Soori, & A. Aminataei, “Numerical solution of space fractional dif-
fusion equation by spline method combined with Richardson extrapola-
tion”, Computational and Applied Mathematics 39 (2020) 1.

[17] A. H. Bhrawy, & M. M. Al-Shomrani, “A shifted Legendre spectral
method for fractional-order multi-point boundary value problems”, Ad-

vances in Difference Equations 2012 (2012) 8.
[18] Y. Luchko1 & M. Yamamoto, “Comparison principles for the time-

fractional diffusion equations with the Robin boundary conditions. Part
I: Linear equations”, Fract. Calc. Appl. Anal., (2023), http://arxiv.
org/abs/2304.07812v1.

[19] L. Qiu, X. Ma & Q. Qin, “A novel meshfree method based on spatio-
temporal homogenization functions for one-dimensional fourth-order
fractional diffusion wave equations”, Applied Math. Letters 142 (2023)
108657.

[20] F. Muhammad, A. Ali, T. T. Merve, M. A. Muhammad, A. Aqeel, E. M.
Emad & S. Y. Ibrahim, “Fractal fractional-order derivative for HIV/AIDS
model with Mittag-Leffler kernel”, Alexandria Engineering Journal 61
(2022) 10965.

[21] F. Muhammad, A. Muhammad, A. Ali & A. Aqeel, “Modeling of
fractional-order COVID-19 epidemic model with quarantine and social
distancing”, Mathematical Method in Applied Sciences 44 (2021) 9334.

[22] Y. E. Aghdam, H. Safdari, Y. Azari, H. Jafari & D. Baleanu, “Numerical
investigation of space fractional order diffusion equation by the Cheby-
shev collocation method of the fourth kind and compact finite difference
scheme”, Discrete & Continuous Dynamical Systems-S 14 (2021) 2025.

[23] J. C. Mason & D. C. Handscomb, Chebyshev polynomials, Chapman and
Hall, CRC Press, 2003.

[24] W. M. Abd-Elhameed1, & Y. H. Youssri, “Fifth-kind orthonormal Cheby-
shev polynomial solutions for fractional differential equations”, Comp.
Appl. Math. 37 (2018) 2897.

9