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: izsakf@cs.elte.hu, tel.: +36 13722500-8428, fax: +36 13812158.
###### 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 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  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

##### Fractional calculus

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

 ¯C(a,b)=C[a,b]|(a,b)and¯C(a,b)/R={f∈¯C(a,b):∫baf=0}.

With these we define

 CI(R)={b−a - periodic extension of f: f∈¯C(a,b)/R}.

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

 −∞Iβxf(x)=1Γ(β)∫x−∞f(s)(x−s)1−βds

and

 xIβ∞f(x)=1Γ(β)∫∞xf(s)(s−x)1−βds.

With this, the left and right-sided Riemann–Liouville derivatives of order are given by

 RL−∞∂αxf(x)=∂nx−∞In−αxf(x)

and

 RLx∂α∞f(x)=(−1)n∂nxxIn−α∞f(x),

where is the integer with .
Accordingly, for a function , where and is a constant function we define the Riesz derivative with

 ∂α|x|f(x)=Cσ(RL−∞∂αxf0(x)+RLx∂α∞f0(x)),

where with a given positive constant .

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

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

 ∫∞xf(s)(s−x)α−1ds<∞. (1)

Obviously, there is a such that and accordingly,

 ∫∞xf(s)(s−x)α−1ds=∫a+(b−a)kxf(s)(s−x)α−1ds+∫∞a+(b−a)kf(s)(s−x)α−1ds. (2)

Here, using the condition we have that the first term is finite. To estimate the second one, we introduce the function with

 F(s)=∫sa+(b−a)kf(s∗)ds∗

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

 ∫∞a+(b−a)kF′(s)(s−x)α−1ds =[F(s)(s−x)α−1]∞a+(b−a)k−∫∞a+(b−a)kF(s)⋅(1−α)(s−x)αds =−∫∞a+(b−a)kF(s)⋅(1−α)(s−x)αds,

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

 (3)

with a given initial function and .

###### Lemma 2

The Cauchy problem in (3) has a unique solution and can be given by

 u(t,x)=(Φα,t∗u0)(x), (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

 {∂tFu(t,s)=−|s|αFu(t,s)t∈R+,s∈RFu(0,s)=Fu0(s)s∈R,

where we have used the identity , see , p. 38. Therefore, such that an inverse Fourier transform implies that

 u(t,x)=F−1(e−t|s|αFu0(s))(x)=F−1(e−t|s|α)∗u0(x).

In this way, the fundamental solution of (3) can be given as

 Φα,t(s)=F−1(e−t|x|α)(s)=F(e−t|x|α)(s).

Using the fact that the function to transform is even and the Fourier transform can be given (see, e.g., , p. 1111) we have

 Φα,t(s) =F(e−t|x|α)(s)=F(e−t|x|e−t|x|α−1)(s)=F(e−t|x|)(s)∗F(e−t|x|α−1)(s) (5) =√2πtt2+s2∗F(e−t|x|α−1)(s).

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

 ∂kΦα,t(s)=∂k√2πtt2+s2∗F(e−t|x|α−1)(s)

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

 Dα,p,h−∞,GLf(x) =1Γ(−α)limM→∞1hαM∑k=0Γ(k−α)Γ(k+1)f(x−(k−p)h) (6) =1hα∞∑k=0gkf(x−(k−p)h)

and

 Dα,p,h∞,GLf(x) =1Γ(−α)limM→∞1hαM∑k=0Γ(k−α)Γ(k+1)f(x+(k−p)h) (7) =1hα∞∑k=0gkf(x+(k−p)h)

depending on the translation parameter , the order of the differentiation and the discretization parameter , where we used the coefficients

 gk=Γ(k−α)Γ(−α)Γ(k+1)=(−1)k(αk).

The principle of the two-sided translated discretizations is depicted in Figure 1.

These coefficients satisfy the following:

 ∞∑k=0gk=0∀α∈(1,2] (8) g1=−α,gj≥0forj≠1.

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

 [Dα,p,h−∞,GLv]j=1hα∞∑k=0(−1)kgkvj+p−kand[Dα,p,h∞,GLv]j=1hα∞∑k=0(−1)kgkvj+k−p, (9)

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

 RL−∞∂αxf(x)=Dα,p,h−∞,GLf(x)+O(h) (10)

and similarly,

 RLx∂α∞f(x)=Dα,p,h∞,GLf(x)+O(h), (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. 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

 GL−∞Dα,p,qx,hu(x)=2q−α2(p−q)Dα,p,h−∞,GLu(x)+2p−α2(p−q)Dα,p,h−∞,GLu(x) (12)

provides a second order accurate approximation of the Riemann–Liouville derivative. 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.

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

 ~⋅:L2(Ω)→CI(R)

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

 {∂t~u(t,x)=∂α|x|~u(t,x)t∈R+,x∈R~u(0,x)=~u0(x)x∈R

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 .

###### Definition 3

We call the 2-periodic extension of the function

 ^uD,1(x)={u(x)forx∈(0,1)−u(−x)forx∈(−1,0)

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

 ^uD(y)=−^uD(−y)and^uD(1+y)=−^uD(1−y)y∈R. (13)

We define similarly the extension of the vector with :

 [^vD]j=⎧⎪⎨⎪⎩vjj=0,1,…,N+1−v2(N+1)−jj=N+2,N+3,…,2N+1vj+2Nj∉{0,1,…,2N+1}. (14)
##### 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

 ^uN,1(x)={u(x)forx∈(0,1)u(−x)forx∈(−1,0)

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 zero-flux, we force a symmetric concentration profile around and .

A simple calculation shows that

 ∂x^fN(1+y)=∂x^fN(1−y)and∂x^fN(y)=∂x^fN(−y)y∈R. (15)

Similar notations are used for the “extended” vector of which is defined as follows:

 [^vN]j=⎧⎪⎨⎪⎩vjj=0,1,…,Nv−j−1j=−N−1,−N,…,−1vj+2N+2j∉{−N−1,…,0,…,N}. (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

 {∂tu=∂α|x|uonR+×Ru(0,⋅)=^uN

with the convolution

 u(t,⋅)=Φα,t∗^uN=^uN∗Φα,t,

such that

 u(t,x)=∫R^uN(x−y)Φα,t(y)dy.

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

 ∂xu(t,1) =limϵn→0−∂xu(t,1−ϵn)=limϵn→0−(∂x^uN∗Φα,t)(1−ϵn) =limϵn→0−∫R∂x^uN(1−ϵn−y)Φα,t(y)dy=−limϵn→0−∫R∂x^uN(1+ϵn+y)Φα,t(y)dy =−limϵn→0−∫R∂x^uN(1+ϵn+y)Φα,t(−y)dy =−limϵn→0−∫R∂x^uN(1+ϵn−y)Φα,t(y)dy=−limϵn→0−(∂x^uN∗Φα,t)(1+ϵn) =−limϵn→0−∂xu(t,1+ϵn)=−∂xu(t,1),

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

 {∂tu=∂α|x|uonR+×Ru(0,⋅)=^uD

with

 u(t,x)=∫R^uD(x−y)Φα,t(y)dy.

Since , the same holds for . Accordingly, the right and left limit of in coincide. Using this and (13) we obtain

 u(t,1) =limϵn→0−u(t,1−ϵn)=limϵn→0−^uD∗Φα,t(1−ϵn) =limϵn→0−∫R^uD(1−ϵn−y)Φα,t(y)dy=−limϵn→0−∫R^uD(1+ϵn+y)Φα,t(y)dy =−limϵn→0−∫R^uD(1+ϵn+y)Φα,t(−y)dy=−limϵn→0−∫R^uD(1+ϵn−y)Φα,t(y)dy =−limϵn→0−^uD∗Φα,t(1+ϵn)=−limϵn→0−u(t,1+ϵn)=−u(t,1),

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

 {∂tu(t,x)=∂α|x|uN(t,⋅)(x)t∈(0,T),x∈(0,1)u(0,x)=u0x∈(0,1). (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 well-posed, and its unique solution satisfies the boundary conditions .

Proof: We first note that the problem

 {∂tu(t,x)=∂α|x|u(t,x)t∈(0,T),x∈Ru(0,x)=^uN0(x)x∈R. (18)

is well-posed and corresponding to Lemma 3 its solution can be given by

 u(t,x)=∫R^uN(x−y)Φα,t(y)dy.

This implies that

 u(t,−x) =∫R^uN(−x−y)Φα,t(y)dy=∫R^uN(x+y)Φα,t(−y)dy =∫R^uN(x−y)Φα,t(y)dy=u(t,x).

and

 u(t,x+2) =∫R^uN(x+2−y)Φα,t(y)dy=∫R^uN(x−y)Φα,t(y)dy=u(t,x).

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

 ∂t^vN(t,−x)=∂t^vN(t,x)x∈R. (19)

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

 Γ(2−α)∞I2−αx^vN0(−x)=∫−x−∞vN0(y∗)(−x−y∗)αdy∗+∫∞−xvN0(y∗)(y∗+x)αdy∗ =∫∞xvN0(−y)(−x+y)αdy+∫x−∞vN0(−y)(−y+x)αdy=∫∞xvN0(y)(−x+y)αdy+∫x−∞vN0(y)(−y+x)αdy =Γ(2−α)∞I2−αx^vN0(x).

Therefore, we also have

 ∂α|x|vN(t,−x)=∂α|x|vN(t,x). (20)

Consequently, (19) and (20) imply

 ∂t^vN(t,−x)=∂t^vN(t,x)=∂α|x|vN(t,x)=∂α|x|vN(t,−x)

and the periodicity obviously gives

 ∂t^vN(t,2+x)=∂t^vN(t,x)=∂α|x|vN(t,x)=∂α|x|vN(t,2+x)

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:

 xj:=πh(j+12)=π(12(N+1)+jN+1). (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

 Aα,hu=Dα,1,h−∞,GLuN+Dα,1,h∞,GLuN. (22)

This is combined with an implicit Euler time stepping to obtain

 un+1j−unjτ=[Aα,hun+1]j. (23)

To make the consecutive formulas more accessible, we expand (22) in a concrete example.
Example We give the first component of for .

 [Aα,hv]1 =g0v0+g1v1+g2v2+g3v3+g4v3+g5v2+g6v1+g7v0+g8v0+… +g0v2+g1v1+g2v0+g3v0+g4v1+g5v2+g6v3+g7v3+g8v2+….

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

 [Aα,hv]j=Cσhα(j+1∑k=0gkvj−k+1+N−j+1∑k=0gkvk+j−1 (24) +∞∑l=0⎛⎝j+N+2∑k=j+2gk+2l(N+1)vk−j−2+2N−j+2∑k=N−j+2gk+2l(N+1)v2N−k−j+2⎞⎠ +∞∑l=0⎛⎝j+2N+3∑k=j+N+3gk+2l(N+1)v2N+j−k+3+3N−j+3∑k=2N−j+3gk+2l(N+1)vk−2N+j−3⎞⎠⎞⎠ j=1,2,…,N−1
 [Av]0=2Cσhα(g0v1+g1v0 (25) +∞∑l=0(N+2∑k=2gk+2l(N+1)vN+2−k+2N+3∑k=N+3gk+2l(N+1)v2N−k+3)),
 [Av]N=2Cσhα(g0vN−1+g1vN (26) +∞∑l=0(N+2∑k=2gk+2l(N+1)vN−k+2+2N+3∑k=N+3gk+2l(N+1