# Stable Numerical Solutions Preserving Qualitative Properties of Nonlocal Biological Dynamic Problems

Authors: M.-A. Piqueras; R. Company; L. Jódar
Journal: Abstract and Applied Analysis
 (2019)

---

## Abstract

This paper deals with solving numerically partial integrodifferential equations appearing in biological dynamics models when nonlocal interaction phenomenon is considered. An explicit finite difference scheme is proposed to get a numerical solution preserving qualitative properties of the solution. Gauss quadrature rules are used for the computation of the integral part of the equation taking advantage of its accuracy and low computational cost. Numerical analysis including consistency, stability, and positivity is included as well as numerical examples illustrating the efficiency of the proposed method.


---

## Body

## 1. 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 integro-differential reaction-diffusion problem (PIDE); see [19]:(1)∂U∂t=DΔU+βUx,t1-aUx,t-b∫Ωψx-yUy,tdy,x∈Ω,t∈0,T,where Ω is a bounded or unbounded domain in R2 and ψ(x) is a nonnegative kernel function satisfying(2)∫R2ψxdx=1,and β, a, and b are some positive constants and D is a positive dispersal rate, together with the initial conditions(3)Ux,0=fx,x∈Ωand the boundary conditions(4)Ux,t=0,x∈∂Ω,where f(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 pseudo-spectral 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 h(x,t) 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. Section2 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.


## 2. 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 bounded numerical domain. Let us consider the numerical domain[-A,A]2×[0,T], with A>0 large enough so that outside of this area the population is negligible and T>0 denoting the time horizon. Let M and N be positive integers, so that the domain [-A,A]2×[0,T] is partitioned in (2M+1)2×(N+1) mesh points denoted by (x1,i,x2,j,tn), where x1,i=ih, -M≤i≤M, x2,j=jh, -M≤j≤M, and tn=nk, 0≤n≤N. The step sizes discretizations h and k verify hM=A and kN=T, respectively. The numerical approximation of the unknown variable at the mesh point (x1,i,x2,j,tn) is denoted by ui,jn≈U(x1,i,x2,j,tn), while for the integral term in (1), we designate(5)gi,jn≈Gx1,i,x2,j,tn=∫Ωψxi,j-yUy,tndy,-M≤i,j≤M,n≥0.where xi,j=(x1,i,x2,j).The approximationgi,jn of the integral term G(x1,i,x2,j,tn) 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 ψ(x) 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 gi,jn; see [1], page 882.According to the expression for the Gauss-Hermite quadrature, we have(6)∫-∞∞∫-∞∞ψxi,j-yUy,tndy≈∑l=1L∑m=1Lwlwmexl2+xm2ψx1,i-y^1,l,x2,j-y^2,mUy^1,l,y^2,m,tn,where wl and wm are the weights and y^1,l and y^2,m, 1≤l,m≤L, are the nodes of the Gauss-Hermite quadrature, respectively.Given a node(y^1,l,y^2,m), 1≤l,m≤L, let us consider the indexes il and jm such that the grid point (y1,il,y2,jm) verifies(7)y1,il≤y^1,l≤y1,il+1,y2,jm≤y^2,m≤y2,jm+1.The bilinear interpolation approximation is denoted by(8)u¯y^1,l,y^2,m,tn=ρlρmuil,jmn+ρl1-ρmuil,jm+1n+1-ρlρmuil+1,jmn+1-ρl1-ρmuil+1,jm+1n,where(9)ρl=y1,il+1-y^1,lh,ρm=y2,jm+1-y^2,mh.With previous notation, the approximationgi,jn of the integral term G(x1,i,x2,j,tn) takes the form(10)gi,jn=∑l=1L∑m=1Lwlwmexl2+xm2ψx1,i-y^1,l,x2,j-y^2,mu¯y^1,l,y^2,m,tn.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,(11)∫-∞∞∫-∞∞ψxi,j-yUy,tndy=∫x2,j-δx2,j+δ∫x1,i-δx1,i+δψxi,j-yUy,tndy.By using the change of variables(12)z=y-xi,jδ,in the right-hand side of (11) it follows that(13)Gx1,i,x2,j,tn=δ2∫-11∫-11ψ-δzUδz+xi,j,tndz.Applying Gauss-Legendre quadrature together with bilinear interpolation in analogously way to (8), the numerical approximation gi,jn of (12) takes the form(14)gi,jn=δ2∑l=1L∑m=1Lwlwmψ-δz1,l,-δz2,mu¯δz1,l+x1,i,δz2,m+x2,j,tn,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,(15)ui,jn+1-ui,jnk≈∂U∂tx1,i,x2,j,tn,ui+1,jn-ui-1,jn2h≈∂U∂x1x1,i,x2,j,tn,ui+1,jn-2ui,jn+ui-1,jnh2≈∂2U∂x12x1,i,x2,j,tnand 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:(16)ui,jn+1=Dkh2ui-1,jn+ui+1,jn+ui,j-1n+ui,j+1n+1-4Dkh2ui,jn+kβui,jn1-aui,jn-bgi,jn,-M+1≤i,j≤M-1,0≤n≤N-1, with initial and transferred boundary conditions of our numerical domain(17)ui,j0=fxi,j,u-M,jn=uM,jn=ui,-Mn=ui,Mn=0,1≤i,j≤M.


## 3. 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 U(x,t) for the PIDE problem (1), (3), and (4) under the assumptions that U(x,0)∈C0(Ω¯) and 0≤U(x,0)≤1/a. Furthermore, this solution verifies the constraints(18)0≤Ux,t≤1a,x∈Ω¯,t≥0.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{ui,jn} 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≤ui,jn≤1/a, we study whether 0≤ui,jn+1≤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(19)gi,jn≤1+ϵa,-M≤i,j≤M,0≤n≤N. Adopting the arbitrary value ϵ=1 and assuming the boundedness of the solution at the temporal level n, 0≤ui,jn≤1/a, we can write, according to expression (16) for the n+1 time level,(20)ui,jn+1≥Dkh2ui-1,jn+ui+1,jn+ui,j-1n+ui,j+1n+1-2k2Dh2+bβaui,jn,resulting that ui,jn+1≥0, -M≤i,j≤M, under the condition for the temporal step size:(21)k<h222D+b/aβh2.Regarding the boundedness of the numerical solution, since the termgi,jn is nonnegative, from the scheme (16) we can write(22)ui,jn+1≤Dkh2ui-1,jn+ui+1,jn+ui,j-1n+ui,j+1n+1-4Dkh2ui,jn+kβui,jn1-aui,jn=ϕui,jn,ui-1,jn,ui+1,jn,ui,j-1n,ui,j+1n,by introducing the function ϕ of several real variables, assumed to be differentiable with respect to all its arguments.Taking partial derivatives with respect toui,jn, and ui,jn≤1/a,(23)∂ϕ∂ui,jn=1-k4Dh2+2βaui,jn-β≥1-k4Dh2+β.Then, under the assumption(24)k<h24D+βh2,the function ϕ is increasing over the range 0≤ui,jn≤1/a, and consequently the following inequality holds:(25)ui,jn+1≤Dkh2ui-1,jn+ui+1,jn+ui,j-1n+ui,j+1n+1-4Dkh21a≤Dkh24a+1-4Dkh21a=1a.In conclusion, assuming that0≤ui,j0≤1/a and taking a temporal step size k such that(26)k<h24D+βαh2,α=max⁡1,2ba,it is guaranteed that 0≤ui,jn≤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(27)k<h22D+βαh2,α=max⁡1,2ba.


## 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 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.Let us consider the (1), in a compact form as L(U)=L(1)(U)+L(2)(U)=0, where(28)L1U=∂U∂t-DΔU,L2U=-βUx,t1-aUx,t-b∫Ωψx-yUy,tdy,x∈Ω,t∈0,T,and the finite difference scheme (16), written as L(u)=L(1)(u)+L(2)(u)=0, where(29)L1u=ui,jn+1-ui,jnk-Dui+1,jn-2ui,jn+ui-1,jnh2+ui,j+1n-2ui,jn+ui,j-1nh2,L2u=-βui,jn1-aui,jn-bgi,jn,-M+1≤i,j≤M-1,0≤n≤N-1.In accordance with [17], scheme L(u) is said to be consistent with problem L(U) if local truncation error Ti,jn(U),(30)Ti,jnU=LUi,jn-LUi,jn,tends to zero as k→0, h→0, where Ui,jn=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 times with respect to x1 and x2, 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)i,jn(U)=L(1)(Ui,jn)-L(1)(Ui,jn), associated to the differential part of (1):(31)T1i,jnU=Ei,jn1k-DEi,jn2+Ei,jn3h2,where(32)Ei,jn1=12∂2U∂t2x1,i,x2,j,τ,tn<τ<tn+1,(33)Ei,jn2=112∂4U∂x14ξ1,x2,j,tn,x1,i-1<ξ1<x1,i+1,(34)Ei,jn3=112∂4U∂x24x1,i,ξ2,tn,x2,j-1<ξ2<x2,j+1.Regarding the local truncation errorT(2)i,jn(U)=L(2)(Ui,jn)-L(2)(Ui,jn), related to the integral term in (1),(35)T2i,jnU=βbUi,jnG¯i,jn-∫Ωψxi,j-yUy,tndy,where G¯i,jn 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(36)U¯y^1,l,y^2,m,tn=ρlρmUil,jmn+ρl1-ρmUil,jm+1n+1-ρlρmUil+1,jmn+1-ρl1-ρmUil+1,jm+1n,making use of relations (9).It can be verified, see [1], that the numerical approximation G¯i,jn satisfies(37)G¯i,jn=G~i,jn+Oh2,being G~i,jn the value of the expression for G¯i,jn replacing the interpolating value U¯(y^1,l,y^2,m,tn) by the exact solution value U(y^1,l,y^2,m,tn).Moreover,(38)G~i,jn-∫Ωψxi,j-yUy,tndy=ϵLis 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 ofT(2)i,jn(U) and (37) and (38) that(39)T2i,jnU=Oh2+ϵL,and the local truncation error Ti,jn(U) satisfies(40)Ti,jnU=Ok+Oh2+ϵL.


## 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(41)fx=14,-4≤x≤4,0,otherwise,(42)ψξ=12,-1≤ξ≤1,0,otherwise.The spatial and temporal step sizes are chosen ash=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 the maximum of gi,jn 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 solution approaches the habitat carrying capacity 1/a and the occupied habitat tends to broaden.Figure 1
Numerical solution in the case of one space dimension and unbounded domain.The convergence rateγ of the solution has also been computed, which is defined here as the quantity γ=log⁡ϵi/ϵj/log⁡(ki/kj), where ϵi holds for the relative root mean squared error of the i-th numerical solution an ki is the step size used to obtain that solution. Table 1 shows the results taking different values for the temporal step size k.Table 1
Relative root mean squared error RRMSE and convergence rateγ, taking k4=0.0040 as the reference step size, and using h=0.05 as spatial step size in all cases.
k 4=0.0005 k 1=0.0040 k 2=0.0020 k 3=0.0010 RRMSE 0.0014733979 0.0006292592 0.0002095394 γ - 1.22742 1.58643Example 2.
Let us consider the same model as in the previous example, with identical parameters values, initial condition, kernel functionψ(ξ), number of nodes L, and time horizon T. If we choose a temporal step size k=2/398≃0.005025, breaking the stability condition (27), it is clear from Figure 2 that the behaviour of the numerical solution U(x,t) becomes unstable and it reaches negatives values.Figure 2
Numerical solution when the positivity and stability condition is broken.Example 3.
In this example, we consider the case of two unbounded spatial dimensions. Here, the parameters are chosen to take the values(D,β,a,b) = (0.2,1,1,1). The initial condition f(x1,x2) and the kernel function ψ(ξ1,ξ2) are taken as(43)fx1,x2=12,x12+x22≤122,0,otherwise,(44)ψξ1,ξ2=1π,ξ12+ξ22≤1,0,otherwise.Note that, with this data, the problem presents radial symmetry. The spatial and temporal step sizes are chosen ash=0.1 and k=0.005, respectively. The number of nodes of the quadrature is L=10 large enough to apply results of Section 3 because the maximum of gi,jn 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 k<0.012195, which is satisfied here. Figure 3 shows the numerical solution U(x1,x2,t) at the time horizon t=T=1.Figure 3
Numerical solution for two space dimension in an unbounded domain fort=T=1.Example 4.
Taking the same data as in Example3, and identical number of nodes L, if we choose a temporal step size k=0.0125, that does not satisfy the condition (27), the numerical solution U(x,t) becomes unstable and it also reaches negatives values, as shown in Figure 4.Figure 4
Numerical solution for two space dimension in an unbounded domain when the positivity and stability condition is broken, fort=T=1.Example 5.
In this example, using the same data and step sizes as in Example3 and identical number of nodes L, we take a Gaussian kernel given by the expression:(45)ψξ1,ξ2=12πe-ξ12+ξ22/2.The numerical solutionU(x,t) is stable and positive, as shown in Figure 5.Figure 5
Numerical solution for two space dimension in an unbounded domain, taking a Gaussian kernel, fort=T=1.



---

*Source: 5787329-2019-07-01.xml*