EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS
Vol. 15, No. 2, 2022, 796-809
ISSN 1307-5543 – ejpam.com
Published by New York Business Global

Toeplitz matrix and Nyström method for solving linear
fractional integro-differential equation

Sammeha Raad1,∗, Khawlah Alqurashi1

1 Mathematical Sciences, Applied Sciences, Umm Al-Qura University, Makkah, Saudi Arabia

Abstract. In this paper, the Volterra-Fredholm integral equation is derived from a linear integro-
differential equation with a fractional order 0 < α < 1 using Riemann–Liouville fractional integral.
The existence and uniqueness of the solution are proved using the Picard method. Popular nu-
merical methods; the Toeplitz matrix, and the product Nyström are used in the solution. These
methods will prove their effective in solving this type of equation. Two examples are solved using
the mentioned methods and the estimation error is calculated. Finally, a comparison between the
numerical results is made.

2020 Mathematics Subject Classifications: 45F15, 45K05, 45E10, 74H15

Key Words and Phrases: Systems of linear Abel equations, Integro-partial differential equa-
tions, Picard method; Toeplitz matrix, Nyström method.

1. Introduction

The wide range of applications for fractional equations has led to increased interest in
recent years. The list of applications has expanded and become more diverse in a short
period of time. One example of such uses is electromagnetic fields. [17] Constructed frac-
tional integro-differential equations (FI-DEs) from electromagnetic waves in a dielectric
material. The existence, uniqueness and the convergence of the solution of FI-DEs by
using the Picard method were discussed in [9]. The existence and uniqueness of the so-
lution of the linear fractional Volterra I-DEs with initial conditions relied on applying
the Picard iteration method to obtain uniformly convergent series for the exact solution
were discussed by [13]. The existence and uniqueness of mild solutions for FI-DEs were
investigated by using Holder’s inequality, p-mean continuity, and Schauder’s fixed point
theorem in Banach spaces in [3]. Several authors are interested in solving FI-DEs by
analytical methods. Furthermore, numerical methods are used to approximate the so-
lutions. In [12] introduced the generalized hat functions and operational matrix of the
fractional integration to solve the FI-DEs of Bratu-type numerically. The decomposition

∗Corresponding author.
DOI: https://doi.org/10.29020/nybg.ejpam.v15i2.4384

Email addresses: saraad@uqu.edu.sa (S. Raad), s44180462@st.uqu.edu.sa (K. Alqurashi)

https://www.ejpam.com 796 © 2022 EJPAM All rights reserved.



S. Raad, K. AlQurashi / Eur. J. Pure Appl. Math, 15 (2) (2022), 796-809 797

method for approximating the solution of FI-DEs systems was implemented by [11]. In
[15] applied the Legendre wavelets method to approximate the solution of FI-DEs. In
[5], Sinc-Collocation method was introduced to solve Volterra-Fredholm FI-DEs. Nu-
merical solutions of linear Fredholm-Volterra FI-DEs using Laguerre polynomials were
investigated by [6]. Chebyshev polynomials were used to solve the Fredholm and Volterra
integro-differential equations in [7]. Using coupled mathematical solutions, the generalized
monotone iterative technique was developed to solve the Caputo FI-DEs of order q[8].

2. The linear fractional integro-differential equation

Consider the linear fractional integro-differential equation (LFI-DE)

∂αu(x, t)

∂tα
= f(x, t) + λ

∫ b

a
k(x, y)u(y, t)dy , (0 < α < 1) (1)

with the initial condition
u(x, 0) = u0(x).

In (1) the unknown function appears on one side of the equation under the fractional order
derivative and appears on the other side under integration.

Definition 1. [16] For all t ∈ [a, b] the left Riemann-Liouville fractional integral of order
α > 0, of the function ϕ : (0,∞) → R is defined as

aI
α
t ϕ(t) =

1

Γ(α)

∫ t

a
(t− τ)α−1ϕ(τ)dτ, t > a. (2)

Applying relation (2) to equation (1), we obtain

u(x, t) = η(x, t) +
λ

Γ(α)

∫ t

0

∫ b

a
(t− s)α−1k(x, y)u(y, s)dyds, (3)

where

η(x, t) = u0(x) +
1

Γ(α)

∫ t

0
(t− s)α−1f(x, s)ds, (4)

equation (3) is a singular Volterra-Fredholm integral equation (V-FIE) of Abel’s type. It
is clear that equation (3) is equivalent to equation (1).

3. Picard method for singularV-FIE [1]

In this section, the Picard method is used to prove the existence of a unique solution
to singular V-FIE (3) where

F (t, s) = (t− s)α−1.

For this goal, consider the following conditions:



S. Raad, K. AlQurashi / Eur. J. Pure Appl. Math, 15 (2) (2022), 796-809 798

(i) Volterra integral’s kernel F (t, s) satisfies the discontinuity condition[ ∫ T

0

∫ T

0
|F (t, s)|2dsdt

] 1
2
= Q, (Q is a constant).

(ii) The kernel of Fredholm k(x, y) belongs to the class C[a, b], and it is bounded, i.e.

|k(x, y)| ≤ B, ∀x, y ∈ [a, b], (B is a constant).

(iii) The function f(x, t), and its partial derivatives with respect to x and t, are continuous
in the Banach space L2[0, T ]× C[a, b], T < 1. Its norm is defined as

||f(x, t)|| = max
a≤x≤b

∣∣∣ ∫ x

a

[ ∫ t

0
|f(x, s)|2ds

] 1
2
dy

∣∣∣ = H, (H is a constant).

(iv) The function u0(x) belongs to the space C[a, b] and has the norm

||u0(x)|| = max
a≤x≤b

|u0(x)| = A, (A is a constant).

Theorem 1. Equation (3) has a unique solution in the Banach space L2[0, T ] × C[a, b],
under the condition

λ(b− a)BQ < Γ(α).

Proof.
The existence of a unique solution of equation (3) can be proved by using the method of

successive approximations which is also called ”Picard’s method ”. This method consists
of the following simple iteration.

un(x, t) = η(x, t) +
λ

Γ(α)

∫ t

0

∫ b

a
F (t, s)k(x, y)un−1(y, s)dyds, (n ≥ 1), (5)

with
u0(x, t) = η(x, t).

For ease manipulation it is convenient to introduce

ξn = un(x, t)− un−1(x, t) (6)

where

un(x, t) =

n∑
i=0

ξi(x, t), ξ0(x, t) = η(x, t), (7)

subtracting from equation (5) a similar equation with n replaced by n− 1, we get

un(x, t)− un−1(x, t) =
λ

Γ(α)

∫ t

0

∫ b

a
F (t, s)k(x, y)[un−1(y, s)− un−2(y, s)]dyds, (8)



S. Raad, K. AlQurashi / Eur. J. Pure Appl. Math, 15 (2) (2022), 796-809 799

inserting (6) in (8), we have

|ξn(x, t)| =
|λ|
Γ(α)

∫ t

0

∫ b

a
|F (t, s)||k(x, y)||ξn−1(y, s)|dyds, (9)

using condition (ii), then applying Cauchy-Schwarz inequality, we obtain

||ξn(x, t)|| ≤
|λ|B
Γ(α)

(∫ T

0

∫ T

0
|F (t, s)|2dsdt

) 1
2
max
a≤x≤b

∣∣∣ ∫ x

a

(∫ t

0
|ξn−1(y, s)|2ds

) 1
2
dy

∣∣∣ ∫ b

a
dy,

in view of condition (i), we get

||ξn(x, t)|| ≤
|λ|BQ(b− a)

Γ(α)
||ξn−1(x, t)||. (10)

Inequality (10) for n = 1, yields

||ξ1(x, t)|| ≤ σ||ξ0(x, t)|| ;
(
σ =

|λ|BQ(b− a)

Γ(α)

)
, (11)

from (4) and (7), we have

|ξ0(x, t)| ≤ |u0(x)|+
1

Γ(α)

∫ t

0
|F (t, s)||f(x, s)|ds,

after applying Cauchy-Schwarz inequality, we get

||ξ0(x, t)|| ≤ max
a≤x≤b

|u0(x)|+
1

Γ(α)

(∫ T

0

∫ T

0
|F (t, s)|2dsdt

) 1
2
max
a≤x≤b

∣∣∣ ∫ x

a

(∫ t

0
|f(x, s)|2ds

) 1
2
dy

∣∣∣,
using conditions (i), (iii) and (iv), the above inequality takes the form

||ξ0(x, t)|| ≤ G ;
(
G = A+

QH

Γ(α)
is a constant

)
. (12)

Introducing (11) in (12), we have

||ξ1(x, t)|| ≤ σG,

by induction, we can prove that

||ξn(x, t)|| ≤ σnG ; n = 0, 1, 2, .... (13)

Since (13) is obviously true for n = 0, 1, then it holds for all n. This bound makes
the sequence {ξn(x, t)} converges under the condition σ < 1, and therefore the sequence
{un(x, t)} in (7) converges. Hence we can write

u(x, t) =
∞∑
i=0

ξi(x, t), (14)



S. Raad, K. AlQurashi / Eur. J. Pure Appl. Math, 15 (2) (2022), 796-809 800

the series (14) is uniformly convergent since the terms ξi(x, t) are limited by σi.
To prove that u(x, t) defined by (14) satisfies equation (3), set

u(x, t) = un(x, t) + ∆n(x, t), (∆n(x, t) → 0 as n→ ∞), (15)

from equation (3), we get

u(x, t)−∆n(x, t) = η(x, t) +
λ

Γ(α)

∫ t

0

∫ b

a
F (t, s)k(x, y)[u(y, s)−∆n−1(y, s)]dyds,

therefore, we have∣∣∣u(x, t)− η(x, t)− λ

Γ(α)

∫ t

0

∫ b

a
F (t, s)k(x, y)u(y, s)dyds

∣∣∣
≤

∣∣∆n(y, s)
∣∣+ |λ|

Γ(α)

∫ t

0

∫ b

a
|F (t, s)||k(x, y)||∆n−1(y, s)|dyds,

using condition (ii), then applying Cauchy-Schwarz inequality to the integral term in the
right-hand side and in view of condition (i), we get

||u(x, t)− η(x, t)− λ

Γ(α)

∫ t

0

∫ b

a
F (t, s)k(x, y)u(y, s)dyds|| ≤ ||∆n(x, t)||+ σ||∆n−1(x, t)||,

(16)
by taking n large enough, the right-hand side of (16) can be made as small as desired.
Consequently, the function u(x, t) defined by (15) satisfies

u(x, t) = η(x, t) +
λ

Γ(α)

∫ t

0

∫ b

a
F (t, s)k(x, y)u(y, s)dyds,

and is therefore the solution to equation (3).
To show that u(x, t) is the only solution of equation (3), we assume the existence of another
solution ũ(x, t), then

|u(x, t)− ũ(x, t)| ≤ |λ|
Γ(α)

∫ t

0

∫ b

a
|F (t, s)||k(x, y)||u(y, s)− ũ(y, s)|dyds,

using condition (ii), then applying Cauchy-Schwarz inequality and finally in view of con-
dition (i), we deduce that

||u(x, t)− ũ(x, t)|| ≤ σ||u(x, t)− ũ(x, t)||, (17)

because of σ < 1, this can be true if u(x, t) = ũ(x, t); that is, the solution of equation (3)
is unique.



S. Raad, K. AlQurashi / Eur. J. Pure Appl. Math, 15 (2) (2022), 796-809 801

4. A system of Volterra integral equations [1, 2]

In this section, the V-FIE will be converted to a system of Volterra integral equation
by dividing the interval [a, b], into M subintervals, such that a = x0 < x1 < x2 < ...... <
xm < ...... < xM = b where x = xm,m = 0, 1, 2...,M. Equation (3) becomes

um(t) = ηm(t) +
λ

Γ(α)

∫ t

0
(t− s)α−1

∫ b

a
k(xm, y)u(y, s)dyds (18)

The Fredholm integral part of equation (18), after using the quadrature formula, takes
the form ∫ b

a
k(xm, y)u(y, s)dy =

M∑
l=0

wlk(xm, xl)u(xl, s), l = 0, 1, .....M (19)

Using (19) in (18), we have the following Volterra integral equations system

um(t) = ηm(t) +
λ

Γ(α)

M∑
l=0

wlkm,l

∫ t

0
(t− s)α−1ul(s)ds. (20)

where um(t) = u(xm, t) , ηm(t) = η(xm, t) , km,l = k(xm, xl), w0 = wm = 1
2h0, wr =

hr, (r ̸= 0,m).

5. Numerical methods to solve singular Volterra integral equation

5.1. The Toeplitz matrix method [2, 4]

Consider the Volterra integral equation:

u∗(t)− λ

∫ t

0
(t− s)(α−1)u∗(s)ds = f∗(t), (21)

write the integral term in the form∫ t

0
(t− s)(α−1)u∗(s)ds =

N−1∑
n=0

∫ nh+h

nh
(t− s)(α−1)u∗(s)ds, (h =

T

N
), (22)

approximate the integral in the right hand side of equation (22) by∫ nh+h

nh
(t− s)(α−1)u∗(s)ds = An(t)u

∗(nh) +Bn(t)u
∗(nh+ h), (23)

where

An(t) =
1

h
[(nh+ h)I(t)− J(t)], Bn(t) =

1

h
[J(t)− nhI(t)], (24)



S. Raad, K. AlQurashi / Eur. J. Pure Appl. Math, 15 (2) (2022), 796-809 802

I(x) =

∫ nh+h

nh
(t− s)(α−1)ds, J(x) =

∫ nh+h

nh
s(t− s)(α−1)ds. (25)

If we set t = ph, the integral equation (21) becomes

u∗(ph)− λDn,pu
∗(nh) = f∗(ph), (26)

or by the matrix expression
[I − λD]U∗ = F ∗ (27)

The function F ∗ is a vector of N +1 elements, but D is a matrix whose elements are given
by

Dn,p=an,p + a′n,p

The matrix an,p = An(ph)+Bn−1(ph) is the Toeplitz matrix of order N+1, 0 ≤ n, p ≤ N .
Whereas, the elements of the second matrix a′n,p are zeros except for a0,p = B1(ph) and
ap−1,p = Bp(ph). Therefore, the solution of the formula (27) will take the form

U∗ = [I − λD]−1F ∗, |I − λD| ≠ 0. (I is the identity matrix.) (28)

5.2. The product Nyström method [1, 10, 14]

Assume the Volterra integral equation can be written in the form

u∗(t)− λ

∫ t

0
p(x, y)k(x, y)u∗(s)ds = f∗(t), (29)

where p(x, y) = (t − s)α−1 and k(x, y) are ‘badly behaved‘ and ‘well behaved‘ of their
arguments, respectively. Equation (29) can be written in the form

u∗(ti)− λ
i∑

j=0

wij(ti − s)α−1u∗(si) = f∗(ti), (30)

here xi = yi = a + ih, i = 0, 1, .., N with h =
T

N
, N must be even, and wij are the

weights to be determined. If we approximate the nonsingular part by the second degree
of Lagrange interpolation polynomial, therefore,

ωi,0 = β1(yi), ωi,2j+1 = 2γj+1(yi),

ωi,2j = ζj(yi) + βj+1(yi), ωi,N (yi) = ζN
2
(yi)

(31)

Assume v = i− 2j + 2, then

ζj(yi) =
h

2

∫ 2

0
ρ(ρ− 1)(vh− ρh)(α−1)dρ,

βj(yi) =
h

2

∫ 2

0
(ρ− 1)(ρ− 2)(vh− ρh)(α−1)dρ,

γj(yi) =
h

2

∫ 2

0
ρ(2− ρ)(vh− ρh)(α−1)dρ.

(32)



S. Raad, K. AlQurashi / Eur. J. Pure Appl. Math, 15 (2) (2022), 796-809 803

Now,

ψi =

∫ 2

0
ρi(vh− ρh)(α−1)dρ,

equation (31) becomes

ωi,0 =
h

2
[2ψ0(v)− 3ψ1(v) + ψ2(v)], v = i

ωi,2j =
h

2
[ψ2(v)− ψ1(v) + 2ψ0(v − 2)− 3ψ1(v − 2) + ψ2(v − 2)], v = i− 2j + 2

ωi,2j+1 = h[2ψ1(v)− ψ2(v)], v = i− 2j

ωi,N =
h

2
[ψ2(v)− ψ1(v)], v = i−N + 2

(33)

Therefore, the integral equation (29) is reduced to a system of linear algebraic equations

(I − λW )U∗ = F , F = f∗i = f∗(ti) ; i = 0, 1, 2, ....., N (34)

which has the solution

U∗ = [I − λW ]−1F , |I − λW | ≠ 0 (I is the identity matrix.) (35)

6. Numerical results

In this section, The fractional integro-differential equation will be solved using the
Toeplitz matrix and the product Nyström methods. The results are calculated at T =
0.009, 0.05, 0.2, and α = 0.6, 0.75, 0.9. Maple 18 software will be used in programming.

Example 1. Consider the FI-DE

∂αu(x, t)

∂tα
= x

t1−α

Γ(2− α)
− 2t2−α

Γ(3− α)
−λx(t−t2e+t2)+λ

∫ 1

0
xeyu(y, t)dy, u0(x) = 0, (36)

after using the Riemann-Liouville fractional integral to the equation (36), we obtain

u(x, t) = η(x, t) +
λ

Γ(α)

∫ t

0

∫ 1

0
(t− s)α−1xeyu(y, s)dyds, (37)

where

η(x, t) = xt− t2 − λx

[
t1+α

Γ(2 + α)
+

2t2+α

Γ(3 + α)
(1− e)

]
. (38)

Applying the Toeplitz matrix and the product Nyström methods to equation (37), the nu-
merical results are shown in the following tables.



S. Raad, K. AlQurashi / Eur. J. Pure Appl. Math, 15 (2) (2022), 796-809 804

Table 1: The numerical solution at α = 0.6 and T = 0.009.

t Exact solution Toeplitz solution Error Nyström solution Error

0 0 0 0 0 0
0.0018 1.7968× 10−3 1.7969× 10−3 1.0585× 10−7 1.7961× 10−3 1.9918× 10−7

0.0036 3.5870× 10−3 3.5880× 10−3 9.6681× 10−7 3.5882× 10−3 1.1977× 10−6

0.0054 5.3708× 10−3 5.3731× 10−3 3.1263× 10−6 5.3743× 10−3 3.4511× 10−6

0.0072 7.1482× 10−3 7.1553× 10−3 7.1190× 10−6 7.1555× 10−3 7.3083× 10−6

0.0090 8.9190× 10−3 8.9346× 10−3 1.5609× 10−5 8.9321× 10−3 1.3067× 10−5

Table 2: The numerical solution at α = 0.9 and T = 0.009.

t Exact solution Toeplitz solution Error Nyström solution Error

0 0 0 0 0 0
0.0018 1.7968× 10−3 1.7968× 10−3 2.5004× 10−8 1.7968× 10−3 2.8887× 10−8

0.0036 3.5870× 10−3 3.5872× 10−3 1.8726× 10−7 3.5872× 10−3 1.9739× 10−7

0.0054 5.3708× 10−3 5.3715× 10−3 6.1311× 10−7 5.3715× 10−3 6.2659× 10−7

0.0072 7.1482× 10−3 7.1496× 10−3 1.4257× 10−6 7.1496× 10−3 1.4294× 10−6

0.0090 8.9190× 10−3 8.9218× 10−3 2.7961× 10−6 8.9217× 10−3 2.7133× 10−6

Table 3: The numerical solution at α = 0.6 and T = 0.2.

t Exact solution Toeplitz solution Error Nyström solution Error

0 0 0 0 0 0
0.04 3.84× 10−2 3.8413× 10−2 1.3325× 10−5 3.8426× 10−2 2.6277× 10−5

0.08 7.36× 10−2 7.3717× 10−2 1.1738× 10−4 7.3748× 10−2 1.4837× 10−4

0.12 1.056× 10−1 1.0596× 10−1 3.5658× 10−4 1.0591× 10−1 3.9846× 10−4

0.16 1.344× 10−1 1.3516× 10−1 7.5928× 10−4 1.3518× 10−1 7.8131× 10−4

0.20 1.60× 10−1 1.6162× 10−1 1.6062× 10−3 1.6128× 10−1 1.2842× 10−3

Table 4: The numerical solution at α = 0.9 and T = 0.2.

t Exact solution Toeplitz solution Error Nyström solution Error

0 0 0 0 0 0
0.04 3.84× 10−2 3.8408× 10−2 8.2129× 10−6 3.8401× 10−2 9.5921× 10−6

0.08 7.36× 10−2 7.3658× 10−2 5.7747× 10−5 7.3661× 10−2 6.1222× 10−5

0.12 1.056× 10−1 1.0578× 10−1 1.7594× 10−4 1.0578× 10−1 1.8036× 10−4

0.16 1.344× 10−1 1.3478× 10−1 3.7816× 10−4 1.3478× 10−1 3.7916× 10−4

0.20 1.60× 10−1 1.6068× 10−1 6.8461× 10−4 1.6066× 10−1 6.5791× 10−4



S. Raad, K. AlQurashi / Eur. J. Pure Appl. Math, 15 (2) (2022), 796-809 805

Figure 1: The error by Toeplitz matrix and Nyström product methods at various values for α = 0.6, 0.75, 0.9
and at T = 0.009, 0.05, 0.2.

Tables and the figures show that the maximum error by Toeplitz matrix is 1.6062×10−3

and by Nyström method is 1.2842× 10−3 both at α = 0.6 and T = 0.2

Example 2. Consider the FI-DE

∂αu(x, t)

∂tα
= x2

2t2−α

Γ(3− α)
− λt2(

x

4
+
x2

5
) + λ

∫ 1

0
(xy + x2y2)u(y, t)dy, u0(x) = 0 (39)

Applying the Riemann-Liouville fractional integral to equation (39), we get

u(x, t) = η(x, t) +
λ

Γ(α)

∫ t

0

∫ 1

0
(t− s)α−1(xy + x2y2)u(y, s)dyds, (40)

and

η(x, t) = x2t2 − λx2
2t2+α

5Γ(3 + α)
− λx

2t2+α

4Γ(3 + α)
. (41)

The numerical results after using the Toeplitz matrix and the product Nyström methods to
equation (40) are demonstrated in the following tables.

Table 5: The numerical results at α = 0.6 and T = 0.009

t Exact solution Toeplitz solution Error Nyström solution Error

0 0 0 0 0 0
0.0018 3.24× 10−6 3.2400× 10−6 3.0434× 10−11 3.2401× 10−6 1.2663× 10−10

0.0036 1.296× 10−5 1.2961× 10−5 7.1041× 10−10 1.2961× 10−5 1.2172× 10−9

0.0054 2.916× 10−5 2.9164× 10−5 3.7224× 10−9 2.915× 10−5 4.8778× 10−9

0.0072 5.184× 10−5 5.1852× 10−5 1.1978× 10−8 5.1853× 10−5 1.3282× 10−8

0.0090 8.1× 10−5 8.1039× 10−5 3.9302× 10−8 8.1029× 10−5 2.9079× 10−8



S. Raad, K. AlQurashi / Eur. J. Pure Appl. Math, 15 (2) (2022), 796-809 806

Table 6: The numerical results at α = 0.9 and T = 0.009

t Exact solution Toeplitz solution Error Nyström solution Error

0 0 0 0 0 0
0.0018 3.24× 10−6 3.2400× 10−6 1.5934× 10−11 3.2400× 10−6 1.8594× 10−11

0.0036 1.296× 10−5 1.2962× 10−5 1.8965× 10−10 1.2960× 10−5 2.0850× 10−10

0.0054 2.916× 10−5 2.9161× 10−5 8.7897× 10−10 2.9161× 10−5 9.2193× 10−10

0.0072 5.184× 10−5 5.1843× 10−5 2.6682× 10−9 5.1843× 10−5 2.7020× 10−9

0.0090 8.1× 10−5 8.1007× 10−5 6.6027× 10−9 8.1006× 10−5 6.2751× 10−9

Table 7: The numerical results at α = 0.6 and T = 0.2

t Exact solution Toeplitz solution Error Nyström solution Error

0 0 0 0 0 0
0.04 1.6× 10−3 1.6001× 10−3 9.6612× 10−8 1.6004× 10−3 4.0211× 10−7

0.08 6.4× 10−3 6.4023× 10−3 2.2562× 10−6 6.4039× 10−3 3.8664× 10−6

0.12 1.44× 10−2 1.4412× 10−2 1.1827× 10−5 1.4416× 10−2 1.5501× 10−5

0.16 2.56× 10−2 2.5638× 10−2 3.8076× 10−5 2.5642× 10−2 4.2227× 10−5

0.20 4× 10−2 4.0125× 10−2 1.2499× 10−4 4.0093× 10−2 9.2503× 10−5

Table 8: The numerical results at α = 0.9 and T = 0.2

t Exact solution Toeplitz solution Error Nyström solution Error

0 0 0 0 0 0
0.04 1.6× 10−3 1.6001× 10−3 1.2824× 10−7 1.6001× 10−3 1.4965× 10−7

0.08 6.4× 10−3 6.4015× 10−3 1.5266× 10−6 6.4017× 10−3 1.6784× 10−6

0.12 1.44× 10−2 1.4407× 10−2 7.0765× 10−6 1.4407× 10−2 7.4226× 10−6

0.16 2.56× 10−2 2.5621× 10−2 2.1487× 10−5 2.5622× 10−2 2.1759× 10−5

0.20 4× 10−2 4.0053× 10−2 5.3186× 10−5 4.0051× 10−2 5.0548× 10−5



S. Raad, K. AlQurashi / Eur. J. Pure Appl. Math, 15 (2) (2022), 796-809 807

Figure 2: The error by Toeplitz matrix and Nyström product methods at various values for α = 0.6, 0.75, 0.9
and at T = 0.009, 0.05, 0.2.

It is clear the largest error by Toeplitz matrix and Nyström product methods is at
α = 0.6 and T = 0.2 . The largest error by Toeplitz matrix is 1.2499 × 10−4, where as it
by Nyström method is 9.2503× 10−5

Notes

Through the results given in the previous tables, we note that both the Toeplitz matrix
and the Nyström product methods are effective in solving our problem. The tables and
the figures show that the numerical solutions were obtained using Toeplitz matrix and the
product Nyström are very close to the exact solution. However, the results using Nyström
are better than corvesponding results by Toeplitz matrix . Also, for any choice of the time
value T , there is a negative proportionality between α and the error, i.e. the higher the
value of α, the smaller the error. We also note that with increasing time values, the error
for any values of α increases.

7. Conclusion

We have assumed a LFI-DE with a fractional-order 0 < α < 1. Then, the Riemann-
Liouville fractional integral of order α > 0 was applied to the FI-DE to convert it to a
V-FIE with Abel’s kernel. After that, the uniqueness of the solution of the V-FIE, which
is equivalent to the LFI-DE, has been proved by using the Picard method. Finally, two
numerical methods; the Toeplitz matrix and product Nyström methods were used to find
the numerical solution. The results showed the efficiency and accuracy of the two methods
through the multi-values of α and T .



REFERENCES 808

References

[1] Abdou M. A., El-Kojak, M. K., and S. A. Raad. Analytic and numeric solution of
linear partial differential equation of fractional order. Global J. and Decision science.
Ins.(USA), 13(3/10):57–71., 2013.

[2] Raad S. A., Abdou M. A., and M. M. El-Kojok. On the solution of fredholm-volterra
integral equation with discontinuous kernel in time. JOURNAL OF ADVANCES IN
MATHEMATICS, 9(4):2553–2563, 2014.

[3] Mohamed I Abbas. Existence and uniqueness of mittage-leffler-ulam stable solution
for fractional integrodifferential equations with nonlocal initial conditions. European
Journal Of Pure And Applied Mathematics, 8(4):478–498, 2015.

[4] MA Abdou, MM El-Borai, and MM El-Kojok. Toeplitz matrix method and nonlin-
ear integral equation of hammerstein type. Journal of Computational and Applied
Mathematics, 223(2):765–776, 2009.

[5] Sertan Alkan and Veysel Fuat Hatipoglu. Approximate solutions of volterra-fredholm
integro-differential equations of fractional order. Tbilisi Mathematical Journal,
10(2):1–13, 2017.

[6] Ayşegül Daşcıoğlu and Dilek Varol. Laguerre polynomial solutions of linear fractional
integro-differential equations. Mathematical Sciences, 15(1):47–54, 2021.

[7] Mishra VN Deepmala, HR Marasi, H Shabanian, and M Nosraty. Solution of voltra-
fredholm integro-differential equations using chebyshev collocation method. Global
Journal of Technology and Optimization, 8(210):1–4, 2017.

[8] J Vasundhara Devi and Chaduvula Venkata Sreedhar. Generalized monotone iterative
method for caputo fractional integro-differential equation. European Journal of Pure
and Applied Mathematics, 9(4):346–359, 2016.

[9] AMA El-Sayed, HHG Hashem, and EAA Ziada. Picard and adomian methods for
quadratic integral equation. Computational & Applied Mathematics, 29(3):447–463,
2010.

[10] Han Guoqiang and Zhang Liqing. Asymptotic expansion for the trapezoidal nyström
method of linear volterra—fredholm equations. Journal of computational and applied
mathematics, 51(3):339–348, 1994.

[11] Shaher Momani and Rami Qaralleh. An efficient method for solving systems of frac-
tional integro-differential equations. Computers & Mathematics with Applications,
52(3-4):459–470, 2006.

[12] Li Mulin, Wang Lifeng, and Yu Liu. Generalized hat functions method for solving
fractional integro-differential equations of bratu-type. IAENG International Journal
of Computer Science, 44(1):105–111, 2017.



REFERENCES 809

[13] Hasan N. N. and Hussien D. A. Existence and uniqueness solution of linear fractional
volterra integro-differential equations. Journal of Engineering & Applied Sciences,
14:418–419, 2019.

[14] Sameeha Ali Raad and Mariam Mohammed Al-Atawi. Nyström method to solve two-
dimensional volterra integral equation with discontinuous kernel. Journal of Compu-
tational and Theoretical Nanoscience, 18(4):1177–1184, 2021.

[15] EA Rawashdeh. Legendre wavelets method for fractional integro-differential equa-
tions. Applied Mathematical Sciences, 5(2):2467–2474, 2011.

[16] N Singha and C Nahak. Solutions of the generalized abel’s integral equation us-
ing laguerre orthogonal approximation. Applications and Applied Mathematics: An
International Journal (AAM), 14(2):27, 2019.

[17] Vasily E Tarasov. Fractional integro-differential equations for electromagnetic waves
in dielectric media. Theoretical and Mathematical Physics, 158(3):355–359, 2009.