Numerical solution of fractional order diffusion problems with Neumann boundary conditions
A finite difference numerical method is investigated for fractional
order diffusion problems in one space dimension.
For this, a mathematical model is developed to incorporate homogeneous
Dirichlet and Neumann type boundary conditions.
The models are based on an appropriate extension
of the initial values. The well-posedness of the obtained initial
value problems is proved and it is pointed out that the extensions
are compatible with the above boundary conditions.
Accordingly, a finite difference scheme is constructed for the Neumann
problem using the shifted Grünwald–Letnikov approximation of the
fractional order derivatives, which is based on infinite many basis points.
The corresponding matrix is expressed in a closed form and the
convergence of an appropriate implicit Euler scheme is proved.
Keywords: fractional order diffusion; Grünwald–Letnikov formula; non-local derivative; Neumann boundary conditions; implicit Euler scheme.
Department of Applied Analysis and Computational Mathematics,
Eötvös Loránd University H-1117, Budapest, Pázmány P. s. 1/C, Hungary.
A widely accepted constitutive relation, the first Fick’s law leads to the standard diffusion model. At the same time, more observations confirmed the presence of super- and subdiffusive dynamics in several phenomena. They range from the plasma physics  through population dynamics  and groundwater flows  to the anomalous diffusion of some chemical compounds. These observations inspired the application of the fractional calculus, which has a long history . An alternative approach for modeling superdiffusion can be given in the framework of the nonlocal calculus, which has been recently developed, see  and  for an up-to-date overview with further references.
Many attempts were made for the numerical solution of the corresponding space-fractional PDE’s. Most of them based on finite difference discretization. It was pointed out that a non-trivial discretization of the one sided fractional derivatives lead to a stable method . This result has been generalized in many aspects: higher-order methods and multi-dimensional schemes ,,,, were constructed and analyzed. Note that recently the finite element (Galerkin) discretization was initiated  for the fractional diffusion equations and a composite approach was analyzed .
The problem which inspired the present research is the following. If a superdiffusive evolution of some density is observed on a physical volume then we have information on the density only on the closure of this. On the other hand, in the mathematical model, nonlocal operators are used, which require the value of also outside of the domain. For an accurate finite difference approximation we also need values outside of the domain. In this way, it is natural to look for an extension of the density .
The majority of the authors consider homogeneous Dirichlet boundary conditions and assume zero values outside of the domain. A similar approach is applied in the non-local calculus .
Regarding Neumann type boundary conditions, at our best knowledge, there is only one attempt , dealing the boundary condition at the operator level and proposing a matrix transformation technique.
The aim of this paper is to
develop a meaningful mathematical approach to model homogeneous Neumann (and Dirichlet) boundary conditions
analyze the well-posedness of the corresponding PDE’s
construct a corresponding finite difference scheme with a full error analysis.
2 Mathematical preliminaries
We summarize some basic notions and properties of the fractional
calculus. For more details and examples we refer to the monographs
To define an appropriate function space for the fractional order derivatives on the real axis we introduce for arbitrary the function spaces
With these we define
1. Functions in are bounded and they are continuous except of the possible discontinuity points .
2. In the points we define to be .
For the exponent the fractional order integral operators and on the space are defined with
With this, the left and right-sided Riemann–Liouville derivatives of order are given by
where is the integer with .
Accordingly, for a function , where and is a constant function we define the Riesz derivative with
where with a given positive constant .
1. For the simplicity, the constant – which corresponds to the intensity of the superdiffusive process in the real-life phenomena – does not appear in the notation .
2. In the original definition, the Riemann–Liouville derivatives are given on a bounded interval such that in the above definition and are substituted with and , respectively. In this case and can be defined for the exponent or even for with , in case of complex valued functions, see Section 2.1 in . Moreover, the definition can be extended to be a bounded operator on the function space , see Theorem 2.6 in . For an overview of the alternating notations and definitions, we refer to the review paper .
3. An advantage of the above approach is that alternative definitions on real axis coincide. For instance, fractional derivatives can be interpreted using the fractional power of the negative Laplacian , see Lemma 1 in . This can be introduced via Fourier transform, see Section 2.6 in . For more information and multidimensional extension of the Riesz derivative see also Section 2.10 in . Note that finite difference discretizations for Riesz fractional derivatives has been studied also in  highlighting its connection with probabilistic models.
4. We have left open the question for which functions does Definition 1 make sense. The general answer requires the discussion of the Triebel–Lizorkin spaces , , which is beyond the scope of this paper. Some sufficient conditions on a bounded interval are also discussed in Lemma 2.2 in . At the same time, since we will approximate the Riesz derivative with finite differences of the fractional integrals and verify that the fractional integrals make sense on the function space .
For each function and any exponent the fractional order integral operators and make sense.
Proof: We prove the statement for the right-sided approximation, the left sided can be handled in a similar way. In concrete terms, we prove that
Obviously, there is a such that and accordingly,
Here, using the condition we have that the first term is finite. To estimate the second one, we introduce the function with
such that on . Also, since
, we have that such that is bounded.
For the forthcoming analysis we analyze the Cauchy problem
with a given initial function and .
The Cauchy problem in (3) has a unique solution and can be given by
where denotes the fundamental solution corresponding to the Riesz fractional differential operator , furthermore, for all .
Proof: Applying the spatial Fourier transform to the equations in (3) we have
where we have used the identity , see , p. 38. Therefore, such that an inverse Fourier transform implies that
In this way, the fundamental solution of (3) can be given as
Using the fact that the function to transform is even and the Fourier transform can be given (see, e.g., , p. 1111) we have
Here for any fixed the real function given by is in and also all of its derivatives are in . On the other hand, for the real function given with is in , therefore is bounded and continuous. Consequently, using (5) the right hand side of the equality
makes sense, which gives statement in the lemma.
The finite difference approximation of the fractional order derivatives is not straightforward. It turns out that an obvious one-sided finite difference approximation of the one-sided Riemann–Liouville derivatives results in an unstable method even if an implicit Euler method is applied for the time marching scheme . To stabilize these schemes, we need to use the translated Grünwald–Letnikov formula, which is given for with
depending on the translation parameter , the order of the differentiation and the discretization parameter , where we used the coefficients
The principle of the two-sided translated discretizations is depicted in Figure 1.
These coefficients satisfy the following:
We use the same notation for the discrete differential (or difference) operators, i.e. for each we write
provided that both of the Fourier transform of and that of
are in . We will point out that in the present framework
no such assumption is necessary.
2. Higher-order finite difference approximations can be obtained as a linear combination of first order ones using different translation parameters. For example, the sum defined by
provides a second order accurate approximation of the Riemann–Liouville
Similar statement holds if is switched to .
For the details, see .
3. The nonlocal effect of the differential operators result in full matrices. At the same time, one can save some computing efforts with an appropriate decomposition of the above matrix .
4. An alternative approximation of fractional order elliptic operators (which can be applied in multidimensional cases) was proposed in , which can be a basis also for finite element discretizations.
Extensions are not only necessary to have well-posed problems involving nonlocal diffusion operators, but also essential at the discrete level. In order to have sufficient accuracy in the finite difference approximation near to the boundary and at the boundary of a nonlocal differential operator, it is necessary to have (virtual) gridpoints outside of the original computational domain . This is clearly shown in (6), (7) and (12).
Summarized, the sketch of our approach is the following:
we extend the problem to to get rid of the boundary conditions
we solve the corresponding Cauchy problem (3) (we will approximate this with finite differences)
we verify that the desired homogeneous Neumann (or homogeneous Dirichlet) boundary conditions are satisfied for the restriction of the solution.
We say that the extension
is compatible with the homogeneous Neumann (no-flux) or Dirichlet boundary conditions and the operator if the function is the unique solution of the following problem
and or for and .
In rough terms, one can say that a correct extension of the solution from is the function which solves the extended problem on the whole real axis such that the boundary condition on is still satisfied.
Extension corresponding to homogeneous Dirichlet boundary condition
As a motivation, we use the idea in  which is generalized to the case of two absorbing walls. For the simplicity, the definition is given for functions .
We call the 2-periodic extension of the function
the Dirichlet type extension of and we denote it with .
1. This is sometimes called the odd extension of and can be obtained first with reflecting the graph of to and to extend it to and reflecting this graph further to and to extend it to and continuing this process.
2. A natural physical interpretation of this extension is the following: To ensure zero-concentration at the end points in the consecutive time steps, we have to force a skew-symmetric concentration profile around and .
3. We may define for but as we will see this is not essential.
A simple calculation shows that we have
We define similarly the extension of the vector with :
Extension corresponding to homogeneous Neumann boundary condition
Similarly to the previous case, the definition is given for functions .
We call the 2-periodic extension of the function
the Neumann type extension of and we denote it with .
1. This is sometimes called the even extension of and can be obtained first with reflecting the graph of to the vertical line given with to extend it to then to the vertical line given with to extend it to and continuing this process.
2. A natural physical interpretation of this extension is the following: To ensure zero-flux, we force a symmetric concentration profile around and .
A simple calculation shows that
Similar notations are used for the “extended” vector of which is defined as follows:
A natural physical interpretation of this extension is that particles are reflected at the boundaries to ensure zero flux. In this way, we also reflect the concentration profile in the model. The principle of the Neumann extension for a vector is visualized in Figure 2.
We verify that the above extensions meet the requirements in Definition 2.
The extensions and are compatible with the Dirichlet and the Neumann boundary conditions, respectively and with the differential operator .
Proof: According to Lemma (2) we can express the solution of
with the convolution
Since , the same holds for . Accordingly, the right and left limit of in coincide. Using this, (15) and the fact that is even we obtain
which gives that the homogeneous Neumann boundary condition is satisfied in . With an obvious modification, using (15), we can also verify the homogeneous Neumann boundary condition .
Similarly, we can express the solution of
Since , the same holds for . Accordingly, the right and left limit of in coincide. Using this and (13) we obtain
which gives that the homogeneous Dirichlet boundary condition is satisfied in . With an obvious modification, using (13), we can verify homogeneous Dirichlet boundary condition also in .
In both cases, the equality is satisfied in which gives the statement in the lemma.
Using the above extension, we will investigate the numerical solution of the problem
The following theorem is the basis of our numerical method, which again confirms the favor of the Neumann type extension.
For any the problem in (17) is well-posed, and its unique solution satisfies the boundary conditions .
Proof: We first note that the problem
is well-posed and corresponding to Lemma 3 its solution can be given by
This implies that
To compute the fractional order integral we introduce the function such that we have
Therefore, we also have
and the periodicity obviously gives
3.2 Numerical methods
Following the classical method of lines technique we first discretize the spatial variables in the extended problem and choose a time stepping scheme for the full discretization.
To streamline the forthcoming computations, the interval will be transformed to , where we use the following grid points:
denotes the numerical approximation at time in the grid point and the values of the analytic solution at time in the grid points.
3.2.1 Analysis of a finite difference scheme
For the spatial discretization we use the Grünwald–Letnikov approximations in (9) introducing with
This is combined with an implicit Euler time stepping to obtain
To make the consecutive formulas more accessible, we
expand (22) in a concrete example.
Example We give the first component of for .