APPLICATIONS OF THE FRACTIONAL CALCULUS: ON A DISCRETIZATION OF FRACTIONAL DIFFUSION EQUATION IN ONE DIMENSION APPLICATIONS OF THE FRACTIONAL CALCULUS: ON A DISCRETIZATION OF FRACTIONAL DIFFUSION EQUATION IN ONE DIMENSION

The paper discusses the problem of classical and fractional diffusion models. It is known that the classical model fails in heterogeneous structures with locations where particles move at a large speed over a long distance. If we replace the second derivative in the space variable in the classical diffusion equation by a fractional derivative of order less than two, we obtain the fractional diffusion equation (FDE) which is more useful in this case. In this paper we introduce a discretization of FDE based on the theory of the difference fractional calculus and we sketch a basic numerical scheme of its solution. Finally, we present some examples comparing classical and fractional diffusion models.


Introduction
Fractional calculus is a mathematical discipline dealing with derivatives and integrals of non-integer orders. Its story started around 1695 when Leibniz first suggested the idea of half-derivative in his correspondence with l'Hospital. During centuries many mathematicians have contributed to the discussion on the "right" definition of a fractional derivative and various approaches have been proposed (see [13]). However, this concept of fractional derivatives and integrals has been still considered only as a pure theoretical matter because no real applications were known at that time.
The situation changed in last decades when many physical and engineering problems turned out to be described more precisely using mathematical tools of the fractional calculus. We focus on one of such phenomena -anomalous diffusion which is faster than the classical model predicts.
The main purpose of this paper is to obtain a numerical solution of the fractional diffusion equation (FDE). In general, considering the discrete approximations of the fractional differential equations, a lot of questions still remain open. One of them is the problem of nonzero initial or boundary conditions of the Riemann-Liouville type (i.e. when the values of fractional derivatives and integrals are given at point zero). Although we present only the example of the FDE with boundary conditions, we believe that the numerical scheme proposed in this paper allows to involve initial or boundary conditions of any type and any values.
The paper is organized as follows. In Section 2 we recall some basic notions of the continuous fractional calculus. Its discrete analogue is considered in the next section. Section 4 summarizes some known facts about the FDE and outlines the background of this model. In Section 5 we introduce a discretization of the FDE based on the theory presented in Section 3. Finally, we illustrate the method with several examples.

Fundamentals of the Fractional Calculus
As we foreshowed above, various definitions of a fractional derivative and integral appeared during the history of the fractional calculus and many of them are used up to this day. More information can be found in relevant monographs, e.g. [10], [11], [13], [14]. We are going to recall and employ the Riemann-Liouville approach which is the most frequent one.
Fractional derivatives are then usually constructed by the repeated differentiation of an appropriate fractional integral.
On this account we formulate the following definition.
Further, let m ʦ Z + be such that m Ϫ 1 Ͻ ␣ Յ m and let f(x) be m-times differentiable on [a, b) except for a set of the zero measure. Then we define the left fractional derivative of f(x) by The left fractional derivative and integral (jointly called the left differintegral) are obviously global operators employing only points from the left of the point x (if the variable x has the meaning of time, the operator depends on past). If x has the meaning of a space variable, there is no reason to prefer the left-side definition to the following right-side one.
Further, let m ʦ Z + , be such that m Ϫ 1 Ͻ ␣ Յ m and let f(x) be m-times differentiable on (a, b] except for a set of the zero measure. Then we define the right fractional derivative of f(x) by , x ʦ (a, b]. If we need to consider all points of the interval (a, b), we have to utilize both left and right differintegral. Sometimes it may be useful to employ the so-called Riesz fractional integral.
Properties of the Riemann-Liouville and Riesz differintegrals are discussed in more details in [10], [11], [13], [14]. Among numerous interesting properties of these operators we emphasize at least their linearity.

Discrete Fractional Calculus
Contrary to the continuous case, the theory of discrete fractional calculus is much less developed. Its concept starts from the time-scale theory, where the notions of delta and nabla differences are generalized to other discrete settings ( [7], [12]). We are not going to recall here related definitions of the time-scale theory (see [4]) because we work only with a special setting anyway, namely with a finite set of equidistant points with a stepsize h ʦ R + . Therefore we introduce the set of mesh points T ϭ {x n ϭ nh, n ϭ ϭ 0,1, … , N}, where the research of fractional differences is advanced (for the methods allowing to solve fractional difference equations see [1]).
In our conception, the notion of the h-backward (nabla) difference is a starting point for the left fractional difference operators. Note that some authors study the h-forward (delta) difference calculus. However, this approach causes that the fractional differences are defined at points which may not belong to the original mesh (see [1], [5]) and therefore it is not suitable from the numerical point of view.
Before we proceed to definitions of fractional differences and sums, we remind the definitions of nabla and delta differences on the mesh T. In the following definition we use the binomial coefficient where Further, let m ʦ Z + be such that m Ϫ 1 Ͻ ␣ Յ m and m Ͻ N. Then we define the left fractional difference of f at x n by , m Յ n Յ N.
Finally, we put 0 ٌ 0 n f n ϭ f n .
Since our mesh T is finite, it is meaningful to define the right fractional sum and difference. Because of the consistency with the previous Definition 3.2 we have to utilize the notation based on delta differences.
We define the right fractional sum of f at x n by , 0 Յ n Յ N.
Further, let m ʦ Z + be such that m Ϫ 1 Ͻ ␣ Յ m and m Ͻ N. Then we define the right fractional difference of f at x n by Finally, we put n ⌬ 0 N f n ϭ f n . Last we introduce a discrete analogy of the Riesz fractional integral.

The Origin of the Fractional Diffusion Equation
Roughly speaking, the diffusion is a process when particles spread from areas of a high concentration to areas of a low concentration. The classical diffusion model in one dimension is described by the well-known partial differential equation , where u(x,t) is an unknown concentration of particles and f(x,t) is a known source term. Besides, we need an initial condition and, moreover if the region is bounded, then also boundary conditions. However, this model is insufficient in the case of fast diffusion, i.e. in the process describing the spreading of particles according to different concentration, but with a significantly greater number of very fast particles. The fast diffusion occurs, e.g., in porous media and it is traditionally described by adding a nonlinearity into the equation (2).
The fractional calculus allows us to model fast diffusion too, but the equation remains still linear. The corresponding FDE has the form (see [2], [3], [14]) where we write the second derivative and the Riesz integral separately due to the discretization process which appears later. The order of differintegration ␣ allows us to control the speed of spreading. Note that the Riesz integral is associated with the space variable x and from the computational viewpoint the time variable t plays a role of a parameter.
From the statistical point of view the derivation of the FDE (3) is quite analogical to the classical case. We start from the usual random walk but in our case we remove the assumption concern-, , ing a finite variance of particle jumps. In other words, we accept the power-law behaviour for long jumps (sometimes called Levy flights). The process of statistical derivation directly implies the fundamental solution of the corresponding fractional Cauchy problem for (3), which coincides with Levy ␣-stable distributions (we emphasize that this solution cannot be expressed analytically and it is defined via its inverse Fourier transform). This is another interesting analogy because Levy ␣-stable distributions play in the generalized central limit theorem the same role as the gaussian distribution in the classical central limit theorem. For more information see e.g. [3], [14] or more detailed [2]. An alternative derivation utilizing the eulerian model based on finite volumes approach was proposed in [17].
Interpreting the fractional term in (3), we note that in the classical model we assume that the particle flux J(x, t) depends only on the neighbourhood of x and, therefore, is local. On the other hand, the fractional model deals with situations when particles move so fast that we cannot neglect even the effect of particles being far from a particular point. Hence, the flux has a global character and is proportional to the derivative of the Riesz (2 Ϫ ␣)-integral of the concentration u(x, t), i.e. .
There are other ways of the generalization of equation (2). For example, we can introduce weight coefficients of the left and the right differintegral for modeling anisotropic media, or we may consider a fractional derivative with respect to time on the left side of the equation. However, in this paper we are focused just on the case described above, in particular, we intend to discuss the numerical solution of the equation (3) on a bounded domain which is not so well established topic like e.g. the Cauchy problem. Now, we introduce the problem we are going to discuss. First, the studied equation should be slightly changed because the formula (3) was derived for an infinite line while we are going to deal with the compact interval [0, 1]. Hence, we assume that the concentration u(x, t) is zero everywhere except for x ʦ [0, 1] where the equation (3) holds. Consequently, we can change the bounds of the differintegration from Ϫϱ and ϱ to 0 and 1 respectively.
The initial condition is standard: we suppose that the distribution of the concentration u(x, t) is known for all x at the time beginning t ϭ 0. Considering the boundary conditions, the situation turns out to be much more complicated. From the mathematical viewpoint there are two possible types of boundary conditions for this FDE. First, we may consider a generalization of the Dirichlet condition, i.e. to prescribe a value for (2 Ϫ ␣)-integral, or a generalization of the Neumann condition which deals with a value of (␣ Ϫ 1)-derivative.
However, we wish to model some real situation so that the question of physical interpretation of the boundary conditions arises. The solution of this issue is not widely accepted so far. In [8] the authors suggest the idea of inseparable twins which essentially means that we have to find another physical quantity which a --ĥ h corresponds to the needed fractional derivative and then we interpret the boundary condition via this new quantity. In our case, the inseparable twins are settled by the relation (4), which gives a clear physical meaning to the Neumann condition (flux). The sense of the Dirichlet condition still remains doubtful.
Considering these ideas, we arrive at the initial value problem with Neumann boundary conditions in the form , t ʦ R + , i.e. we prescribe the flux on both ends of the interval [0, 1].

The Discretization of the Fractional Diffusion Equation
Before we start to deal with the problem (5), we should clarify our discrete fractional tools and perhaps mention other approaches. Some authors use, probably for the sake of simplicity, local approximation of the differintegral (see [2]). However, this approach is slightly rough because it suppresses the whole idea of fractional derivative as a global operator.
As a much better way it turns out to use fractional differences. This approach starts from the Grunwald-Letnikov definition of a differintegral which coincides under some assumptions with the Riemann-Liouville one (see [14]). As one can observe in many papers (e.g. [6], [9], [15] and [16]), the fractional differences are the most frequently used approximations of differintegrals of the Riemann-Liouville and the Caputo types (more information about the Caputo differintegrals in [14]). On the other hand, the fractional differences do not naturally work with initial and boundary conditions which are needed in the Riemann-Liouville approach. Therefore, all papers mentioned above deal either with the Caputo differintegral or with the Riemann-Liouville differintegral and zero initial conditions. If nonzero initial conditions occur in the Riemann-Liouville case, Podlubny introduces a change of variable which transforms the initial conditions into the zero ones (see [15]).
Being inspired by the fractional calculus on time-scales (Definitions 3.3 and 3.2), we consider a fractional difference to be an ordinary difference of an appropriate fractional sum. That is why the numerical scheme presented in this paper removes the problems with nonzero Riemann-Liouville initial or boundary conditions as we illustrate with problem (5).
We introduce an equidistant discretization of the space variable x n ϭ nh, where n ϭ 0, 1, …, N, h ϭ 1/N and of the time variable t i ϭ i, where i ϭ 0, 1, …, M, ϭ T/M. N is a number of space  Fig. 1). The discretized problem (5) then looks like where the symbol G 2 n denotes the classical second central difference . The Riesz fractional difference has not been defined previously, because it is not clear which type of a difference should be applied to the corresponding Riesz integral. Considering our problem, we suggest the above introduction of G 2 n with respect to its symmetry. This choice implies that the symbol 0 G ␣Ϫ1 N,n used in the boundary conditions is not defined uniquely because we do not determine whether the nabla or delta difference is supposed to be utilized.
This uncertainty allows to employ the delta difference for n ϭ 1 and the nabla one for n ϭ N Ϫ 1 which is necessary for involving the boundary conditions to the system of equations. Rewrite all integer-order differences and study the equations originated from (6).
At the inner points (black circles in Fig. 1), i.e for 1 Յ i Յ M and 2 Յ n Յ N Ϫ 2 we arrive to the expected relation (7) where for i ϭ 1 the initial condition has to be used. At the points adjacent to the boundary points (white circles in Fig. 1) the situ- Similarly at the right boundary (n ϭ N Ϫ 1) we utilize the second boundary condition b 1 i (9) For an easier implementation it is necessary to rewrite the equations (7)-(9) into a compact matrix form.
Obviously, it is possible to write the Riesz fractional sum of a function f (without the first and the last point) in the form of matrix product where the matrix of the generalized binomial coefficients with the factor h ␣ /2 will be denoted by 0 G Ϫ␣ N and its dimension will always correspond with the dimension of the vector f.
After an analysis of the relations (7)-(9) we arrive at the matrix of coefficients . Now we denote E the unit matrix of the dimension N Ϫ 1 and define the vectors of all the involved quantities as u 0 ϭ (g 1 … g NϪ1 ) T (the initial concentration), Therefore we can calculate the solutions in the space region recursively.

Examples
We present some examples of the problem (5), where we consider the parameters ␣ ϭ 2, 1.9, 1.7, 1.4. The initial state is depicted in all the figures with the red line. The space step is h ϭ 0.01 (i.e. N ϭ 100), the time step ϭ 0.01 and we follow the system until the time is T ϭ 0.2 (i.e. M ϭ 20).
First, we consider the insulated system with no source term, hence the total mass of the system is conserved. Translated to the mathematical language it means , , The initial distribution of particles (i.e. the initial condition) is The solutions obtained via our discretization scheme are drawn in Fig. 2. The conservation of mass is evident but there is another interesting effect. With the decreasing value of ␣, the particles are accumulated close to the boundaries. Roughly speaking, the lower value of ␣ means greater probability of very long particle jumps, hence due to the boundedness of the region more and more particles are forced to stop on the boundaries although their jump was supposed to be longer.
The second example describes the constant inflow to the system from the right side. We assume that the region is empty at the beginning and there is no source term: The numerical solutions of this problem are depicted in Fig.  3. The amount of the influent particles is the same for all ␣, hence due to the boundary effect described in the previous example, there are fewer particles in the middle of the region for the decreasing ␣.
Last we present the system with insulated boundaries, but with nonzero source term. The source works for the first five timesteps (blue lines in Fig. 4) and then it is cut off: The solutions drawn in Fig. 4 confirm the observation from the first example, i.e. the lower values of ␣ emphasize sudden changes of the particle concentration.

Summary
We develop a numerical scheme which provides numerical solutions for both initial and boundary value problems involving the Riemann-Liouville fractional differential equations. In particular, we can take also the case when these conditions are nonzero.
The key idea is to consider the fractional difference strictly as a composition of the integer-order difference and the fractional , , ,  sum. The choice of the type of this integer-order difference (delta versus nabla) depends on the particular equation. We illustrate this choice by the discretization of the fractional diffusion equation.