Numerical analysis of finite difference schemes arising from time-memory partial integro-differential equations

This paper investigates the partial integro-differential equation of memory type numerically. The differential operator is discretized based on θ-finite difference schemes, while the integral operator is approximated using Simpson's rule. The mesh points of an integral part are adapted to coincide with the nodes of the computational domain using the Heaviside function. The stability of the proposed numerical methods is established based on Gerschgoren's theorems. Also, its consistency is investigated to guarantee the numerical solutions' convergence. Several examples are provided to discuss the efficiency of the used finite difference schemes and compare them with previous studies.


. Introduction
The realm of differential equations is full of various categories of differential equations. When the differential equation has one independent variable is called an ordinary differential equation (ODE). In comparison, it has more than one called a partial differential equation (PDE). An equation contains an unknown function in the integral operator, known as an integral equation (IE). An integro-differential equation (IDE) is an equation that includes the derivatives and integration of the unknown function with one independent variable. For an unknown function with several independent variables, it is called a partial integro-differential equation (PIDE). Besides that, there is a classification based on the linearity of the unknown function. Partial differential equations (PDEs) are used to model several complicated phenomena in physics, chemistry, biology, and finance. They exhibit versatility in various models. However, there are so many phenomena in these fields that PDEs fail to model them. Because they need to consider their hierarchical status, we can split them into two main categories; the first requires considering its status in past time, and the second demand considering the accumulation behavior for one of its variables (spatial variable). The first category is represented by a convolution on the time variable by Yanik and Fairweather [1], while the second uses the convolution in the spatial variable [2].
Moreover, some phenomena use nonlocal operators in nonconvolution forms, such as options under the Lévy process and Bates models [3][4][5]. Here comes the importance of integrodifferential equations (IDEs) and partial integro-differential equations (PIDEs), which use the nonlocal operator to capture the behaviors of these phenomena. There are enormous physical, economic, and biological phenomena modeled based on these equations. For instance, neural IDEs are used to model brain dynamics, and it has been investigated based on the deep learning framework, see Zappala et al. [6]. The population dynamics are modeled based on nonlinear PIDEs and solved using the pseudo-spectral technique [2]. Various categories of IDEs are studied using different numerical techniques, see Mirzaee and Rafei [7], Mirzaee and Piroozfar [8], and Sakran [9].
Many researchers solved PIDEs using different mathematical techniques, which can be classified into two main categories: numerical [10][11][12] and spectral techniques [2,13,14]. The collocation method has been adapted to solve a broad class of PIDEs, see Yüzbaşı and Yıldırım [15] and Xu [11]. The Pell-Lucas collocation method is used to solve the parabolic PIDEs [15]. A space-time Sinc collocation technique is applied to obtain numerical solutions of a PIDE with a weak-singular integral-kernel [11]. In Kumar and Vijesh [16], the authors solved a two-spatial dimension PIDE of the memory type based on the Legendre wavelet iterative method. Also, the numerical methods are extended to solve fractional PIDEs [14,17]. For instance, Caputo-Fabrizio fractional Volterra PIDEs are solved using Legendre-Gauss-Lobatto quadrature incorporated with shifted Hahn polynomials, see Dehestani and Ordokhani [14]. Here, we consider a one spatial dimension parabolic partial differential equation with a source term and incorporated with a nonlocal integral part, given by subjected to initial and boundary conditions where, F(x, τ ) represents the source function, K(x, τ − s) is the kernel of integration, and ∪ Ŵ = [0, L]. The problem (Equations 1, 2) is used to describe the behavior of a continuous medium nuclear reactor such that the unknown function represents the total reactor power and also a wide class of models in thermoelasticity, viscoelasticity, refer to Lakshmikantham [18, Chap. 5]. Most mathematical techniques that are used to solve the problem (Equations 1, 2) treat the discretization of the unknowns in an integral part by applying additional expansion to make them coincide with the unknowns of the differential part consequently, increasing the computational time cost. Here, we introduce a mathematical treatment to overcome this issue.
The outline of the paper is as follows. Section 2 is dedicated to constructing the corresponding finite difference scheme of the problem (Equations 1, 2). The differential part is discretized using θ -finite difference approximations, while the integral part is discretized using Simpson's rule incorporated with the Heaviside step function. It is used to convert the variable upper bound of an integral part into a constant without adding any numerical computational cost, which is discussed in detail next section. To the best of our knowledge, we are the first authors who introduced this idea. The stability analysis of the proposed scheme is studied based on the second norm of matrices and Gerschgoren's theorems in Section 3. The consistency is discussed in Section 4. In Section 5, several examples are implemented to reveal the versatility and reliability of the proposed technique.
. The finite di erence scheme establishment We start this section with the construction of the numerical grid of the computational domain (x, τ ) ∈ [0, L] × [0, T], where T is the total elapsed time, the time variable τ is discretized using uniform mesh-points τ n = n τ , τ = T/N τ , and the spatial variable x is divided into subintervals of equal length using the spatial nodes x k = k x, x = L/N x , N τ , and N x is the total number of subintervals. Let V n k be the approximation of V at the point (x k , τ n ). The first-time partial derivative is approximated by and the second spatial derivative is discretized based on the θmethod, which incorporates a combination of explicit-implicit discretization, given by The θ -parameter is responsible for the ratio between explicit and implicit finite difference schemes, and its value θ ∈ [0, 1]. When θ = 0, then we have an explicit scheme, while it is a fully implicit scheme for θ = 1. For θ ∈ (0, 1), it represents a combination of explicit and implicit schemes. From Equations (3) and (4) into Equation (1), we have where, A and B are difference operators, given by To obtain the first time-level, we put n = 0 into Equation (5), then the lower and upper bounds become equal; consequently, the integral part vanishes, and the first time-level becomes The integral part of Equation (5) has two drawbacks; the first one is the upper bound of the integration changing as the index n changes, so, for each time level, we have a different integral interval. The second drawback appears when we try to apply a suitable substitution that converts the interval of variable length into an interval with a fixed size. It leads to the discretization of V k (s) does not coincide with the unknowns V n k at the constructing mesh points; in other words, it gives new unknowns other than V n k . To overcome the first drawback, we use the Heaviside step function H(τ − a), which is defined by consequently, the integral part of Equation (5) is rewritten as Next, we approximate the integration based on Simpson's rule Here, we discretize the interval of integration as the same as the discretization of the time variable, i.e., s = τ , and s j = τ j , j = 0, 1, . . . , N τ to guarantee the resultant unknowns V j k coincide with the original unknowns V n k , consequently, Equation (8) becomes The discretization of initial and boundary conditions is given by Hence, the finite difference problem becomes (10) Note that from Equation (10), the values of the unknown function V at time level n + 1 are obtained as a linear combination of all its values from 0 to n.
We mention Lax's equivalence theorem [19, p. 72] before commencing investigating the stability and consistency of the proposed finite difference scheme (Equation 10). Theorem 1. For a linear finite-difference scheme that yields from a discretization of a given well-posed linear initialvalue problem, when the numerical approximation satisfies the consistency condition and stability are the necessary and sufficient conditions for convergence.
In light of Lax's equivalence theorem, when stability and consistency are established for a linear finite difference scheme, the resultant approximations are guaranteed to converge to the solution of the given initial value problem regardless of whether the exact solution is known or not. Moreover, this fact is valid for the discretization of the integral part [20]. Based on that fact, we will investigate the stability and consistency of the proposed finite difference scheme (Equation 10).

. Stability analysis
To investigate the stability of the proposed finite difference scheme, we rewrite the problem (Equation 10) in the following matrix form Frontiers in Applied Mathematics and Statistics frontiersin.org . /fams. .
the vectors V n , F n ∈ R (N x −1)×1 are vectors of the unknowns and the discretization of the source term combined with the boundary conditions, given by and { (n,j) } n j=0 are a sequence of diagonal matrices such that To investigate the stability of Equation (11), first, we show that the inverse matrix of A exists. Based on the following theorem, we use a generic matrix M with entries like the entries of matrix A.
Theorem 2. For a tridiagonal symmetric matrix M ∈ R N×N and its entries are where α ∈ R + . We introduce the sequence {µ p } N p=1 , N ∈ N, such that µ p = det(M), M ∈ R p×p . Then the following properties hold 3. The ratio between any two successive elements is greater than 1 + α, i.e., µ p+1 Proof. For M ∈ R p×p , p = 1, 2, we have µ 1 = 1 + 2α, and For p > 2, we use the following recurrence relation [21] to obtain the determinant µ p Using the mathematical induction, we assume µ p µ p−1 > 1 + α holds. For the next ratio, we substitute p + 1 into Equation (16), we have µ p+1 Based on Theorem 2, the inverse of matrix A of the system (Equation 11) exists. Next, we recall some useful definitions and properties that relate the norm of symmetric tridiagonal and its inverse with its eigenvalues. First, we start with the eigenvalues (denoted bym) of a symmetric tridiagonal matrix M with entries given in Theorem 2 is given by Smith [19].m Since which means that all its eigenvalues are positive and greater than 1. Moreover, the maximum and minimum eigenvalues are obtained when j = 1, and N, respectively. Let us denote the maximum eigenvalue bym and the minimum bym, then The second norm of a matrix is defined by the square root of the spectral radius of the matrix M T M. Let us denote the spectral radius of a matrix M by ρ(M). Based on the properties of symmetric tridiagonal matrices [22,23], we have The stability analysis of scheme (Equation 11) is investigated based on Gerschgoren's theorems [19,24]. First, we rewrite the matrix B in the following form Frontiers in Applied Mathematics and Statistics frontiersin.org . /fams. . (19) into Equation (11), multiplying both sides by A −1 , and rearranging the elements, one gets

Next, by substitution Equation
where is the identity matrix. The necessary condition that guarantees the stability of the scheme (Equation 20) is that the eigenvalues of I − rA −1 C, and , must be less than or equal to 1.
The sequence { (n,j) } n j=0 are diagonal matrices with entries representing the discretization of the integral kernel. Here, we deal with differentiable kernels; consequently, they are continuous and bounded over closed intervals. Let K = max x,τ |K(x, τ )| then any eigenvalue of (n,j) lies in the interval (−2K( τ ) 2 , 2K( τ ) 2 ). For the first time level, we have by applying the limit as N x → +∞, yields while for θ ∈ [0.5, 1], the inequality (Equation 21) unconditionally holds. For all n ≥ 1, we have (23) When we say a given function is bounded, there is a constant such that all its absolute values over the given domain are less than that constant. However, this constant may be of the order less than O(10), or greater than this value. For instance, when K(x, τ ) = e xτ , or e x 2 τ 2 , and x ∈ [0, 5], τ ∈ [0, 1], then the maximum values become greater than or equal to e 5 ≈ 148.41, and e 25 ≈ 7.2 × 10 10 , which means the magnitude K( τ ) 2 cannot be neglected, while for bounded kernels with K < O(10), it can be neglected. For the case when K < O(10), inequality (Equation 23) reduces to the same condition as inequality (Equation 21). We investigate inequality (Equation 23) when K >> 1, it can be written as by simplifying inequality (Equation 24), and applying the limit as N x → +∞ one gets In order to discuss the inequality (Equation 26), we use the deterministicˆ of the quadrature equation such that Whenˆ ≤ 0, the inequality (Equation 26) unconditionally holds. Forˆ > 0, i.e., ( x) 2 < (1−2θ) √ K , then we have two real roots say δ 1 , and δ 2 , between them the inequality (Equation 26) is broken, while outside this interval the inequality holds. Consequently, we have two choices for τ , either τ ∈ (0, δ 1 ) or τ ∈ (δ 2 , +∞). In other words, we can find a suitable τ that is not necessarily to be so small to satisfy the inequality (Equation 26). Finally, for the remaining terms of an integral part, i.e., {A −1 (n,j) } n−1 j=0 , to investigate its stability, it is sufficient to show that the norm or A −1 is less than 1. Based on Equation (18), we have We summarize the stability of scheme (Equation 10) in the following theorem.

. Consistency
Here, we study the consistency of the proposed finite difference scheme to guarantee the convergence of the obtained numerical solutions. It is said that a finite difference scheme is consistent with a PIDE problem if the absolute difference between the exact theoretical solution and its corresponding numerical solution obtained from the used finite difference scheme tends to zeros as the discretization stepsizes tend to zero [19,20]. Let v n k = V(x k , τ n ) represent the exact solution at the point (x k , τ n ) and it is smooth enough, the local truncation error T n k (V) is given by where L n k (V), and J n k (V) are the local truncated errors of the differential, and integral operators, respectively. Our aim is to show that T n k (V) tends to 0, as x, and τ approach 0. We start with the differential part, by applying Taylor's expansion of a single valued function on v n k+1 , v n k−1 , and v n+1 k around the point (x k , τ n ), and the expansion of a function of two variables [25] on v n+1 k−1 , and v n+1 k+1 around the same point, one gets and the local truncated error of an integral part is given by By applying the modulus and taking the limits on Equations (31) and (32), as x, τ → 0, then the principal part of the local truncated error is Based on Equation (33), we obtain a local truncated error or order ( τ + ( x) 2 ), for θ ∈ [0, 1]\{ 1 2 }, while for θ = 1 2 , its order (( τ ) 2 + ( x) 2 ). In other words, the proposed θ -finite difference scheme is unconditionally consistent.

. Numerical experiments
Here, we present three examples that solve three different problems; two have exact solutions, while the exact solution of the last problem is unknown. In Example 1, we plot the numerical solutions for different values of the parameter θ since it represents the ratio between explicit and implicit finite difference schemes to explore the behavior of the numerical solutions comparable to the exact solution. Next, we investigate three aspects of the error associated with the numerical computations; first, the related maximum absolute errors are computed for each value of θ , and their orders are examined in view of Equation (33). Second, the orders of absolute errors are plotted. Third, the rate of convergence due to the change of spatial and time stepsizes is calculated, and the elapsed time to obtain the numerical solutions is tabulated for various values of mesh grids and θ . After that, the behavior of the spatial and time stepsizes are investigated when they violate the stability conditions. Example 2 is dedicated to comparing the performance of the proposed finite difference scheme with other works. To reveal the importance of Lax's equivalent theorem 1, we select a PIDE problem (Example 3) that its exact solution is unknown and obtain the numerical approximations when the stepsizes satisfy the stability conditions, and also when the stability conditions are broken. These examples are done using a 2.7GHz Xeon processor and Matlab simulation.
We examine the properties of the maximum absolute error associated with the numerical computations from two aspects; .
/fams. . its consistency based on Equation (33) and its convergence rates as the computational domain get finer. First, the order of the consistency of E max is calculated when (N x , N τ ) = (32, 256) for each value of the θ -parameter as shown in Table 1 which matches with Equation (33). The maximum absolute error is computed for various mesh grids of (N x , N τ ); after that, the ratio between two successive maximum errors is calculated to reveal the associated convergence rates as illustrated in Table 2. On the other hand, the order of the absolute error (ln E(x, τ = T)) is plotted as a function of the spatial variable x for various values of θ , as shown in Figure 3. The behavior of its order under the explicit, Crank-Nicolson, and fully implicit are grouped in Figure 3A. Figure 3B represents its behavior under the mixed schemes θ = 1 6 , 1 3 , 2 3 . It has been observed that the maximum error is of order 10 −3 when x = x p around the value 1. The computational interval [0, 2π ], can be written as the union of two intervalsˆ p = [0, x p ], and c p = [x p , 2π ]. In the first intervalˆ p , the order of error swiftly increases, while it gradually decreases in the second intervalˆ c p , and it reaches to 10 −7 , 10 −9 , 10 −13 at some points for θ = 0, 1 2 , and 1, respectively as shown in Figure 3A. Moreover, the order of the absolute error significantly decreases as θ increases.
Next, we discuss the behavior of the associated absolute error (AE) of the proposed finite difference scheme with respect to the change of time-discretization, while the spatial discretization is fixed and vice versa. The absolute errors at x = π/8, 5π/8, 5π/4, and 27π/16 are listed in Table 3 for various values of N τ together with the corresponding convergence rates β, which is defined by , and its computational cost in seconds, while the value of N x is selected as 16. Table 4 is dedicated to revealing the convergence rates and computational cost of the absolute error at τ = 1, 2, 3, 4, due to the change of the number of spatial mesh points N x , while N τ remains constant and equals 1,024 for the mixed finite difference scheme θ = 1 6 , 1 3 , 2 3 . Finally, we investigate the importance of the stability conditions (Equations 27-29). The computational domain is consequently, we investigate the last condition (Equation 29) for explicit scheme θ = 0, and a mixed scheme θ = 1 3 . The number of mesh points is selected to be N x = 16, N τ = 32, which implies r = 0.81 violates the stability condition r ≤ 0.5, θ = 0.5, 0.75, θ = 1 3 .
On the one hand, Figure 4 represents the surface solution as a function in x, and τ , when the stability condition is broken. On the other hand, Figure 5 shows a cross-section of the solution when τ = T.
Example 2. The goal of this example is to compare the behavior of the proposed finite difference scheme with other techniques, such as the Legendre-collocation method [26]. Problem 1 is chosen [26], where the PIDE problem has kernel K(x, τ ) = e xτ , boundary conditions; V(0, τ ) = sin τ , V(1, τ ) = 0, and initial condition V(x, 0) = 0, and source function The problem has the exact solution Here, the numerical solutions are obtained for explicit, Crank-Nicolson, and fully implicit cases, also their corresponding absolute errors are calculated for mesh grid points (N x , N τ ) = (16, 512) as illustrated in Figure 6.
For more investigation of the error behavior associated with the computational solutions, we calculate the maximum errors .
/fams. .  for several grids, the ratios of consecutive E max , and the CPU time as reported in Table 5. It has been observed that for a given mesh grid, the absolute maximum error can be reduced by increasing the value of the θ -parameter. Furthermore, for a given value of θ , the associated error can be reduced by increasing the number of space and time mesh points. However, their increments have different effects; for instance, a modest increase in spatial node number enhances the related error but with a higher computational cost than the same increment in time nodes. Decreasing the space and time stepsizes together leads to second-order convergence. To wrap up, Figure 6 and Table 5 illustrate the performance of the proposed θ -finite   difference scheme considering the stability and consistency in solving this problem. On the other hand, the maximum absolute error associated with solving the same problem based on the Legendre-collocation method is available in Table 1 [26, p. 7].
Example 3. This example aims to obtain numerical approximations of a PIDE problem based on the proposed finite difference scheme whose exact solution is unknown. The kernel is K(x, τ ) = e x 2 τ 2 , initial and boundary conditions are in the absence of the source function, i.e., F(x, τ ) = 0. The computational domain is selected such that (x, τ ) ∈ [0,     numerical solutions at these points are calculated for various computational grids. Therefore, we calculate so-called the root mean square relative error (RMSRE) such that where V ref i (τ = T), and V N τ i are the numerical solutions at the refine grid and the corresponding values of other grids, respectively. Based on these calculations, we present the RMSRE, its ratio of convergence and CPU time in Table 6.
Finally, we show that the stability conditions (Equations 27-29) cannot be neglected. The computational domain (x, τ ) ∈ [0, √ 6] × [0, 1], consequently, the maximum value of the integral-kernel is K = e 6 ≈ 403.429 >> 1. So that we investigate the stability condition (Equation 28). Note that it consists of two inequalities; one for the time stepsize and the other for the spatial stepsize. In light of this fact, we consider two cases. The first case is when τ satisfies its condition, while the inequality of x does not hold. In the second case, when τ breaks its inequality, while x satisfies its inequality. For the . /fams. .   Figure 8 shows the surface of numerical solutions, and a cross-section when τ = T of case one, i.e., x, does not satisfy its stability condition. On the other hand, Figure 9 describes the numerical solution behavior under the second case when (N x , N τ ) = (10, 20) (i.e., τ = 0.05 > 1
Therefore, stability conditions have a crucial role in obtaining reliable approximations.

. Conclusion
The PIDE of memory type has been investigated based on a suitable finite difference scheme. Here the differential part is discretized using the θ -approximation that combines explicitimplicit discretization of the second-order spatial derivative, while the time derivative is approximated using explicit forward discretization. The integral part has two obstacles; first, the . /fams. . upper boundary of the integration is variable, and second, the discretization of an integral part produces more unknowns other than the unknowns of the differential part. We implement the Heaviside step function to convert the upper variable boundary into a constant without changing the independent variables of the unknown function. Consequently, it enables us to discretize it on the same mesh points of the computational domain. The stability of the resulting finite difference scheme has been investigated based on the properties of the matrix norm and Gerschgoren's theorems. Also, consistency has been considered to guarantee the convergence of the numerical approximation. In light of consistency and stability analysis, we have a conditional stable finite difference scheme for θ ∈ [0, 0.5), while it is unconditionally stable for θ ∈ [0.5, 1]. Moreover, it is unconditional consistency. Finally, three examples have been studied to examine the accuracy and stability of the used scheme. For future research, we will extend the technique used to solve this category of partial integro-differential equations in higher dimensions. Also, we will introduce nonstandard finite difference schemes and investigate their numerical properties.

Data availability statement
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.