# Approximation Techniques for Solving Linear Systems of Volterra Integro-Differential Equations

Authors: Ahmad Issa; Naji Qatanani; Adnan Daraghmeh
Journal: Journal of Applied Mathematics
 (2020)

---

## Abstract

In this paper, a collocation method using sinc functions and Chebyshev wavelet method is implemented to solve linear systems of Volterra integro-differential equations. To test the validity of these methods, two numerical examples with known exact solution are presented. Numerical results indicate that the convergence and accuracy of these methods are in good a agreement with the analytical solution. However, according to comparison of these methods, we conclude that the Chebyshev wavelet method provides more accurate results.


---

## Body

## 1. Introduction

Systems of integro-differential equations have motivated huge amounts of research in recent years. They arise in many physical phenomena like wind ripple in the desert, nano-hydrodynamics, population growth model, glass-forming process, and oceanography [1–3]. Various numerical methods for solving systems of linear integro-differential equations have been developed by many researchers. Hesameddini and Rahimi [4] used the reconstruction of variational iteration method (RVIM) for solving systems of Volterra integro-differential equations. In [5], Hesameddini and Asadolhifard, implemented the sinc-collocation method to approximate the solution of systems of linear Volterra integro-differential equations with initial conditions. Aminikhah and Hosseini [6] applied the wavelet method for the numerical solution of systems of integro-differential equations. They used the operational matrix of integration to solve these systems. Draidi and Qatanani [7] emplemented product Nystrom and sinc-collocation methods to solve Volterra integral equation with Carleman kernel. Hamaydi and Qatanani [8] used the Taylor expansion and the variational iteration methods to give approximate solution of Volterra integral equation of the second kind. In addition, Issa [9] has employed several numerical techniques for solving systems of Volterra integro-differential equations. Other numerical methods for systems of integro-differential equations are (power) functions and Chebyshev polynomials [10], single term Walsh series [11], Chebyshev collocation [12], rationalized Haar functions [13], differential transform [14], homotopy perturbation [15], power series [16], and finite difference approximation [17]. Regarding the stability of a system of Volterra integro-differential equations, some stability results are proposed for the linear system VIDEs in the 1980s, those of Burton are worthy to mention. His work [18, 19] laid the foundation for a systematic treatment of the basic structure and stability properties of VIDEs via the direct method of Lyapunov. A more recent result is by Elaydi [20], who proposed a type of Lyapunov functional that is also applicable to delay equations. Moreover, Zhang [21] proposed recently a stability result from which certain well-known result could be derived. Also, Vanualailai and Nakagiri [22] have proposed a new stability criteria based on new and known forms of Lyapunov functionals for a system of Volterra integro-differential equations. In this article, we propose two numerical methods, namely, a collocation method using sinc functions and Chebyshev wavelet method to approximate the solution of a system of linear Volterra integro-differential equations given by(1)
u
i
n
x
=
f
i
x
+
∫
a
x
∑
j
=
1
N
k
i
j
x
,
t
u
j
t
d
t
,
a
≤
x
≤
b
,
1
≤
i
≤
Nsubject to the initial conditions(2)
u
i
s
=
a
i
s
,
i
=
1
,
2
,
3
,
…
,
N
,
s
=
0
,
1
,
2
,
…
,
n
−
1
.The kernelsk
i
j
x
,
t and the function f
i
x are given real valued functions and the unknown functions u
i
x are to be determined. A comparison between these methods is carried out by solving some numerical examples.The paper is organized as follows: In Section2 the sinc-collocation method based on sinc functions is presented. The Chebyshev wavelet method is addressed in Section 3. In Section 5, the proposed methods are implemented using two numerical examples with known analytical solution by applying MAPLE software. Conclusions are followed in Section 6.


## 2. Sinc Collocation Method Based on Sinc Functions

The sinc collocation method based on sinc functions is widely used for obtaining the approximate solution of ordinary and partial differential equations and integral equations [5]. It is well-known that the sinc approximate solution converges exponentially to the exact solution.Definition 1 (see [23]).
The sinc function is defined on the whole real line−
∞
<
x
<
∞ by(3)
sinc
x
=
sin
π
x
π
x
,
x
≠
0
1
,
x
=
0
as shown in Figure1.Figure 1
Sinc function.Definition 2 (see [23]).
Letk
=
0
,
±
1
,
±
2
,
±
3
,
… then the translated sinc basis functions are defined as(4)
s
k
,
r
x
=
sinc
x
−
k
r
r
=
sin
π
/
r
x
−
k
r
π
/
r
x
−
k
r
,
x
≠
k
r
1
,
x
=
k
r
which are called thek
th sinc functions.Corollary 1 (see [5]).
The sinc function for the interpolating pointsx
d
=
d
r, is given by(5)
s
k
,
r
d
r
=
μ
k
d
=
0
,
k
≠
d
1
,
k
=
dCorollary 2 (see [5]).
Ifp
x is defined on the real axis and r is a positive integer, then the series(6)
a
p
,
r
x
=
∑
k
=
−
∞
∞
p
k
r
s
k
,
r
x
is called the Whittaker Cardinal expansion ofp
x.The properties of the Whittaker Cardinal expansion have been extensively studied in [22]. These properties are derived in the infinite strip D
s of the complex w–planes where for any g
>
0,(7)
G
s
=
w
=
t
+
i
s
:
s
<
g
≤
π
2
.To construct an approximation on the intervala
,
b, we use the conformal map:(8)
w
=
φ
z
=
ln
z
−
a
b
−
z
.This map carries the eye-shaped region(9)
G
E
=
z
∈
C
:
arg
z
−
a
b
−
z
<
g
≤
π
2
.For the sinc method, the basis functions ona
,
b for z
∈
G
E are derived from the composite translated sinc functions,(10)
s
k
,
r
∘
φ
x
=
sinc
φ
x
−
k
r
r
,wheres
k
,
r
∘
φ
x
=
s
k
,
r
φ
x.The inverse map ofw
=
φ
z is(11)
z
=
φ
−
1
w
=
a
+
b
e
w
1
+
e
w
.Also we define the range ofφ
−
1 on the real line as(12)
τ
=
y
<
u
<
∞
φ
−
1
u
∈
G
E
:
−
∞
<
u
<
∞and the interpolation points{
x
d
} are then given by:(13)
x
d
=
φ
−
1
d
r
=
a
+
b
e
d
r
1
+
e
d
r
,
d
=
0
±
1
,
±
2
,
…
.Definition 3 (see [23]).
LetL
α
G
E be the set of all analytic functions. Then, there exists a constant c, such that:(14)
u
z
≤
c
β
z
α
1
+
β
z
2
α
,
z
∈
G
E
,
0
<
α
≤
1
,
whereβ
z
=
e
φ
z.Theorem 1 (see [23]).
Letu
∈
L
α
G
E and N is a natural number, and r be selected by the formula(15)
r
=
π
g
α
N
1
/
2
,
where(16)
0
<
α
≤
1
,
g
≤
π
2
.
Then, there exists a positive constantc
1, independent of N, such that(17)
sup
z
∈
τ
u
z
−
∑
d
=
−
N
N
u
z
d
s
d
,
r
∘
φ
z
≤
c
1
e
−
π
g
α
N
1
/
2
.Theorem 2 (see [5]).
Letu
/
φ
∈
L
α
G
E and N is a natural number, and r be given as(18)
r
=
π
g
α
N
1
/
2
,
where(19)
0
<
α
≤
1
,
g
≤
π
2
.
Moreover, letμ
k
d
−
1 be defined as(20)
μ
k
d
−
1
=
1
2
+
∫
0
k
−
d
sin
π
t
π
t
d
t
.
Then, there exists a positive constantc
2, independent of N, such that(21)
∫
a
z
k
u
t
d
t
−
r
∑
d
=
−
N
N
μ
k
d
−
1
u
z
d
φ
z
d
≤
c
2
e
−
π
g
α
N
1
/
2
.Theorem 3 (see [5]).
Letφ be a conformal injective map of the simply connected domain G
E
onto
G
E. Then(22)
μ
k
d
0
=
s
k
,
d
∘
φ
x
x
=
x
d
=
1
,
k
=
d
0
,
k
≠
d
μ
k
d
1
=
r
d
d
φ
s
k
,
d
∘
φ
x
x
=
x
d
=
0
,
k
=
d
−
1
d
−
k
d
−
k
,
k
≠
d
μ
k
d
2
=
r
2
d
2
d
φ
2
s
k
,
d
∘
φ
x
x
=
x
d
=
−
π
2
3
,
k
=
d
−
2
−
1
d
−
k
d
−
k
2
,
k
≠
dWe consider the system of linear Volterra integro-differential equations of the form:(23)
u
i
n
x
=
f
i
x
+
∫
a
x
∑
j
=
1
N
k
i
j
x
,
t
u
j
t
d
t
,
1
≤
i
≤
N
.Subject to the initial conditions:(24)
u
i
s
0
=
a
i
s
,
i
=
1
,
2
,
…
,
N
,
s
=
0
,
1
,
2
,
…
,
n
−
1
,in the domain0
,
1, and let u
i
x
∈
L
α
G
E. By using Theorem 1, u
i
x is approximated as follows:(25)
u
i
x
=
R
i
x
+
A
i
x
,where(26)
R
i
x
=
∑
k
=
−
N
N
c
k
i
w
x
sinc
φ
x
−
k
r
r
,
A
i
x
=
∑
j
=
0
n
a
j
i
x
j
,wherec
k are unknown coefficients and w
x
=
x
n
x
−
1
n.Integrating both sides of Equation (25) from 0 to x we get(27)
∫
0
x
u
i
t
d
t
=
∫
0
x
R
i
t
d
t
+
∫
0
x
A
i
t
d
tand by differentiating both sides of Equation (25) with respect to x we get(28)
u
i
n
x
=
R
i
n
x
+
A
i
n
x
,where(29)
R
i
n
x
=
∑
k
=
−
N
N
c
k
d
n
w
x
sinc
φ
x
−
k
r
/
r
d
x
n
.Substituting Equations (27) and (29) into Equation (23), and by evaluating the result at the sinc points(30)
x
j
=
e
j
r
1
+
e
j
r
,wherej
=
−
N
−
1
,
…
,
N, and using Theorems 2 and 3, we obtain a system of algebraic equations. Solving this system we obtain the unknown coefficients(31)
c
k
i
k
=
−
N
Nand(32)
a
j
i
j
=
0
n
.


## 3. Chebyshev Wavelets Method (CWM)

The main idea of using Chebyshev basis is that the problem under study reduces to a system of linear or nonlinear algebraic equations. This may be done by truncated series of orthogonal basis functions for the solution of the problem and using the operational matrices [6]. Wavelets constitute a family of functions constructed from dilation and translation of a single function called the mother wavelet [24–26]. When dilation p and translation q vary continuously, we have the following family of continuous wavelets as(33)
ψ
p
,
q
x
=
1
p
ψ
x
−
q
p
,
p
,
q
∈
R
,
p
≠
0
.If we choose the dilation and translationp
−
a, and b
q
p
−
a, respectively, where p
>
1, q
>
0, then we have the following family of continuous wavelets as(34)
ψ
a
,
b
x
=
p
a
ψ
p
−
a
x
−
b
q
,
a
,
b
∈
Z
+
,whereψ
a
,
b forms a wavelet basis for(35)
L
2
R
=
g
:
R
→
C
∫
−
∞
∞
g
x
2
d
x
<
∞
,whereL
2
R: set of all square integrable functions equipped with norm(36)
f
L
2
a
,
b
=
∫
a
b
f
x
2
d
x
1
/
2
.For the particular case, whenp
=
2 and q
=
1, then ψ
a
,
b
x forms an orthogonal basis.Chebyshev waveletsψ
b
,
c
x
=
ψ
a
,
b
,
c
x have four parameters, b
=
1
,
2
,
3
,
…
,
2
a
−
1, a
∈
Z
+ and c is the degree of Chebyshev polynomials of the first kind. They are defined on the interval 0
≤
x
≤
1 by:(37)
ψ
b
,
c
x
=
ψ
a
,
b
,
c
x
=
2
a
/
2
T
~
c
2
a
x
−
2
b
+
1
,
for
b
−
1
2
a
−
1
≤
x
≤
b
2
a
−
1
,
0
,
otherwise,where(38)
T
~
c
x
=
1
π
,
c
=
0
1
π
T
c
x
,
c
>
0andc
=
0
,
1
,
2
,
…
,
C
−
1, and b
=
1
,
2
,
3
,
…
,
2
a
−
1.T
c
x are the famous Chebyshev polynomials of the first kind of degree c which are orthogonal with respect to the weight function(39)
w
x
=
1
−
x
2
,
−
1
≤
x
≤
1and satisfy the following recurrence relation(40)
T
0
x
=
1
,
T
1
x
=
x
,
T
c
+
1
x
=
2
x
T
c
x
−
T
c
−
1
x
,
c
=
1
,
2
,
3
,
…Corollary 3 (see [6]).
The set of Chebyshev wavelets is an orthogonal set with respect to the weight function(41)
w
b
x
=
w
2
a
x
−
2
b
+
1
.Definition 4 (see [6]).
A functionh
x defined on the interval 0
≤
x
≤
1 is called the wavelet series if this function is written in the following form(42)
h
x
=
∑
b
=
1
∞
∑
c
=
0
∞
d
b
c
ψ
b
c
x
,
where(43)
d
b
c
=
h
x
,
ψ
b
c
x
w
b
x
is the inner product inL
w
b
2
0
,
1.Corollary 4 (see [6]).
The wavelet series inL
2
0
,
1 is convergent if(44)
lim
u
1
,
u
2
→
∞
h
x
−
∑
b
=
1
u
1
∑
c
=
0
u
2
d
b
c
ψ
b
c
x
=
0
.
If the wavelet series is truncated, then it can be written as(45)
h
x
≅
∑
b
=
1
2
a
−
1
∑
c
=
0
c
−
1
d
b
c
ψ
b
c
x
=
D
T
ψ
x
,
whereD and ψ
x are 2
a
−
1
c
×
1 matrices given by(46)
D
=
d
1
,
0
,
d
1
,
1
,
…
,
d
1
,
c
−
1
,
d
2
,
0
,
d
2
,
1
,
…
,
d
2
,
c
−
1
,
…
,
d
2
a
−
1
,
0
,
…
,
d
2
a
−
1
,
c
−
1
T
=
d
1
,
d
2
,
…
,
d
c
,
d
c
+
1
,
…
,
d
2
a
−
1
,
c
T
,
ψ
x
=
ψ
1
,
0
x
,
ψ
1
,
1
x
,
…
,
ψ
1
,
c
−
1
x
,
ψ
2
,
0
x
,
ψ
2
,
1
x
,
…
,
ψ
2
,
c
−
1
x
,
…
,
ψ
2
a
−
1
,
0
x
,
…
,
ψ
2
a
−
1
,
c
−
1
x
T
=
ψ
1
x
,
ψ
2
x
,
…
,
ψ
c
x
,
ψ
c
+
1
x
,
…
,
ψ
2
a
−
1
,
c
x
T
.Corollary 5 (see [27]).
The integral of the multiple of two Chebyshev wavelets vector functions with respect tow
b
x from 0 to 1 is an identity matrix. Moreover, a function h
x
,
y defined on 0
,
1
×
0
,
1 can be approximated as:(47)
h
x
,
y
≅
∑
i
=
1
2
a
−
1
c
∑
j
=
1
2
a
−
1
c
a
i
j
ψ
i
x
ψ
j
y
=
ψ
T
x
A
ψ
y
,
whereA
=
a
i
j is a matrix of the entries 2
a
−
1
c
×
2
a
−
1
c, that can be determined by:(48)
a
i
j
=
ψ
i
x
,
h
x
,
y
,
ψ
j
y
w
b
y
w
b
x
,
i
=
1
,
2
,
3
,
…
,
2
a
−
1
c and j
=
1
,
2
,
3
,
…
,
2
a
−
1
c.The integral of the vectorψ
x defined in Equation (46), is given as:(49)
∫
0
x
ψ
t
d
t
=
B
ψ
x
,whereB is the 2
a
−
1
c
×
2
a
−
1
c operational matrix of integration [2].This matrix has the form:(50)
B
=
2
−
a
M
E
E
⋯
E
O
M
E
⋱
⋮
O
O
M
⋱
⋮
⋮
⋱
⋱
⋱
E
O
…
O
O
M
,whereM, E, and O are C
×
C matrices given by(51)
M
=
1
1
2
0
0
0
⋯
0
−
2
4
0
1
4
0
0
⋯
0
−
2
3
−
1
2
0
1
6
0
⋯
0
⋮
⋱
⋱
⋱
⋱
⋱
⋮
2
−
1
λ
2
1
λ
−
1
−
1
λ
⋯
−
1
2
λ
−
1
0
1
2
λ
⋯
0
⋮
⋱
⋱
⋱
⋱
⋱
⋮
2
−
1
c
2
1
c
−
1
−
1
c
0
0
0
⋯
−
1
2
c
−
1
0
(52)
E
=
2
0
0
⋯
0
0
0
0
⋯
0
−
2
2
3
0
0
⋯
0
⋮
⋮
⋮
⋱
⋮
2
2
1
−
−
1
λ
λ
−
1
−
−
1
λ
−
2
λ
−
2
0
0
⋯
0
⋮
⋮
:
⋱
⋮
2
2
1
−
−
1
C
C
−
1
−
−
1
C
−
2
C
−
2
0
0
⋯
0
and
O
=
0
⋯
0
⋮
⋱
⋮
0
⋯
0
.The product characteristic of two Chebyshev wavelets vector functions are given as(53)
ψ
x
ψ
T
Z
≈
Z
~
ψ
x
,whereZ is a given vector and Z
∼
=
Z
∼
i
j
2
a
−
1
c
×
2
a
−
1
c is an operational matrix. We consider the system of linear Volterra integro-differential equations of the form:(54)
u
i
n
x
=
f
i
x
+
∑
j
=
1
c
∫
0
x
k
i
j
x
,
t
u
j
t
d
twith the following conditions(55)
u
i
r
0
=
a
i
r
,
i
=
1
,
2
,
3
,
…
,
n
,
r
=
0
,
1
,
2
,
…
,
n
−
1
,
c
=
1
,
2
,
3
,
…
.Now we approximateu
i
n
x by using Chebyshev wavelet space as follows(56)
u
i
n
x
=
D
i
T
ψ
x
,
i
=
1
,
2
,
3
,
…
,
n
.Therefore we have(57)
u
i
r
x
=
D
i
T
B
n
−
r
ψ
x
+
∑
j
=
0
n
−
r
−
1
a
i
r
x
j
j
!
,
i
=
1
,
2
,
3
,
…
,
n
,whereD and ψ
x are 2
a
−
1
c
×
1 matrices given by(58)
D
=
d
1
,
0
i
,
d
1
,
1
i
,
…
,
d
1
,
c
−
1
i
,
d
2
,
0
i
,
d
2
,
1
i
,
…
,
d
2
,
c
−
1
i
,
…
,
d
2
a
−
1
,
0
i
,
…
,
d
2
a
−
1
,
c
−
1
i
T
=
d
i
,
1
,
d
i
,
2
,
…
,
d
i
,
c
,
d
i
,
c
+
1
,
…
,
d
i
,
2
a
−
1
,
c
T
,
ψ
x
=
ψ
1
,
0
x
,
ψ
1
,
1
x
,
…
,
ψ
1
,
c
−
1
x
,
ψ
2
,
0
x
,
ψ
2
,
1
x
,
…
,
ψ
2
,
c
−
1
x
,
…
,
ψ
2
a
−
1
,
0
x
,
…
,
ψ
2
a
−
1
,
c
−
1
x
T
=
ψ
1
x
,
ψ
2
x
,
…
,
ψ
c
x
,
ψ
c
+
1
x
,
…
,
ψ
2
a
−
1
,
c
x
T
.In virtue of Equations (55) and (56) we have the following approximations:(59)
f
i
x
≅
F
i
T
ψ
x
,
u
i
x
≅
D
i
T
B
ψ
x
+
E
i
T
ψ
x
,
k
i
j
x
,
t
≅
ψ
T
k
i
j
ψ
t
,wherek
i
j and F
i are known matrices for i
=
1
,
2
,
3
,
…
,
n
,
j
=
1
,
2
,
3
,
…
,
c.By substituting the approximations Equations (56) and (59) into the system equation (54), we obtain:(60)
D
i
T
ψ
x
≅
F
i
T
ψ
x
+
∑
j
=
1
c
∫
0
x
ψ
T
x
k
i
j
ψ
t
D
i
T
B
ψ
t
+
E
i
T
ψ
t
d
t
=
F
i
T
ψ
x
+
∑
j
=
1
c
ψ
T
x
k
i
j
∫
0
x
ψ
t
D
i
T
B
ψ
t
+
E
i
T
ψ
t
d
t
=
F
i
T
ψ
x
+
∑
j
=
1
c
ψ
T
x
k
i
j
∫
0
x
ψ
t
D
i
T
B
+
E
i
T
ψ
t
d
t
.Therefore,(61)
D
i
T
ψ
x
≅
F
i
T
ψ
x
+
∑
j
=
1
c
ψ
T
x
k
i
j
Z
~
i
B
ψ
t
,
i
=
1
,
2
,
3
,
…
,
nc
=
1
,
2
,
3
,
…, where B is the 2
a
−
1
c
×
2
a
−
1
c operational matrix of integration and Z
∼
i are 2
a
−
1
c
×
1 matrices.We multiply both sides of Equation (61) by w
n
x
ψ
T
x and integrating with respect to x from 0 to 1, we obtain a linear system in terms of input D
i
,
i
=
1
,
2
,
3
,
…
,
n. Consequently, the vector functions D
i elements are calculated by solving this system.


## 4. Stability of Systems of Volterra Integro-Differential Equations (VIDEs)

In this section, we present some important results on the stability of VIDEs (1) (for more details see [22]).Definition 5 (see [22]).
Ifμ
:
0
,
x
0
→
R
n is a continuous initial function, then u
x
,
x
0
,
μ will denote the solution of (1) on x
0
,
∞
. Frequently, it is sufficient to write u
x. If f
i
0
=
u
i
0
=
0, then u
x
≡
0 is a solution of (1) called the zero solution. The norm on the initial function μ
x
=
μ
1
x
,
μ
2
x
,
⋯
,
μ
n
x is given by(62)
μ
=
sup
μ
x
=
∑
i
=
1
n
μ
i
:
0
≤
x
≤
x
0
.
The definition of stability of the zero solution is given in Burton [18] and is restated below.Definition 6 (see [22]).
The zero solution of (1) is stable if for each ϵ
>
0 and each x
0
≥
0, there exists δ such that μ
x
≤
δ on 0
,
x
0 and x
≥
x
0 imply u
x
,
x
0
,
μ
≤
ϵ.
We next define the statement “a Lyapunov functional for system (1)”. Let V
x
,
Ψ
. be defined for x
≥
0 and Ψ
∈
0
,
x
;
R
n and let V be locally Lipschitz in Ψ. For each x
≥
0 and every Ψ
∈
0
,
x
;
R
n, we define the derivative V along a solution of (1) by(63)
V
1
′
x
,
Ψ
.
=
lim
Δ
x
→
0
+
sup
V
x
+
Δ
x
,
u
·
,
x
,
Ψ
−
V
x
,
Ψ
.
Δ
x
,
whereu
ξ
;
x
,
Ψ is the unique solution of (1) with initial conditions x and Ψ. Then the following result by Driver ([28]) gives a definition of the Lyapunov functional.Theorem 4 (see [28]).
IfV
x
,
Ψ
. is defined for x
≥
0 and every Ψ
∈
0
,
x
;
R
n with(1)
V
x
,
0
≡
0.
(2)
V continuous in x and Lipschitz in Ψ.
(3)
V
x
,
Ψ
.
≥
W
Ψ
., where W
:
0
,
∞
→
0
,
∞ is a continuous function with(64)
W
0
=
0
,
W
r
>
0
,
r
>
0
andW strictly increasing (positive definiteness).
(4)
V
1
′
x
,
Ψ
.
≤
0.
then the zero solution of (1) is stable, and(65)
V
x
,
Ψ
.
=
V
x
,
Ψ
t
:
0
≤
t
≤
x
is called a Lyapunov functional for system (1).Finally, we assumed that the functions in (1) are well behaved, that continuous initial functions generate solutions, and that solutions which remain bounded can be continued.


## 5. Numerical Examples and Results

In this section, some numerical examples are presented to show the validity of the proposed methods. In addition, the numerical results are compared with exact solution.Example 1.
Consider the system of Volterra integro-differential equations:(66)
u
1
′
x
=
−
2
+
x
2
−
x
4
+
3
x
5
20
+
2
x
6
+
x
7
5
−
x
8
8
+
∫
0
x
t
3
−
x
2
u
1
+
12
t
2
−
x
u
2
d
t
u
2
′
x
=
4
−
8
x
−
x
4
−
x
3
3
+
2
x
4
−
8
x
5
5
+
x
6
30
−
4
e
x
+
∫
0
x
t
−
x
u
1
+
8
1
−
t
u
2
+
2
u
3
d
t
u
3
′
x
=
3
−
7
x
2
2
+
4
x
3
3
+
6
x
5
5
−
7
x
6
30
+
∫
0
x
2
x
−
t
u
1
+
6
t
u
2
+
u
3
d
t
together with the initial conditions(67)
u
1
0
=
0
,
u
2
0
=
1
,
u
3
0
=
2
.The exact solution of system (66) is(68)
u
1
x
=
x
4
−
2
x
,
u
2
x
=
1
−
x
3
,
u
3
x
=
x
+
2
e
x
.We start by implementing Algorithm1 to solve system (66) using the Sinc collocation method based on sinc functions.Algorithm 1: Numerical Realization Using the Sinc Collocation Method.
1:
Input:(i)
N
,
M
,
α
,
d
,
n
(ii)
f
i
x
f
o
r
i
=
1
,
2
,
…
,
N
(iii)
k
i
j
x
f
o
r
i
=
1
,
2
,
…
,
N
,
j
=
1
,
2
,
…
,
N
(iv)
Initial ConditionU
i
0
f
o
r
i
=
1
,
2
,
…
,
N
(v)
Initial ConditionU
i
d
0
f
o
r
i
=
1
,
2
,
…
,
N
,
d
=
1
,
2
,
…
,
n
−
1
2:
Define:(i)
Sinc FunctionS
k
,
h
,
x
(ii)
P
i
x
=
∑
j
=
0
n
a
i
j
x
j
f
o
r
i
=
1
,
2
,
…
,
n
(iii)
Y
i
x
=
∑
j
=
−
M
M
C
i
j
S
k
,
h
,
x
f
o
r
i
=
1
,
2
,
…
,
N
3:
Replaceu
i
x
=
Y
i
x
+
P
i
x
f
o
r
i
=
1
,
2
,
…
,
N, in each equation
4:
Calculate Sinc pointx
p
f
o
r
p
=
−
M
−
1
,
−
M
,
−
M
+
1
,
…
,
0
,
1
,
…
,
M
5:
Evaluateu
i
x
p
=
f
i
x
p
+
∫
0
x
p
∑
j
=
1
N
k
i
j
x
p
,
t
∗
u
j
t
d
t
f
o
r
p
=
−
M
−
1
,
−
M
,
−
M
+
1
,
…
,
0
,
1
,
…
,
M, for
i
=
1
,
2
,
…
,
N
From this step, we get2
M
+
1
∗
N equation
6:
Evaluateu
i
0
=
U
i
0
f
o
r
i
=
1
,
2
,
…
,
N. From this step, we get N equation
7:
Evaluateu
i
d
0
=
U
i
d
0
f
o
r
i
=
1
,
2
,
…
,
N
,
d
=
1
,
2
,
…
,
n
−
1. From this step, we get N
∗
n
−
1 equation
8:
Solving the algebraic system to getc
i
j and a
i
k
f
o
r
i
=
1
,
2
,
…
,
N
,
j
=
1
,
2
,
…
,
N
,
k
=
1
,
2
,
…
,
n
9:
Setu
a
p
r
r
o
x
i
x
=
sup
c
i
j and a
i
k in u
i
x
10:
Inputu
e
x
a
c
t
i
x
11:
Plotu
a
p
r
r
o
x
i
x
;
u
e
x
a
c
t
i
x
12:
Define the error|
u
e
x
a
c
t
i
x
−
u
a
p
r
r
o
x
i
x
|
13:
Plot the errorTables1 and 2 contain the exact and numerical solutions together with the resulting error with M
=
8.Table 1
The exact and numerical solutions of applying Algorithm1 for system (66) with M
=
8.
x
Exact solutionu
1
x
=
x
4
−
2
x
Numerical solutionu
1
a
p
p
Exact solutionu
2
x
=
1
−
x
3
Numerical solutionu
2
a
p
p
Exact solutionu
3
x
=
x
+
2
e
x
Numerical solutionu
1
a
p
p
0
0
0
1
1
2
2
0.1
−0.1999
−0.1998999986410037
0.999
0.9989999996164405
2.3103418361512955
2.310341836911406
0.2
−0.3984
−0.39840000469382025
0.992
0.992000001268595
2.64280551632034
2.642805513660997
0.3
−0.5919
−0.5919000795392859
0.973
0.9730000216099085
2.999717615152006
2.9997175706314225
0.4
−0.7744
−0.77440022803092
0.936
0.9360000615945837
3.3836493952825406
3.3836492678242998
0.5
−0.9375
−0.9375003156846542
0.875
0.8750000858865024
3.7974425414002564
3.7974423632300143
0.6
−1.0704
−1.0704001918522101
0.784
0.784000056734032
4.244237600781018
4.244237484652817
0.7
−1.1599
−1.1598998894031196
0.657
0.6569999842789827
4.727505414940953
4.72750545433764
0.8
−1.1904
−1.190399643757478
0.488
0.4879999282244343
5.251081856984936
5.251082019406921
0.9
−1.1439
−1.1438996142131692
0.271
0.2709999270259329
5.8192062223139
5.819206393371196Table 2
The resulting error for the numerical solution.
x
Absolute erroru
1
−
u
1
a
p
p
Absolute erroru
2
−
u
2
a
p
p
Absolute erroru
3
−
u
3
a
p
p
0
0
−1
0
0.1
1.358996248868038e
−
9
3.835595174805917e
−
10
7.601106410959346e
−
10
0.2
4.693820221390865e
−
9
1.268595006820305e
−
9
2.659342968058808e
−
9
0.3
7.953928582438152e
−
8
2.160990852928535e
−
8
4.452058366410938e
−
8
0.4
2.280309290281224e
−
7
6.159458376675531e
−
8
1.274582408505864e
−
7
0.5
3.156846541951807e
−
7
8.588650235452633e
−
8
1.781702421155273e
−
7
0.6
1.91852210118526e
−
7
5.673403202788307e
−
8
1.161282012773767e
−
7
0.7
1.10596880320557e
−
7
1.572101726576846e
−
8
3.939668680175146e
−
8
0.8
3.562425212599862e
−
7
7.177556560211684e
−
8
1.624219851947828e
−
7
0.9
3.857868307033385e
−
7
7.297406701134435e
−
8
1.710572963276035e
−
7Figure2 shows a comparison between the exact and numerical solutions for system (66). The maximum error corresponding to u
1
,
u
2, and u
3 is E
1
≈
3.9
e
−
7, E
2
≈
8.5
e
−
8, and E
3
≈
1.79
e
−
7
, respectively.Figure 2
A comparison between exact and numerical solutions by applying Algorithm1 for system (66) with M
=
8. (a) The exact and numerical solutions of u
1. (b) The exact and numerical solutions of u
2. (c) The exact and numerical solutions of u
3.
(a)
(b)
(c)Algorithm 2: Numerical Realization Using the Chebyshev Wavelets Method.
1:
Input:(i)
N
,
M
,
k
,
b
,
n
(ii)
f
i
x
f
o
r
i
=
1
,
2
,
…
,
N
(iii)
k
i
j
x
f
o
r
i
=
1
,
2
,
…
,
N
,
j
=
1
,
2
,
…
,
N
(iv)
Initial ConditionU
i
0
f
o
r
i
=
1
,
2
,
…
,
N
(v)
Initial ConditionU
i
r
0
f
o
r
i
=
1
,
2
,
…
,
N
,
r
=
1
,
2
,
…
,
n
−
1
2:
Define:(i)
Chebyshev FunctionT
k
,
b
,
m
,
x
(ii)
weight functionW
x
(iii)
D
i
x
=
d
i
1
,
d
i
2
,
…
,
d
i
M
f
o
r
i
=
1
,
2
,
…
,
N
(iv)
ψ
x
=
T
k
,
0
,
m
,
x
,
…
,
T
k
,
b
,
M
−
1
,
x
(v)
D
M
h
x as definition (4)
(vi)
operatorD
M
h
x
,
y as remark (5)
3:
CalculateF
i
=
D
M
F
i
x
,
f
o
r
i
=
1
,
2
,
…
,
N
4:
CalculateQ
i
=
D
M
U
i
0
,
f
o
r
i
=
1
,
2
,
…
,
N
5:
Calculatek
i
=
D
M
k
i
j
x
,
f
o
r
i
=
1
,
2
,
…
,
N
,
j
=
1
,
2
,
…
,
N
6:
Define operation matrixB
7:
Defineu
i
n
x
=
D
i
T
·
ψ
x
,
f
o
r
i
=
1
,
2
,
…
,
N
8:
Defineu
i
r
x
=
D
i
·
B
n
−
r
·
ψ
x
+
∑
j
=
0
n
−
r
−
1
a
i
r
x
j
/
j
!, for i
=
1
,
2
,
…
,
N
,
r
=
1
,
2
,
…
,
n
−
1
9:
Substitutingu
i
r
x
,
f
o
r
i
=
1
,
2
,
…
,
N
f
o
r
r
=
1
,
2
,
…
,
n in the system
10:
Multiplying each equation byW
x
·
x
ψ
x
11:
Applying∫
0
1
·
d
x for all equations
12:
From this step, we getM
∗
N equation
13:
Solving the algebraic system to getd
i
j
,
f
o
r
i
=
1
,
2
,
…
,
N
,
j
=
1
,
2
,
…
,
N
14:
Setu
a
p
r
r
o
x
i
x
=
sup
d
i
j in u
i
x
15:
Inputu
e
x
a
c
t
i
x
16:
Plotu
a
p
r
r
o
x
i
x
;
u
e
x
a
c
t
i
x
17:
Define the error|
u
e
x
a
c
t
i
x
−
u
a
p
r
r
o
x
i
x
|
18:
Plot the errorNext, we implement Algorithm2 to solve system (66) using the Chebyshev wavelets method. Tables 3 and 4 contain the exact and numerical solutions for system (66) together with the resulting error with M
=
8.Table 3
The exact and numerical solutions of applying Algorithm2 for system (66).
x
Exact solutionu
1
x
=
x
4
−
2
x
Numerical solutionu
1
a
p
p
Exact solutionu
2
x
=
1
−
x
3
Numerical solutionu
2
a
p
p
Exact solutionu
3
x
=
x
+
2
e
x
Numerical solutionu
1
a
p
p
0
0
2.690833666996184e
−
13
1
0.9999999999543171
2
1.9999999974820484
0.1
−0.1999
−0.19990000000022942
0.999
0.9990000000727636
2.3103418361512955
2.3103418351947798
0.2
−0.3984
−0.39839999999984066
0.992
0.991999999990548
2.64280551632034
2.6428055151638277
0.4
−0.7744
−0.7744000000009819
0.936
0.9360000000694738
3.3836493952825406
3.3836493955712172
0.5
−0.9375
−0.9374999999997702
0.875
0.875000000131186
3.7974425414002564
3.7974425389296056
0.6
−1.0704
−1.0703999999956941
0.784
0.7840000001191981
4.244237600781018
4.244237600772698
0.7
−1.1599
−1.1598999999883794
0.657
0.6570000000729469
4.727505414940953
4.727505417496077
0.8
−1.1904
−1.190399999977122
0.488
0.4880000000954549
5.251081856984936
5.2510818561222665
0.9
−1.1439
−1.1438999999585011
0.271
0.27100000021849296
5.8192062223139
5.819206221276965Table 4
The resulting error for the numerical solution.
x
Absolute erroru
1
−
u
1
a
p
p
Absolute erroru
2
−
u
2
a
p
p
Absolute erroru
3
−
u
3
a
p
p
0
2.690833666996184e
−
13
4.568290190576363e
−
11
2.517951624980696e
−
9
0.1
2.293998324631729e
−
13
7.276357294472291e
−
11
9.565157554902726e
−
10
0.2
1.593725151849412e
−
13
9.451994742448733e
−
12
1.156512219324668e
−
9
0.3
4.383160501220118e
−
13
9.397482791939638e
−
12
2.494912720862885e
−
9
0.4
9.818812429784884e
−
13
6.947387110045611e
−
11
2.886766381493544e
−
10
0.5
2.298161660974074e
−
13
1.311859509911528e
−
10
2.470650795061146e
−
9
0.6
4.305888978706207e
−
12
1.191982068604602e
−
10
8.319567257331073e
−
12
0.7
1.162048235414658e
−
11
7.294698178839099e
−
11
2.555123224112776e
−
9
0.8
2.287792177924075e
−
11
9.545503276697787e
−
11
8.62669047307918e
−
10
0.9
4.14988043928588e
−
11
2.184930569804066e
−
10
1.0369349823236e
−
9The maximum error corresponding tou
1
,
u
2 and u
3 is E
1
≈
7.7
e
−
11, E
2
≈
2.5
e
−
10, and E
3
≈
2.6
e
−
9, respectively. Figure 3 displays a comparison between the exact and numerical solutions for systems (66) using Algorithm 2.Figure 3
A comparison between exact and numerical solutions by applying Algorithm2 for system (66) with M
=
8. (a) The exact and numerical solutions of u
1. (b) The exact and numerical solutions of u
2. (c) The exact and numerical solutions of u
3.
(a)
(b)
(c)Example 2.
Consider the system of Volterra integro-differential equations:(69)
u
1
″
x
=
−
1
−
x
+
cosh
x
−
sin
3
x
3
−
sinh
x
+
e
x
+
∫
0
x
e
−
t
u
1
+
sin
2
t
u
2
d
t
u
2
″
x
=
−
3
+
x
2
−
2
x
3
3
−
2
e
x
x
−
1
+
∫
0
x
x
2
−
t
2
u
1
+
x
−
t
u
2
d
t
together with the initial conditions(70)
u
1
0
=
2
,
u
1
′
0
=
1
,
u
2
0
=
1
,
u
2
′
0
=
0
.
The exact solution of system (69) is(71)
u
1
x
=
e
x
+
1
,
u
2
x
=
cos
x
.We implement Algorithm1 to solve system (69) using the sinc collocation method based on sinc functions. Table 5 contains the exact and numerical solutions using Algorithm 1 for system (69) together with the resulting error with M
=
8.Table 5
The exact and numerical solutions of applying Algorithm1 for system (69).
x
Exact solutionu
1
x
=
e
x
+
1
Numerical solutionu
1
a
p
p
Absolute erroru
1
−
u
1
a
p
p
Exact solutionu
2
x
=
cos
x
Numerical solutionu
2
a
p
p
Absolute erroru
2
−
u
2
a
p
p
0
2
2
0
1
1
0
0.1
2.1051709180756477
2.1051397636262448
0.00003115444940293699
0.9950041652780258
0.9949953813521188
0.00000878392590697441
0.2
2.2214027581601696
2.2212621771873557
0.00014058097281388626
0.9800665778412416
0.9800269511507578
0.00003962669048385159
0.3
2.349858807576003
2.353150577162296
0.0032917695862928475
0.955336489125606
0.9562638272781463
0.0009273381525403135
0.4
2.4918246976412703
2.5101191909529916
0.018294493311721283
0.9210609940028851
0.9262143139859589
0.005153319983073779
0.5
2.648721270700128
2.694898115680843
0.046176844980714726
0.8775825618903728
0.890587846895059
0.013005285004686229
0.6
2.822118800390509
2.8955876162197285
0.07346881582921938
0.8253356149096782
0.8460218395509642
0.020686224641286
0.7
3.0137527074704766
3.094820636647796
0.08106792917731953
0.7648421872844884
0.7876576255740133
0.022815438289524925
0.8
3.225540928492468
3.289093529586907
0.06355260109443917
0.6967067093471654
0.7145821371224699
0.017875427775304487
0.9
3.45960311115695
3.4945194299151083
0.03491631875815848
0.6216099682706644
0.631444245918761
0.009834277648096634Figure4 compares the exact and numerical solutions u
1
x
=
e
x
+
1 and the approximate solution with M
=
8. The maximum error corresponding to u
1 and u
2 is E
1
≈
0.081 and E
2
≈
0.229.Figure 4
The exact and numerical solutions ofu
1 using Algorithm 1 for system (69).Figure5 compares the exact solution u
2
x
=
cos
x and the approximate solution with M
=
8.Figure 5
The exact and numerical solutions ofu
2 using Algorithm 1 for system (69).Next, we implement Algorithm2 to solve system (69) using the Chebyshev wavelets method. Table 6 contains the exact and numerical solutions with M
=
8.Table 6
The exact and numerical solutions of applying Algorithm2 for system (69).
x
Exact solutionu
1
x
=
e
x
+
1
Numerical solutionu
1
a
p
p
Absolute erroru
1
−
u
1
a
p
p
Exact solutionu
2
x
=
cos
x
Numerical solutionu
2
a
p
p
Absolute erroru
2
−
u
2
a
p
p
0
2
1.999999998740645
1.25935506467556e
−
9
1
0.9999999993413167
6.586833301014394e
−
10
0.1
2.1051709180756477
2.105170917580036
4.956115517984472e
−
10
0.9950041652780258
0.9950041649575281
3.204977394588582e
−
10
0.2
2.2214027581601696
2.22140275754089
6.192797386006532e
−
10
0.9800665778412416
0.9800665775316376
3.096040090966312e
−
10
0.3
2.349858807576003
2.349858808759274
1.183270814664183e
−
9
0.955336489125606
0.9553364897343352
6.087291781753379e
−
10
0.4
2.4918246976412703
2.491824697676863
3.559286199106282e
−
11
0.9210609940028851
0.9210609939384683
6.441680522328852e
−
10
0.6
2.822118800390509
2.822118800114823
2.756861405828203e
−
10
0.8253356149096782
0.8253356147660844
1.435938035143635e
−
10
0.7
3.0137527074704766
3.01375270838393
9.134533129895317e
−
10
0.7648421872844884
0.764842187697923
4.134346198725325e
−
10
0.8
3.225540928492468
3.2255409275988933
8.935745476890133e
−
10
0.6967067093471654
0.6967067087823426
5.64822744131277e
−
10
0.9
3.45960311115695
3.459603110051884
1.105065816631167e
−
9
0.6216099682706644
0.6216099676731845
5.974798433783235e
−
10Figures6 and 7 compare the exact solutions using Chebyshev wavlet method with M
=
8. The maximum error corresponding to u
1 and u
2 is E
1
≈
1.7
e
−
9 and E
2
≈
9.57
e
−
9 respectively.Figure 6
The exact and numerical solutions ofu
1 using Algorithm 2 for system (69).Figure7 compares the exact solution u
2
x
=
cos
x and the approximate solution with M
=
8.Figure 7
The exact and numerical solutions ofu
2 using Algorithm 2 for system (69).


## 6. Conclusions

In this article, a collocation method using sinc functions and Chebyshev wavelets method is proposed to solve linear system of integro-differential equations. The numerical results show that the convergence and accuracy of both methods were in good agreement with the analytical solution. Comparison of numerical results mentioned in tables and figures shows clearly that the Chebyshev wavelets method provides more accurate results and is therefore more effective than any other methods for solving systems of integro-differential equations.



---

*Source: 2360487-2020-02-04.xml*