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

Numerical simulation of initial value problem of
integro-differential equation models

F. M. Alharbi1, S. S. Althubiti2,∗

1 Mathematics Department, Faculty of Sciences, Umm Al-Quraa University, Makkah, Saudi
Arabia

Abstract. In this article, the Volterra-Fredholm integral equation is derived from an initial
value problem of kind integro-differential equation. We discuss the existence and uniqueness of
the solution to the problem in Hilbert space. A numerical method is used to reduce this type
of equation to the system of Fredholm integral equations of the second kind. In light of this,
the Collocation method and the Galerkin method are used to solve the system of second-order
Fredholm integral equations and calculate the error in each case. Finally, the approximate and
exact solutions are plotted on the same coordinate plane using MATLAB code (2022).

2020 Mathematics Subject Classifications: 34K05,45E20
Key Words and Phrases: Integro-differential equations IDE, Volterra-Fredholm integral equa-
tion V-FIE, system of second-order Fredholm integral equations SFIEs, Collocation Method and
Galerkin method

1. Introduction

Integro-differential equations IDE have garnered growing interest from the math-
ematical and physics communities. These equations appear often in a wide range of
application domains including engineering, mechanics, elastic theory, probability theory,
and mathematical physics. Also, arise in fluid dynamics, such as the glass-forming pro-
cess and nano-hydrodynamics, falling condensation, biological models, chemical kinetics,
ecology, and control theory in financial mathematics, space systems, and industrial math-
ematics. see [12][10][15]. Recently, the authors have used various methods to display the
numerical or analytical solutions of IDEs. Therefore, many authors worked on semi-
analytical methods such as sequential Taylor expansion method, see [2][13], Bessel collo-
cation method, see[17], the variational iteration method, see[14], Haar functions method,
see[5][3], Legendre-Spectral method, see[16] [11], multi-wave Legendre method, see [9],
Legendre matrix method, see [18] and differential transform method, see [4].

In this article, we study the numerical solutions for IDEs of order two.

∗Corresponding author.
DOI: https://doi.org/10.29020/nybg.ejpam.v17i1.5024

Email addresses: fmharbi@uqu.edu.sa (F. M. Alharbi),
sharifah.althubiti10@gmail.com (S. S. Althubiti)

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



F. M. Alharbi, S. S. Althubiti / Eur. J. Pure Appl. Math, 17 (1) (2024), 286-299 287

2. Formulation of the problem

Consider the IDE,

µy′′(t) +B1y
′(t) +B2y(t) +

∫ a

0
k(t, τ)y(τ)dτ = f(t) (1)

with initial conditions:
y′(0) = q1 y(0) = q0 (2)

where : y′′(t) = d2y
dt2

and y(t) is the unknown function in the Hilbert space and it is
continuous with its derivatives. The parameters µ may have a physical meaning. The
known function f(t) is continuous. the k(t, τ) is the IDE’s kernel which is a continuous
function or at least satisfies Fredholm’s condition. B1 and B2 are known continuous
functions in the class space L2[0, a] with their derivatives.

To solve the IDE, we must transform it into a Volterra-Fredholm integral equation
V-FIE. Then the existence of a unique solution of the equation (2.1) under the given
conditions (2.2) is provided. the essential requirements that ensure the result is unique.
see[1].

After that, by using algebraic techniques, we reduce the V-FIE to an aliner system
of Fredholm integral equation SFIEs. The Collocation and Galerkin methods are used
to generate a numerical solution of a linear system of algebraic equations. which are then
solved using numerical methods. Furthermore, in each method, the estimating error is
computed and plotted. see[7]

Assume that,
y′′(t) = g(t) (3)

Integrating equation (2.3) two times. to follow,

y(t) =

∫ t

0
(t− x)g(x)dx+ tq1 + q0 (4)

Using the preceding results from equation (2.1),

µg(t) +

∫ t

0
ϕ(t, x)g(x)dx+

∫ a

0
ψ(t, τ)g(x)dx = F (t) (5)

where,
ϕ(t, x) = B1 + (t− x)B2 (6)

ψ(t, τ) =

∫ t

0
k(t, τ)(τ − x)dτ (7)

F (x) = f(x)− [B1q1 +B2tq1 +B2q0]− [

∫ a

0
k(t, τ)(tq1 + q0)dτ ] (8)

In the space L2[0, a] × C[0, T ], T < ∞ t, τ ∈ [0, a], equation (2.5) is known as the
V-FIE. where the Fredholm integral term is positive and continuous kernel ψ(t, τ). while
the Volterra integral term is considered in time with a positive continuous kernel ϕ(t, x)
For all t, x ∈ [0, T ], T < ∞. The free term F (x, t) is the surface integral equation (2.5)
and it is a known continuous function in the space L2[a, b]



F. M. Alharbi, S. S. Althubiti / Eur. J. Pure Appl. Math, 17 (1) (2024), 286-299 288

2.1. The main conditions

In order to ensure a unique solution to equation(2.1), we assume the following
conditions must be satisfies :
(i) The kernel of the integral term K(t, τ) must be continuous or at least satisfy the
Fredholm condition.

|
∫ a

0

∫ a

0
k2(t, τ)dtdτ |

1
2≤ α α is constant.

(ii) For the constants A1, A2 the given continuous functions B1(t), B2(t) satisfies the fol-
lowing conditions:

| B1(t) |≤ A1, | B2(t) |≤ A2

(iii) The norm of surtace function f(x) is defined as:

∥ f(x) ∥= [

∫ a

0
f2(x)dx]

1
2 ≤ γ, γ is constant.

(iv) The unknown function y(t) in the Hilbert space L2[0, a] behaves as the given function
f(t)

2.2. The normality and continuity of the integral operator:

Theorem 1. The integral equation (2.1) under the previous conditions (i)-(iv) has a
distinct solution.

proof :
To prove the existence of a unique solution of equation (2.5), we use the normality

and continuity of the mixed integral equation. For this, the integral equation (2.5) can be
written in the integral operator form:

W̄g(t) = F (t) +Wg(t) (9)

Wg(t) = ψg(x) + ϕg(x) (10)

ψg(t) = − 1

µ

∫ a

0
| (t− x) | g(x)dx (11)

ϕg(t) = − 1

µ

∫ t

0
ϕ(t, x)g(x)dx (12)

If conditions (i)-(iii) are satisfied, then we have respectively found:
(a) There exists a constant Ωψ, such that:[∫ a

0

∫ a

0
ψ2(|t− τ |)dτdt

] 1
2

≤ Ωψ (13)



F. M. Alharbi, S. S. Althubiti / Eur. J. Pure Appl. Math, 17 (1) (2024), 286-299 289

To obtain this, we take the norm of equation (2.7) and apply the Cauchy-Schwarz
inequality.

∥ψ(t, τ)∥ ≤ 1

µ

∥∥∥∥∥
∫ a

0

(
k2(t, τ)dτ

) 1
2

(∫ a

0
(τ − x)2dτ

) 1
2

∥∥∥∥∥ (14)

By using condition (i), we get:

∥ψ(t, τ)∥ ≤ α

µ

(∫ a

0

[
−x3 − 3ax2 + 3a2x− a3

3

]
dt

) 1
2

(15)

Hence

∥ψ(|t− τ |)∥ ≤ α

µ

(
a4

12

) 1
2

= Ωψ < 1 (16)

(b) There exists a constant Ωϕ such that,

∥ϕ(t, x)∥ ≤ Ωϕ, (17)

Taking the norm of equation (2.6) and applying condition (ii), we get:

∥ϕ(t, x)∥ ≤ 1

µ
[|B1|+ |B2|∥(t− x)∥]

≤ 1

µ
(A1 +A2∥t∥) ≤

1

µ
(A1 +A2(

a3

3
)
1
2 ) ≤ Ωϕ < 1

(18)

(c) There exists a constant ΩF such that

∥F (x)∥ ≤ ΩF (19)

Where
1

µ
∥f(x)∥ ≤ γ < 1 (20)

1

µ
∥(B1q1 +B2q0 +B2q1t)∥ ≤ ϵ < 1 (21)

1

µ
∥[
∫ a

0
k(t, τ)(tq1 + q0)dτ ]∥ ≤ α

µ
(q0 + q1(

a3

3
)
1
2 ) = σ < 1 (22)

ΩF = (γ + ϵ+ σ) (23)

Theorem 2. If conditions (a)-(c) are satisfied, and the integration factor (2.9) is normal
and continuous, then equation (2.5) has a unique solution in Banach space L2 [0, a], under
the condition,

|λ| ≤
1− Ωϕ
Ωψ

(24)



F. M. Alharbi, S. S. Althubiti / Eur. J. Pure Appl. Math, 17 (1) (2024), 286-299 290

Proof :
We establish the normality and continuity of the integral operator (2.9).

(a) For the normality of the integral operator Wg, we write:

∥Wg∥ ≤ −1

µ
[

∥∥∥∥∫ t

0
ϕ(t, x)g(x)dx

∥∥∥∥+

∥∥∥∥∫ a

0
ψ(|t− x|)g(x)dx

∥∥∥∥]. (25)

Then,

∥ϕg(x)∥ ≤ 1

µ
(A1 +A2∥t∥)∥g(x)∥ ≤ 1

µ
(A1 +A2(

a3

3
)
1
2 )∥g(x)∥, (26)

which can be adapted as:
∥ϕg(x)∥ ≤ Ωϕ∥g(x)∥, . (27)

also,

∥ψg(x)∥ ≤
∥∥∥∥∫ a

0
ψ(|t− x|)g(x)dx

∥∥∥∥ . (28)

By using condition (a), we get:

∥ψg(x) ∥≤ Ωψ∥ g(x)∥, (29)

and them,
∥Wg(x)∥ ≤ χ∥g(x)∥, χ = (Ωϕ +Ωψ) . (30)

So, W is a norm operator that leads straight to the normality of the operator W̄ after
applying condition (c).
(b) We assume that the two potential functions g1(x), g2(x) in the Helbert space, then,

∥∥W̄ (g1 − g2)
∥∥ ≤

∥∥∥∥∫ t

0
ϕ(t, x) (g1(x)− g2(x)) dx

∥∥∥∥
+

∥∥∥∥∫ a

0
ψ(|t− x|) (g1(x)− g2(x)) dx

∥∥∥∥ . (31)

Using conditions (a) and (b), we get:∥∥W̄ (g1 − g2)
∥∥ ≤ χ ∥g1(x)− g2(x)∥ . (32)

That proves W̄ is a continuous operator,then, by using the condition χ < 1, we deduce
that W̄ is a contraction operator, and it has a unique solution.

2.3. System of Fredholm integral equations SFIEs

The quadratic method has wide application in mathematical and physics problems,
where the eigenvalues and eigenfunctions of integral equations are often studied and dis-
cussed. It also has wide applications in applied sciences, especially in the theory of elas-
ticity, mixed problems in the fluid of mechanics, and communication problems. This
numerical technique will be applied in this section to reduce V-FIEs to linear SFIEs.



F. M. Alharbi, S. S. Althubiti / Eur. J. Pure Appl. Math, 17 (1) (2024), 286-299 291

Consider,

µg(t) =F (t) +

∫ t

0
ϕ(t, x)g(x)dx+

∫ a

0
ψ(t, τ)g(x)dx (33)

Divide the interval [0, T ] as 0 = t0 ≤ t1 ≤ · · · ≤ tN = T . Using the quadrature formula,
equation (4.1) becomes:∫ t

0
ϕ(t, x)g(x)dx =

n∑
m=0

umϕ(tn, xm)g(xm) (34)

Where n = 0, 1, 2, · · · , N − 1, and u0 =
1
2h0, un = 1

2hn, ui = hi, (i ̸= 0, n).
Using (2.34) in (2.33), we have:

g(t) =F (t) +

∫ a

0
ψ(t, τ)g(x)dx+

n∑
m=0

umϕ (tn, xm) g(xm) (35)

Then:

gn = Yn(t) +

∫ a

0
ψ(t, τ)gn(x)dx (36)

Where Yn(t) = Fn(t) +
∑n

m=0 umϕ (tn, xm) g(xm), n = 0, 1, · · · , N .
Formula (2.36) represents a SFIEs , and we have N unknown functions gn(t) correspond-
ing to time interval [0, T ].

3. Numerical Methods

In this section, the Collocation method and Galerkin method are used to solve
V-FIE of the second kind.

3.1. Collocation Method

To find the solution to equation (2.5), use the Collocation method. We approxi-
mate the unidentified function g(x) by using the function Li(x).

S(t, xi) =
n∑
i=1

ciLi(x) (37)

Given a set of n linearly independent functions L1(x), L2(x), . . . , LN (x) defined on the
interval (0, a). Therefore, we have:

µSn(t) ≈ F (t, xi) +

∫ a

0
ψ(t, x)Sn(x)dx+

n−1∑
m=0

umϕjmSn(x)

+ ϵ
(
x, c1(t), c2(t), . . . , cN (t) +R

(
hp+1
i

))
.j = 0, 1, . . . , n

(38)



F. M. Alharbi, S. S. Althubiti / Eur. J. Pure Appl. Math, 17 (1) (2024), 286-299 292

Of course, if the approximate solution (3.1) is substituted with (2.36) for the function
g(t, x), there would inevitably be an error denoted as ε (x, c1(t), c2(t), . . . , cN (t)). The
extent of this error is dependent on the selection of coefficients in the formula. The value
of x is selected. Given that x is equal to xi for i ranging from 0 to N , we may express
this as:

µ
n∑
i=1

ciLi(x) ≈ F (t, xi) +

∫ a

0
ψ(t, x)

n∑
i=1

ciLi(x)dx

+
n−1∑
m=0

umϕjm

n∑
i=1

ciLi(x)

+ ϵ
(
x, c1(t), c2(t), . . . , cN (t) +R

(
hp+1
i

))
, j = 1, . . . , n

(39)

For determining the coefficients c1 (ti) , c2 (ti) , . . . , cn (ti) of the approximate solution
Sn (xj), as given in equation (3.1), using n linearly independent functions L1(x), L2(x), . . . , Ln(x)
defined on the numerical interval [0, a]. Hence, we need to carry out the process of in-
tegration and subsequently replace x with x1, x2, . . . , xN to identify the points at which
the error ϵ (x, c1(t), c2(t), . . . , cN (t)) becomes zero. By substituting the equation (3.1) into
equation (3.2), we obtain:

µ

n∑
i=1

ciLi(x) ≈ F (t, xi) +

∫ a

0
ψ(t, x)

n∑
i=1

ciLi(x)dx+

n−1∑
m=0

umϕjm

n∑
i=1

ciLi(x)

+ ϵ
(
x, c1(t), c2(t), . . . , cN (t) +R

(
hp+1
i

))
, j = 1, . . . , n

This gives us:

µ
n∑
i=1

ciLi(x) ≈ F (t, xi) +
n∑
i=1

ci

∫ a

0
ψ(t, x)Li(x)dx

+
n−1∑
m=0

n∑
i=1

umϕjmciLi(x)

+ ϵ
(
x, c1(t), c2(t), . . . , cN (t) +R

(
hp+1
i

))
, j = 0, . . . , n

see[6, 7]

3.2. Galerkin method

The Galerkin method is used to obtain an approximate solution to equation (2.36).
The approach sets the necessary conditions for calculating a set of n coefficients, as stated
in Equation (3.1). By introducing the error ε (x, c1(t), c2(t), . . . , cN (t)) in equation (3.2)
which is perpendicular to n linearly independent functions L1(x), L2(x), . . . , LN (x) on the



F. M. Alharbi, S. S. Althubiti / Eur. J. Pure Appl. Math, 17 (1) (2024), 286-299 293

interval (0, a), i.e. ∫ a

0
Lj(x)ϵ(x, c1(t), c2(t)), . . . , cN (t)dx = 0 (40)

Then from (2.34), we have

µ
n∑
i=1

ciLi(x)−
∫ a

0
ψ(t, x)

n∑
i=1

ciLi(x)dx−
n−1∑
m=0

umϕjm

n∑
i=1

ciLi(x)

+ ϵ
(
x, c1(t), c2(t), . . . , cN (t) +R

(
hp+1
i

))
= F (t)

(41)

Then the Galerkin equations are obtained by multiplying both sides of (3.5) by Lj(x)
and then integrating with respect to x from 0 to a, we obtain

n∑
i=1

ci

∫ a

0

[
µLi(x)−

∫ a

0
ψ(t, x)Li(x)dx−

n−1∑
m=0

umϕjmLi(x)

]
Li(x)dx (42)

=

∫ a

0
F (t)Li(x)dx (43)

The unknown parameters ci are determined by solving the system of equations men-
tioned above and inserting these values of parameters in trial functions. We obtain an
approximate solution, denoted as gh, of the V-FIE .see [7, 8]

4. Applications.

Consider the following applications:

Application (1)

Consider the IVP,

y′′(t) = 2− cos t+

∫ π

0
ty(t)dt (44)

under the initial conditions
y′(0) = 0 y(0) = 0 (45)

The exact solution is y(t) = | cos t|.
After converting it to V-FIE we get:

ϕ(t)−
∫ π

0
ψ(x)ϕ(x)dx = F (t)

where,

ψ(x) =

∫ t

0
t(t− x)dt

F (x) = 2− cos t.

(46)



F. M. Alharbi, S. S. Althubiti / Eur. J. Pure Appl. Math, 17 (1) (2024), 286-299 294

• Using Collocation method

Let us consider the approximate answer the equation (4.1), as the three independent
functions,
L0(t) = 1, L1 = t, L2(t) = |cos(t)|. By substituting these functions into equation (4.3) and
then solving the resulting equation when t = 0, π2 , π, we get:

c0 = −1.34404356

c1 = 0.382981562

c2 = 2.34404356

Therefore, the approximate solution is

S1(t) = −1.34404356 + 0.382981562t+ 2.34404356|cos(t)|.

• Using Galerkin method

As the same consideration independent approximate solution and the same arbitrary
points in the collocation method, we get,

c0 = −1.32446935

c1 = 0.332647197

c2 = 2.18043309

Therefore, the approximate solution :

S2(t) = −1.32446935 + 0.332647197t+ 2.18043309|cos(t)|.

The following figures discuss the shape of the numerical solution of the two methods and
the relation between the estimating errors that were obtained.



F. M. Alharbi, S. S. Althubiti / Eur. J. Pure Appl. Math, 17 (1) (2024), 286-299 295

Figure 1: The relation between the exact solution and numerical solution in Collocation approximate.

Figure 2: relation between the exact solution and numerical solution in Galerkin approximate.

Figure 3: The relation between estimating error of Collocation approximate and Galerkin approximate.

Application (2)

Consider the IVP,

y′′(t) = 1− e+ et +

∫ 1

0
y(t)dt (47)

under the initial conditions
y′(0) = 1 y(0) = 1 (48)

The exact solution is y(t) = et

After converting it to F-VIE We get:



F. M. Alharbi, S. S. Althubiti / Eur. J. Pure Appl. Math, 17 (1) (2024), 286-299 296

ϕ(t)−
∫ 1

0
ψ(x)ϕ(x)dx = F (x)

where,

ψ(x) =

∫ t

0
(t− x)dtF (x) =

5

2
− e+ et. (49)

• Using Collocation method

Let us consider the approximate answer the equation (4.1), as the three independent
functions,
L0(t) = 1, L1 = t, L2(t) = t2. By substituting these functions into equation (4.6) and then
solving the resulting equation when t = 0, π2 , π, we get:

c0 = 0.7817181715

c1 = −0.267748300

c2 = 1.98603013

Therefore, the approximate solution is

S2(t) = 0.7817181715− 0.267748300t+ 1.98603013t2.

• Using Galerkin method

As the same consideration independent approximate solution and the same arbitrary
points in the collocation method, we get,

c0 = 0.844299425

c1 = −0.370315407

c2 = 1.90695998

Therefore, the approximate solution:

S2(t) =0.844299425− 0.370315407t+ 1.90695998t2.

The results are shown in figure 4 and 5 for n = 6.



F. M. Alharbi, S. S. Althubiti / Eur. J. Pure Appl. Math, 17 (1) (2024), 286-299 297

Figure 4: The relation between the exact solution and numerical solution in Collocation approximate.

Figure 5: relation between the exact solution and numerical solution in Galerkin approximate.

Figure 6: The relation between estimating error of Collocation approximate and Galerkin approximate.

5. Conclusion

In this study, we focus on the initial value problem of a linear integro-differentail
equation with a continuous kernel, or at least satisfy the Fredholm condition. By convert-
ing the studding equation to a Volterra-Fredholm integral equation of the second kind it
transforms into a linear system of Fredholm integral equations. The collocation method
and the Galerkin method are used to efficiently discretize the integral equations using a
set of points. We discuss the error estimates associated with both methods.



REFERENCES 298

References

[1] FM Alharbi and MA Abdou. Boundary and initial value problems and integral op-
erator. Adv. Differ. Equ. Control Process, 19, 2018.

[2] Azizallah Alvandi and Mahmoud Paripour. Reproducing kernel method with tay-
lor expansion for linear volterra integro-differential equations. Communications in
Numerical Analysis, 1, 2017.

[3] Rohul Amin, Ibrahim Mahariq, Kamal Shah, Muhammad Awais, and Fahmi Elsayed.
Numerical solution of the second order linear and nonlinear integro-differential equa-
tions using haar wavelet method. Arab Journal of Basic and Applied Sciences, 28,
2021.

[4] Aytac Arikoglu and Ibrahim Ozkol. Solutions of integral and integro-differential
equation systems by using differential transform method. Computers & Mathematics
with Applications, 56, 2008.

[5] Imran Aziz et al. New algorithms for the numerical solution of nonlinear fredholm
and volterra integral equations using haar wavelets. Journal of Computational and
Applied Mathematics, 239, 2013.

[6] Angelamaria Cardone, Dajana Conte, Raffaele D’Ambrosio, and Beatrice Paternoster.
Collocation methods for volterra integral and integro-differential equations: A review.
axioms, 7, 2018.

[7] Fatheah Ahmad Hendi and Abeer Majed Albugami. Numerical solution for fredholm–
volterra integral equation of the second kind by using collocation and galerkin meth-
ods. Journal of King Saud University-Science, 22, 2010.

[8] K Issa and F Salehi. Approximate solution of perturbed volterra-fredholm integrod-
ifferential equations by chebyshev-galerkin method. Journal of Mathematics, 2017,
2017.

[9] Imran Khan, Muhammad Asif, Rohul Amin, Qasem Al-Mdallal, and Fahd Jarad. On
a new method for finding numerical solutions to integro-differential equations based
on legendre multi-wavelets collocation. Alexandria Engineering Journal, 61, 2022.

[10] Vangipuram Lakshmikantham. Theory of integro-differential equations, volume 1.
CRC press, 1995.

[11] Mahmoud Lotfi and Amjad Alipanah. Legendre spectral element method for solving
volterra-integro differential equations. Results in Applied Mathematics, 7, 2020.

[12] Allen C Pipkin. A course on integral equations. Number 9. Springer Science &
Business Media, 1991.



REFERENCES 299

[13] J Rashidinia and A Tahmasebi. Taylor series method for the system of linear volterra
integro-differential equations. System, 50, 2012.

[14] Jafar Saberi-Nadjafi and Mohamadreza Tamamgar. The variational iteration method:
A highly promising method for solving the system of integro-differential equations.
Computers Mathematics with Applications, 56, 2008.

[15] Abdul-Majid Wazwaz. First Course In Integral Equations, A. World Scientific Pub-
lishing Company, 2015.

[16] Yunxia Wei and Yanping Chen. Convergence analysis of the legendre spectral col-
location methods for second order volterra integro-differential equations. Numerical
Mathematics: Theory, Methods and Applications, 4, 2011.

[17] Şuayip Yüzbaşı. Improved bessel collocation method for linear volterra integro-
differential equations with piecewise intervals and application of a volterra population
model. Applied Mathematical Modelling, 40, 2016.

[18] Şuayip Yüzbaşı, Mehmet Sezer, and Bayram Kemancı. Numerical solutions of integro-
differential equations and application of a population model with an improved legen-
dre method. Applied Mathematical Modelling, 37, 2013.