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

Journal of the
Nigerian Society

of Physical
Sciences

Collocation Approach for the Computational Solution Of
Fredholm-Volterra Fractional Order of Integro-Differential

Equations

G. Ajileyea,∗, A.A. Jamesb, A. M. Ayindec, T. Oyedepod

aDepartment of Mathematics and Statistics, Federal University Wukari, Taraba State, Nigeria.
bDepartment of Mathematics, American University of Nigeria, Yola, Adamawa State, Nigeria.

cDepartment of Mathematics, Moddibo-Adama University, Yola Nigeria.
dFederal College of Dental Technology and Therapy, Enugu, Nigeria.

Abstract

In this work, a collocation technique is used to determine the computational solution to fractional order Fredholm-Volterra integro-differential
equations with boundary conditions using Caputo sense. We obtained the linear integral form of the problem and transformed it into a system
of linear algebraic equations using standard collocation points. The matrix inversion approach is adopted to solve the algebraic equation and
substituted it into the approximate solution. We established the uniqueness and convergence of the method and some modelled numerical examples
are provided to demonstrate the method’s correctness and efficiency. It is observed that the results obtained by our new method are accurate and
performed better than the results obtained in the literature. The study will be useful to engineers and scientists. It is advantageous because
it addresses the difficulty in tackling fractional order Fredholm-Volterra integro-differential problems using a simple collocation strategy. The
approach has the advantage of being more accurate and reducing computer running time.

DOI:10.46481/jnsps.2022.834

Keywords: Collocation, Fredholm-Volterra, Integro-differential, Fractional derivatives, Approximate Solution.

Article History :
Received: 27 May 2022
Received in revised form: 13 August 2022
Accepted for publication: 18 August 2022
Published: 03 October 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: J. Ndam

1. Introduction

Fractional calculus is an aspect in mathematics that studies
the properties of integrals combined with noninteger order
derivatives.The concept and method of solving differential
equations containing fractional derivatives of unknown
functions are covered in this field (fractional differential

∗Corresponding author tel. no: +2348077092831
Email address: adewale.james@aun.edu.ng (G. Ajileye )

equations). Many prominent mathematicians, such as
Liouville, Grunwald, Riemann, Euler, Langrange, Heaviside,
Fourier, Abel, and others, established fractional calculus on
a formal foundation[1]. Very recently, scholars developed
huge interest in fractional calculus due to its relevance and
application to many areas of scientific endeavors. Fractional
differential equations, or those containing real and complex
order derivatives have been more significant in describing the
most broad fields of science and technology with peculiar
dynamics of various processes involving complex systems [2].

1



Ajileye et al. / J. Nig. Soc. Phys. Sci. 4 (2022) 834 2

Many different approaches have been adopted to investigate
the solution of fractional integrodifferential equations, such as
Adomian decompositions method [3-5], collocation method
[6, 7], Laplace decomposition method [8, 9], Taylor
expansion method [10], Least square method [9], differential
transform method [11], homotopy perturbation method [12-18],
sinc-collocation method [15-17] and variational iteration
method [12, 13]. Linna et al. [20] considered a
numerical method to solve fractional variational problem.
They simplified the fractional variational problems by the
operational matrices. The operational matrices are based
on the Chelyshkov polynomials. The fractional variational
problem was transformed into a set of algebraic equations.
The unknown coefficients were solved using the Lagrange
multiplier techniques. Avcı & Mahmudov [19] proposed
a numerical method based on the fractional Taylor vector
for solving multi-term fractional differential equations. The
main idea of this method was to reduce the given problems
to a set of algebraic equations by utilizing the fractional
Taylor operational matrix of fractional integration. Some
numerical examples were given to demonstrate the accuracy
and applicability. The results show that the presented method
was efficient and applicable. Fadugba [21] presented the
Mellin transform for the solution of the fractional order
equations. The Mellin transform approach occurs in many
areas of applied mathematics and technology. The Mellin
transform of fractional calculus of different flavours; namely
the Riemann-Liouville fractionalderivative, Riemann-Liouville
fractional integral, Caputo fractional derivative and the
Miller-Ross sequential fractional derivative were obtained.

In this study, we consider Fredholm-Volterra
Integrodifferential equation of fractional order of the the
form:

D
α

xy(x) − P1y
′

(x) − P2y
′′

(x) − P0y(x)

= g (x) + λ1

∫ x
0

k1 (x, t) y (t) dt

+λ2

∫ 1
0

k2 (x, t) y (t) dt, (1)

with the given boundary conditions

y (a) = 0, y(b) = 0, a < x < b, (2)

where y(x) is to be determined, Dαx is the Caputo’s derivative,
k1 (x, t) and
k2 (x, t) are the Fredholm and Volterra integral kanel function

respectively. P j, Pα,λ j are known constants. g(x) is the known
function

2. Basic Definitions

Under this section, we present some definitions and basic
concepts of fractional calculus for the formulation of the given
problem

Definition 2.1: The Caputo derivative with order α > 0 of the
given function f (x), x ∈ (a, b) is defined as [ Litfi, Dehghan
and Yousefi, 2011]

C
x D

α
a y(x) =

1
Γ(m −α)

∫ x
a

(x − s)m−α−1y(m)(s)d s, (3)

where m − 1 ≤ α ≤ m, m ∈ N, x > 0
Definition 2.2: Let (an) , n ≥ 0 be a sequence of real numbers.
The power series in x with coefficients an is an expression [
Edward, Ford and Simpson, 2002]

y(x) = a0 + a1 x + a2 x
2

+ a3 x
3

+ · · ·aN x
N

=

N∑
n=0

an x
n

= φ(x) A, (4)

where φ(x) = [1 x x2 · · · xN ], A = [a0 a1 · · · aN ]T ,
then y(x, n) = xnA, n = 0(1)N, n ∈ Z+.
Definition 2.3: Standard Collocation Method (SCM). This
approach is used to find the collocation points that are desired
within a certain interval. i.e [a,b] and is given by

xi = a +
(b − a)i

N
, i = 1, 2, 3, ...N. (5)

Definition 2.4: A metric on a set M is a function d : M×M −→
R with the following properties. For all x, y ∈ M
(a) d(x, y) ≥ 0;
(b) d(x, y) = 0 ⇐⇒ x = y
(c) d(x, y) = d(y, x)
(d) d(x, y) ≤ d(x, z) + d(x, y)
If d is a metric on M, then the pair (M, d) is called a metric
space.
Definition 2.5: Let (X, d) be a metric space, A mapping T :
X −→ X is Lipschitzian if ∃ a constant L > 0 such that
d(T x, T y) ≤ Ld(x, y) ∀ x, y ∈ X.
Definition 2.6: Let y(x) be a continuous function, then

0 I
β
x

(
C
0 D

β
xy(x)

)
= y(x) −

N∑
k=0

y(k)(0)
k!

xk, (6)

where m − 1 < β < 1.
Let p(s) be an integrable function, then

0 I
β
x ( p(s)) =

1
Γ(β)

∫ x
0

(x − s)β−1 p(s)d s. (7)

3. Mathematical Background

In this section, combination of collocation method and
power series approximation is employed for the computational
solution of (FVID) of fractional order.
Theorem 3.0: Banach’s fixed point theorem Let (X, d) be
a complete metric space, then each contraction mapping T :
X −→ X has a unique fixed point x of T in X, such that T (x) = x
Lemma (3.1). Let y(x) be the solution to (1) subject to (2) , the
integral form

y (x) = W(x) +
1

Γ(α)

∫ x
0

(x − s)α−1
(
P2y

′′

(s)
)

d s (8)

2



Ajileye et al. / J. Nig. Soc. Phys. Sci. 4 (2022) 834 3

+
1

Γ(α)

∫ x
0

(x − s)α−1
(
P1y

′

(s)
)

d s

+
1

Γ(α)

∫ x
0

(x − s)α−1 (P0y(s)) d s

−
1

Γ(α)

∫ x
0

(x − s)α−1
[
λ1

∫ s
0

k1 (s, t) y (t) dt

+ λ2

∫ 1
0

k2 (s, t) y (t) dt
]

where

W(x) =
N∑

k=0

y(k)(0)
k!

xk −
1

Γ(α)

∫ x
0

(x − s)α−1 g(s)d s.

Multiply equation (1) by 0 I
β
x (.) gives

0 I
α
x (D

αy(x)) = 0 I
β
x

(
P2y

′′

(x)
)

+ 0 I
β
x

(
P1y

′

(s)
)

+0 I
β
x (P0y(s)) −0 I

β
x (g(x))

− 0 I
β
x

[
λ1

∫ s
0

k1 (s, t) y (t) dt

+ λ2

∫ 1
0

k2 (s, t) y (t) dt
]
. (9)

Using(6) on equation (9)gives

y(x) =
N∑

k=0

y(k)(0)
k!

xk + 0 I
β
x

(
P2y

′′

(x)
)

+ 0 I
β
x

(
P1y

′

(s)
)

+0 I
β
x (P0y(s)) (10)

0 I
β
x (g(x)) + 0 I

β
x

(∫ b
0

k1(x, t)y(t)dt
)

+ 0 I
β
x

(∫ s
0

k2(s, t)y(t)dt
)

Applying equation (4) and (7) to (10) gives

y(x) =
N∑

k=0

y(k)(0)
k!

xk

+
1

Γ(α)

∫ x
0

(x − s)α−1 P2
(
φ(s)

′′
)

d sA

+
1

Γ(α)

∫ x
0

(x − s)α−1 P1
(
φ(s)

′
)

d sA

+
1

Γ(α)

∫ x
0

(x − s)α−1 P0 (φ(s)) d sA

−
1

Γ(α)

∫ x
0

(x − s)α−1 (g(s)) d s

−
1

Γ(α)

∫ x
0

(x − s)α−1
[
λ1

∫ s
0

k1 (s, t) φ(t) dt (11)

+λ2

∫ 1
0

k2 (s, t) φ(t) dt
]

d sA .

3.1. Method of Solution
Collocatiing at xi in equation (11) gives

y(xi) =
N∑

k=0

y(k)(0)
k!

xk +

1
Γ(α)

∫ xi
0

(xi − s)
α−1 P2

(
φ(s)

′′
)

d sA

+
1

Γ(α)

∫ xi
0

(xi − s)
α−1 P1

(
φ(s)

′
)

d sA

+
1

Γ(α)

∫ xi
0

(xi − s)
α−1 P0 (φ(s)) d sA

−
1

Γ(α)

∫ xi
0

(xi − s)
α−1 (g(s)) d s −

1
Γ(α)

∫ xi
0

(xi − s)
α−1

[
λ1

∫ s
0

k1 (s, t) φ(t) dt (12)

+λ2

∫ 1
0

k2 (s, t) φ(t) dt
]

d sA

Applying (4) on (12) gives

φ(xi)A = W(xi)

+

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

1
Γ(α)

∫ xi
0

(xi − s)
α−1 P2

(
φ(s)

′′
)

d s+
1

Γ(α)

∫ xi
0

(xi − s)
α−1 P1

(
φ(s)

′
)

d s+
1

Γ(α)

∫ xi
0

(xi − s)
α−1 P0 (φ(s)) d s

−
1

Γ(α)

∫ xi
0

(xi − s)
α−1

(
λ1

∫ xi
0

k1 (xi, t) φ(t) dt

+λ2
∫ 1

0
k2 (xi, t) φ(t) dt

)
d s


A, (13)

where

W(xi) =
N∑

k=0

y(k)(0)
k!

xk −
1

Γ(α)

∫ xi
0

(xi − s)
α−1 g(s)d s.

Factorise the values of A from equation (13) gives

W(xi) =

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

φ(xi) −
1

Γ(α)

∫ xi
0

(xi − s)
α−1 P2

(
φ(s)

′′
)

d s−
1

Γ(α)

∫ xi
0

(xi − s)
α−1 P1

(
φ(s)

′
)

d s−
1

Γ(α)

∫ xi
0

(xi − s)
α−1 P0 (φ(s)) d s

+ 1
Γ(α)

∫ xi
0

(xi − s)
α−1(

λ1
∫ s

0
k1 (s, t) φ(t) dt

+ λ2
∫ 1

0
k2 (s, t) φ(t) dt

)
d s


A. (14)

Equation 14 can be in the form

τn(xi)A =W(xi), (15)

where

τn(xi) =


φ(xi) −

1
Γ(α)

∫ xi
0

(xi − s)
α−1 P2

(
φ(s)

′′
)

d s−
1

Γ(α)

∫ xi
0

(xi − s)
α−1 P1

(
φ(s)

′
)

d s−
1

Γ(α)

∫ xi
0

(xi − s)
α−1 P0 (φ(s)) d s +

1
Γ(α)

∫ xi
0

(xi − s)
α−1(

λ1
∫ s

0
k1 (s, t) φ(t) dt + λ2

∫ 1
0

k2 (s, t) φ(t) dt
)

d s,


and

A = [a0 a1 · · · aN ]T .
Multiply both sides of equation(15) byτ−1n (xi) gives

A =τ−1n (xi)W(xi). (16)

Substituting A into the approximate solution (4) gives

y(x) = φ(xi)τ
−1
n (xi) W(xi) (17)

3



Ajileye et al. / J. Nig. Soc. Phys. Sci. 4 (2022) 834 4

4. Uniqueness of the Method

In order to establish the uniqueness of the method, we introduce
the following hypothesis

H1 : k
∗
1 = maxx∈[0,1]

∫ x
0
λ1 |k1(x, t)|dt

H2 : k
∗
2 = maxx∈[0,1]

∫ 1
0
λ2 |k2(x, t)|dt

H3 :
∣∣∣y(m)1 − y(m)2 ∣∣∣ ≤ Lm |y1 − y2| .

Lemma (4.1) (q-contraction) Let T : X −→ X be a mapping
defined by theorem (3.0) for y1, y2 ∈ X. T is q-contraction if
and only if 

∑2
n=0 Ln −

(
k∗1 + k

∗
2

)
Γ(α + 1)

 < 1,
then there exist a unique solution of T.

(T y1) (x) = W(x) +
1

Γ(α)

∫ x
0

(x − s)α−1
(
P2y

′′

1 (s)
)

d s

+
1

Γ(α)

∫ x
0

(x − s)α−1
(
P1y

′

1(s)
)

d s

+
1

Γ(α)

∫ x
0

(x − s)α−1 (P0y1(s)) d s

−
1

Γ(α)

∫ x
0

(x − s)α−1
[
λ1

∫ s
0

k1 (s, t) y1 (t) dt

+λ2

∫ 1
0

k2 (s, t) y1 (t) dt
]

d s

and

(T y2) (x) = W(x) +
1

Γ(α)

∫ x
0

(x − s)α−1
(
P2y

′′

2 (s)
)

d s

+
1

Γ(α)

∫ x
0

(x − s)α−1
(
P1y

′

2(s)
)

d s

+
1

Γ(α)

∫ x
0

(x − s)α−1 (P0y2(s)) d s

−
1

Γ(α)

∫ x
0

(x − s)α−1
[
λ1

∫ s
0

k1 (s, t) y2 (t) dt

+λ2

∫ 1
0

k2 (s, t) y2 (t) dt
]

d s.

Thus,

|(T y1) (x) − (T y2) (x)|

=
1

Γ(α)

∫ x
0

(x − s)α−1 P2
∣∣∣y′′1 (s) − y′′2 (s)∣∣∣ d s

+
1

Γ(α)

∫ x
0

(x − s)α−1 P1
∣∣∣y′1(s) − y′2(s)∣∣∣ d s

+
1

Γ(α)

∫ x
0

(x − s)α−1 P0 |y1(s) − y2(s)|d s

−
1

Γ(α)

∫ x
0

(x − s)α−1


∫ s
0
λ1 |k1(s, t)| |y1(t) − y2(t)|dt

+
∫ 1

0
λ1 |k2(s, t)| |(y1(t) − y2(t))|dt

 d s.

Taking maximum of both sides and using H1 to H3 gives

d (T y1(x), T y2(x)) ≤

 L2 + L1 + L0 −
(
k∗1 + k

∗
2

)
Γ(α + 1)

 d(y1, y2)
Since T is a contraction

∑2
n=0 Ln −

(
k∗1 + k

∗
2

)
Γ(α + 1)

 < 1.
We can conclude that T has a unique solution.

5. Convergence Analysis

We establish the convergence of the method

yN (x) = W(x) +
1

Γ(α)

∫ x
0

(x − s)α−1
(
P2y

′′

N (s)
)

d s

+
1

Γ(α)

∫ x
0

(x − s)α−1
(
P1y

′

N (s)
)

d s

+
1

Γ(α)

∫ x
0

(x − s)α−1 (P0yN (s)) d s (18)

−
1

Γ(α)

∫ x
0

(x − s)α−1
[
λ1

∫ s
0

k1 (s, t) yN (t) dt

+λ2

∫ 1
0

k2 (s, t) yN (t) dt
]

d s.

Subtract (8) from (18) gives

EN (x) = yN (x) − y(x).

Hence

|EN (x)| =
1

Γ(α)

∫ x
0

(x − s)α−1 |EN (s) (P2 + P1)|d s

−
1

Γ(α)

∫ x
0

(x − s)α−1[
|λ1|

∣∣∣∣∣
∫ s

0
k1 (s, t) EN (t)dt

∣∣∣∣∣
+ |λ2|

∣∣∣∣∣∣
∫ 1

0
k2 (s, t) EN (t)dt

∣∣∣∣∣∣d s.
(20)

Therefore

‖EN (x)‖∞
‖EN (t)‖∞

≤
1

Γ(α)

∫ x
0

(x − s)α−1

×


|En(s) (P2 + P1)|
−

1
Γ(α)

∫ x
0

(x − s)α−1(
|λ1|

∣∣∣∫ s
0

k1 (s, t) En(t)dt
∣∣∣

+ |λ2|
∣∣∣∣∫ 10 k2 (s, t) En(t)dt∣∣∣∣

 d s.
The method converges

4



Ajileye et al. / J. Nig. Soc. Phys. Sci. 4 (2022) 834 5

6. Numerical Examples

In this section, two numerical problems with boundary
conditions are presented to text the efficiency and simplicity of
the method. All computations are done with the aid of Maple
18. Let yn(x) and y(x) be the approximate and exact solution
respectively. ErrorN = |yn(x) − y(x)|

Example 1:

Considering fractional integro-differential equation [17]

y′′ (x) +
1
x

D0.5x y(x) +
1
x2

y(x) = f (x)

+

∫ x
0

sin (x − t) y(t)dt +
∫ 1

0
cos (x − t) y (t) dt,

subject to the boundary conditions

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

where

f (x) = 5 + 1.50451x0.5 − 13x − 1.80541x1.5

−x2 + x3 − 2.067x cos(x) + 5.95385 sin(x).

Exact solution y (x) = x2 (1 − x)

Solution 1

We solve this problem at N = 3 and 5 but we use N = 3 for
demonstration. Integral form of example 1 is

y′′ (x) +
1
x

(
1

Γ(1 − 0.5)

∫ x
0

(x − t)−0.5 y(1)(t)dt
)

+
1
x2

y(x)

= f (x) +
∫ x

0
sin (x − t) y(t)dt +

∫ 1
0

cos (x − t) y (t) dt.(21)

Using approximate solution (4) on (14) gives
φ
′′

(x) + 1x
(

1
Γ(1−0.5)

∫ x
0

(x − t)−0.5 dd x (φ(t))dt
)

+ 1
x2
φ
′′

(x) −
∫ x

0
sin (x − t) φ(t)dt

−
∫ 1

0
cos (x − t) φ (t) dt

 A = f (x) (22)
Equation (15) gives

τ(x)A = f (x), (23)

where

τ(x) =
φ
′′

(x) + 1x
(

1
Γ(1−0.5)

∫ x
0

(x − t)−0.5 dd x (φ(t))dt
)

+ 1x2 φ
′′

(x) −
∫ x

0
sin (x − t) φ(t)dt −

∫ 1
0

cos (x − t) φ (t) dt.

Collocating at x3 =
[

1
3

2
3 1

]
and substituting the

boundary conditions gives

τ(x)∗A = f (x)∗, (24)

where

τi (x) =


8.1837497830 4.5298367530 3.5647166670 2.4321976680
1.0664602830 2.3357941530 3.8888778550 5.4288519160
-0.3011686789 1.5101524580 4.1068429140 8.5147669310

1 0 0 0
1 1 1 1

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

Table 1. Exact and approximate values, Example 1
x Exact N = 3 N = 5

0.2 0.032000000000 0.031999821270 0.031999851830
0.4 0.096000000000 0.095999891010 0.095999961860
0.6 0.144000000000 0.143999795600 0.143999892900
0.8 0.128000000000 0.127999546700 0.127999648400
1.0 0.000000000000 -8.444000000000e-7 -7.382000000000e-7

Table 2. Absolute Error for Example 1
x error3 error5 error [17]=32

0.2 -1.787300000000e-7 -1.481700000000e-7 2.048e-5
0.4 -1.089900000000e-7 -3.814000000000e-8 2.503e-5
0.6 -2.044000000000e-7 -1.071000000000e-8 1.789e-5
0.8 -1.071000000000e-7 -3.516000000000e-8 7.682e-6
1 -8.444000000000e-7 -7.382000000000e-8 3.034e-6

f (x) =
[

1.1325152410 -1.5399784720 -4.4079289640 0 0
]
.

We now solve for the unknown values A (17) making use of
matrix inversion results in;

y3 =
(
−4.2497275388 × 10−7 + 0.16917e − 5x+

0.9999976497x2 − 0.9999997608x3

)
.

Applying the same procedure for N = 5 gives

y5 =
(
−3.4773620430 × 10−7 + 8.9350760390 × 107 x + 1.0000017198x2

−1.0000070682x3 + 0.57092e − 5x4 − 0.16454e − 5x5

)

Example 2

Considering fractional integro-differential equation [17]

y′′ (x)+ Dαx y(x) = f (x)+2
∫ x

0
(x − t) y(t)dt +

∫ 1
0

(
x2 − t

)
y (t) dt

with the given boundary conditions

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

where

f (x) = −
1
30
− 6x +

181x2

20
+ 4x3 −

x5

10
+

x6

15
and α = 1

Exact solution y (x) = x3 (x − 1)

Solution 2

We solve this problem at N = 4 and 5, however, make use of
N = 4 for demonstration.
Integral form of example 2 is

y′′ (x) +
(

1
Γ(2 − 1)

∫ x
0

(x − t)0 y(2)(t)dt
)

(25)

= f (x) + 2
∫ x

0
(x − t) y(t)dt +

∫ 1
0

(
x2 − t

)
y (t) dt

5



Ajileye et al. / J. Nig. Soc. Phys. Sci. 4 (2022) 834 6

Table 3. Exact and approximate values, Example 2
x Exact N = 4 N = 5

0.2 -0.6400000000e-2 -0.6400000000e-2 -0.6400000000e-2
0.4 -0.3840000000e-1 -0.3840000000e-1 -0.3840000000e-1
0.6 -0.864000000e-1 -0.8640000000e-1 -0.8640000000e-1
0.8 -0.1024000000 -0.1024000000 -0.1024000000
1.0 -0.9000000000e-3 -0.900000000e-3 -0.900000000e-3

Table 4. Absolute Error for Example 2
x error4 error5 error [17]=64

0.2 0.0 0.0 2.931e-7
0.4 0.0 0.0 3.915e-7
0.6 0.0 0.0 2.696e-7
0.8 0.0 0.0 1.011e-7
1 0.0 0.0 3.596e-8

Using approximate solution (4) on (18) gives φ′′(x) +
(

1
Γ(2−1)

∫ x
0

(x − t)0 d
2

d x2
(φ(t))dt

)
−2

∫ x
0

(x − t) φ(t)dt −
∫ 1

0

(
x2 − t

)
φ (t) dt

 A = f (x). (26)
Equation (6.6) gives

τ(x)A = f (x) (27)

where

τ(x) =
φ
′′

(x) +
(

1
Γ(2−1)

∫ x
0

(x − t)0 d
2

d x2 (φ(t))dt
)

−2
∫ x

0
(x − t) φ(t)dt −

∫ 1
0

(
x2 − t

)
φ (t) dt.

Collocating at x4 =
[
0 14

2
4

3
4 1

]
and

substituting the boundary conditions gives

τ(x)∗A = f (x)∗ (28)

where

τi (x) =

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

1
2 0.3333333333 2.250000000 0.2000000000 0.1666666667
1
2

59
192

4193
1536

19169
10240

59393
61440

1
2

1
4

305
96

249
64

3473
960

1
2

37
192

1851
512

64211
10240

522457
61440

1
2

1
6

49
12

181
20

481
30

1 0 0 0 0
1 1 1 1 1


f (x) =

[
- 130 -

55621
61440 -

131
480

137191
61440

419
60 0 0

]
We solve for the unknown values A using matrix inversion,
results;

y4 =


−2.2211399390 × 10−13

+3.0783153800 × 10−12 x
−1.4125589590 × 10−11 x2

−1.0000000000x3 + 1.0000000000x4


Applying the same procedure for N = 5 gives

y5 =



5.329070518000 × 10−15

+1.781685910000 × 10−12 x
−1.580957587000 × 10−11 x2

−1.000000000000x3

+1.000000000000x4

+1.062971933000 × 10−11 x5



7. Conclusion

Collocation approach is utilized to solving the fractional order
Fredholm-Volterra integro-differential problems in this paper.
The applied method is consistent, efficient and easy to compute.
The results of the numerical example 1 as shown in table 1
shows that the approximate solution at N = 3 gives y3(x) =
−4.2497275388 × 10−7 + 0.16917e − 5x + 0.9999976497x2 −
0.9999997608x3 and solving at N = 5 indicates that as the value
of N increases, the error becomes smaller. We also compare our
absolute errors with Alkan et al (2017) as shown in table 2, this
clearly shows that our method performs better. The results of
the numerical example 2 as shown in table 3 shows that the
approximate solution obtained at N = 4 and 5 converges to
the exact solution. Table 4 compares the errors in the method
proposed by Alkan et al (2017) at N = 64 and errors in our new
method which converges to the exact solution at N = 4.

References

[1] P. Kamal, Higher order numerical methods for fractional order
differential equations, Doctoral dissertation, University of Chester,
United Kingdom, 2015. http://hdl.handle.net/10034/613354.

[2] D. Baleanu, K. Diethelm & J. J. Trujillo, Fractional calculus: Models and
numerical methods, Series on complexity nonlinearity and chaos, World
Scientific 2012.

[3] R. H. Khan & H. O. Bakodah, “Adomian decomposition method and
its modification for nonlinear Abel’s integral equations”, Computers and
Mathematics with Applications 7 (2013) 2349.

[4] R. C. Mittal & R. Nigam, “Solution of fractional integro-differential
equations by Adomiandecomposition method”, The International Journal
of Applied Mathematics and Mechanics 2 (2008) 87.

[5] S. Momani & M. A. Noor, “Numerical methods for fourth-order fractional
integro-differential equations”, Applied Mathematics and Computation 1
(2006) 754.

[6] A. O. Agbolade & T. A. Anake, “Solution of first order volterra linear
integro differential equations by collocation method”, J. Appl. Math.,
Article ID. 1510267 (2017). doi:10.1155/2017/1510267

[7] O. Cem, S. Mehmet & D. O. Arzu, “Chelyshkov collocation approach
to solving the systems of linear functional differential equations”, New
Trend in Mathematical Sciences (NTMSCI) 4 (2015) 83.

[8] D. Nazari & M. Jahanshahi, “Some notes on multi-order frational
integro-differential equations”, 43rd Annual Iranian Mathematics
Conference (2012) pp 27.

[9] D. S. Mohammed, ”Numerical solution of fractional integro-differential
equations by least square method and shifted Chebyshev polynomial”,
Mathematical Problems in Engineering 2014 (2014) 431965.
https://doi.org/10.1155/2014/431965

[10] L. Huang, X. F. Li, Y. Zhao & X. Y. Duan, “Approximate solution of
fractional integro-differential equations by Taylor expansion method”,
Computers & Mathematics with Applications 3 (2011) 1127.

[11] A. Arikoglu & I. Ozkol, ”Solution of fractional integro-differential
equations by using fractional differential transform method”, Chaos,
Solitons & Fractals 2 (2009) 521.

[12] A. A. Elbeleze, A. Klman & B. M. Taib, “Approximate solution of
integro-differential equation of fractional (arbitrary) order”, Journal of
King Saud University - Science 28 (2016) 61.

[13] Y. Nawaz, “Variational iteration method and homotopy perturbation
method for fourth-order fractional integro-differential equations”,
Computers & Mathematics with Applications 8 (2011) 2330.

[14] R. K. Saeed & H. M. Sdeq, “Solving a system of linear fredholm
fractional integro-differential equations using homotopy perturbation
method”, Australian Journal of Basic and Applied Sciences 4 (2010) 633.

[15] S. Alkan, “A new solution method for nonlinear fractional
integro-differential equations”, Discrete and Continuous Dynamical
Systems - Series S 6 (2015) 1065.

6



Ajileye et al. / J. Nig. Soc. Phys. Sci. 4 (2022) 834 7

[16] S. Alkan & A. Secer, “Solution of nonlinear fractional boundary value
problems with nonho-mogeneous boundary conditions”, Applied and
Computational Mathematics 3 (2015).

[17] S. Alkan & V. F. Hatipoglu, “Approximate solutions of Volterra-Fredholm
integro-differential equations of fractional order”, Tbilisi Mathematical
Journal 2 (2017) 1.

[18] O. A. Uwaheren, A. F. Adebisi & O. A. Taiwo, “Perturbed Collocation
Method For Solving Singular Multi-order Fractional Differential
Equations of Lane-Emden Type”, Journal of the Nigerian Society of
Physical Sciences 3 (2020) 141. https://doi.org/10.46481/jnsps.2020.69.

[19] I. Avcı & N. I. Mahmudov, “Numerical Solutions for Multi-Term

Fractional Order Differential Equations with Fractional Taylor
Operational Matrix of Fractional Integration”, Mathematics 1 (2020) 96.
https://doi.org/10.3390/math8010096.

[20] L. Linna, W. Zhirou & H. Qiongdan, “A numerical method for
solving feactional variational problems by operational matrix based on
Chelyshekov polynomial”, Enginerring Letters 28 (2020) 486.

[21] S. E. Fadugba, “Solution of Fractional Order Equations in the Domain
of the Mellin Transform”, Journal of the Nigerian Society of Physical
Sciences 4 (2019) 138. https://doi.org/10.46481/jnsps.2019.31

7