**Preface to "Differential/Difference Equations: Mathematical Modeling, Oscillation and Applications"**

The study of oscillatory phenomena is an important part of the theory of differential equations. Oscillations naturally occur in virtually every area of applied science including, e.g., mechanics, electrical, radio engineering, and vibrotechnics.

This Special Issue includes 19 high-quality papers with original research results in theoretical research, and recent progress in the study of applied problems in science and technology.

This Special Issue brought together mathematicians with physicists, engineers, as well as other scientists. Topics covered in this issue:

Oscillation theory; Differential/difference equations; Partial differential equations; Dynamical systems; Fractional calculus; Delays; Mathematical modeling and oscillations.

> **Ioannis Dassios, Omar Bazighifan, Osama Moaaz** *Editors*

### **On the Approximate Solution of Partial Integro-Differential Equations Using the Pseudospectral Method Based on Chebyshev Cardinal Functions**

**Fairouz Tchier 1,\* ,† , Ioannis Dassios 2,†, Ferdous Tawfiq 3,† and Lakhdar Ragoub 4,†**


**Abstract:** In this paper, we apply the pseudospectral method based on the Chebyshev cardinal function to solve the parabolic partial integro-differential equations (PIDEs). Since these equations play a key role in mathematics, physics, and engineering, finding an appropriate solution is important. We use an efficient method to solve PIDEs, especially for the integral part. Unlike when using Chebyshev functions, when using Chebyshev cardinal functions it is no longer necessary to integrate to find expansion coefficients of a given function. This reduces the computation. The convergence analysis is investigated and some numerical examples guarantee our theoretical results. We compare the presented method with others. The results confirm the efficiency and accuracy of the method.

**Keywords:** interpolating scaling functions; hyperbolic equation; Galerkin method

#### **1. Introduction**

In this paper, we apply the pseudospectral method based on Chebyshev cardinal functions to solve one-dimensional partial integro-differential equations (PIDEs)

$$aw\_l(\mathbf{x},t) + aw\_{\mathbf{x}\mathbf{x}}(\mathbf{x},t) = \mathfrak{F} \int\_0^t k(\mathbf{x},t,\mathbf{s}, w(\mathbf{x},\mathbf{s}))d\mathbf{s} + f(\mathbf{x},t), \quad \mathbf{x} \in [a,b], \quad t \in [0,T], \tag{1}$$

with initial and boundary conditions

$$w(\mathbf{x},0) = \mathbf{g}(\mathbf{x}), \quad \mathbf{x} \in [a,b]. \tag{2}$$

$$w(0,t) = h\_0(t), \qquad w(1,t) = h\_1(t), \quad t \in [0,T]. \tag{3}$$

where *α* and *β* are constants and the functions *f*(*x*, *t*) and *k*(*x*, *t*,*s*, *w*) are assumed to be sufficiently smooth on D := [0, 1] × [0, *T*] and S with S := {(*x*, *t*,*s*) : *x* ∈ [0, 1],*s*, *t* ∈ [0, *T*]}, respectively, as prescribed before and such that (1) has a unique solution *w*(*x*, *t*) ∈ *C*(*D*). In addition, we assume that the kernel function is of diffusion type which is given by

$$k(\mathbf{x}, t, \mathbf{s}, w(\mathbf{x}, \mathbf{s})) := k\_1(\mathbf{x}, t - s)w(\mathbf{x}, \mathbf{s}), \tag{4}$$

and satisfies the Lipschitz condition as follows

$$|k(\mathbf{x}, t, \mathbf{s}, w(\mathbf{x}, \mathbf{s})) - k(\mathbf{x}, t, \mathbf{s}, v(\mathbf{x}, \mathbf{s})) \le \mathcal{A} |w(\mathbf{x}, \mathbf{s}) - v(\mathbf{x}, \mathbf{s})|,\tag{5}$$

where A ≥ 0 is referred to as a Lipschitz constant.

**Citation:** Tchier, F.; Dassios, I.; Tawfiq, F.; Ragoub, L. On the Approximate Solution of Partial Integro-Differential Equations Using the Pseudospectral Method Based on Chebyshev Cardinal Functions. *Mathematics* **2021**, *9*, 286. https://dx.doi.org/10.3390/ math9030286

Academic Editor: Konstantinas Pileckas Received: 4 January 2021 Accepted: 28 January 2021 Published: 1 February 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

In various fields of physics and engineering, systems are often functions of space and time and are described by partial differential equations. But in some cases, such a formulation can not accurately model this system. Because we can not take into account the effect of a past time when the system is a function of a given time. Such systems appear in heat transfer, thermoelasticity and nuclear reactor dynamics. This phenomenon has resulted in the inclusion of an integral term in the basic partial differential equation that leads to a PIDEs [1]. The existence, uniqueness, and asymptotic behavior of the solution of this equation are discussed in [2]. In this paper, we can find the physical situation that leads to Equation (1). A Simple example that refers to a PIDEs is considered by Habetler and Schiffman [3] where the compression of viscoelastic media is studied. For more applications, we refer readers to [4–7].

Spectral methods are schemes to discretize the PDEs. To this end, they utilize the polynomials to approximate the exact solution. Since any analytic function can be exponentially approximated by polynomials. In contrast to other methods such as finite elements and finite differences, these methods can achieve an infinite degree of accuracy. That's mean the order of the convergence of the approximate solution is limited only by the regularity of the exact solution. In other words, for numerical simulations, fewer degrees of freedom are necessary to obtain a given accuracy. The Galerkin method is a class of spectral techniques that convert a continuous operator problem to a discrete problem. In other words, this scheme applies the method of variation of parameters to function space by transforming the equation to a weak formulation. To implement this method, we can not compute the integrals analytically. That's why we can't use this method in most cases [8,9]. Another method that is closely related to spectral methods is the pseudospectral method. The pseudospectral methods are a special type of numerical method that used scientific computing and applied mathematics to solve partial differential equations. These methods allow the representation of functions on a quadrature grid and cause simplification of the calculations [10,11].

Several techniques have been used to solve one-dimensional partial differential equations, such as the finite difference method, finite element method, and spectral method. In [12], the Legendre-collocation method is used to solve the parabolic Volterra integrodifferential equation. For an infinite domain, Dehghan et al. [12] used the algebraic mapping to obtain a finite domain and then they utilized their proposed method. The Legendre multiwavelets collocation method is used to find the numerical solution of PIDEs [13]. To find the approximate solution of PIDEs, Avazzadeh et al. [14] applied the radial basis functions (RBFs) and finite difference method (FDM). To solve nonlinear parabolic PIDEs in one space variable, Douglas and Jones [15] proposed backward difference and Crank-Nicolson type methods. Han et al. [16] approximated the solution of (1) with kernel function of diffusion type and on unbounded spatial domains using artificial boundary method. In [17], a finite difference scheme is considered to solve PIDEs with a weakly singular kernel.

According to the above, considerable attention has been devoted to solving PIDEs numerically. In this paper, we introduce a simple numerical method with high accuracy. To this end, while introducing the Chebyshev cardinal functions, the pseudospectral method applies to obtain the approximate solution of PIDEs (1). Generally, cardinal functions {*Ci*} are polynomials of a given degree that *C<sup>i</sup>* vanishes at all interpolation grids except *x<sup>i</sup>* . These bases are also called the shape functions, Lagrange basis, and so on. One of the advantages of using such bases is the reduction of calculations to find the expansion coefficients of a given function. In other words, to find the expansion coefficients based on these bases, there is no need to integrate, and this is due to the cardinality, which makes these bases superior to other functions. Laksetani and Dehghan [18] is used Chebyshev cardinal functions to solve a PDE with an unknown time-dependent coefficient. In [19], these functions are used to solve the fractional differential equation. Heydari [20] described a new direct scheme for solving variable-order fractional optimal control problem via

Chebyshev cardinal functions. For more details about the Chebyshev cardinal functions and their applications, we refer the reader to [21,22].

This paper is organized as follows, Section 2 is devoted to a brief introduction to Chebyshev cardinal functions. In Section 3, we presented an efficient and applicable method based on Chebyshev cardinal functions to solve PIDEs (1). In Section 4, the convergence analysis is investigated and we proved that the proposed method is convergence. Section 5 is devoted to some numerical tests to show the ability ad accuracy of the method. Finally, Section 6 contains a few concluding remarks.

#### **2. Chebyshev Cardinal Functions**

Given *<sup>M</sup>* <sup>∈</sup> <sup>N</sup>, assume that <sup>M</sup> :<sup>=</sup> {1, 2, . . . , *<sup>M</sup>* <sup>+</sup> <sup>1</sup>} and <sup>X</sup> :<sup>=</sup> {*x<sup>i</sup>* : *TM*+1(*xi*) = 0, *i* ∈ M} where *TM*+<sup>1</sup> is the first kind Chebyshev function of order *M* + 1 on [−1, 1]. Recall that the Chebyshev grid is obtained by

$$\mu\_i := \cos\left(\frac{(2i - 1)\pi}{2M + 2}\right), \quad \forall i \in \mathcal{M}. \tag{6}$$

To utilize the Chebyshev functions of any arbitrary interval [*a*, *b*], one can apply the change the variable *x* = 2(*t*−*a*) *<sup>b</sup>*−*<sup>a</sup>* <sup>−</sup> <sup>1</sup> to obtain the shifted Chebyshev functions, viz

$$T\_{M+1}^\*(t) := T\_{M+1} \left( \frac{2(t-a)}{b-a} - 1 \right), \quad t \in [a,b]. \tag{7}$$

Note that it is easy to show that the grids of shifted Chebyshev function *T* ∗ *M*+1 is equal to *t<sup>i</sup>* = (*x*+1)(*b*−*a*) <sup>2</sup> + *a*.

A significant example of the cardinal functions for orthogonal polynomials is the Chebyshev cardinal functions. The cardinal Chebyshev functions of order *M* + 1 are defined as

$$\mathcal{C}\_{i}(\mathbf{x}) = \frac{T\_{M+1}(\mathbf{x})}{T\_{M+1,\mathbf{x}}(\mathbf{x}\_{i})(\mathbf{x} - \mathbf{x}\_{i})}, \quad i \in \mathcal{M}, \tag{8}$$

where the subscript *x* denotes *x*-differentiation. It is obvious that the functions *Ci*(*x*) are polynomials of degree *M* which satisfy the condition

$$\mathsf{C}\_{i}(\mathsf{x}\_{l}) = \delta\_{\mathsf{il}} \tag{9}$$

where *δil* is the Kronecker *δ*-function.

In view of (9), the cardinal functions are nonzero at one and only one of the points *x<sup>i</sup>* ∈ X implies that for arbitrary function *p*(*t*), the function can be approximated by

$$p(t) \approx \sum\_{i=1}^{M+1} p(t\_i)\mathbb{C}\_i(t). \tag{10}$$

Assume that *H<sup>n</sup>* ([*a*, *<sup>b</sup>*]), *<sup>n</sup>* <sup>∈</sup> <sup>N</sup> (Sobolev spaces) denotes the space of all functions *p* ∈ *C n* ([*a*, *<sup>b</sup>*]) such that *<sup>D</sup><sup>α</sup> <sup>p</sup>* <sup>∈</sup> *<sup>L</sup>* 2 ([*a*, *b*]) for all *α* ≤ *n*, where *α* is a nonnegative integer and *D* is the derivative operator. Sobolov space *H<sup>n</sup>* ([*a*, *b*]) is equipped with a norm defined by

$$\|\|p\|\|\_{H^n([a,b])}^2 = \sum\_{l=0}^n \|\|p^{(l)}(t)\|\|\_{L^2([a,b])}^2. \tag{11}$$

There exista a semi-norm that is defined as follows

$$\left\| p \right\|\_{H^{n,M}([a,b])}^2 = \sum\_{l=\min n,M}^M \left\| p^{(l)}(t) \right\|\_{L^2([a,b])}^2. \tag{12}$$

It follows from [23] that the error of expansion (10) can be bounded by the following lemma.

**Lemma 1.** *Let* {*ti*}*i*∈M ∈ X <sup>∗</sup> *denotes shifted Gauss-Chebyshev points where* X <sup>∗</sup> := {*t<sup>i</sup>* : *T* ∗ *M*+1 (*ti*) = 0, *<sup>i</sup>* ∈ M} *and that p*(*t*) ∈ H*<sup>n</sup>* ([*a*, *b*]) *can be approximated by p<sup>M</sup> via*

$$p\_M(t) = \sum\_{i=1}^{M+1} p(t\_i) \mathcal{C}\_i(t).$$

*Then one can prove that*

$$\|\|p - p\_M\|\|\_{L^2([a, b])} \le \mathcal{C}M^{-n} |p|\_{H^{n, M}([a, b])}\tag{13}$$

*where C is a constant and independent of M.*

#### **3. Pseudospectral Method**

In this section, we apply the pseudospectral method to solve PIDEs (1) based on Chebyshev cardinal functions. Let us consider the partial integro-differential Equation (1) on the region <sup>Ω</sup> × T. We introduce differential operator

$$\mathcal{L} := \frac{\partial}{\partial t} + \mathfrak{a} \frac{\partial^2}{\partial \mathfrak{x}^2} \,. \tag{14}$$

and integral operator

$$\mathcal{T} := \beta \int\_0^t k(\mathbf{x}, t, \mathbf{s}, \cdot) ds. \tag{15}$$

Applying these operators, PIDEs (1) can be rewritten in the operator form

$$(\mathcal{L} + \mathcal{L})(w) = f. \tag{16}$$

Let the solution of (1) is approximated by the polynomial *w*˜(*x*, *t*), via

$$\mathfrak{w}(\mathbf{x},t) = \sum\_{i=1}^{M+1} \sum\_{j=1}^{M+1} w^n(t\_i, t\_j) \mathbb{C}\_i(\mathbf{x}) \mathbb{C}\_j(t). \tag{17}$$

If we define a matrix *W* of dimension (*M* + 1) × (*M* + 1) whose (*i*, *j*)-th element is *w*(*t<sup>i</sup>* , *tj*), then Equation (17) becomes the matrix problem

$$\mathfrak{w}(\mathbf{x},t) = \mathcal{C}^{T}(\mathbf{x})\mathcal{W}\mathcal{C}(t),\tag{18}$$

where the vector elements of C(*x*) are the Chebyshev cardinal functions {*Ci*(*x*)}.

Inasmuch as the Chebyshev cardinal functions are polynomial, it is easy to evaluate their derivatives. In view of (17), one can write

$$\mathfrak{w}\_{\mathbf{x}}(\mathbf{x},t) = \sum\_{i=1}^{M+1} \sum\_{i=1}^{M+1} w(t\_i, t\_j) \mathbb{C}\_{i, \mathbf{x}}(\mathbf{x}) \mathbb{C}\_j(t) = \mathcal{C}\_{\mathbf{x}}^T(\mathbf{x}) \mathsf{W} \mathcal{C}(t), \tag{19}$$

where C*x*(*x*) is a vector of dimension (*M* + 1) whose *i*-th element is *Ci*,*x*(*x*). Similarly we have

$$
\tilde{w}\_l(\mathbf{x}, t) = \sum\_{i=1}^{M+1} \sum\_{i=1}^{M+1} w(t\_i, t\_j) \mathbb{C}\_{i, \mathbf{x}}(\mathbf{x}) \mathbb{C}\_l(t) = \mathcal{C}^T(\mathbf{x}) \mathsf{W} \mathcal{C}\_l(t), \tag{20}
$$

where C*t*(*t*) is a vector of dimension (*M* + 1) whose *i*-th element is *Ci*,*t*(*t*). Suppose that D ∈ <sup>R</sup>*M*+1,*M*+<sup>1</sup> is the operational matrix of derivative whose (*i*, *j*)-th element is D*i*,*<sup>j</sup>* = *Ci*,*t*(*tj*). Thus, it follows from C*x*(*x*) = DC(*x*) that

$$\tilde{w}\_{\mathbf{x}}(\mathbf{x},t) = \mathcal{C}^{T}(\mathbf{x})\mathcal{D}^{T}\mathcal{W}\mathcal{C}(t),\tag{21}$$

and

$$
\mathfrak{w}\_t(\mathfrak{x}, t) = \mathcal{C}^T(\mathfrak{x}) \mathcal{W} \mathcal{D} \mathcal{C}(t). \tag{22}
$$

It can easily be shown that *w*˜ *xx*(*x*, *t*) is approximated as follows

$$\mathfrak{w}\_{\text{xx}}(\mathfrak{x},t) = \mathcal{C}^{T}(\mathfrak{x}) \mathcal{D}^{T^{2}} \mathcal{W} \mathcal{C}(t). \tag{23}$$

Thus, by substituting (22) and (23) into the differential part of desired Equation (16), we can approximate the differential operator L (14), via

$$\mathcal{L}(w)(\mathbf{x}, t) \approx \mathcal{C}^{T}(\mathbf{x}) \mathsf{W} \mathcal{D} \mathcal{C}(t) + \mathsf{u} \mathcal{C}^{T}(\mathbf{x}) \mathsf{D}^{T^{2}} \mathsf{W} \mathcal{C}(t), \tag{24}$$

To approximate the integral part, we assume that

$$\int\_{0}^{t} \mathcal{C}(\mathbf{x}) d\mathbf{x} = I \mathcal{C}(t), \tag{25}$$

where *<sup>I</sup>* <sup>∈</sup> <sup>R</sup>*M*+1,*M*+<sup>1</sup> is the operational matrix of integral. It follows from (15) that

$$\mathcal{L}(w)(\mathbf{x},t) = \beta \int\_0^t k(\mathbf{x}, t, \mathbf{s}, w(\mathbf{x}, \mathbf{s}))d\mathbf{s}.\tag{26}$$

If we replace *w* with *w*˜, then one can write

$$\mathcal{Z}(w)(\mathbf{x},t) \approx \beta \int\_0^t k(\mathbf{x},t,\mathbf{s},\mathfrak{w}(\mathbf{x},\mathbf{s}))d\mathbf{s}.\tag{27}$$

Assume that *k*(*x*, *t*,*s*, *w*˜(*x*,*s*)) can be approximated by C *T* (*x*)*K*C(*t*) where *K* is a matrix whose elements depend on *t* and unknown coefficients *W*. Replacing C *T* (*x*)*K*C(*t*) into (27), and using the operational matrix of integration *I*, we get

$$\begin{split} \mathcal{I}(w)(\mathbf{x},t) &\approx \beta \int\_0^t \mathcal{C}^T(\mathbf{x}) K \mathcal{C}(s) ds \\ &= \beta \mathcal{C}^T(\mathbf{x}) K \int\_0^t \mathcal{C}(s) ds \\ &= \beta \mathcal{C}^T(\mathbf{x}) K \mathcal{C}(t) \\ &= q(\mathbf{x},t) = \mathcal{C}^T(\mathbf{x}) Q \mathcal{C}(t). \end{split} \tag{28}$$

where (*i*, *j*)-th element of matrix *Q* is *q*(*t<sup>i</sup>* , *tj*). Substituting (25) and (28) into (16), one can write

$$\mathcal{C}^{T}(\mathbf{x})(\mathcal{W}\mathcal{D} + \mathcal{a}\mathcal{D}^{T^{2}}\mathcal{W} + \mathcal{Q})\mathcal{C}(t) = \mathcal{C}^{T}(\mathbf{x})\mathcal{F}\mathcal{C}(t). \tag{29}$$

The Chebyshev cardinal functions {*Ci*(*x*)} are orthogonal with respect to weighted inner product on [−1, 1]

$$\langle \mathbf{C}\_{i}(\mathbf{x}), \mathbf{C}\_{j}(\mathbf{x}) \rangle\_{\omega(\mathbf{x})} = \begin{cases} \frac{\pi}{M+1} & i = j, \\ 0, & i \neq j. \end{cases}$$

where *<sup>ω</sup>*(*x*) = 1/<sup>√</sup> 1 − *x* 2 . This gives rise to equation

$$
\alpha \mathcal{W} \mathcal{D} + \mathcal{a} \mathcal{D}^{\mathsf{T}^2} \mathcal{W} + Q = F. \tag{30}
$$

Let us rewrite this system as

$$\mathcal{F}(\mathcal{W}) := \mathcal{W}\mathcal{D} + \mathfrak{a}\mathcal{D}^{T^2}\mathcal{W} + \mathcal{Q} - F = 0. \tag{31}$$

We Replace the first column of (31) with the initial condition (2) and the first and last rows of (31) with the boundary conditions (3), i.e.,

$$\begin{aligned} [\mathcal{F}(\mathcal{W})]\_{i,1} &= [\mathcal{W}\mathcal{C}(\mathbf{0})]\_i - g(t\_i)\_{\prime} \\ [\mathcal{F}(\mathcal{W})]\_{1,i} &= [\mathcal{C}^T(\mathbf{0})\mathcal{W}]\_i - h\_0(t\_i)\_{\prime} \\ [\mathcal{F}(\mathcal{W})]\_{M+1,i} &= [\mathcal{C}^T(\mathbf{1})\mathcal{W}]\_i - h\_1(t\_i)\_{\prime} \\ i &= 1, \dots, M+1. \end{aligned}$$

Using the matrix to vector conversion, this system is changed to a new system by (*M* + 1) 2 equations with (*M* + 1) <sup>2</sup> unknowns

$$\begin{cases} \quad \bar{W}\Gamma = \mathfrak{F}\_{\prime} & \text{if } \mathbf{k} \text{ is a nonlinear function of } \mathbf{w}\_{\prime} \\ \quad \mathcal{F} = \mathfrak{F}\_{\prime} & \text{if } \mathbf{k} \text{ is a linear function of } \mathbf{w}\_{\prime} \end{cases} \tag{32}$$

where *<sup>W</sup>*¯ , <sup>F</sup>, and <sup>F</sup>¯ are obtained using the matrix to vector conversion of *<sup>W</sup>*, *<sup>F</sup>*, and <sup>F</sup> respectively.

After solving the linear or nonlinear system (32) using the generalized minimal residual method (GMRES) [24] and Newton-Raphson method, respectively, the unknowns *W* are found, and then the approximate solution can be obtained using (18).

#### **4. Convergence Analysis**

Because the function *f*(*x*, *t*) is a continuous function on *D*, the approximate error by comparing the function *f* with ˜ *f* may be bounded, established by the following theorem.

**Theorem 1.** *Let <sup>f</sup>* : *<sup>D</sup>* <sup>→</sup> <sup>R</sup><sup>2</sup> *be a sufficiently smooth function. Thus Chebyshev cardinal approximation to function f can be written as*

$$\|f - f\| \approx \mathcal{O}(2^{-2M}).\tag{33}$$

**Proof.** Let *PM*+1(*x*) denote that polynomial of degree *M* + 1 which interpolates to the function *f* at the *M* + 1 zeros of the first kind Chebyshev polynomials. It follows from [25] that

$$\begin{split} |f(\mathbf{x},t) - \mathcal{P}\_{M+1}(\mathbf{x},t)| &= \frac{\partial^{M+1}}{\partial \mathbf{x}^{M+1}} f(\boldsymbol{\xi},t) \frac{\Pi\_{i=1}^{M+1}(\mathbf{x}-t\_{i})}{(M+1)!} + \frac{\partial^{M+1}}{\partial t^{M+1}} f(\mathbf{x},\boldsymbol{\eta}) \frac{\Pi\_{j=1}^{M+1}(t-t\_{j})}{(M+1)!} \\ &- \frac{\partial^{2M+2}}{\partial \mathbf{x}^{M+1}t^{M+1}} f(\boldsymbol{\xi}',\boldsymbol{\eta}') \frac{\Pi\_{i=1}^{M+1}(\mathbf{x}-t\_{i})\Pi\_{j=1}^{M+1}(t-t\_{j})}{(M+1)!(M+1)!}. \end{split}$$

Since the leading coefficient of the first kind Chebyshev functions is 2*M*, and <sup>|</sup>*Ti*(*x*)| ≤ 1, ∀*i* ∈ M. It is possible to write

$$\begin{split} |f(\mathbf{x},t) - \mathcal{P}\_{M+1}(\mathbf{x},t)| &\leq \left(\frac{b-a}{2}\right)^{M+1} \frac{1}{2^{M}(M+1)!} \left(\sup\_{\tilde{\xi}\in[a,b]} |\frac{\partial^{M+1}}{\partial \mathbf{x}^{M+1}} f(\tilde{\xi},t)| + \sup\_{\eta\in[0,T]} |\frac{\partial^{M+1}}{\partial t^{\prime}} f(\mathbf{x},\eta)|\right) \\ &+ \left(\frac{b-a}{2}\right)^{2M+2} \frac{1}{4^{M}((M+1)!)^{2}} \sup\_{(\tilde{\xi}',\eta')\in D} |\frac{\partial^{2M+2}}{\partial \mathbf{x}^{\prime}\partial t^{M+1}} f(\tilde{\xi}',\eta')|. \end{split}$$

Since ˜ *f* is approximated by Chebyshev cardinal functions and these bases are polynomials, thus one can obtain

$$\begin{split} \|\|f - f\|\|^{2} &= \iint\_{D} |f(\mathbf{x}, t) - f(\mathbf{x}, t)|^{2} dt d\mathbf{x} \\ &\leq \iint\_{D} |f(\mathbf{x}, t) - P\_{\mathbf{M} + 1}(\mathbf{x}, t)|^{2} dt d\mathbf{x} \\ &\leq \iint\_{D} \left(\frac{b - a}{2}\right)^{2M + 1} \frac{1}{2^{M}(M + 1)!} \left(\sup\_{\xi \in [a, b]} |\frac{\partial^{M + 1} f(\xi, t)|}{\partial x^{M + 1}} f(\xi, t)| + \sup\_{\eta \in [0, T]} |\frac{\partial^{M + 1} f(\mathbf{x}, \eta)|}{\partial t'} f(\mathbf{x}, \eta)|\right) dt d\mathbf{x} \\ &+ \iint\_{D} \left(\frac{b - a}{2}\right)^{2M + 2} \frac{1}{4^{M}((M + 1)!)^{2}} \sup\_{\xi, \eta' \in D} |\frac{\partial^{2M + 2}}{\partial x^{\eta'} \partial t'} f(\xi', \eta')| dt d\mathbf{x} \\ &\leq 2^{-2M} \frac{(b - a)^{2M}}{(M + 1)!} C\_{\max}(1/2 + 2^{-2M - 2}/(M + 1)!) \iint\_{D} dt d\mathbf{x} \\ &\leq C\_{1} 2^{-2M} \end{split}$$
 
$$\text{where } \mathcal{C}\_{1} := \frac{(b - a)^{2M}}{(M + 1)!} C\_{\max}(1/2 + 2^{-2M - 2}/(M + 1)!) |D| \text{ and}$$

$$\mathcal{C}\_{\max} := \max \{ \sup\_{\substack{\mathfrak{T} \in [a,b]}} | \frac{\partial^{M+1}}{\partial \mathbf{x}^{M+1}} f(\mathfrak{f}, t) |\_{\prime} \sup\_{\eta \in [0, T]} | \frac{\partial^{M+1}}{\partial t^{r}} |\_{\prime} \sup\_{(\mathfrak{f}^{\prime}, \eta^{\prime}) \in D} | \frac{\partial^{2M+2}}{\partial \mathbf{x}^{r} \partial t^{M+1}} | \}.$$

**Proof.** Let *w*˜ denotes the approximate solution of (1) for which *e* = *w* − *w*˜. We subtract Equation (1) from

$$
\mathfrak{w}\_t(\mathbf{x},t) + a\mathfrak{w}\_{\mathbf{x}\mathbf{x}}(\mathbf{x},t) = \mathfrak{z} \int\_0^t k(\mathbf{x},t,\mathbf{s},\mathfrak{w}(\mathbf{x},\mathbf{s}))d\mathbf{s} + \mathfrak{f}(\mathbf{x},t),\tag{34}
$$

to obtain the following equation

$$e\_t(\mathbf{x},t) + ae\_{\mathbf{x}t}(\mathbf{x},t) = \beta \int\_0^t k(\mathbf{x},t,\mathbf{s},e(\mathbf{x},\mathbf{s}))d\mathbf{s} + f(\mathbf{x},t) - \tilde{f}(\mathbf{x},t). \tag{35}$$

Now, Assume that we can approximate the error function *e*(*x*, *t*) as follows

$$e(\mathbf{x}, t) \approx \mathcal{C}^{T}(\mathbf{x}) E \mathcal{C}(t),\tag{36}$$

where *E* is a matrix whose (*i*, *j*)-th element is *e*(*t<sup>i</sup>* , *tj*). Using this approximation and Lipschitz condition (5), Equation (35) may be written as

$$\mathcal{C}^{T}(\mathbf{x})\mathsf{E}\mathcal{D}\mathcal{C}(t) + \mathfrak{a}\mathcal{C}^{T}(\mathbf{x})\mathsf{D}^{T^{2}}\mathsf{E}\mathcal{C}(t) \leq \mathfrak{f}\mathcal{A}\mathcal{C}^{T}(\mathbf{x})\mathsf{E}\mathcal{C}(t) + \mathcal{C}^{T}(\mathbf{x})\eta\mathcal{C}(t),\tag{37}$$

where <sup>|</sup> *<sup>f</sup>* <sup>−</sup> ˜ *<sup>f</sup>* | ≈ C*<sup>T</sup>* (*x*)*η*C(*t*). By dropping the second term in the left to the other side of the inequality and taking norm from both sides, we have

$$\|\|E\mathcal{D}\|\| \le \mathcal{A} |\beta| \|\|E\|\| + |\alpha| \|\|\mathcal{D}^T E\|\| + \|\|\eta\|\|.\tag{38}$$

Because {*Ci*} are orthogonal functions, we removed kCk from both sides. Multiplying the right side of (38) by kDk, it follows that

$$\begin{aligned} \|\mathcal{ED}\| &\leq \mathcal{A} |\mathcal{B}| \|\mathcal{E}\| \|\mathcal{D}\| + |\mathfrak{a}| \|\mathcal{D}^{T^2} \mathcal{E}\| \|\mathcal{D}\| + \|\mathfrak{q}\| \|\mathcal{D}\| \\ &\leq \mathcal{A} |\mathcal{B}| \|\mathcal{E}\| \|\mathcal{I}\| \|\mathcal{D}\| + |\mathfrak{a}| \|\mathcal{D}^{T^2}\| \|\mathcal{E}\| \|\mathcal{D}\| + \|\mathfrak{q}\| \|\mathcal{D}\|\_{\prime} \end{aligned}$$

and then

$$\begin{aligned} &\|\|E\|\|\mathcal{D}\| \leq \mathcal{A}|\beta| \|\|E\|\| \|\mathcal{D}\| + |\mathfrak{a}| \|\|\mathcal{D}^{T^2} E\| \|\|\mathcal{D}\| + \|\|\eta\|\| \|\mathcal{D}\| \\ &\Rightarrow \|\|E\|\| \leq \mathcal{A}|\beta| \|\|E\|\| \|\|\mathcal{I}\| + |\mathfrak{a}| \|\|\mathcal{D}^{T^2}\| \|\|E\|\| + \|\eta\|\|. \end{aligned}$$

So, it is obvious that we shall have

$$\|\|E\|\|\mathbf{1} - \mathcal{A}|\beta||\|I\| - |\alpha| \|\mathcal{D}^2\|\| \le \|\|\eta\|\|. \tag{39}$$

Consequently, we obtain

$$\|E\| \le \left| 1 - \mathcal{A} |\mathcal{B}| \|I\| - |\mathfrak{a}| \|\mathcal{D}^2\| \right|^{-1} \|\eta\|. \tag{40}$$

If *<sup>f</sup>* be a sufficiently smooth function, then k*η*k → 0 as *<sup>M</sup>* → <sup>∞</sup>. Thus, we have

$$||e|| \to 0, \text{ as } \quad M \to \infty.$$

Therefore, the proposed method is convergent.

#### **5. Test Problems**

**Example 1.** *Let us dedicate the first example to the case that the desired Equation (1) is of form*

$$w\_t(\mathbf{x}, t) - w\_{\mathbf{x}\mathbf{x}}(\mathbf{x}, t) = f(\mathbf{x}, t) - \int\_0^t e^{\mathbf{x}(t-s)} w(\mathbf{x}, s) ds,$$

*with initial and boundary conditions*

$$\begin{aligned} w(\mathbf{x},0)&=0, & \mathbf{x} \in [0,1],\\ w(0,t)&=\sin(t), & w(1,t)&=0, & t \in [0,1].\end{aligned}$$

*and also f*(*x*, *t*) := (−*<sup>x</sup>* <sup>2</sup>+1)e *xt*+(*x* <sup>3</sup>+2 *x* <sup>2</sup>−*x*+2) sin(*t*)+(−*<sup>x</sup>* <sup>4</sup>+*x* 2 ) cos(*t*) *x* <sup>2</sup>+1 *. The exact solution for this example is given by [13] w*(*x*, *t*) = (1 − *x* 2 ) sin(*t*).

Table 1 shows a comparison between the proposed method and Legendre multiwavelets collocation method [13]. As you can see, our proposed method gives better results than [13]. According to Table 1, we can see that with fewer bases, we have achieved much better accuracy than the method in [13]. For different values of *M*, the errors in Table 2 are given with *L* <sup>∞</sup>, *L* <sup>2</sup> norms applying pseudospectral method based on Chebyshev cardinal functions. In Figure 1, the approximate solution, and absolute value of error are depicted.


**Table 1.** Comparison of the maximum absolute errors at different times for Example 1.

**Table 2.** The *L* <sup>∞</sup>, *L* 2 errors and CPU time for Example 1.


**Figure 1.** Plot of the approximate solution and absolute value of the error for Example 1.

**Example 2.** *Consider the following PIDEs [14]*

$$w\_t(\mathbf{x},t) + w\_{\text{xx}}(\mathbf{x},t) = \frac{\left(-\mathbf{x}^3 + (t^2+1)\mathbf{x}^2 - (t+1)^2\mathbf{x} + 2t\right)\mathbf{e}^{-\mathbf{x}t} + \mathbf{e}^{-t}\mathbf{x}}{\mathbf{x}-1} - \int\_0^t e^{\mathbf{s}-t} w(\mathbf{x}, \mathbf{s})d\mathbf{s}$$

*with initial and boundary conditions*

$$\begin{aligned} w(\mathbf{x},0) &= \mathbf{x}, & \mathbf{x} &\in [0,1], \\ w(0,t) &= 0, & w(1,t) &= e^{-t}, & t &\in [0,1]. \end{aligned}$$

*The exact solution for this example is w*(*x*, *t*) = *xe*−*xt .*

In Table 3, we report the *L* <sup>∞</sup>, *L* 2 errors and CPU time for different values of *M*. These results guarantee our convergence investigation in Section 4. When *M* increases, the error decreases, and approaches zero. The *L* <sup>∞</sup>, *L* 2 errors obtained by presented method are compared with Hermite-Taylor matrix method [26] and radial basis functions [14] in Table 4. According to Table 4, we can see that our presented method is better than

Hermite-Taylor matrix method [26] and radial basis functions [14]. Finally, we illustrate the approximate solution and absolute error in Figure 2.


**Table 3.** The *L* <sup>∞</sup>, *L* 2 errors and CPU time for Example 2.

**Table 4.** Comparison of the *L* <sup>∞</sup> and Ł<sup>2</sup> errors at different times for Example 2.


**Figure 2.** Plot of the approximate solution and absolute value of the error for Example 2.

**Example 3.** *To show the ability of the proposed method for solving nonlinear PIDEs (1), we consider the following equation.*

$$w\_t(\mathbf{x},t) + w\_{\mathbf{x}\mathbf{x}}(\mathbf{x},t) = \int\_0^t e^{\mathbf{x}+t+s} w^2(\mathbf{x},s) + f(\mathbf{x},t)\,\mu$$

*where*

$$f(\mathbf{x},t) = \frac{\left(\mathbf{x}\left(\left(\cos(t)\right)^2 + 2\cos(t)\sin(t) + 2\right)\mathbf{e}^{\mathbf{x}+2\cdot t} - 3\mathbf{e}^{\mathbf{x}+t}\mathbf{x} - 5\sin(t)\right)\mathbf{x}}{5},$$

*with the boundary and initial conditions*

$$\begin{aligned} w(x,0) &= x, & x &\in [0,1], \\ w(0,t) &= 0, & w(1,t) &= \cos(t), \quad t \in [0,1]. \end{aligned}$$

*The exact solution for this Example is given by w*(*x*, *t*) := *x* cos(*t*)*. Thus, we can easily judge the accuracy and convergency of the method.*

Figure 3 illustrates the log(*L* 2 *errors*), taking different values for *M*. To show the order of convergence, we also plotted the linear regression. The slope of this line is equal to the order of convergence (1.03248915355714). The numerical values with associated *L* 2 error and *L* <sup>∞</sup> error are tabulated in Table 5. Finally, we illustrate the approximate solution and absolute error, taking *M* = 8 in Figure 4.

**Table 5.** The *L* <sup>∞</sup> and *L* 2 errors for Example 3.


**Figure 3.** Plot of the log(*L* 2 *errors*) and the linear regression for Example 3.

**Figure 4.** Plot of the approximate solution and absolute value of the error for Example 3.

**Example 4.** *The last example is dedicated to equation*

$$w\_t(\mathbf{x}, t) - w\_{\mathbf{x}\mathbf{x}}(\mathbf{x}, t) = f(z, t) + \int\_0^t \mathfrak{Z}\mathbf{x}ste^{w(\mathbf{x}, s)}ds,$$

*where*

$$f(\mathbf{x},t) := \frac{-3t^2 \mathbf{x} \cos(\sin(\mathbf{x})t) \sin(\mathbf{x}) + 3\tan \mathbf{sin}(\sin(\mathbf{x})t) - \sin(\mathbf{x})(\cos(\mathbf{x}) - 1)(\cos(\mathbf{x}) + 1)(t+1)}{(\sin(\mathbf{x}))^2},$$

*and*

$$\begin{aligned} w(x,0)&=0, & x \in [0,1],\\ w(0,t)&=0, & w(1,t)&=\sin(1)t, & t \in [0,1].\end{aligned}$$

Since the closed form of the exact solution to the problem is unavailable, we compute a reference solution by picking a large *M* = 12. The *L* <sup>∞</sup>, *L* 2 errors, CPU time and order of convergence are tabulated in Table 6 for different values of *M*. Figure 5 illustrates the approximate solution and absolute error, taking *M* = 9. Table 7 shows the *L* <sup>∞</sup>, *L* 2 errors at the different times, taking different *M*.

**Table 6.** The *L* <sup>∞</sup>, *L* 2 errors, CPU time and order of convergence for Example 4.


**Table 7.** Comparison of the *L* <sup>∞</sup> and *L* 2 errors at different times for Example 4.


**Figure 5.** Plot of the approximate solution and absolute value of the error for Example 4.

#### **6. Conclusions**

In this paper, an efficient and novel numerical method is applied to solve partial integro-differential equations using the pseudospectral method based on Chebyshev cardinal functions. Due to the simplicity of using cardinal functions, the presented method is good for solving PIDEs. The convergence analysis is investigated and we can show when the number of bases increases, the accuracy is also increased. The presented method was applied to solve some numerical tests and the results guarantee our convergence investigation and application of the proposed method to this problem shows that it performs extremely well in terms of accuracy.

**Author Contributions:** Conceptualization, F.T. (Fairouz Tchier), and I.D.; methodology, software, F.T. (Fairouz Tchier) and F.T. ( Ferdous Tawfiq) and F.B.; validation, formal analysis, F.T. (Fairouz Tchier) and I.D. and L.R.; writing—original draft preparation, investigation, funding acquisition, F.T. (Fairouz Tchier) and F.T. ( Ferdous Tawfiq) and L.R.; writing—review and editing, F.T. (Fairouz Tchier) and F.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** The authors extend their appreciation to the Deanship of Scientific Research at King Saud University for funding this work through Research Group no RG-1440-010.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:


#### **References**


*Article*
