Stable Numerical Solutions Preserving Qualitative Properties of Nonlocal Biological Dynamic Problems

and Applied Analysis 3 By using the change of variables z = y − xi,j δ , (12) in the right-hand side of (11) it follows that G(x1,i, x2,j, tn) = δ2 ∫1 −1 ∫1 −1 ψ (−δz) U (δz + xi,j, tn) dz. (13) ApplyingGauss-Legendre quadrature togetherwith bilinear interpolation in analogously way to (8), the numerical approximation gn i,j of (12) takes the form gn i,j = δ2 L ∑ l=1 L ∑ m=1 wlwmψ (−δz1,l, −δz2,m) ⋅ u (δz1,l + x1,i, δz2,m + x2,j, tn) , (14) using the weights and nodes of the Gauss-Legendre formula. Regarding the differential part of PIDE (1), let us consider the forward approximation of the time derivatives and central approximation of the spatial derivatives, un+1 i,j − un i,j k ≈ ∂U ∂t (x1,i, x2,j, tn) , un i+1,j − un i−1,j 2h ≈ ∂U ∂x1 (x1,i, x2,j, t n) , un i+1,j − 2un i,j + un i−1,j h2 ≈ ∂2U ∂x2 1 (x1,i, x2,j, t n) (15) and analogous expressions for derivatives with respect to the second spatial variable x2. From (5) and (15) the following explicit numerical scheme for (1),(3), and(4) has been constructed: un+1 i,j = Dk h2 (un i−1,j + un i+1,j + un i,j−1 + un i,j+1) + (1 − 4Dk h2 )un i,j + kβun i,j (1 − aun i,j − bgn i,j) , −M + 1 ≤ i, j ≤ M − 1, 0 ≤ n ≤ N − 1, (16) with initial and transferred boundary conditions of our numerical domain u0 i,j = f (xi,j) , un −M,j = un M,j = un i,−M = un i,M = 0, 1 ≤ i, j ≤ M. (17) 3. Positivity, Stability, and Carrying Capacity Bound According to Lemma 3.1 and Theorem 3.2 of [15] it holds true that there exists a global solution U(x, t) for the PIDE problem (1), (3), and (4) under the assumptions thatU(x, 0) ∈ C0(Ω) and 0 ≤ U(x, 0) ≤ 1/a. Furthermore, this solution verifies the constraints 0 ≤ U (x, t) ≤ 1 a , x ∈ Ω, t ≥ 0. (18) This section is devoted to the numerical analysis of the proposed scheme, guaranteeing the preservation of the qualitative properties of the theoretical solution. In the followingwe show that under appropriate step size conditions the numerical solution {un i,j} is nonnegative and is bounded by the carrying capacity 1/a in agreement with (18). Thus the stability of the numerical solution is granted because it is bounded. Using the induction principle on the temporal index n and assuming the induction hypothesis 0 ≤ un i,j ≤ 1/a, we study whether 0 ≤ un+1 i,j ≤ 1/a also remains true. Taking into account the equality (2) for the argument xi,j− y, it holds, for a large enough value of L, that gn i,j ≤ 1 + ε a , − M ≤ i, j ≤ M, 0 ≤ n ≤ N. (19) Adopting the arbitrary value ε = 1 and assuming the boundedness of the solution at the temporal level n, 0 ≤ un i,j ≤ 1/a, we can write, according to expression (16) for the n + 1 time level, un+1 i,j ≥ Dk h2 (un i−1,j + un i+1,j + un i,j−1 + un i,j+1) + (1 − 2k(2D h2 + bβ a ))un i,j, (20) resulting that un+1 i,j ≥ 0, −M ≤ i, j ≤ M, under the condition for the temporal step size: k < h2 2 (2D + (b/a) βh2) . (21) Regarding the boundedness of the numerical solution, since the term gn i,j is nonnegative, from the scheme (16) we can write un+1 i,j ≤ Dk h2 (un i−1,j + un i+1,j + un i,j−1 + un i,j+1) + (1 − 4Dk h2 )un i,j + kβun i,j (1 − aun i,j) = φ (un i,j, un i−1,j, un i+1,j, un i,j−1, un i,j+1) , (22) by introducing the function φ of several real variables, assumed to be differentiable with respect to all its arguments. Taking partial derivatives with respect to un i,j, and un i,j ≤ 1/a, ∂φ ∂un i,j = 1 − k ( 4D h2 + 2βaun i,j − β) ≥ 1 − k (4D h2 + β) . (23) 4 Abstract and Applied Analysis Then, under the assumption k < h2 4D + βh2 , (24) the function φ is increasing over the range 0 ≤ un i,j ≤ 1/a, and consequently the following inequality holds: un+1 i,j ≤ Dk h2 (un i−1,j + un i+1,j + un i,j−1 + un i,j+1) + (1 − 4Dk h2 ) 1 a ≤ Dk h2 4 a + (1 − 4Dk h2 ) 1 a = 1 a . (25) In conclusion, assuming that 0 ≤ u0 i,j ≤ 1/a and taking a temporal step size k such that k < h2 4D + βαh2 , α = max{1, 2b a } , (26) it is guaranteed that 0 ≤ un i,j ≤ 1/a, 1 ≤ n ≤ N. Remark 1. Note that stability and positivity step size condition is linked to the problem dimension. In particular, for the one dimensional case, the condition becomes k < h2 2D + βαh2 , α = max{1, 2b a } . (27) 4. Consistency In this section we study the consistency of the numerical solution, given by the scheme (16), with the problem (1)(4). Consistency of a numerical scheme with the respective continuous problemmeans that the theoretical solution of the problem approximates well the numerical scheme when the step size discretizations tend to zero. So, a numerical scheme can be consistent with an equation and not with another one; see [20], Chapter 2. Thus, it is important to address the consistency of a numerical scheme with a problem. Let us consider the (1), in a compact form as L(U) = L(1)(U) +L(2)(U) = 0, where L (1) (U) = ∂U ∂t − DΔU, L (2) (U) = −βU (x, t) ⋅ (1 − aU (x, t) − b∫ Ω ψ (x − y) U (y, t) dy) , x ∈ Ω, t ∈ [0, T] , (28) and the finite difference scheme (16), written as L(u) = L(1)(u) + L(2)(u) = 0, where L (1) (u) = u n+1 i,j − un i,j k − D( un i+1,j − 2un i,j + un i−1,j h2 + u n i,j+1 − 2un i,j + un i,j−1 h2 ) , L (2) (u) = −βun i,j (1 − aun i,j − bgn i,j) , −M + 1 ≤ i, j ≤ M − 1, 0 ≤ n ≤ N − 1. (29) In accordance with [20], scheme L(u) is said to be consistent with problemL(U) if local truncation error Tn i,j(U), Tn i,j (U) = L (Un i,j) −L (Un i,j) , (30) tends to zero as k 󳨀→ 0, h 󳨀→ 0, where Un i,j = U(x1,i, x2,j, tn) is the value of the exact solution of problem (1)-(4) of the PIDE at the point (x1,i, x2,j, tn). Now let us assume that the exact solution U(x1, x2, t) is continuously partial differentiable four timeswith respect tox1 andx2, and, two times with respect to t. By using Taylor’s expansion about (x1,i, x2,j, tn), one gets the expression of the local truncation error T(1)ni,j(U) = L(1)(Un i,j) − L(1)(Un i,j), associated to the differential part of (1): T (1)ni,j (U) = En i,j (1) k − D (En i,j (2) + En i,j (3)) h2, (31) where En i,j (1) = 1 2 ∂2U ∂t2 (x1,i, x2,j, τ) , tn < τ < tn+1, (32) En i,j (2) = 1 12 ∂4U ∂x14 (ξ1, x2,j, t n) , x1,i−1 < ξ1 < x1,i+1, (33) En i,j (3) = 1 12 ∂4U ∂x24 (x1,i, ξ2, t n) , x2,j−1 < ξ2 < x2,j+1. (34) Regarding the local truncation error T(2)ni,j(U) = L(2)(Un i,j) −L(2)(Un i,j), related to the integral term in (1), T (2)ni,j (U) = βbUn i,j (Gi,j − ∫ Ω ψ (xi,j − y)U (y, tn) dy) , (35) where Gi,j is the numerical approximation of ∫ Ω ψ(xi,j − y)U(y, tn)dy by means of the corresponding Gauss-quadrature rule together with the use of a linear interpolation U (y1,l, y2,m, tn) = ρlρmUn il,jm + ρl (1 − ρm) Un il,jm+1 + (1 − ρl) ρmUn il+1,jm + (1 − ρl) (1 − ρm) Un il+1,jm+1, (36) making use of relations (9). Abstract and Applied Analysis 5and Applied Analysis 5 It can be verified, see [19], that the numerical approximation Gi,j satisfies Gi,j = Gn i,j + O (h2) , (37) being Gn i,j the value of the expression for Gi,j replacing the interpolating valueU(y1,l, y2,m, tn) by the exact solution value U(y1,l, y2,m, tn). Moreover, Gn i,j − ∫ Ω ψ (xi,j − y)U (y, tn) dy = ε (L) (38) is the associated quadrature error of the two-dimensional corresponding Gauss quadrature formula. An estimation of the error bound for Gaussian quadrature rules in two dimensions can be found in [21] using divided differences, assuming the integrand to be differentiable in the region of integration. It can be verified from the expression of T(2)ni,j(U) and (37) and (38) that T (2)ni,j (U) = O (h2) + ε (L) , (39) and the local truncation error Tn i,j(U) satisfies Tn i,j (U) = O (k) + O (h2) + ε (L) . (40) 5. Numerical Examples This section illustrates the behaviour of the numerical solution of the problem (1)-(4) with some numerical experiments, making use of the proposed scheme (16), in both cases of one and two space dimensions. Example 1. Let us consider the nonlocal logistic diffusion model (1)-(4) in an unbounded one space dimension, with parameters values (D, β, a, b) = (0.25, 5, 1, 1). The initial condition f(x) and the function ψ(ξ) in the integral term are taken, respectively, as f (x) = {{ 1 4 , −4 ≤ x ≤ 4, 0, otherwise, (41) ψ (ξ) = {{ 1 2 , −1 ≤ ξ ≤ 1, 0, otherwise. (42) The spatial and temporal step sizes are chosen as h = 0.05 and k = 0.004, respectively. The number of nodes of the quadrature is taken as L = 10 and it is large enough to apply results of stability Section 3 because themaximumof gn i,j does not exceed value 0.675 and assumption (19) is verified for ε = 1. According to the expression (27), if k < 0.004762, which is fulfilled in this case, the positivity and stability of the solution are guaranteed. Figure 1 shows the behaviour of the numerical solution U(x, t) from t = 0 to the time horizon T = 2. Note that as time increases, the numerical 1.5


Introduction
Since the seminal papers by Fisher [8] and Kolmogorov-Petrovsky-Piskunov (KPP) [13], the diffusive logistic models related to local interaction dynamic biological models have been successfully developed [3,14,15,20,21]. In ecological context, there is no real justification for assuming that the interactions are local [9]. Also, in evolutionary theory, where interactions are due not only to intraspecific competition but also random mutations, the nonlocal interaction approach is necessary [6,11].
Important theoretical results about existence of solutions and qualitative properties of nonlocal biological dynamic problems have been treated in [4,10,12,18,19].
In this paper we consider the nonlocal interaction biological dynamic model described by the partial integrodifferential reaction-diffusion problem (PIDE); see [19] where Ω is a bounded or unbounded domain in R 2 and (x) is a nonnegative kernel function satisfying and , , and are some positive constants and is a positive dispersal rate, together with the initial conditions and the boundary conditions where (x) represents an arbitrary continuous function.
For the sake of interest in consistency issues that we study in Section 4, we assume that the kernel function (x) is differentiable in the region of integration. From the biological point of view, the first term of the right-hand side models the diffusion, and the second includes the pure logistic quadratic term and the consumption of resources in some area around the average location. Note that if the kernel (x) is the Dirac delta function centered at the origin, (1) recovers the Fisher-KPP equation.
Numerical analysis of the problem is suitable because the best model may be wasted by a disregarded numerical treatment. Numerical methods dealing with problem (1)- (4) are not extensive in the literature. Some relevant and easy to implement worthy exceptions are [7] using a pseudospectral Galerkin method and [16] that uses a meshless local radial point interpolation method. Both papers consider the bounded domain case and add a source term ℎ(x, ) in the right-hand side of the PIDE that is useful to check accuracy simulations.
In this paper we develop an explicit finite difference scheme for the numerical computation and analysis of qualitative properties preserving numerical solutions of problem (1)- (4). For the best of our knowledge about the publications in the field this analysis seems not easily reachable with other techniques such as finite elements, meshless, and even implicit difference methods. It is relevant to point out that the integral term of the PIDE is treated using Gauss quadrature rules having the versatility advantage of including both the bounded and unbounded domain cases, just adapting the quadrature rule.
Positivity of the numerical solutions is essential dealing with a population problem and needs to be guaranteed. It is also important to check that numerical solutions are bounded by the carrying capacity of the problem, in agreement with the behaviour of the theoretical solution [19].
This paper is organized as follows. Section 2 presents the discretization of the continuous problem, reaching an explicit finite difference scheme. Section 3 deals with the positivity, stability, and qualitative properties of the numerical solution. The consistency of the scheme with the PIDE is studied in Section 4. Finally, Section 5 features some numerical experiments, showing that if the step size requirements for stability are not satisfied results are wrong.

Discretization and Numerical Scheme Construction
In this section, we perform the discretization of the continuous problem, with the goal to reach an explicit finite difference scheme. Hereafter, we will work in a suitable  The approximation , of the integral term ( 1, , 2, , ) is performed by means of the accurate and computationally cheap Gauss quadrature rule; see [5], Chapters 2 and 3. Gauss-Hermite or Gauss-Legendre quadratures are used depending on whether the support of the kernel function ( ) is unbounded or compact, respectively. As the nodes of the quadrature rule are not necessarily mesh points of the grid, a bilinear interpolation is used for the computation of the terms , ; see [1], page 882.
According to the expression for the Gauss-Hermite quadrature, we have where and are the weights and̂1 , and̂2 , , 1 ≤ , ≤ , are the nodes of the Gauss-Hermite quadrature, respectively.
Given a node (̂1 , ,̂2 , ), 1 ≤ , ≤ , let us consider the indexes and such that the grid point ( 1, , 2, ) verifies The bilinear interpolation approximation is denoted by where With previous notation, the approximation , of the integral term ( 1, , 2, , ) takes the form The Gauss-Legendre quadrature rule is appropriate in the case that the kernel function has a compact support. For instance, let us consider (x) with square support [− , ] 2 . Then, Abstract and Applied Analysis 3 By using the change of variables in the right-hand side of (11) it follows that Applying Gauss-Legendre quadrature together with bilinear interpolation in analogously way to (8), the numerical approximation , of (12) takes the form using the weights and nodes of the Gauss-Legendre formula.
Regarding the differential part of PIDE (1), let us consider the forward approximation of the time derivatives and central approximation of the spatial derivatives, and analogous expressions for derivatives with respect to the second spatial variable 2 . From (5) and (15) the following explicit numerical scheme for (1), (3), and(4) has been constructed: with initial and transferred boundary conditions of our numerical domain 0 , = (x , ) ,

Positivity, Stability, and Carrying Capacity Bound
According to Lemma 3.1 and Theorem 3.2 of [19] it holds true that there exists a global solution (x, ) for the PIDE problem (1), (3), and (4) under the assumptions that (x, 0) ∈ 0 (Ω) and 0 ≤ (x, 0) ≤ 1/ . Furthermore, this solution verifies the constraints This section is devoted to the numerical analysis of the proposed scheme, guaranteeing the preservation of the qualitative properties of the theoretical solution. In the following we show that under appropriate step size conditions the numerical solution { , } is nonnegative and is bounded by the carrying capacity 1/ in agreement with (18). Thus the stability of the numerical solution is granted because it is bounded. Using the induction principle on the temporal index and assuming the induction hypothesis 0 ≤ , ≤ 1/ , we study whether 0 ≤ +1 , ≤ 1/ also remains true. Taking into account the equality (2) for the argument x , − y, it holds, for a large enough value of , that Adopting the arbitrary value = 1 and assuming the boundedness of the solution at the temporal level , 0 ≤ , ≤ 1/ , we can write, according to expression (16) for the + 1 time level, +1 , ≥ ℎ 2 ( −1, + i+1, + , −1 + , +1 ) resulting that +1 , ≥ 0, − ≤ , ≤ , under the condition for the temporal step size: Regarding the boundedness of the numerical solution, since the term , is nonnegative, from the scheme (16) we can write by introducing the function of several real variables, assumed to be differentiable with respect to all its arguments.

Abstract and Applied Analysis
Taking partial derivatives with respect to , , and , ≤ 1/ , Then, under the assumption the function is increasing over the range 0 ≤ , ≤ 1/ , and consequently the following inequality holds: In conclusion, assuming that 0 ≤ 0 , ≤ 1/ and taking a temporal step size such that it is guaranteed that 0 ≤ , ≤ 1/ , 1 ≤ ≤ .

Remark 1.
Note that stability and positivity step size condition is linked to the problem dimension. In particular, for the one dimensional case, the condition becomes

Consistency
In this section we study the consistency of the numerical solution, given by the scheme (16), with the problem (1)-(4). Consistency of a numerical scheme with the respective continuous problem means that the theoretical solution of the problem approximates well the numerical scheme when the step size discretizations tend to zero. So, a numerical scheme can be consistent with an equation and not with another one; see [17], Chapter 2. Thus, it is important to address the consistency of a numerical scheme with a problem.
In accordance with [17], scheme ( ) is said to be consistent with problem L( ) if local truncation error , ( ), tends to zero as → 0, ℎ → 0, where , = ( 1, , 2, , ) is the value of the exact solution of problem (1)-(4) of the PIDE at the point ( 1, , 2, , ). Now let us assume that the exact solution ( 1 , 2 , ) is continuously partial differentiable four times with respect to 1 and 2 , and, two times with respect to . By using Taylor's expansion about ( 1, , 2, , ), one gets the expression of the local truncation error (1) , ( ) = (1)( , ) − L(1)( , ), associated to the differential part of (1): where Abstract and Applied Analysis 5 where , is the numerical approximation of ∫ Ω (x , − y) (y, ) y by means of the corresponding Gauss-quadrature rule together with the use of a linear interpolation (̂1 , ,̂2 , , ) = , + (1 − ) , +1 making use of relations (9). It can be verified, see [1], that the numerical approximation , satisfies being̃, the value of the expression for , replacing the interpolating value (̂1 , ,̂2 , , ) by the exact solution value (̂1 , ,̂2 , , ). Moreover, is the associated quadrature error of the two-dimensional corresponding Gauss quadrature formula. An estimation of the error bound for Gaussian quadrature rules in two dimensions can be found in [2] using divided differences, assuming the integrand to be differentiable in the region of integration. It can be verified from the expression of (2) , ( ) and (37) and (38) that and the local truncation error , ( ) satisfies

Numerical Examples
This section illustrates the behaviour of the numerical solution of the problem (1)-(4) with some numerical experiments, making use of the proposed scheme (16), in both cases of one and two space dimensions.  The spatial and temporal step sizes are chosen as ℎ = 0.05 and = 0.004, respectively. The number of nodes of the quadrature is taken as = 10 and it is large enough to apply results of stability Section 3 because the maximum of , does not exceed value 0.675 and assumption (19) is verified for = 1. According to the expression (27), if < 0.004762, which is fulfilled in this case, the positivity and stability of the solution are guaranteed. Figure 1 shows the behaviour of the numerical solution ( , ) from = 0 to the time horizon = 2. Note that as time increases, the numerical solution approaches the habitat carrying capacity 1/ and the occupied habitat tends to broaden.
The convergence rate of the solution has also been computed, which is defined here as the quantity = log( / )/ log( / ), where holds for the relative root mean squared error of the − ℎ numerical solution an is the step size used to obtain that solution. Table 1 shows the results taking different values for the temporal step size .   that the behaviour of the numerical solution ( , ) becomes unstable and it reaches negatives values.
Example 3. In this example, we consider the case of two unbounded spatial dimensions. Here, the parameters are chosen to take the values ( , , , ) = (0.2, 1, 1, 1). The initial condition ( 1 , 2 ) and the kernel function ( 1 , 2 ) are taken as Note that, with this data, the problem presents radial symmetry. The spatial and temporal step sizes are chosen as ℎ = 0.1 and = 0.005, respectively. The number of nodes of the quadrature is = 10 large enough to apply results of Section 3 because the maximum of , does not exceed value 0.022 and assumption (19) is verified for = 1. Applying (26), it results that the positivity and stability are guaranteed when < 0.012195, which is satisfied here. Figure 3 shows the numerical solution ( 1 , 2 , ) at the time horizon = = 1.

Example 4.
Taking the same data as in Example 3, and identical number of nodes , if we choose a temporal step size = 0.0125, that does not satisfy the condition (27), the numerical solution ( , ) becomes unstable and it also reaches negatives values, as shown in Figure 4.  Example 5. In this example, using the same data and step sizes as in Example 3 and identical number of nodes , we take a Gaussian kernel given by the expression: The numerical solution ( , ) is stable and positive, as shown in Figure 5.

Data Availability
No data were used to support this study.

Conflicts of Interest
The authors of this work declare that there are no conflicts of interest regarding the publication of this paper.