Numerical solution of fractional order diffusion problems with Neumann boundary conditions

Numerical solution of fractional order diffusion problems with Neumann boundary conditions

Béla J. Szekeres and Ferenc Izsák 111Corresponding author. E-mail address:, tel.: +36 13722500-8428, fax: +36 13812158.

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.

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 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 [14]. This result has been generalized in many aspects: higher-order methods and multi-dimensional 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 non-local 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 well-posedness 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

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


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 [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 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.

With this the second term on the right hand side of (2) can be rewritten as

which is also finite since is bounded and . Therefore, (1) is finite, which proves the lemma.  

Fundamental solutions

For the forthcoming analysis we analyze the Cauchy problem


with a given initial function and .

Lemma 2

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 [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


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 [14]. 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.


Figure 1: Basis points for the the left-sided () and the right-sided () translated Grünwald–Letnikov formula (given in (6) and (7)) applied in with the translation parameter .

These coefficients satisfy the following:


We use the same notation for the discrete differential (or difference) operators, i.e. for each we write


where the superscript denotes the th component.
1. One can prove [14] that the integrals in (6) and (7) approximate the corresponding Riemann–Liouville derivatives in the following sense:


and similarly,


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. 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 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 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.

Definition 2

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 [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 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 .

Definition 4

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.

Figure 2: Entries of the Neumann extension of the vector .

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


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.

Theorem 1

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


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


To compute the fractional order integral we introduce the function such that we have

Therefore, we also have


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:


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 .

In the general case, using (22), (9) and (16) we obtain that