Existence and Uniqueness of Positive and Bounded Solutions of a Discrete Population Model with Fractional Dynamics

We depart from the well-known one-dimensional Fisher’s equation from population dynamics and consider an extension of this model using Riesz fractional derivatives in space. Positive and bounded initial-boundary data are imposed on a closed and bounded domain, and a fully discrete form of this fractional initial-boundary-value problem is provided next using fractional centered differences. The fully discrete population model is implicit and linear, so a convenient vector representation is readily derived. Under suitable conditions, the matrix representing the implicit problem is an inverse-positive matrix. Using this fact, we establish that the discrete population model is capable of preserving the positivity and the boundedness of the discrete initial-boundary conditions. Moreover, the computational solubility of the discrete model is tackled in the closing remarks.


Introduction
The development in recent decades of fractional calculus has led to important discoveries in many scientific areas [1,2].For example, research in the physical sciences has developed toward the construction of a physically meaningful calculus of variations for fractional systems [3,4].Here, the problem lies in the fact that the notions of "energy" and "Hamiltonian" still lack a concrete physical connotation, so the determination of the physical significance of those concepts is still a fruitful area of study [5].In the field of differential/difference equations, the determination of theorems on the existence and the uniqueness of solutions of fractional systems is nowadays an area of analytic importance [6].In this context, the determination of the properties of the relevant solutions of fractional systems is also a transited avenue of research, with the conditions of positivity and boundedness being of particular interest [7].From a numerical perspective, the design of computational techniques with desirable numerical properties (convergence, stability, consistency, etc.) is also a scientific problem that has attracted the attention of researchers in the area [8].However, the development of numerical methods that are capable of preserving structural properties of their continuous counterparts (the positivity, the boundedness, the monotonicity, or the preservation of energy, mass, momentum, etc.) is still scarce.
One of the most studied partial differential equations in the literature is the well-known diffusion-reaction model investigated simultaneously and independently in 1937 by Fisher [9] and Kolmogorov et al. [10].That model possesses positive, bounded, and traveling-wave solutions, and it was used initially to describe the propagation of mutant genes that are advantageous to the survival of populations distributed in linear habitats [11].Fisher's equation has also found applications in the investigation of the neutron flux and temperature in prompt feedback nuclear reactors [12], where the governing law is the Pearl-Verhulst equation.As many other equations from the physical sciences, that model has also been extended to the fractional scenario in order to describe more accurately the phenomena of interest (see [13] and the references therein).A vast number of criteria for fractional differentiation have been employed to that end, but some of these extensions are capable of resembling the properties of the solutions of the model investigated in 1937, like the existence and uniqueness of traveling-wave solutions that are positive and bounded [14].
On the other hand, the exact determination of solutions of fractional extensions of Fisher's equation is a difficult task.
To solve this problem, various computational approaches have been designed in order to approximate the solutions of continuous (fractional or integral) Fisher's equations.Some approximation techniques have been proposed in the literature to that effect, including finite-difference schemes [15], differential and integral quadrature techniques [16,17], Legendre spectral collocation methods [18], discrete local discontinuous Galerkin methods [19], finite volume schemes with preconditioned Lanczos techniques [20], and homotopy perturbation methods [21].In most of the cases, the methods proposed are capable of approximating the solutions of the fractional Fisher's equation with a high degree of precision, but the preservation of the most important features of the solutions of interest is neglected.The present manuscript is motivated by a Riesz space-fractional extension of Fisher's equation for which traveling-wave solutions that are positive and bounded exist.Full discretization that uses fractional centered differences is considered in this work and we show that, as its continuous counterpart, the discrete model possesses solutions that are positive and bounded.More precisely, we provide a dynamically consistent fully discrete form of our fractional model in the spirit of many of the works by Mickens and coworkers [22].
The present note is sectioned as follows.In Section 2, we introduce the fractional partial differential equation that motivates our investigation.The model is an extension of Fisher's equation from population dynamics in which the diffusion term is expressed as a Riesz fractional derivative.Both the continuous and the discrete Riesz fractional operators are presented in that section, together with some useful lemmas quoted from the standard literature.The discrete fractional Fisher's model is presented in Section 3. The discrete nomenclature and pertinent initial-boundary conditions are given therein as well as a vector form of the discrete model.In Section 4, we establish the most important properties of the discrete fractional Fisher's equation, namely, the existence and the uniqueness of positive and bounded solutions.The cornerstone of the proof is the concept of Minkowski matrices, and we recall this notion and its main properties in that section.This note closes with a section of concluding remarks and directions of future research.
For the sake of convenience, we will define (, ) = 0 for each  ∈ (−∞, ) ∪ (, ∞) and each  ∈ [0, ].Using this convention, the space-fractional operator in (1) is the Riesz fractional derivative of order , which is given by for each (, ) ∈ Ω.Here, Γ is the gamma function, which is defined in R \ { : − ∈ N or  = 0} by Note that the partial differential equation of ( 1) is actually a space-fractional extension of the model from population dynamics investigated simultaneously and independently by Fisher [9] and Kolmogorov et al. [10] in 1937.
In this work, we follow a fractional difference approach to approximate the solutions of (1) and use fractional centered differences to approximate Riesz space-fractional derivatives.For any function  : R → R, any ℎ > 0, and any  > −1, the fractional centered difference of order  of  at the point  is defined as where , ∀ ∈ Z. (5) In the case that 1 <  ≤ 2, the fractional centered differences satisfy (see [23]) Lemma 1 (see [24]).If 1 <  < 2, then the coefficients Lemma 2 (see [24]).Let  ∈ C 5 (R) and assume that all its derivatives up to order five belong to
(ii) Also, note that the fully discrete fractional model (11) is computationally an implicit two-step discretization of the continuous system (1).
Note that the discrete model (11) has a convenient vector representation.To describe it, let us consider the mesh and let V ℎ be the real vector space of all grid functions on R ℎ , that is, real functions V defined on R ℎ .For each V ∈ V ℎ and each  ∈ {0, 1, . . ., }, we convey that V  = V(  ).For each  ∈ {0, 1, . . ., }, let   ∈ V ℎ be the vector whose components are    with  ∈ {0, 1, . . ., }.In what follows, we let  = /ℎ  and define For the sake of briefness, in the following, we will obviate the dependence of ĝ,  0, on  and .Recall from Lemma 1 that    =   − for each  ≥ 1.For each  ∈ {0, 1, . . .,  − 1}, define the matrix    of size ( + 1) × ( + 1) and the ( + 1)dimensional real vectors û and Φ by Under these circumstances, the fully discrete population model (11) with the set of initial-boundary data (Φ, , ) may be expressed equivalently as the vector problem

Properties
For convenience, we denote the ranges of , , and  by R  , R  , and R  , respectively.Recall that a matrix  is nonnegative if every entry of  is a nonnegative number, a fact that is denoted by  ≥ 0. If  is any real number, we will say that  is bounded from above by  if every entry of  is less than or equal to , a fact that will be represented by  ≤ .
Obviously, an -dimensional real vector V satisfies V ≤  if and only if  − V ≥ 0, where  is the -dimensional vector whose components are all equal to 1.
A real square matrix  is a -matrix if all its off-diagonal entries are less than or equal to zero.We say that  is a Minkowski matrix if the following are satisfied: (i)  is a -matrix.
(ii) All the diagonal entries of  are positive.
(iii) There exists a diagonal matrix  with positive diagonal elements, such that  is strictly diagonally dominant.
Obviously, if  is strictly diagonally dominant, then condition (iii) is trivially satisfied with  equal to the identity.
The following result is part of the standard literature.
We provide an additional algebraic property of the matrices    .To that end, consider  = R +1 as a real vector space with the usual operations and denote the vector space of linear operators on R +1 by hom R (, ).Lemma 6.Let  < 1 and let  :  → hom R (, ) be given by  () () =   , ∀,  ∈ . ( Then,  is a linear transformation such that, for each , V,  ∈ , if 0 ≤  ≤ V ≤  and  ≥ 0, then 0 ≤   () ≤  V ().
The following is the most important result of this work.
Proof.Note that  0 satisfies the conclusion by hypothesis, so assume that it is true for some  ∈ {0, 1, . . .,  − 1}.By Lemmas 4 and 5, the matrix    is nonsingular and the entries of its inverse are positive numbers.Moreover, û ≥ 0 by the induction hypothesis and the assumptions on the boundary conditions.It follows that and we only need to show that  +1 ≤  or, equivalently, that V +1 =  −  +1 ≥ 0, where  is the vector of the same dimension as  +1 , all of whose entries are equal to 1.Note that the recursive identity in ( 15) is equivalent to    V +1 =   , where Clearly, the first and the last components of   are nonnegative.Meanwhile, after a simplification, one may readily check that the ( + 1)st component of   is of the form for each  ∈ {1, . . .,  − 1}.It follows that   ≥ 0, and the fact that    is a Minkowski matrix yields that V +1 ≥ 0 or, equivalently, that  +1 ≤ .The result follows now by induction.
From a practical point of view, Theorem 7 establishes that the fully discrete model (11) possesses unique solutions which may be interpreted as population densities when  = 1.Indeed, the theorem guarantees that, for any initial-boundary data which are nonnegative and bounded from above by , there exists a unique solution of (11) with the same properties of nonnegativity and boundedness.In that sense, our fully discrete model may be an ideal paradigm to describe the dynamics of variables measured in absolute scales which are normalized with respect to some maximum value.These features are in agreement also with respect to nonfractional models of science and engineering such as Fisher's equation and various of its generalizations.
The following example is presented to illustrate the capability of the fully discrete system (11) to preserve the nonnegativity and the boundedness of solutions.
Computationally, let ℎ = 1,  = 0.01.Under these circumstances, Figure 2 provides snapshots of the solutions of the fully discrete model (11) at various times for several values of , namely, 2 (solid), 1.9 (dashed), 1.8 (dasheddotted), and 1.7 (dotted).The graphs show that the diffusion increases when  decreases.This is a typical behavior of partial differential equations with fractional diffusion and is in full agreement with results obtained in the literature [26].
Before we close this section, it is important to note that the matrices   are Hermitian and positive-definite under the hypotheses of the theorem, so they have Cholesky factorizations of the form  =  * , where  is a lower triangular matrix of the same size as  and  * is its conjugate transpose.This means in particular that the discrete population model ( 11) may be efficiently implemented using the Cholesky decomposition method.It is important to point out that the fully discrete model (11) has been implemented computationally using this technique.Some of the results are shown in Example 8. We have conducted more numerical experiments which are not presented here in view that they also confirm the preservation of the properties of positivity and boundedness of the solutions

Conclusions and Perspectives
In this note, we considered a fully discrete model associated with a diffusion-reaction partial differential equation with Riesz space-fractional derivatives which generalizes Fisher's equation from population dynamics and nuclear physics.The spatial derivative orders are in the interval (1, 2), and the spatial discretization of the model uses fractional centered differences, which are consistent discretizations of the Riesz fractional derivatives.The fully discrete population model is implicit and linear, and it can be represented in vector form using an -matrix under suitable parameter conditions.As a consequence, the scheme is always solvable, and the properties of -matrices guarantee the positivity and the boundedness of the solutions.We calculated the solution of some initial-boundary-value problems governed by the discrete fractional Fisher's equation through an implementation of the Cholesky decomposition algorithm.The results are not provided in this work, but they confirm the capacity of the model to preserve the positivity and the boundedness when the initial-boundary data satisfy those properties.
Various avenues of research open after the conclusion of this manuscript.For example, it would be interesting to determine whether theorems on the existence and the uniqueness of positive and bounded solutions for more general discrete fractional systems can be established.There are many generalizations of Fisher's equation in the context of population dynamics, and various approaches may be followed in order to provide fractional extensions for those models.An interesting question that arises is whether those discrete systems accept solutions with the characteristics of interest.In that sense, the present note is an affirmative step toward the structural analysis of discrete fractional systems arising in population dynamics.Once we have established relevant theorems on the solutions of the discrete fractional Fisher's equation, another interesting avenue of research is the numerical analysis of structure-preserving methods to approximate continuous fractional models.The literature already has some relevant properties on the consistency of fractional centered differences, but the design of structurepreserving methods for fractional partial differential equations is still an area that is yet to be explored.The topic is outside the scope of this journal, but it is a problem that deserves a thorough numerical investigation.

Figure 1 :
Figure 1: Stencil of the fully discrete model(11).The black circles represent the known solutions at the time   and at the boundary, while the crosses denote the unknown approximations at the time  +1 .