EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS
Vol. 17, No. 3, 2024, 1429-1448
ISSN 1307-5543 – ejpam.com
Published by New York Business Global

A Numerical Method for Investigating Fractional
Volterra-Fredholm Integro-Differential Model

Muhammed I. Syam1,∗, Mwaffag Sharadga2, Ishak Hashim2,3

1 Department of Mathematical Sciences, United Arab Emirates University, Al-Ain,
United Arab Emirates
2 Department of Mathematical Sciences, Faculty of Science and Technology,
Universiti Kebangsaan Malaysia, 43600 UKM Bangi Selangor, Malaysia
3 Nonlinear Dynamics Research Center (NDRC), Ajman University, Ajman PO Box 346,
United Arab Emirates

Abstract. In this article, we investigate the fractional Volterra-Fredholm integro-differential equa-
tions. These equations appear in several applications such as control theory, biology, and particle
dynamics in physics. We derive a numerical method based on the operational matrix method to
solve this class of integro-differential equations. We prove the existence and uniqueness of the ex-
act solution. Additionally, we demonstrate the uniform convergence of the numerical solutions to
the exact solution. We present several numerical examples to show the numerical efficiency of the
proposed method. In the first example, we choose a linear problem and find that the approximate
solution converges to the exact solution when the number of block pulse functions is very large.
In the next two examples, we consider the nonlinear case and compute the L2-local truncation
error since exact solutions are not available. The error was of order 10−12. Furthermore, we sketch
the graph of the approximate solutions for different values of the fractional derivative to observe
the influence of the fractional derivative on the profile of the solutions. Theoretical and numerical
results show that the proposed method is accurate and can be applied to other nonlinear problems
in science.

2020 Mathematics Subject Classifications: 74S40, 65D15

Key Words and Phrases: Block pulse function, Fractional derivative, Nonlinear dynamics,
Volterra-Fredholm integro-differential

1. Introduction

Fractional calculus, a branch of mathematical analysis, extends the concepts of differ-
entiation and integration to non-integer orders. The numerous applications it has found
in physics, chemistry, biology, and control theory have garnered considerable attention.

∗Corresponding author.
DOI: https://doi.org/10.29020/nybg.ejpam.v17i3.5280

Email addresses: m.syam@uaeu.ac.ae (M. Syam),
P101953@siswa.ukm.edu.my (M. Sharadga), ishak h@ukm.edu.my (I. Hashim)

https://www.ejpam.com 1429 © 2024 EJPAM All rights reserved.



M.I. Syam, M. Sharadga, I. Hashim / Eur. J. Pure Appl. Math, 17 (3) (2024), 1429-1448 1430

For example, fractional calculus is applied in anomalous diffusion, viscoelasticity, and elec-
tromagnetism in physics, as well as control systems and signal processing in engineering.
For more details, see [9, 10, 14, 16]. Incorporating fractional derivatives and integrals
is one of the valuable tools FC provides for simulating complex systems with memory
and long-term interactions. This theory has revolutionized the approach to many previ-
ously challenging problems and has found applications in modeling anomalous diffusion,
viscoelasticity, fractional control systems, and more.

The Caputo derivative, introduced by Michele Caputo in 1967, stands out for its
practicality. The Caputo derivative is more suitable for real-world simulations because it
doesn’t require the function to be different throughout the domain, unlike the Riemann-
Liouville fractional derivative.

Phenomena like fractional control systems, viscoelasticity, and atypical diffusion are
all modeled using the Caputo derivative in engineering and research. It is used for mathe-
matical modeling of physical and biological systems such as heat transfer and groundwater
flow. The following references are recommended for further reading [6, 11, 12].

A powerful technique for solving fractional differential equations is the operational
matrix method. A linear combination of basis functions, usually starting with block pulse
functions, is used to approximate the solution. These functions are constructed using
operational matrixes, which represent differentiation and integration operations, and the
coefficients of the linear combination are determined by resolving algebraic equations. The
OMM has been successfully applied to various issues, including fractional delay equations
[15], Ricatti equations [5], systems of differential equations [8], and nonlinear differential
equations [1]. Researchers investigated alternative basis functions to enhance computa-
tional effectiveness and precision [2–4, 7].

In this article, we will study the following class of integro-differential equations:

DµΩ(t) = G1(t,Ω(t)) +G2(t) +

∫ t

0
Π1(t, s)G3(Ω(s))ds+

∫ 1

0
Π2(t, s)G4(Ω(s))ds, (1)

Ω(0) = ω0, (2)

where t ∈ [0, 1], Π1,Π2 ∈ C([0, 1]×[0, 1]),G1 ∈ C([0, 1]×ℜ), G3, G4 ∈ C(ℜ), G2 ∈ C([0, 1]),
and 0 < µ ≤ 1. The derivative her is in the Caputo sense.

The fractional Volterra-Fredholm integro-differential model (FVFIDM) represents a
significant extension of classical Volterra and Fredholm integral equations by incorporat-
ing fractional calculus concepts. This model plays a crucial role in various fields such
as physics, engineering, biology, and finance, where systems exhibit memory effects, long-
range interactions, and complex dynamics. Volterra and Fredholm integrodifferential equa-
tions are two types of integral equations with important applications in various fields such
as physics, engineering, and applied mathematics. Although they are related, they differ in
their formulation and properties. Volterra integrodifferential equations involve an integral
that extends over a variable limit of integration, typically from a fixed starting point to
the independent variable, while Fredholm integrodifferential equations involve an integral
with fixed limits of integration that do not depend on the independent variable.



M.I. Syam, M. Sharadga, I. Hashim / Eur. J. Pure Appl. Math, 17 (3) (2024), 1429-1448 1431

Applications of FVFIDMs are widespread, ranging from modeling viscoelastic mate-
rials and anomalous diffusion processes to analyzing population dynamics in ecology and
pricing financial derivatives. Understanding the behavior of FVFIDMs is essential for
predicting and controlling complex systems in various scientific and engineering domains.
Fore more details, see [6, 11, 12].

The structure of our paper is listed below. Several definitions and lemmas will be
introduced in the next section. Section three will develop a version of the OMM to solve
the proposed problem. Some theoretical outcomes, such as existence and uniqueness,
estimations of errors, and convergence of the solution, will be demonstrated in section
4. Three illustrative examples will be presented in section 5. Numerical validation of
the proposed method’s convergence to the unique solution of our problems provided by
these examples. Finally, we will draw conclusions and provide closing remarks in the final
section.

2. Foundational Concepts

This section presents various fundamental concepts and results employed within this
paper.

Definition 1. [6, 11, 12] For µ ∈ (0, 1) and t > 0, the Caputo derivative of Ω(t) is defined
as follows

DµΩ(t) =
1

Γ(1− µ)

∫ t

0
(t− τ)−µΩ′(τ)dτ, (3)

and the fractional integral operator is given by

IµΩ(t) =
1

Γ(µ)

∫ t

0
(t− τ)µ−1Ω(τ)dτ. (4)

The fractional fundamental theorem of calculus is given as

Lemma 1. [13, 20] For µ > 0 ∈ (0, 1), and Ω(t) ∈ C[0, 1], the following equations hold

IµDµΩ(t) = Ω(t)− Ω(0), (5)

DµIµΩ(t) = Ω(t). (6)

Another important concept in this paper is the block pulse function (BPF) which is
defined as follows.

Definition 2. [15, 18, 19] Let M be a postive integer. tr = jr, j ∈ 0, 1, 2, ...,M − 1, r = 1
M .

Then, the j-block pulse function is give by

γj(t) =

{
1, t ∈ [tj , tj+1)

0, [0, 1]− [tj , tj+1)
, 0 ≤ s < M. (7)

In the next theorem, we present two properties of the BPFs which are the disjoint and
orthogonal properties.



M.I. Syam, M. Sharadga, I. Hashim / Eur. J. Pure Appl. Math, 17 (3) (2024), 1429-1448 1432

Theorem 1. [5, 17, 21] Then, for any i, j ∈ {0, 1, . . . ,M − 1}, we have

γi(t)γj(t) =

{
γi(t), if i = j

0, if i ̸= j
, (8)

and ∫ 1

0
γi(t)γj(t)dt =

{
r, if i = j

0, if i ̸= j
. (9)

For numerical purposes, we present the following important lemma.

Lemma 2. [5, 17, 21] If Ω ∈ L2[0, 1], then

Ω(t) = lim
M→∞

M−1∑
j=0

ωjγj(t), (10)

where

ωj =
1

r

∫ (j+1)r

jr
Ω(t)dt. (11)

For the approximation purposes, we choose M large enough. For such m, we can
rewrite Ω in the matrix form as

Ω(t) ≈ Ω
T
γ(t) (12)

where

Ω =


ω0

ω1
...

ωM−1

 , γ(t)


γ0(t)
γ1(t)
...

γM−1(t)

 . (13)

3. Method of Solution

In this section, we present the method of solution for problem (1)-(2). Let us assume
that

I1(t) =

∫ t

0
Π1(t, s)G3(Ω(s))ds (14)

and

I2(t) =

∫ 1

0
Π2(t, s)G4(Ω(s))ds. (15)

Let us approximate Ω(t) by

Ω(t) ≈
M−1∑
i=0

ωiγi(t). (16)



M.I. Syam, M. Sharadga, I. Hashim / Eur. J. Pure Appl. Math, 17 (3) (2024), 1429-1448 1433

Let

G3(Ω(t)) =

n3∑
i=0

a3,iΩ
i(t), G4(Ω(t)) =

n4∑
i=0

a4,iΩ
i(t), (17)

Π1(t, s) =

M−1∑
i=0

M−1∑
j=0

b1,i,jγi(t)γj(s), Π2(t, s) =

M−1∑
i=0

M−1∑
j=0

b2,i,jγi(t)γj(s). (18)

Using Equation (8) and by substituting Equation (16) into Equation (17), we get

G3(Ω(t)) =

n3∑
i=0

a3,iΩ
i(t)

=

n3∑
i=0

a3,i

M−1∑
j=0

ωjγj(t)

i

=

n3∑
i=0

a3,i

M−1∑
j=0

ωi
jγj(t)


=

n3∑
i=0

M−1∑
j=0

a3,iω
i
jγj(t) (19)

and

G4(Ω(t)) =

n4∑
i=0

a4,iΩ
i(t)

=

n4∑
i=0

a4,i

M−1∑
j=0

ωjγj(t)

i

=

n4∑
i=0

a4,i

M−1∑
j=0

ωi
jγj(t)


=

n4∑
i=0

M−1∑
j=0

a4,iω
i
jγj(t). (20)

Substitute Equation (19) and (18) into Equation (14) and using Equation (8) to get

I1(t) =

∫ t

0
Π1(t, s)G3(Ω(s))ds

=

∫ t

0

M−1∑
i=0

M−1∑
j=0

b1,i,jγi(t)γj(s)ds

 n3∑
i=0

M−1∑
j=0

a3,iω
i
jγj(s)





M.I. Syam, M. Sharadga, I. Hashim / Eur. J. Pure Appl. Math, 17 (3) (2024), 1429-1448 1434

=
M−1∑
i=0

M−1∑
j=0

n3∑
l=0

M−1∑
k=0

b1,i,ja3,lω
l
kγj(t)

∫ t

0
γj(s)γk(s)ds

=

M−1∑
i=0

M−1∑
j=0

n3∑
l=0

b1,i,ja3,lω
l
jγj(t)

∫ t

0
γj(s)ds. (21)

Using the definition of the BPFs, we can see that

∫ t

0
γj(s)ds =


∫ t
0 0ds, t < jr∫ t
jr 1ds, r < t < (j + 1)r∫ (j+1)r
jr 1ds, t ≥ (j + 1)r

=


0, t < jr

(t− jr), jr < t < (j + 1)r

1, t ≥ (j + 1)r

. (22)

Using Lemma (2), we have ∫ t

0
γj(s)ds =

M−1∑
k=0

θkγk(t) (23)

where

θk =


0, k < j
r
2 , k = j

r, j < k < M

. (24)

Thus, Equation (21) becomes

I1(t) =
M−1∑
i=0

M−1∑
j=0

n3∑
l=0

b1,i,ja3,lω
l
jγj(t)

(
M−1∑
k=0

θkγk(t)

)

=

M−1∑
i=0

M−1∑
j=0

n3∑
l=0

b1,i,ja3,lω
l
jθjγj(t)

= P1(Ω)
Tγ(t) (25)

where

(P1(Ω))j =
M−1∑
i=0

M−1∑
j=0

n3∑
l=0

b1,i,ja3,lω
l
jθj . (26)

Similarly, Substitute Equation (19) and (20) into Equation (15) and using Equation (8)
to get

I2(t) =

∫ 1

0
Π2(t, s)G4(Ω(s))ds



M.I. Syam, M. Sharadga, I. Hashim / Eur. J. Pure Appl. Math, 17 (3) (2024), 1429-1448 1435

=

∫ 1

0

M−1∑
i=0

M−1∑
j=0

b2,i,jγi(t)γj(s)ds

 n4∑
i=0

M−1∑
j=0

a4,iω
i
jγj(s)


=

M−1∑
i=0

M−1∑
j=0

n4∑
l=0

M−1∑
k=0

b2,i,ja4,lω
l
kγj(t)

∫ 1

0
γj(s)γk(s)ds

=

M−1∑
i=0

M−1∑
j=0

n4∑
l=0

b2,i,ja4,lω
l
jγj(t)

∫ 1

0
γj(s)ds. (27)

Using the definition of the BPFs, we can see that∫ 1

0
γj(s)ds = r. (28)

Using Lemma (2), we have ∫ 1

0
γj(s)ds =

M−1∑
k=0

rγk(t). (29)

Thus, Equation (27) becomes

I2(t) =

M−1∑
i=0

M−1∑
j=0

n3∑
l=0

b1,i,ja3,lω
l
jγj(t)

(
M−1∑
k=0

rγk(t)

)

=
M−1∑
i=0

M−1∑
j=0

n3∑
l=0

b1,i,ja3,lω
l
jrγj(t)

= P2(Ω)
Tγ(t) (30)

where

(P2(Ω))j =
M−1∑
i=0

M−1∑
j=0

n3∑
l=0

b1,i,ja3,lω
l
jr. (31)

Let

G1(t,Ω(t)) =

n1∑
i=0

b1,i(t)Ω
i(t)

=

n1∑
i=0

M−1∑
j=0

a1,i,jγj(t)

M−1∑
j=0

ωjγj(t)

i

=

n1∑
i=0

M−1∑
j=0

a1,i,jγj(t)

M−1∑
j=0

ωi
jγj(t)


=

n1∑
i=0

M−1∑
j=0

a1,i,jω
i
jγj(t)



M.I. Syam, M. Sharadga, I. Hashim / Eur. J. Pure Appl. Math, 17 (3) (2024), 1429-1448 1436

= P3(Ω)
Tγ(t) (32)

where

(P3(Ω))j =

n1∑
i=0

a1,i,jω
i
j . (33)

Finally, we write G2(t) as

G2(t) =
M−1∑
i=0

a2,iγi(t),

= P T
4 γ(t) (34)

where
(P4)i = a2,i. (35)

Now, we will derive the the Riemann-Liouville fractional integral operator which plays an
important role in our derivation.

Theorem 2. The operational matrix of Iµ is given by

Oµ =
rµ

Γ(µ+ 2)



1 σ1 σ2 . . . σm−2 σM−1

0 1 σ1 . . . σM−3 σM−2

0 0 1
. . . σM−4 σM−3

...
...

...
. . .

. . .
...

0 0 0 . . . 1 σ1
0 0 0 . . . 0 1


(36)

where σς = (ς + 1)µ+1 − 2ςµ+1 + (ς − 1)µ+1, ς = 1, 2, . . . ,M − 1.

Proof. Let l ∈ {0, 1, . . . ,M − 1}. Then,

Iµγl(t) =
1

Γ(µ)

∫ t

0
(t− s)µ−1γl(s)ds

=


0, t < lr
(t−lr)µ

Γ(µ+1) , lr ≤ t < (l + 1)r
(t−lr)µ−(t−lr−r)µ

Γ(µ+1) , (l + 1)r ≤ t < 1

. (37)

Let

Iµγl(t) =

M−1∑
i=0

ci,lγi(t). (38)

Then,

ci,l =
1

r

∫ 1

0
(Iµγl(t)) γi(t)dt



M.I. Syam, M. Sharadga, I. Hashim / Eur. J. Pure Appl. Math, 17 (3) (2024), 1429-1448 1437

=
1

r

∫ (i+1)r

ir
(Iµγl(t)) dt

=


rµ

Γ(µ+2) , 0 ≤ i = l ≤ M − 1
rµ((l−i+1)µ+1−2(l−i)µ+1+(l−i−1)µ+1)

Γ(µ+2) , 0 ≤ i < l ≤ M − 1

0, 0 ≤ l < i ≤ M − 1

. (39)

Let ς = j − i and σς = (ς + 1)µ+1 − 2ςµ+1 + (ς − 1)µ+1. Then, the operational matrix of
Iµ is

Oµ =
rµ

Γ(µ+ 2)



1 σ1 σ2 . . . σm−2 σM−1

0 1 σ1 . . . σM−3 σM−2

0 0 1
. . . σM−4 σM−3

...
...

...
. . .

. . .
...

0 0 0 . . . 1 σ1
0 0 0 . . . 0 1


. (40)

Now, we can rewrite Equation (1) as

DµΩ(t) =
(
P3(Ω) + P4 + P1(Ω) + P2(Ω)

)T
γ(t). (41)

Using Lemma (1), we have

Ω(t) = Ω(0) +
(
P3(Ω) + P4 + P1(Ω) + P2(Ω)

)T
Iµγ(t)

= ω0 +
(
P3(Ω) + P4 + P1(Ω) + P2(Ω)

)T
Iµγ(t)

=
(
ω0P5 +

(
P3(Ω) + P4 + P1(Ω) + P2(Ω)

)T
Oµ

)
γ(t) (42)

where
(P5)i = 1. (43)

Hence,

Ω
T
γ(t) =

(
ω0P5 +

(
P3(Ω) + P4 + P1(Ω) + P2(Ω)

)T
Oµ

)
γ(t). (44)

Using the orthogonality property of the BPFs, Equation (44) becomes

ω0P5 +
(
P3(Ω) + P4 + P1(Ω) + P2(Ω)

)T
Oµ − Ω = 0. (45)

Then, we solve the algebraic System (45) to find the coefficients Ω of the approximate
solution.

We can summarize (OMM) as follows:
Algorithm 1

(i) Find the operational matrices for integral operators.



M.I. Syam, M. Sharadga, I. Hashim / Eur. J. Pure Appl. Math, 17 (3) (2024), 1429-1448 1438

(ii) Approximate the solution in terms of the Block Pulse Functions (BPF).

(iii) Take the integral of both sides of the proposed problem.

(iv) Substitute the operational matrices to generate an algebraic system.

(v) Solve the algebraic system to obtain the coefficients of the solution.

4. Theoretical Results

In the first theorem, we want to prove that Problem (1)-(2) has a unique solution.

Theorem 3. Let Π1,Π2 : [0, 1] × [0, 1] → ℜ and G2 : [0, 1] → ℜ be continuous functions
which are bounded by Q1 and Q2, respectively, G3, G4 : ℜ → ℜ be continuous Lipschitz
functions with Lipschitz constants L3 and L4, respectively, and G1 : [0, 1] × ℜ → ℜ be
continuous Lipschitz function with respect to the second component with Lipschitz constant
L1. Then, the the following problem

DµΩ(t) = G1(t,Ω(t)) +G2(t) +

∫ t

0
Π1(t, s)G3(Ω(s))ds

+

∫ 1

0
Π2(t, s)G4(Ω(s))ds, (46)

Ω(0) = ω0, (47)

has a unique solution if

L1(µ+ 1) + (L3Q1 + L4Q2)

Γ(µ+ 2)
< 1. (48)

Proof. Take the fractional integral operator for both sides of Equation (46) to get

Ω(t)− ω0 =
1

Γ(µ)

∫ t

0

(
G1(z,Ω(z)) +G2(z) +

∫ z

0
Π1(z, s)G3(Ω(s))ds

)
(t− z)µ−1dz

+
1

Γ(µ)

∫ t

0

(∫ 1

0
Π2(z, s)G4(Ω(s))ds

)
(t− z)µ−1dz. (49)

Let

℘(Ω) = ω0 +
1

Γ(µ)

∫ t

0

(
G1(z,Ω(z)) +G2(z) +

∫ z

0
Π1(z, s)G3(Ω(s))ds

)
(t− z)µ−1dz

+
1

Γ(µ)

∫ t

0

(∫ 1

0
Π2(z, s)G4(Ω(s))ds

)
(t− z)µ−1dz. (50)

Therefor,

|℘(Ω1)− ℘(Ω2)| ≤ 1

Γ(µ)
|
∫ t

0
(G1(z,Ω1(z))−G1(z,Ω2(z))) (t− z)µ−1dz |



M.I. Syam, M. Sharadga, I. Hashim / Eur. J. Pure Appl. Math, 17 (3) (2024), 1429-1448 1439

+
1

Γ(µ)
|
∫ t

0

∫ z

0
Π1(z, s) (G3(Ω1(s))−G3(Ω2(s))) ds(t− z)µ−1dz |

+
1

Γ(µ)
|
∫ t

0

∫ 1

0
Π2(z, s) (G4(Ω1(s))−G4(Ω2(s))) ds(t− z)µ−1dz | .

Since G1, G3, G4 are Lipschitz functions, then

|℘(Ω1)− ℘(Ω2)| ≤ L1 ∥ Ω1 − Ω2 ∥
Γ(µ)

|
∫ t

0
(t− z)µ−1dz |

+
L3 ∥ Ω1 − Ω2 ∥

Γ(µ)
|
∫ t

0

∫ z

0
Π1(z, s)ds(t− z)µ−1dz |

+
L4 ∥ Ω1 − Ω2 ∥

Γ(µ)
|
∫ t

0

∫ 1

0
Π2(z, s)ds(t− z)µ−1dz | . (51)

Since Π1 and Π2 are continuous on a compact set [0, 1]× [0, 1], then they are bounded by
Q1 and Q2, respectively. Thus,

|℘(Ω1)− ℘(Ω2)| ≤ L1 ∥ Ω1 − Ω2 ∥
Γ(µ)

|
∫ t

0
(t− z)µ−1dz |

+
L3Q1 ∥ Ω1 − Ω2 ∥

Γ(µ)
|
∫ t

0

∫ z

0
ds(t− z)µ−1dz |

+
L4Q2 ∥ Ω1 − Ω2 ∥

Γ(µ)
|
∫ t

0

∫ 1

0
ds(t− z)µ−1dz |

≤
(

L1t
µ

Γ(µ+ 1)
+

(L3Q1 + L4Q2)t
µ+1

Γ(µ+ 2)

)
∥ Ω1 − Ω2 ∥

≤
(
L1(µ+ 1) + (L3Q1 + L4Q2)

Γ(µ+ 2)

)
∥ Ω1 − Ω2 ∥ . (52)

Since L1(µ+1)+(L3Q1+L4Q2)
Γ(µ+2) < 1, then ℘ is contraction on Ω. By Banach fixed point theorem,

Problem (46)-(47) has unique solution.

Next, our target is to prove that the sequence
{∑M−1

i=0 ωiγi(t)
}∞

M=1
is uniformly con-

vergent to the unique solution of Problem (1)-(2) on [0, 1].

Theorem 4. Let Π1,Π2 : [0, 1] × [0, 1] → ℜ and G2 : [0, 1] → ℜ be continuous functions
which are bounded by Q1 and Q2, respectively, G3, G4 : ℜ → ℜ be continuous Lipschitz
functions with Lipschitz constants L3 and L4, respectively, and G1 : [0, 1] × ℜ → ℜ be
continuous Lipschitz function with respect to the second component with Lipschitz constant
L1, then the sequence {

M−1∑
i=0

ωiγi(t)

}∞

M=0

(53)

converges to the unique solution of Problem (1)-(2) if

L1(µ+ 1) + (L3Q1 + L4Q2)

Γ(µ+ 2)
< 1. (54)



M.I. Syam, M. Sharadga, I. Hashim / Eur. J. Pure Appl. Math, 17 (3) (2024), 1429-1448 1440

Proof. Let Ω(t) be the exact solution of Problem (1)-(2). Using Lemma (2), Ω can be
written as

Ω(t) =

∞∑
i=0

ωiγi(t). (55)

For any 0 ≤ i, j ≤ M − 1, we have

γi(tj) =

{
1, i = j

0, i ̸= j
. (56)

Then,

Ω(tj) =

∞∑
i=0

ωiγi(tj) = ωj . (57)

Using Equation (49) and (57), we have

ωj = ω0 +
1

Γ(µ)

∫ tj

0

(
G1(z,Ω(z)) +G2(z) +

∫ z

0
Π1(z, s)G3(Ω(s))ds

)
(t− z)µ−1dz

+
1

Γ(µ)

∫ tj

0

(∫ 1

0
Π2(z, s)G4(Ω(s))ds

)
(t− z)µ−1dz. (58)

Let ΩM (t) be the approximate solution of Problem (1)-(2). Then,

ΩM (t) =
M−1∑
i=0

ωiγi(t). (59)

Using previous argument, we have

ωj = ω0 +
1

Γ(µ)

∫ tj

0

(
G1(z,ΩM (z)) +G2(z) +

∫ z

0
Π1(z, s)G3(ΩM (s))ds

)
(t− z)µ−1dz

+
1

Γ(µ)

∫ tj

0

(∫ 1

0
Π2(z, s)G4(ΩM (s))ds

)
(t− z)µ−1dz. (60)

From Equations (58) and (60), we have

| ωj − ωj | ≤ 1

Γ(µ)
|
∫ tj

0
(G1(z,Ω(z))−G1(z,ΩM (z))) (t− z)µ−1dz |

+
1

Γ(µ)
|
∫ tj

0

∫ z

0
Π1(z, s) (G3(Ω(s))−G3(ΩM (s))) ds(t− z)µ−1dz |

+
1

Γ(µ)
|
∫ tj

0

∫ 1

0
Π2(z, s) (G4(Ω(s))−G4(ΩM (s))) ds(t− z)µ−1dz | .

Thus,

| ωj − ωj | ≤ L1 ∥ Ω− ΩM ∥
Γ(µ)

|
∫ tj

0
(tj − z)µ−1dz |



M.I. Syam, M. Sharadga, I. Hashim / Eur. J. Pure Appl. Math, 17 (3) (2024), 1429-1448 1441

+
L3 ∥ Ω− ΩM ∥

Γ(µ)
|
∫ tj

0

∫ z

0
Π1(z, s)ds(tj − z)µ−1dz |

+
L4 ∥ Ω− ΩM ∥

Γ(µ)
|
∫ tj

0

∫ 1

0
Π2(z, s)ds(tj − z)µ−1dz | . (61)

Since Π1 and Π2 are continuous on a compact set [0, 1]× [0, 1], then

| ωj − ωj | ≤ L1 ∥ Ω− ΩM ∥
Γ(µ)

|
∫ tj

0
(tj − z)µ−1dz |

+
L3Q1 ∥ Ω− ΩM ∥

Γ(µ)
|
∫ tj

0

∫ z

0
ds(tj − z)µ−1dz |

+
L4Q2 ∥ Ω− ΩM ∥

Γ(µ)
|
∫ tj

0

∫ 1

0
ds(tj − z)µ−1dz |

≤
(

L1tj
µ

Γ(µ+ 1)
+

(L3Q1 + L4Q2)tj
µ+1

Γ(µ+ 2)

)
∥ Ω− ΩM ∥

≤
(
L1(µ+ 1) + (L3Q1 + L4Q2)

Γ(µ+ 2)

)
∥ Ω− ΩM ∥ . (62)

Therefore,

| Ω(t)− ΩM (t) | = ∥
M−1∑
i=0

(ωi − ωi)γi(t) +

∞∑
j=M

ωjγj(t)∥

≤ ∥
M−1∑
i=0

(ωi − ωi)γi(t)∥+ ∥
∞∑

j=M

ωjγj(t)∥

≤
M−1∑
i=0

(
L1(µ+ 1) + (L3Q1 + L4Q2)

Γ(µ+ 2)

)
∥ Ω− ΩM ∥ +∥

∞∑
j=M

ωjγj(t) ∥ .

Since the series in Equation (55) converges uniformly in [0, 1], then for ϵ = 1, there exist
positive integer N1 such that

∥
∞∑

j=M

ωjγj(t) ∥≤ 1, M ≥ N1, t ∈ [0, 1]. (63)

Therefore,

| Ω(t)− ΩM (t) | ≤
(
L1(µ+ 1) + (L3Q1 + L4Q2)M

Γ(µ+ 2)

)
∥ Ω− ΩM ∥ +1. (64)

Take the supreme over [0, 1] to get

∥ Ω− ΩM ∥ ≤
(
L1(µ+ 1) + (L3Q1 + L4Q2)M

Γ(µ+ 2)

)
∥ Ω− ΩM ∥ +1 (65)



M.I. Syam, M. Sharadga, I. Hashim / Eur. J. Pure Appl. Math, 17 (3) (2024), 1429-1448 1442

which can simplify as

∥ Ω− ΩM ∥≤ Γ(µ+ 2)

Γ(µ+ 2)− L1(µ+ 1) + (L3Q1 + L4Q2)M
→ 0 (66)

as M approaches to infinity. Hence, {ΩM (t)}∞M=1 converges uniformly to the unique
solution of Problem (1)-(2) on [0, 1].

5. Numerical Results

In this section, we will present three examples to demonstrate the effectiveness of the
proposed method.

Example 1. Consider the following class of integro-differential equations:

DµΩ(t) =
−t2et

3
Ω(t) +

t1−µ

Γ(2− µ)
− t2

2
+

∫ t

0
setΩ(s)ds+

∫ 1

0
t2Ω(s)ds,

Ω(0) = 0,

where t ∈ [0, 1], and 0 < µ ≤ 1. Then, the exact solution is Q(t) = t. We approximate
Ω(t) as follows

Ω(t) =
M−1∑
i=0

ωjµj(t).

Then,

G1(t,Ω(t)) =
−t2et

3
Ω(t), G2(t) =

t1−µ

Γ(2− µ)
− t2

2
,

G3(Ω(s)) = Ω(s), G4(Ω(s)) = Ω(s),

Π1(t, s) = set, Π2(t, s) = t2, ω0 = 0.

Direct calculations implies that

P1(Ω) = Ω
T
A1, P2(Ω) = Ω

T
A2, P3(Ω) = Ω

T
A3,

where

A1 =
r

2


f(0) f(1) . . . f(M − 1)
3f(0) 3f(1) . . . 3f(M − 1)

...
... . . .

...
(2M + 1)f(0) (2M + 1)f(1) . . . (2M + 1)f(M − 1)



A2 =
r2

3


1 23 − 13 . . . M3 − (M − 1)3

1 23 − 13 . . . M3 − (M − 1)3

...
... . . .

...
1 23 − 13 . . . M3 − (M − 1)3

 , A3 =
1

r


g(0) 0 . . . 0
0 g(1) . . . 0
...

... . . .
...

0 0 . . . g(M − 1)

 ,



M.I. Syam, M. Sharadga, I. Hashim / Eur. J. Pure Appl. Math, 17 (3) (2024), 1429-1448 1443

P4 =
(
h(0)− h(1) h(1)− h(2) . . . h(M − 1)− h(M)

)
f(j) = er(j+1) − erj , g(j) =

−1

3
(er(j+1)(2− 2(j + 1)r + r2(j + 1)2)− erj(2− 2rj + r2j2)),

h(j) =
1

r

(
(rj)2−µ

Γ(3− µ)
− (rj)3

3

)
.

Then,

Ω
T
((A1 +A2 +A3)Oµ − IM ) = P4Oµ

which can be written as

((A1 +A2 +A3)Oµ − IM )T Ω = (P4Oµ)
T . (67)

One can see that System (67) is linear system. For M = 500, We solve the System (67)
in terms of r and we found that

ωj =
2j + 1

2
r, j = 0, 1, . . . ,M − 1.

Thus, the approximate solution is given as

r
M−1∑
j=0

2j + 1

2
µj(t) → x

as M approaches to ∞. Thus, our approximate solution converges to the exact solution.
This result can be obtained using Mathematical for linear cases. Let us study a nonlinear
problem in the next example.

Example 2. Consider the following class of integro-differential equations:

DµΩ(t) = Ω2(t) +G2(t) +

∫ t

0
(sµ + 1)(tµ + 1)Ω2(s)ds+

∫ 1

0
(sµ + 1)tµΩ2(s)ds,

Ω(0) = 1,

where t ∈ [0, 1], 0 < µ ≤ 1, and Then,

G2(t) = −(5µ+ 4)t3µ+2

3µ2 + 5µ+ 2
− 3(3µ+ 2)t5µ+3

20µ2 + 27µ+ 9
− (5µ+ 4)(µ(2µ+ 3)(6µ+ 23) + 21)tµ

(µ+ 1)(3µ+ 2)(4µ+ 3)(5µ+ 3)

− (µ+ 2)tµ+1

µ+ 1
− (2µ+ t+ 3)t2µ+1

µ+ 1
− (3µ+ 4)t4µ+2

3µ+ 2
− t4µ+3

4µ+ 3
− t6µ+3

5µ+ 3

+
Γ(2µ+ 2)tµ+1

Γ(µ+ 2)
− t− 1.

First, let’s examine the influence of the fractional order derivative on the behavior of
the solution. Figures 1 shows the approximate solutions Ω for different values of µ =



M.I. Syam, M. Sharadga, I. Hashim / Eur. J. Pure Appl. Math, 17 (3) (2024), 1429-1448 1444

0.2 0.4 0.6 0.8 1.0

1.2

1.4

1.6

1.8

2.0

μ=0.5

μ=0.6

μ=0.7

μ=0.8

μ=0.9

μ=1

Figure 1: The approximate solution Ω for µ = 0.5, 0.6, 0.7, 0.8, 0.9, 1.

0.5, 0.6, 0.7, 0.8, 0.9, 1 with r = 0.005. Since the exact solution is not given, we assess the
accuracy of our approximation by computing the L2-truncation errors, defined as

ϵ(µ) =

(∫ 1

0

(
DµΩ(t)− Ω2(t)−G2(t) +

∫ t

0
(sµ + 1)(tµ + 1)Ω2(s)ds−

∫ 1

0
(sµ + 1)tµΩ2(s)ds

)
dt

) 1
2

.

The errors are presented in Table 1.

Table 1: The L2-error for µ = 0.5, 0.6, 0.7, 0.8, 0.9, 1..

µ ϵ(µ)

0.5 1.99× 10−12

0.6 1.98× 10−12

0.7 1.74× 10−12

0.8 1.66× 10−12

9 1.23× 10−12

1 1.10× 10−12

Example 3. Consider the following class of integro-differential equations:

DµΩ(t) = Ω3(t) +G2(t) +

∫ t

0
(s+ 1)(t+ 1)Ω3(s)ds+

∫ 1

0
(sµ + 1)t3Ω2(s)ds,

Ω(0) = 0,

where t ∈ [0, 1], 0 < µ ≤ 1, and Then,

G2(t) = Γ(µ+ 1)−
(
9µ2 + 9µ+ 3µt(t+ 1)2 + t(t+ 2)(t+ 1) + 2

)
t3µ

9µ(µ+ 1) + 2

−
3
(
16µ2 + 20µ+ 4µt(t+ 1)2 + t(2t+ 3)(t+ 1) + 6

)
t4µ+1

2(2µ+ 1)(4µ+ 3)

−
3
(
25µ2 + 35µ+ 5µt(t+ 1)2 + t(3t+ 4)(t+ 1) + 12

)
t5µ+2

(5µ+ 3)(5µ+ 4)



M.I. Syam, M. Sharadga, I. Hashim / Eur. J. Pure Appl. Math, 17 (3) (2024), 1429-1448 1445

−
(
36µ2 + 54µ+ 6µt(t+ 1)2 + t(4t+ 5)(t+ 1) + 20

)
t6µ+3

2(3µ+ 2)(6µ+ 5)

− (7µ+ 4)(µ(7µ(18µ+ 31) + 120) + 21)t3

(2µ+ 1)(3µ+ 1)(3µ+ 2)(4µ+ 3)(5µ+ 3)
+

Γ(2µ+ 2)tµ+1

Γ(µ+ 2)
.

First, let’s examine the influence of the fractional order derivative on the behavior of
the solution. Figures 2 shows the approximate solutions Ω for different values of µ =
0.5, 0.6, 0.7, 0.8, 0.9, 1 with r = 0.005. Since the exact solution is not given, we assess the

0.2 0.4 0.6 0.8 1.0

0.5

1.0

1.5

2.0

μ=0.5

μ=0.6

μ=0.7

μ=0.8

μ=0.9

μ=1

Figure 2: The approximate solution Ω for µ = 0.5, 0.6, 0.7, 0.8, 0.9, 1.

accuracy of our approximation by computing the L2-truncation errors, defined as

ϵ(µ) =

(∫ 1

0

(
DµΩ(t)− Ω2(t)−G2(t) +

∫ t

0
(sµ + 1)(tµ + 1)Ω2(s)ds−

∫ 1

0
(sµ + 1)tµΩ2(s)ds

)
dt

) 1
2

.

The errors are presented in Table 2.

Table 2: The L2-error for µ = 0.5, 0.6, 0.7, 0.8, 0.9, 1.

µ ϵ(µ)

0.5 2.34× 10−12

0.6 2.21× 10−12

0.7 2.01× 10−12

0.8 1.98× 10−12

9 1.95× 10−12

1 1.82× 10−12

6. Conclusion

The operational matrix method is a useful approach for solving fractional Volterra-
Fredholm integro-differential equations. It involves converting the differential system into
a system of algebraic equations to determine the coefficients of the approximating solution.



REFERENCES 1446

Typically, these coefficients are computed by establishing operational matrices correspond-
ing to integral, derivative, and product operators. The Block Pulse functions, upon which
our approximation basis, possess three key properties—disjointness, orthogonality, and
completeness—that facilitate the computations. In this paper, we investigate fractional
Volterra-Fredholm integro-differential equations, a fundamental problem in various fields
such as control theory, biology, and particle dynamics in physics. We develop a numerical
method based on the operational matrix method to solve these equations, proving the
existence and uniqueness of the exact solution. Additionally, we demonstrate the uniform
convergence of numerical solutions to the exact solution and present several numerical
examples illustrating the method’s efficiency. In one example, for a linear problem, we
observe convergence of the approximate solution to the exact solution as the number of
Block Pulse functions increases. In nonlinear cases, where exact solutions are unavailable,
we compute the L2-local truncation error, which is on the order of 10−12. We also exam-
ine the influence of the fractional derivative on solution profiles through graph sketches.
Theoretical and numerical results affirm the accuracy and applicability of our proposed
method to nonlinear problems in science.

Concluding this paper, we highlight the following observations:

(i) Example 1 demonstrates that in linear cases, the approximate solution converges
to the exact solution with increasing numbers of Block Pulse functions, yielding
solutions in closed form.

(ii) Example 2 illustrates the decreasing influence of the fractional order on solution
profiles as the fractional derivative increases. The L2-truncation error is of order
10−12, as shown in Table 1 and Figure 1.

(iii) Similar results are obtained in Example 3, as presented in Table 2 and Figure 2.

(iv) Our findings suggest that our proposed method is promising and applicable across
diverse models in physics and engineering, even in the presence of significant non-
linearity.

Acknowledgements

The authors express their gratitude to the editor and the reviewers for their valuable
comments.

References

[1] D Baleanu, K Diethelm, E Scalas, and J.J. Trujillo, editors. Fractional Calculus:
Models and Numerical Methods., volume 3. World Scientific, 2012.

[2] A H Bhrawy and A S Alofi. Operational matrix method for solving fractional dif-
ferential equations: A review. Journal of Computational and Applied Mathematics,
330:1056–1066, 2017.



REFERENCES 1447

[3] A Ghorbani and S A Yousefi. Numerical solution of nonlinear fractional differen-
tial equations by operational matrix method based on difference basis. Journal of
Computational and Applied Mathematics, 347:603–618, 2019.

[4] A Ghorbani and S A Yousefi. Solving fractional differential equations by operational
matrix method based on difference basis. Communications in Nonlinear Science and
Numerical Simulation, 67:416–433, 2019.

[5] M Khashan and M I Syam. An efficient method for solving fractional ricatti equations.
Advances in Difference Equations, 2019(1):1–12, 2019.

[6] A Kilbas, H M Srivastava, and J J Trujillo. Theory and Application of Fractional
Differential Equations., volume 204 of North Holland Mathematics Studies. Elsevier,
Amsterdam, 2006.

[7] P Kumar and D Kumar. Numerical solution of fractional differential equations using
difference operational matrix. Results in Physics, 15:102616, 2019.

[8] S Momani and Z Odibat. Analytical solution of a time-fractional navier–stokes equa-
tion by adomian decomposition method. Applied Mathematics and Computation,
190(1):102–110, 2008.

[9] A-H I Mourad, A M Ghazal, M M Syam, O D Al Qadi, and H Al Jassmi. Utilization of
additive manufacturing in evaluating the performance of internally defected materials.
In IOP Conf. Ser.: Mater. Sci. Eng., volume 362, page 012026, 2018.

[10] S Al Omari, A M Ghazal, M Syam, R Al Najjar, and M Y Selim. An investigation
on the thermal degradation performance of crude glycerol and date seeds blends
using thermogravimetric analysis (tga). In 5th International Conference on Renewable
Energy: Generation and Application, ICREGA 2018, pages 102–106, 2018.

[11] I Podlubny. Fractional Differential Equations. Academic Press, San Diego, 1999.

[12] S G Samko, A A Kilbas, and O I Marichev. Fractional Integrals and Derivatives:
Theory and Applications. Gordon and Breach, Yverdon, 1993.

[13] M I Syam, MN Y Anwar, A Yildirim, and et al. The modified fractional power series
method for solving fractional non-isothermal reaction–diffusion model equations in a
spherical catalyst. International Journal of Applied and Computational Mathematics,
5:38, 2019.

[14] M I Syam, M A Raja, M M Syam, and et al. An accurate method for solving the
undamped duffing equation with cubic nonlinearity. International Journal of Applied
and Computational Mathematics, 4:69, 2018.

[15] M I Syam, M Sharadga, and I Hashim. A numerical method for solving fractional
delay differential equations based on the operational matrix method. Chaos, Solitons
& Fractals, 147:110977, 2021.



REFERENCES 1448

[16] M M Syam, S Cabrera-Calderon, K A Vijayan, V Balaji, P E Phelan, and J R
Villalobos. Mini containers to improve the cold chain energy efficiency and carbon
footprint. Climate, 10:76, 2022.

[17] S M Syam, Z Siri, S H Altoum, M A Aigo, and R Md Kasmani. A novel study for
solving systems of nonlinear fractional integral equations. Applied Mathematics in
Science and Engineering, 31(1), 2023.

[18] S M Syam, Z Siri, S H Altoum, and R M Kasmani. An efficient numerical approach for
solving systems of fractional problems and their applications in science. Mathematics,
11:3132, 2023.

[19] S M Syam, Z Siri, S H Altoum, and R Md Kasmani. Analytical and numerical
methods for solving second-order two-dimensional symmetric sequential fractional
integro-differential equations. Symmetry, 15:1263, 2023.

[20] S M Syam, Z Siri, and R M Kasmani. Operational matrix method for solving frac-
tional system of riccati equations. In 2023 International Conference on Fractional
Differentiation and Its Applications (ICFDA), pages 1–6, Ajman, United Arab Emi-
rates, 2023.

[21] S M Syam, Z Siri, R Md Kasmani, and Kenan Yildirim. A new method for solving se-
quential fractional wave equations. Journal of Mathematics, 2023:Article ID 5888010,
16 pages, 2023.