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 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 .
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 MTM2017-89664-P.