Research Article  Open Access
Stable Numerical Solutions Preserving Qualitative Properties of Nonlocal Biological Dynamic Problems
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.
1. Introduction
Since the seminal papers by Fisher [8] and KolmogorovPetrovskyPiskunov (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 reactiondiffusion problem (PIDE); see [19]:where is a bounded or unbounded domain in and is a nonnegative kernel function satisfyingand , , and are some positive constants and is a positive dispersal rate, together with the initial conditionsand the boundary conditionswhere 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 is differentiable in the region of integration.
From the biological point of view, the first term of the righthand 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 is the Dirac delta function centered at the origin, (1) recovers the FisherKPP 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 in the righthand 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.
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 , with large enough so that outside of this area the population is negligible and denoting the time horizon. Let and be positive integers, so that the domain is partitioned in mesh points denoted by , where , , , , and , . The step sizes discretizations and verify and , respectively. The numerical approximation of the unknown variable at the mesh point is denoted by , while for the integral term in (1), we designatewhere .
The approximation of the integral term is performed by means of the accurate and computationally cheap Gauss quadrature rule; see [5], Chapters 2 and 3. GaussHermite or GaussLegendre 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 GaussHermite quadrature, we havewhere and are the weights and and , , are the nodes of the GaussHermite quadrature, respectively.
Given a node , , let us consider the indexes and such that the grid point verifiesThe bilinear interpolation approximation is denoted bywhere
With previous notation, the approximation of the integral term takes the form
The GaussLegendre quadrature rule is appropriate in the case that the kernel function has a compact support. For instance, let us consider with square support . Then,
By using the change of variablesin the righthand side of (11) it follows that
Applying GaussLegendre quadrature together with bilinear interpolation in analogously way to (8), the numerical approximation of (12) takes the formusing the weights and nodes of the GaussLegendre 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 . 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
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 for the PIDE problem (1), (3), and (4) under the assumptions that and . 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 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 , we study whether also remains true.
Taking into account the equality (2) for the argument , it holds, for a large enough value of , that Adopting the arbitrary value and assuming the boundedness of the solution at the temporal level , , we can write, according to expression (16) for the time level,resulting that , , 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 writeby introducing the function of several real variables, assumed to be differentiable with respect to all its arguments.
Taking partial derivatives with respect to , and ,
Then, under the assumptionthe function is increasing over the range , and consequently the following inequality holds:
In conclusion, assuming that and taking a temporal step size such thatit is guaranteed that , .
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
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 , whereand the finite difference scheme (16), written as , where
In accordance with [17], scheme is said to be consistent with problem if local truncation error ,tends to zero as , , where is the value of the exact solution of problem (1)(4) of the PIDE at the point . Now let us assume that the exact solution is continuously partial differentiable four times with respect to and , and, two times with respect to . By using Taylor’s expansion about , one gets the expression of the local truncation error , associated to the differential part of (1):where
Regarding the local truncation error , related to the integral term in (1),where is the numerical approximation of by means of the corresponding Gaussquadrature rule together with the use of a linear interpolationmaking use of relations (9).
It can be verified, see [1], that the numerical approximation satisfiesbeing the value of the expression for replacing the interpolating value by the exact solution value .
Moreover,is the associated quadrature error of the twodimensional 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 and (37) and (38) thatand the local truncation error satisfies
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 = . The initial condition and the function in the integral term are taken, respectively, as
The spatial and temporal step sizes are chosen as and , respectively. The number of nodes of the quadrature is taken as 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 . According to the expression (27), if , 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 to the time horizon . Note that as time increases, the numerical solution approaches the habitat carrying capacity and the occupied habitat tends to broaden.
The convergence rate of the solution has also been computed, which is defined here as the quantity , 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 .

Example 2. Let us consider the same model as in the previous example, with identical parameters values, initial condition, kernel function , number of nodes , and time horizon . If we choose a temporal step size , breaking the stability condition (27), it is clear from Figure 2 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 = . The initial condition and the kernel function are taken as
Note that, with this data, the problem presents radial symmetry. The spatial and temporal step sizes are chosen as and , respectively. The number of nodes of the quadrature is large enough to apply results of Section 3 because the maximum of does not exceed value 0.022 and assumption (19) is verified for . Applying (26), it results that the positivity and stability are guaranteed when , which is satisfied here. Figure 3 shows the numerical solution at the time horizon .
Example 4. Taking the same data as in Example 3, and identical number of nodes , if we choose a temporal step size , 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.
Acknowledgments
This work has been partially supported by the Ministerio de Economía y Competitividad Spanish grant MTM201789664P.
References
 M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions, with Formulas, Graphs, and Mathematical Tables, Dover, New York, NY, USA, 1972. View at: MathSciNet
 A. C. Ahlin, “On error bounds for Gaussian cubature,” SIAM Review, vol. 4, pp. 25–39, 1962. View at: Publisher Site  Google Scholar  MathSciNet
 D. G. Aronson and H. F. Weinberger, “Multidimensional nonlinear diffusion arising in population genetics,” Advances in Mathematics, vol. 30, no. 1, pp. 33–76, 1978. View at: Publisher Site  Google Scholar  MathSciNet
 H. Berestycki, G. Nadin, B. Perthame, and L. Ryzhik, “The nonlocal FisherKPP equation: travelling waves and steady states,” Nonlinearity, vol. 22, no. 12, pp. 2813–2844, 2009. View at: Publisher Site  Google Scholar  MathSciNet
 P. J. Davis and P. Rabinowitz, Methods of Numerical Integration, Academic Press, New York, NY, USA, 2nd edition, 1984. View at: MathSciNet
 G. M. Edelman and J. A. Gally, “Degeneracy and complexity in biological systems,” Proceedings of the National Acadamy of Sciences of the United States of America, vol. 98, no. 24, pp. 13763–13768, 2001. View at: Publisher Site  Google Scholar
 F. FakharIzadi and M. Dehghan, “An efficient pseudospectral LegendreGalerkin method for solving a nonlinear partial integrodifferential equation arising in population dynamics,” Mathematical Methods in the Applied Sciences, vol. 36, no. 12, pp. 1485–1511, 2013. View at: Publisher Site  Google Scholar  MathSciNet
 R. A. Fisher, “The wave of advance of advantageous genes,” Annals of Eugenics, vol. 7, no. 4, pp. 355–369, 1937. View at: Publisher Site  Google Scholar
 J. Furter and M. Grinfeld, “Local vs. nonlocal interactions in population dynamics,” Journal of Mathematical Biology, vol. 27, no. 1, pp. 65–80, 1989. View at: Publisher Site  Google Scholar  MathSciNet
 S. Genieys, V. Volpert, and P. Auger, “Pattern and waves for a model in population dynamics with nonlocal consumption of resources,” Mathematical Modelling of Natural Phenomena, vol. 1, no. 1, pp. 65–82, 2006. View at: Publisher Site  Google Scholar  MathSciNet
 S. Genieys, N. Bessonov, and V. Volpert, “Mathematical model of evolutionary branching,” Mathematical and Computer Modelling, vol. 49, no. 1112, pp. 2109–2115, 2009. View at: Publisher Site  Google Scholar  MathSciNet
 F. Hamel and L. Ryzhik, “On the nonlocal FisherKPP equation: steady states, spreading speed and global bounds,” Nonlinearity, vol. 27, no. 11, pp. 2735–2753, 2014. View at: Publisher Site  Google Scholar  MathSciNet
 A. Kolmogorov, I. Petrovskii, and N. Piskunov, “Study of a diffusion equation that is related to the growth of a quality of matter and its application to a biological problem,” Moscow University Mathematics Bulletin, vol. 1, pp. 1–26, 1937. View at: Google Scholar
 A. Okubo, Diffusion and Ecological Problem: Mathematical Models, vol. 10 of Biomathematics, Springer, Berlin, Germany, 1980. View at: MathSciNet
 N. Shigesada and K. Kawasaki, Biological Invasions: Theory and Practice, Oxford Series in Ecology and Evolution, Oxford University Press, Oxford, UK, 1997.
 E. Shivanian, “Analysis of Meshless Local Radial Point Interpolation (MLRPI) on a nonlinear partial integrodifferential equation arising in population dynamics,” Engineering Analysis with Boundary Elements, vol. 37, no. 12, pp. 1693–1702, 2013. View at: Publisher Site  Google Scholar  MathSciNet
 G. D. Smith, Numerical Solution of Partial Differential Equations: Finite Difference Methods, Oxford University Press, Oxford, UK, 3rd edition, 1985. View at: MathSciNet
 C. Tian, Z. Ling, and L. Zhang, “Nonlocal interaction driven pattern formation in a preypredator model,” Applied Mathematics and Computation, vol. 308, pp. 73–83, 2017. View at: Publisher Site  Google Scholar  MathSciNet
 N. Apreutesei, N. Bessonov, V. Volpert, and V. Vougalter, “Spatial structures and generalized travelling waves for an integrodifferential equation,” Discrete and Continuous Dynamical Systems  Series B, vol. 13, no. 3, pp. 537–557, 2010. View at: Publisher Site  Google Scholar  MathSciNet
 H. F. Weinberger, “On spreading speeds and traveling waves for growth and migration models in a periodic habitat,” Journal of Mathematical Biology, vol. 45, no. 6, pp. 511–548, 2002. View at: Publisher Site  Google Scholar  MathSciNet
 H. F. Weinberger, M. A. Lewis, and B. Li, “Anomalous spreading speeds of cooperative recursion systems,” Journal of Mathematical Biology, vol. 55, no. 2, pp. 207–222, 2007. View at: Publisher Site  Google Scholar  MathSciNet
Copyright
Copyright © 2019 M.A. Piqueras et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.