Numerical solution of fractional order diffusion problems with Neumann boundary conditions
Abstract
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 wellposedness 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; nonlocal derivative; Neumann boundary conditions;
implicit Euler scheme.
Department of Applied Analysis and Computational Mathematics,
Eötvös Loránd University H1117, Budapest, Pázmány P. s. 1/C, Hungary.
1 Introduction
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 [24] through population dynamics [7] and groundwater flows [2] to the anomalous diffusion of some chemical compounds. These observations inspired the application of the fractional calculus, which has a long history [17]. An alternative approach for modeling superdiffusion can be given in the framework of the nonlocal calculus, which has been recently developed, see [6] and [5] for an uptodate overview with further references.
Many attempts were made for the numerical solution of the corresponding spacefractional PDE’s. Most of them based on finite difference discretization. It was pointed out that a nontrivial discretization of the one sided fractional derivatives lead to a stable method [14]. This result has been generalized in many aspects: higherorder methods and multidimensional schemes [4],[21],[22],[23],[27] were constructed and analyzed. Note that recently the finite element (Galerkin) discretization was initiated [16] for the fractional diffusion equations and a composite approach was analyzed [11].
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 nonlocal calculus [5].
Regarding Neumann type boundary conditions, at our best knowledge, there is only one attempt [12], 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 wellposedness of the corresponding PDE’s

construct a corresponding finite difference scheme with a full error analysis.
2 Mathematical preliminaries
Fractional calculus
We summarize some basic notions and properties of the fractional
calculus. For more details and examples we refer to the monographs
[13],[15],[18].
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
Remarks:
1. Functions in are bounded and
they are continuous except of the possible discontinuity points
.
2. In the points we define
to be .
Definition 1
For the exponent the fractional order integral operators and on the space are defined with
and
With this, the left and rightsided Riemann–Liouville derivatives of order are given by
and
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 .
Remarks:
1. For the simplicity, the constant – which corresponds
to the intensity of the superdiffusive process in the reallife
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 [13].
Moreover, the definition can be extended to be a bounded operator
on the function space , see Theorem 2.6 in [18].
For an overview of the alternating notations and definitions, we
refer to the review paper [10].
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 [26]. This can be introduced via
Fourier transform, see Section 2.6 in [18]. For more
information and multidimensional extension of the Riesz derivative
see also Section 2.10 in [13]. Note that finite difference
discretizations for Riesz fractional derivatives has been studied also
in [19] 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 [1], [18], which is
beyond the scope of this paper. Some sufficient conditions on a bounded
interval are also discussed in Lemma 2.2 in [13].
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 .
Lemma 1
For each function and any exponent the fractional order integral operators and make sense.
Proof: We prove the statement for the rightsided approximation, the left sided can be handled in a similar way. In concrete terms, we prove that
(1) 
Obviously, there is a such that and accordingly,
(2) 
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.
Fundamental solutions
For the forthcoming analysis we analyze the Cauchy problem
(3) 
with a given initial function and .
Lemma 2
The Cauchy problem in (3) has a unique solution and can be given by
(4) 
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 [10], 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., [9], p. 1111) we have
(5)  
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.
Discretization
The finite difference approximation of the fractional order derivatives is not straightforward. It turns out that an obvious onesided finite difference approximation of the onesided Riemann–Liouville derivatives results in an unstable method even if an implicit Euler method is applied for the time marching scheme [14]. To stabilize these schemes, we need to use the translated Grünwald–Letnikov formula, which is given for with
(6)  
and
(7)  
depending on the translation parameter , the order of the differentiation and the discretization parameter , where we used the coefficients
The principle of the twosided translated discretizations is depicted in Figure 1.
These coefficients satisfy the following:
(8)  
We use the same notation for the discrete differential (or difference) operators, i.e. for each we write
(9) 
where the superscript denotes the th component.
Remarks:
1. One can prove [14] that the integrals in
(6) and (7) approximate
the corresponding Riemann–Liouville derivatives in the following sense:
(10) 
and similarly,
(11) 
provided that both of the Fourier transform of and that of
and
are in . We will point out that in the present framework
no such assumption is necessary.
2. Higherorder finite difference approximations can be obtained
as a linear combination of first order ones using different translation
parameters.
For example, the sum defined by
(12) 
provides a second order accurate approximation of the Riemann–Liouville
derivative.
Similar statement holds if is switched to .
For the details, see [27].
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 [25].
4. An alternative approximation of fractional order elliptic operators
(which can be applied in multidimensional cases) was proposed in
[3], which can be a basis also for finite element
discretizations.
3 Results
3.1 Extensions
Extensions are not only necessary to have wellposed 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.
Definition 2
We say that the extension
is compatible with the homogeneous Neumann (noflux) 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 [20] which is generalized to the case of two absorbing walls. For the simplicity, the definition is given for functions .
Definition 3
We call the 2periodic extension of the function
the Dirichlet type extension of and we denote it with .
Remarks:
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 zeroconcentration at the end points in the consecutive
time steps, we have to force a skewsymmetric 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
(13) 
We define similarly the extension of the vector with :
(14) 
Extension corresponding to homogeneous Neumann boundary condition
Similarly to the previous case, the definition is given for functions .
Definition 4
We call the 2periodic extension of the function
the Neumann type extension of and we denote it with .
Remarks:
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 zeroflux, we force a symmetric concentration profile
around and .
A simple calculation shows that
(15) 
Similar notations are used for the “extended” vector of which is defined as follows:
(16) 
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.
Lemma 3
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
such that
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
with
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
(17) 
The following theorem is the basis of our numerical method, which again confirms the favor of the Neumann type extension.
Theorem 1
For any the problem in (17) is wellposed, and its unique solution satisfies the boundary conditions .
Proof: We first note that the problem
(18) 
is wellposed and corresponding to Lemma 3 its solution can be given by
This implies that
and
Therefore, the restriction of the solution of (18) solves the problem in (17). Here we have also used throughout that for the Neumann extension: such that the Riesz derivative makes sense.
To prove the uniqueness, we assume that solves (17). According to (15) we obviously have that the Neumann extension satisfies
(19) 
To compute the fractional order integral we introduce the function such that we have
Therefore, we also have
(20) 
Consequently, (19) and (20) imply
and the periodicity obviously gives
such that the Neumann extension solves (18). This, however, has a unique solution, which gives the uniqueness of the solution of (17).
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:
(21) 
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
(22) 
This is combined with an implicit Euler time stepping to obtain
(23) 
To make the consecutive formulas more accessible, we
expand (22) in a concrete example.
Example We give the first component of
for .
In the general case, using (22), (9) and (16) we obtain that
(24)  
(25)  
(26)  