A numerical method for space‐fractional diffusion models with mass‐conserving boundary conditions

The study of the spread of epidemics has gained significant attention in recent years, due to ongoing and recurring outbreaks of diseases such as COVID‐19, dengue, Ebola, and West Nile virus. In particular, modeling the spatial spread of these epidemics is crucial. This article explores the use of fractional diffusion as a means of describing non‐local infection spread. The Grünwald–Letnikov formulation of fractional diffusion is presented, along with several mass‐conserving boundary conditions, that is, we aim to design the boundary conditions in a mass‐conserving way, by not allowing gain or loss of the total population. The stationary points of the model for both sticky and reflecting boundary conditions are discussed, with numerical examples provided to illustrate the results. It is shown that reflecting boundary conditions are more reasonable, as the stationary point for sticky boundary conditions is infinite at the boundaries, while reflecting boundary conditions only have the trivial stationary point, given sufficiently fine discretization. The numerical results were applied to an SI model with fractional diffusion, highlighting the dependence of the system on the value of the fractional derivative. Results indicate that as the order of the derivative increases, the diffusivity also increases, accompanied by a slight increase in the average number of infected individuals. These models have the potential to provide valuable insights into the dynamics of disease spread and aid in the development of effective control strategies.


INTRODUCTION
The authors' aim is to find an apt representation about diffusion patterns in the spreading of infections with special kinds of diffusion patterns and especially to find reasonable boundary conditions for (future) epidemiological investigations. For a short time scale, we can consider countries as isolated, due to, for example, the geographical or political situation; in early 2020, for example, border closures prevented the majority of non-economical traffic in many parts of the world. Therefore, epidemics mostly spread only inside those countries. This allows a country-based investigation of the spread neglecting external factors. Spatial spread of epidemics can be modeled locally and non-locally. Local models can be gained by adding mobility matrices in several compartments (see, e.g., Goel et al. [1] and Schäfer and Götz [2]), traveler matrices (Schäfer et al. [3]), diffusion models or partial differential equation systems (Gai et al. [4] or Kuniya and Wang [5]). This can be described by the equation One issue is that the local context is disrupted by numerous sparsely populated areas, resulting in limited interaction between neighboring villages, but more important exchange between further apart villages. Thus, non-local models that include spread on higher distances can be a helpful tool. While one possibility of attaining non-local spreading effects are integro-differential equations (see, e.g., Kergaßner et al. [6], Kuniya and Wang [5], or Schäfer et al. [7]), the focus of this article is the usage of a fractional spatial derivative, that is, an equation of using a fractional diffusion operator Δ , so that Equation (1) transforms into t (x, t) = Δ (x, t) + ( (x, t)).
Fractional calculus, that is, calculus using non-integer powers of integration and derivation, is nowadays used in many applications, for example, network dynamics (cf. Baleanu and Kumar [8]), viscoelasticity (cf. Baleanu and Kumar [8] and Soczkiewicz [9]), field theory and gravity (cf. Baleanu and Kumar [8] and Calcagni [10]), advection-dispersion flow equations (cf. Meerschaert and Tadjeran [11]), diffusion processes in epidemiology (cf. Baleanu and Kumar [8] and Oldham and Spanier [12]), and price fluctuations in financial markets (cf. Baleanu et al. [13]. Specifically, in modeling of epidemics, both space-and time-fractional models are used, specifically time-fractional models: Rezapour et al. [14] have set up an SEIR model for transmission of COVID-19). Goufo et al. [15] investigated an time-fractional SEIR model concerning measles. Sidi Ammi et al. [16] provided some useful theoretical results concerning time-fractional SIR models. Kuehn and Mölter [17] also investigate transport effects on epidemics using two coupled models, a static epidemic network and a dynamical transport network, also with non-local, fractional transport dynamics. Hamdan and Kilicman [18] use a system of fractional-order derivative systems to model and analyze control strategies for dengue fever epidemics in Malaysia. In this article, we will consider aspects of space-fractional models which have been, for example, considered by Gorenflo and Mainardi [19]. One can show that space-fractional diffusion corresponds to a change from a Markovian random walk to a non-Markovian random walk with Levy flights; compare Oldham and Spanier [12] and Skwara et al. [20].

Structure of the paper
In the methods section, we provide an outline of the definition of fractional derivatives as far as they are relevant for the subsequent model(s). Furthermore, we introduce a finite difference approximation of fractional derivatives that allows numerical simulations. Two different kinds of boundary conditions are introduced, and theoretical results for approximations and convergence of the schemes are given. Further on, numerical results for some one-dimensional examples using both types of boundary conditions are provided. Lastly, both theoretical and numerical results are discussed.

Fractional derivatives
To generalize the concept of derivatives from integer to non-integer order, there exist several concepts. In our work, we rely on the definitions introduced by Riemann and Liouville and the according discretization given by Grünwald and Letnikov [8][9][10][11]. Let ∶ [a, b] → R denote a sufficiently smooth function on an interval [a, b], where a = −∞ and b = ∞ can be included. Let Γ(n) = ∫ ∞ 0 t n−1 e −t dt = (n − 1)! denote the gamma function for n ∈ N + . Then, its right-sided n-fold repeated integral is given by (cf. Tadjeran et al. [21] and Tadjeran and Meerschaert [22]) We can straightforwardly extend this to non-integer order ∈ R∖ {−1, −2, … } by The above integral is called the right-sided fractional Riemann-Liouville integral of order . Analogously, we may define the left-sided fractional integral by Remark 1. It is easy to prove (see Meerschaert and Tadjeran [11]) that for the fractional Riemann-Liouville integral, it holds According to the fundamental theorem of calculus, we identify the derivative of order k ∈ N as the inverse operators to the integrals of order k, that is, To define the fractional derivative of order > 0, let n = ⌈ ⌉ = min k∈N {k > } denote the smallest integer larger than . Formally, taking the n-derivative of the fractional Riemann-Liouville integral of order n − , that is, d n dx n I (n− ) yields the Riemann-Liouville definition of the fractional derivative of order .
Definition 1 (Left-and Right-sided Riemann-Liouville derivatives). Let ∈  n [a, b]. We call the Riemann-Liouville derivative of fractional order with left boundary a and the Riemann-Liouville derivative of fractional order with right boundary b (cf. Soczkiewicz [9]).
Remark 2. For this definition, requiring ∈  n−1 ([a, b]) would actually be too strong; for example, for ∈ [1, 2), the function would only need to be absolutely continuous, that is, ∈ ([a, b]); for this, compare also Salem [23]. However, for the proof of the mass conservation, we require differentiability, so that we will keep this condition.
Remark 3. The above fractional derivatives are linear, that is, for all , g ∈  n [a, b], all , ∈ R, and any > 0, we have Example 1. We consider the monomial x k and compute its fractional derivative of order . In alignment with the formal generalization of the integer-order derivative we obtain and analogously, For a = 0 and = n ∈ N, this corresponds to the well-known nth derivatives of the monomials.

Example 2.
As a second example, we consider the exponential (x) = e cx for c > 0 and x ∈ [0, ∞). Using the series representation of the exponential and the above right-sided fractional derivative of the monomials, we obtain After some algebraic manipulations and by introducing the incomplete gamma function Γ(− , cx) = ∫ ∞ cx t − −1 e −t dt, we arrive at Again, for = n ∈ N, this corresponds to the regular nth derivative of the exponential function.

Sticky boundary conditions
As a next step, we define the fractional diffusion operator Δ as a linear combination of the right-sided and the left-sided fractional derivative. For the sake of simplicity, we restrict ourselves on the spatial one-dimensional case on [a, b] with ≡ 0. In any non-local fractional diffusion process, sticky boundary conditions imply that in a closed domain [a, b], mass that leaves a certain point x and, without the boundaries, would come to rest at a place < a or > b, is assumed to remain in either boundary a or b (cf. Borges et al. [24]). This kind of boundary condition naturally meets the mass conservation rule for ≡ 0 as no mass can leave (or enter) the domain. The fractional diffusion equation for an order of diffusion with 1 < ≤ 2 then reads as the sum of both right-sided and left-sided fractional derivatives: Remark 4. The sticky boundary conditions are already defined over the definition of the fractional derivatives in the domain, that is, In order to prevent confusion in the notation, [a, b] = [−1, 1] will be our standard interval from now onwards. Because of this, we can simplify the notation by D ( ) Example 3. Note that for a symmetrical diffusion operator with > 0, the fractional derivatives of the constant c (cf. Example 1) do not vanish; in fact, they are not even constant: We can show that the PDE (15) is mass conservative: Given the fractional PDE (15) with initial condition 0 (x) and an ∈ [1,2). Define as the mass of the at time t. Then, dM Proof. The result can be easily attained by deriving Equation (17): □

Proposition 2. The only stationary solution of the equation is
Proof. To determine stationary solutions of the fractional diffusion, we aim to solve the equation Using the definition of the fractional derivatives, this reads as This implies that the integral ∫ Symmetry of the solution requires c 1 = 0; hence, the integral has to equal the constant c 0 .
The work of Shinbrot [25, pp. 23-24] provides a solution for equations of the type for an arbitrary ∈ [0, 1] is provided (which, for our problem, is guaranteed by the choice = − 1). Defining coefficients and functions the solution of the above integral equation is given by the series In our case, g(x) ≡ c 0 . Thus, The series term vanishes, since for n ≥ 1 all the terms (1 − x 2 ) n+ +1 2 or its integer derivatives vanish at the boundaries. Hence, where c 0 is a constant defined by the total initial mass. Also, notice that for = 2, the constant reduces to c 0 as . For all 1 < < 2, that is, 0 < = − 1 < 1, it holds that −1 2 < 0 and the solution will tend to infinity at the boundary x = ±1. □

Reflecting boundary conditions
As it can be observed, the sticky boundary conditions satisfy the mass conservation requirement; however, they do not provide a physically or biologically accurate description of diffusion scenarios. To avoid the singularities mentioned above, we instead assume that mass moving into the boundary is reflected and moves back in the opposite direction. To comply with these reflecting boundary conditions, we consider a periodically extended version of , that is, While is defined on whole R, we only consider fractional diffusion in our closed domain, that is, for x ∈ [−1, 1]. This is then denoted by the following equation: .
In fact, the description as "boundary condition" is rather misleading, as the equations have to be updated for all grid points in this case.

Proposition 3. The only stationary solution is the constant, that is
Proof. Stationary equilibria * (x) are solutions of the equation .
Shinbrot [25,Corollary,p. 24] shows that equation has at most one solution given in the form. We show that * (x) = c 0 (and thus (x) ≡ c 0 ) solves the problem; hence, it is the unique solution: To show uniqueness of the (trivial) solution, we assume * (x) ≢ c 0 solves Equation (27). Then, due to similar thoughts as for the sticky boundary conditions (symmetry of the solution), we know that We reconsider Shinbrot [25]: Construct a Fourier transform function n (x). Then, the solution can be explicitly calculated by Fourier transforms using functions n (x), n ∈ N: The definition of all n (x) is unique, so that (x) * ≡ c 0 which contradicts the assumption. Due to mass conservation, (x) ≡ c 0 thus is the unique solution for the problem with reflecting boundary conditions. □ Additionally, due to the definition of the mass transport, the mass conservation result in (26) is conveyed to the case of reflecting boundary conditions.

NUMERICAL METHODS
In order to derive a finite difference approximation to the fractional derivative, we generalize the standard difference approximation to fractional orders using (−1) k . This leads us to the following: Definition 2 ( cf. Letnikov [26]). Let ∈ C n−1 [a, b] and ∈ [n−1, n). The right-shifted Grünwald-Letnikov difference formulation to the fractional derivative of order is given by The left-shifted Grünwald-Letnikov difference formulation to the fractional derivative of order is given by In order to simplify notations, we introduce the Grünwald coefficient for k ∈ N 0 and > 0 (cf. Baeumer et al. [27]).
Remark 5. In Meerschaert and Sikorskii [28], it is shown that the sum of all Grünwald coefficients vanishes, that is, ∑ ∞ k=0 g k = 0, so that the method is conservative. In order to approximate Equation (15), assume that we have found a numerical approximation, that is, using a discrete operator Δ ,h which is yet to be found. Then, we can use a standard ODE solver, for example, the explicit or implicit Euler system or the Runge-Kutta method, to find the numerical solution of this equation for h > 0 and The discrete operator Δ ,h will now be adapted to the two different boundary conditions.

Sticky boundary conditions
All mass that moves out of the boundary is assumed to find a rest at the boundary, so the mass from x that would be moved to x − (k − 1)h < a will be moved to a instead. Analogously, the mass from x that would be moved to x + (k − 1)h > b will be moved to b instead. The same thing happens for the transportation of mass from x = a to a + (k − 1)h and k = 0 for the right-sided fractional derivative and the transportation of mass from x = b to b − (k − 1)h and k = 0 for the left-sided fractional derivative. As (x) = 0 for x < x 0 ∶= a and x > x m ∶= b, the sum in Equations (33) and (34) becomes finite.
, respectively. It is assumed that a part of the mass at value x − (k − 1)h is moved to x, namely, g k h − (x − (k − 1)h). This results into the following discrete formulations: Also, terms 1− result of the boundary conditions and the formulations for the right-sided fractional derivative at x = a and the left-sided fractional derivative at x = b. For the boundary points, we have to update the corresponding Grünwald coefficients according to the absorption of mass: For an equidistant grid with x 0 = a, x m = b and x i = x 0 + ih, i = 0, … , n, h = 1 m+1 , the discrete operator then reads as with transition matrices for the right-sided fractional diffusion, that is, as ∑ ∞ p=0 g p = 0, and the left-sided fractional diffusion, that is, Remark 6. It can be shown (cf. Meerschaert and Tadjeran [29]) that this shifted method used is first-order consistent, so that the error terms r 1 ,

Reflecting boundary conditions
For the numerical scheme, consider the discretization on an equidistant grid consisting of m + 1 points x 0 = a, x m = b and When the mass has moved m steps to the left and k + 1 steps are actually needed to be covered, then the mass is transported from the left boundary to the right for the remaining k + 1 − m steps. Same holds true for the right boundary. As the Grünwald coefficient g k → 0 for k → ∞ but remains nonzero for finite k, we have an infinite process where a decreasing mass is transported infinitely often between the left and right boundaries. Note that there are two possible ways to design the mirror: One would be that if mass reaches the boundary and "has yet one step to cover," then it can either return into the other direction for one step or remain at the boundary and only moves on into the other direction when there is more than one step yet to cover. The first method is used in this paper to achieve a certain symmetry between all grid points. We know that the Grünwald-Letnikov formulation is consistent of order h for reflecting boundary conditions and the sum from 0 to x−a h + 1 or b−x h + 1, respectively, because the function values vanish for points outside the boundary. For h > 0, it holds: For the left-sided fractional derivative and arbitrary k, the mass moving from x i to x is considered. It first moves i + 2 steps to the left boundary a, then one step within the left boundary, then steps to the point x , making it i + + 3 steps in total. Because the process is repeated every 2(m + 1) steps, all elements i + + 3 + 2(m + 1)l, l ∈ N 0 , have to be summed up as well. After the i + + 3 steps, it also moves m − points to the right boundary b, another point inside the right boundary and then m − points back to point x , making it i − + 2m + 4 points in total, similar to above it follows that mass is moved from x i to x for every distance i − + 2m + 4 + 2(m + 1)l, l ∈ N 0 ; compare also Figure 1.
For the right-sided fractional derivative, for an arbitrary (but large enough) k, the mass moving from x i to x is considered. It first moves m − i + 2 steps to the right boundary b, then one step within the right boundary, then m − steps to the point x , making it −i − + 2m + 2 steps in total. Thus, as described above, all elements −i − + 2m + 2 + 2(m + 1)l, l ∈ N 0 , have to be summed up as well. After the 2m + 2 − i − steps, it also moves points to the left boundary a, another point inside the left boundary, and then points back to point x , making it − i + 2m + 3 points in total; similar to above, it follows that mass is moved from i to for every distance − i + 2m + 3 + 2(m + 1)l, l ∈ N 0 ; compare also Figure 2.  Until now, the actual position of x i and x has not been regarded. We have to consider the possibilities x i < x , x i = x , and x i > x . These can easily be condensed if we consider a modified modulo formula: mod(a, n) mod(a, n) ≠ 0, n mod(a, n) = 0.
A problem lies in the the infinite number of summands in each formula. However, as we know that g k → 0 for k → ∞, reconsider the formulation of the Grünwald coefficient as of Equation (35). If for an arbitrary, but fixed N ∈ N, the partial sum ∑ ∞ k=N g k = ∑ N−1 k=0 g k can be estimated, we can find a decent approximation by the apt choice of N after which the sum of the Grünwald coefficients is stopped.

Proposition 4. It holds
Proof. The term Γ(− ) is independent of k and thus can be seen as constant. The term Γ(k + 1) is equal to k! as k is integer. The term Γ(k − ) can be represented by Altogether, this implies For the sum, it then holds If we now set N = (h −p ), it follows ∑ ∞ k=N g k = (h p ). □ We can even improve this estimation by computing the logarithm of the Grünwald coefficient.
Proof. For that, we use the approximation by Rocktäschel [30]: We take the natural logarithm of the non-constant terms of the Grünwald coefficient and apply this approximation: Additionally, we can approximate This way, we have In total, we have found This implies and thus, the better approximation which by a choice of N = (h −p ) yields Additionally, the equations are subject to reflecting boundary conditions, as discussed earlier. As s = 1 − z, the only remaining equation is The results of the SI model can be seen in the visually appealing figures, which depict the evolution of the infected population over time. These figures should then effectively demonstrate the effects of using fractional partial differential equations instead of classical diffusion models.

Sticky boundary conditions
Numerical examples are provided for (x) ≡ 1, = 10 −2 and various . The border points at x = ±1 can be evaluated numerically at any time t < ∞, but as t → ∞, the value will diverge to ∞ as seen in Lemma 2. The analytical solution of the stationary solution as of (5) is compared to the numerical results with applied reflecting boundary conditions. The amount of runs of the algorithm is terminated by the  2 -norm. In Figures 3-5, it is visible that the singularity causes a larger discretization error near the boundaries, such that a comparison of the norms is not useful.
This equilibrium is independent of starting values of the function, given the total mass is the same. Also, the singularity at the boundaries is not what we would expect for our problem (infectives in Sri Lanka), so reflecting boundaries are not a good method for fractional diffusion in our case.

Reflecting boundary conditions
Numerical examples are provided for the domain Ω = [−1, 1] ⊂ R: These numerical solutions support the theoretical findings that the system with mirrored boundary condition runs into the trivial equilibrium (x) ≡ c 0 , where c 0 is the mean value of the initial mass. The "velocity" of diffusion increases with the value of the fractional derivative .

SI model
The results of the fractional SI model as of Equation (62) for = 0.1, = 0.01, and different values of are shown in the following. The initial distribution of infected z 0 (x) is set as Figures 8-12 show the spatiotemporal dynamics of the infectives at time t. As the value of , which represents the fractional order of the derivative, increases, the average number of infected individuals also increases, although only by a slight amount. There is minimal change observed in this aspect. However, the infection spreads more rapidly with increasing value of and therefore prevents a high percentage of infected individuals in the region where the disease first occurred.

DISCUSSION AND CONCLUSION
This article deals with the numerical approximation of space-fractional derivatives, which can be used to solve space-fractional PDE systems which have many applications. We aimed to find a representation with respect to disease spread to be used in future works. Using the finite difference approximation to the Riemann-Liouville definition, the Grünwald-Letnikov formulation provides an easy to discretize the derivative of fractional order. In order to conserve the mass, that is, the total population in a country, and assuming that there is no exchange of population to other countries or regions, we considered two different boundary conditions; using these, we can also handle the non-local nature of the fractional diffusion. For sticky boundary conditions, where all mass moving out of the boundary remains at the boundary, we use the property of the Grünwald coefficients ∑ k g k = 0 and a lemma of Shinbrot to calculated the stationary solution, which tends to infinity at the boundaries even if the integral is finite. This makes the usage of sticky boundary conditions a less realistic choice than reflecting boundary conditions, as by using von Neumann stability analysis, we can show the only stationary point is the trivial equilibrium. For these kind of "boundary" conditions that impose an infinite sum on each grid point, a finite formulation of order has been found. Numerical one-dimensional results for both formulations were applied to illustrate the scheme and the structure of the solutions which support the theoretical findings.
The application of numerical results to an SI model with fractional diffusion was performed to demonstrate the scheme and to reveal the dependence of the system behavior on the value of the fractional derivative . The results indicate that as the value of increases, the rate of diffusion increases accordingly. In the case of a regional epidemic, the fractional model predicts that with a rise in , the disease will spread rapidly to previously unaffected regions. However, it also suggests a decline in the number of cases in the original epicenter of the pandemic. On average, a slight increase in the number of infected individuals is observed as increases.
In the future, these models will be applied to describe specific epidemic (e.g., Dengue fever or COVID-19 outbreaks in a certain country or region) to validate the results of space-fractional or even time-fractional partial differential equation models to epidemiological data and better understand the behavior of disease spread in a specific area. These models have the potential to provide valuable insights into the dynamics of disease spread, which can be useful in developing control strategies. In order to meaningfully include fractional diffusion, we aim to make use of the fractional diffusion model in chapter 5.3 to describe recent infection waves. We aim to approximate the medical data by parameter estimation, which will not only be done for the classical transmission parameters, but also for the diffusion parameters and .