Abstract and Applied Analysis

Abstract and Applied Analysis / 2019 / Article

Research Article | Open Access

Volume 2019 |Article ID 5787329 | https://doi.org/10.1155/2019/5787329

M.-A. Piqueras, R. Company, L. Jódar, "Stable Numerical Solutions Preserving Qualitative Properties of Nonlocal Biological Dynamic Problems", Abstract and Applied Analysis, vol. 2019, Article ID 5787329, 7 pages, 2019. https://doi.org/10.1155/2019/5787329

Stable Numerical Solutions Preserving Qualitative Properties of Nonlocal Biological Dynamic Problems

Academic Editor: Maurizio Grasselli
Received26 Feb 2019
Accepted09 Jun 2019
Published01 Jul 2019


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 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]: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 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 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 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.

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. 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 havewhere and are the weights and and , , are the nodes of the Gauss-Hermite 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 Gauss-Legendre 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 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 formusing 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 . 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 Gauss-quadrature 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 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 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 .

=0.0005 =0.0040 =0.0020 =0.0010

0.0014733979 0.0006292592 0.0002095394
- 1.22742 1.58643

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.


This work has been partially supported by the Ministerio de Economía y Competitividad Spanish grant MTM2017-89664-P.


  1. 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
  2. A. C. Ahlin, “On error bounds for Gaussian cubature,” SIAM Review, vol. 4, pp. 25–39, 1962. View at: Publisher Site | Google Scholar | MathSciNet
  3. 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
  4. H. Berestycki, G. Nadin, B. Perthame, and L. Ryzhik, “The non-local Fisher-KPP equation: travelling waves and steady states,” Nonlinearity, vol. 22, no. 12, pp. 2813–2844, 2009. View at: Publisher Site | Google Scholar | MathSciNet
  5. P. J. Davis and P. Rabinowitz, Methods of Numerical Integration, Academic Press, New York, NY, USA, 2nd edition, 1984. View at: MathSciNet
  6. 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
  7. F. Fakhar-Izadi and M. Dehghan, “An efficient pseudo-spectral Legendre-Galerkin method for solving a nonlinear partial integro-differential 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
  8. 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
  9. J. Furter and M. Grinfeld, “Local vs. non-local interactions in population dynamics,” Journal of Mathematical Biology, vol. 27, no. 1, pp. 65–80, 1989. View at: Publisher Site | Google Scholar | MathSciNet
  10. 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
  11. S. Genieys, N. Bessonov, and V. Volpert, “Mathematical model of evolutionary branching,” Mathematical and Computer Modelling, vol. 49, no. 11-12, pp. 2109–2115, 2009. View at: Publisher Site | Google Scholar | MathSciNet
  12. F. Hamel and L. Ryzhik, “On the nonlocal Fisher-KPP equation: steady states, spreading speed and global bounds,” Nonlinearity, vol. 27, no. 11, pp. 2735–2753, 2014. View at: Publisher Site | Google Scholar | MathSciNet
  13. 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
  14. A. Okubo, Diffusion and Ecological Problem: Mathematical Models, vol. 10 of Biomathematics, Springer, Berlin, Germany, 1980. View at: MathSciNet
  15. N. Shigesada and K. Kawasaki, Biological Invasions: Theory and Practice, Oxford Series in Ecology and Evolution, Oxford University Press, Oxford, UK, 1997.
  16. E. Shivanian, “Analysis of Meshless Local Radial Point Interpolation (MLRPI) on a nonlinear partial integro-differential 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
  17. G. D. Smith, Numerical Solution of Partial Differential Equations: Finite Difference Methods, Oxford University Press, Oxford, UK, 3rd edition, 1985. View at: MathSciNet
  18. C. Tian, Z. Ling, and L. Zhang, “Nonlocal interaction driven pattern formation in a prey-predator model,” Applied Mathematics and Computation, vol. 308, pp. 73–83, 2017. View at: Publisher Site | Google Scholar | MathSciNet
  19. N. Apreutesei, N. Bessonov, V. Volpert, and V. Vougalter, “Spatial structures and generalized travelling waves for an integro-differential equation,” Discrete and Continuous Dynamical Systems - Series B, vol. 13, no. 3, pp. 537–557, 2010. View at: Publisher Site | Google Scholar | MathSciNet
  20. 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
  21. 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 © 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.

More related articles

 PDF Download Citation Citation
 Download other formatsMore
 Order printed copiesOrder

Related articles

Article of the Year Award: Outstanding research contributions of 2020, as selected by our Chief Editors. Read the winning articles.