A finite difference method for fractional diffusion equations with Neumann boundary conditions

Abstract A finite difference numerical method is investigated for fractional order diffusion problems in one space dimension. The basis of the mathematical model and the numerical approximation is an appropriate extension of the initial values, which incorporates homogeneous Dirichlet or Neumann type boundary conditions. The wellposedness of the obtained initial value problem is proved and it is pointed out that each extension is compatible with the original 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.

can be combined with tau methods [14,[17][18][19] or with Gauss-Lobatto collocation methods [20]. A furher alternative to solve fractional diffusion problems numerically is the application of fractional order Laplace transform [21], [22]. The last two approaches can also deliver higher-order spatial accuracy.
The problem which inspired the present research is the following. If a superdiffusive evolution of some density u 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 u 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 u.
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 [6].
At the same time, the homogeneous Neumann boundary conditions have a central importance in the modeling of real-life phenomena, since this corresponds to zero flux at the boundary. To our best knowledge, there is only one attempt [23], dealing with the Neumann type boundary condition at the operator level and proposing a matrix transformation technique.
The main objective of this paper is to develop a finite difference numerical solution for space-fractional diffusion problems with Neumann boundary conditions. In details, the main steps in the articles are 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.

Mathematical preliminaries 2.1 Fractional calculus
We summarize some basic notions and properties of the fractional calculus. For more details and examples we refer to the monographs [24][25][26] and the recent work [27].
To define an appropriate function space for the fractional order derivatives on the real axis we introduce for arbitrary a; b 2 R the function spaces where n is the integer with n 1 <˛< n. Accordingly, for a function f D f 0 C C f , where f 0 2 C I .R/ and C f is a constant function we define the Riesz derivative with where C D 2 cos˛ 2 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 @j xj . 2. In the original definition, the Riemann-Liouville derivatives are given on a bounded interval .a; b/ R such that in the above definition 1 and 1 are substituted with a and b, respectively. In this case, a Ix and x Ib can be defined for the exponentˇ2 R C or even forˇ2 C with Reˇ> 0, in case of complex valued functions, see Section 2.1 in [24]. Moreover, the definition can be extended to be a bounded operator on the function space L 1 .a; b/, see Theorem 2.6 in [26]. For an overview of the alternating notations and definitions, we refer to the review paper [28]. 3. An advantage of the above approach is that alternative definitions on the real axis coincide. For instance, fractional derivatives can be interpreted using the fractional power of the negative Laplacian , see Lemma 1 in [29]. This can be introduced via Fourier transform, see Section 2.6 in [26]. For more information and multidimensional extension of the Riesz derivative see also Section 2.10 in [24]. Note that finite difference discretizations for Riesz fractional derivatives has been studied also in [30] highlighting its connection with probabilistic models.

4.
We have left open the question for which functions does Definition 2.1 make sense. The general answer requires the discussion of the Triebel-Lizorkin spaces [31], [26], which is beyond the scope of this paper. Some sufficient conditions on a bounded interval .a; b/ are also discussed in Lemma 2.2 in [24]. At the same time, 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 C I .R/.

Lemma 2.2.
For each function f 2 C I .R/ and any exponent˛2 .1; 2/ the fractional order integral operators 1 I 2 x f and x I 2 1 f 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 k 2 Z such that a C .b a/k > x and accordingly, Here, using the condition 0 <˛ 1 < 1 we have that the first term is finite. To estimate the second one, we introduce the function F W .a C .b a/k; 1/ ! R with Also, since f 2 C I .R/, we have that 0 D F .a C .b a/k/ D F .a C .b a/.k C 1// D : : : such that F is bounded.
With this the second term on the right hand side of (2) can be rewritten as which is also finite since F is bounded and 2 ˇ> 1. Therefore, (1) is finite, which proves the lemma.

Fundamental solutions
For the forthcoming analysis, we analyze the Cauchy problem with a given initial function u 0 2 C I .R/ and˛2 .1; 2.
Lemma 2.3. The Cauchy problem in (3) has a unique solution and can be given by whereˆ˛; t denotes the fundamental solution corresponding to the Riesz fractional differential operator @j xj , furthermore, u.t; / 2 C 1 .R/ for all t 2 R C .
Proof. Applying the spatial Fourier transform F to the equations in (3) we have where we have used the identity F h @j xj u.t; x/ i .s/ D jsj˛Fu.t; s/, see [28], p. 38. Therefore, Fu.t; s/ D e t jsj˛F u 0 .s/ such that an inverse Fourier transform F 1 implies that In this way, the fundamental solution of (3) can be given aŝ˛; t .s/ D F 1 .e t jxj˛/ .s/ D F.e tjxj˛/ .s/: Using the fact that the function to transform is even and the Fourier transform F expf aj jg.s/ can be given (see, e.g., [32], p. 1111) we havê˛; Here for any fixed t the real function given by s ! q 2 t t 2 Cs 2 is in C 1 .R/ and also all of its derivatives are in L 1 .R/. On the other hand, for˛> 1 the real function given with e t jxj˛ 1 is in L 1 .R/, therefore F.e t jxj˛ 1 / is bounded and continuous. Consequently, using (5) the right hand side of the equality @ kˆ˛; t .s/ D @ k q 2 t t 2 C s 2 F.e tjxj˛ 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 [7]. To stabilize these schemes, we need to use the translated Grünwald-Letnikov formula, which is given for and depending on the translation parameter p 2 N, the order of the differentiation˛2 .1; 2 and the discretization parameter h 2 R C , where we used the coefficients The principle of the two-sided translated discretizations is depicted in Figure 1.  (7)) applied in x with the translation parameter p D 2.
These coefficients satisfy the following: We use the same notation for the discrete differential (or difference) operators, i.e. for each v D .: : : ; v 1 ; v 0 ; v 1 ; : : : where the superscript j denotes the j th component.

Remarks.
1. One can prove [7] that the integrals in (6) and (7) approximate the corresponding Riemann-Liouville derivatives in the following sense: RL and similarly, provided that both of the Fourier transform of f and that of RL 1 @x f .x/ and RL x @1f .x/ are in L 1 .R/. 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 1 is switched to 1. For the details, see [12]. 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 [33]. 4. An alternative approximation of fractional order elliptic operators (which can be applied in multidimensional cases) was proposed in [34], which can be a basis also for finite element discretizations.

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 D .a; b/. This is clearly shown in (6), (7) and (12). To summarize, the sketch of our approach is the following: -we extend the problem to R 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 3.1. We say that the extension is compatible with the homogeneous Neumann (no-flux) or Dirichlet boundary conditions and the operator @j xj if the function Q u is the unique solution of the following problem 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 [35] which is generalized to the case of two absorbing walls. For the simplicity, the definition is given for functions u W .0; 1/ ! R.
Definition 3.2. We call the 2-periodic extension of the function the Dirichlet type extension of u and we denote it with O u D .

Remarks.
1. This is sometimes called the odd extension of u and can be obtained first with reflecting the graph of u to .0; 0/ and .1; 0/ 2 R 2 to extend it to . 1; 2/ and reflecting this graph further to . 1; 0/ and .2; 0/ 2 R 2 to extend it to . 2; 3/ 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 0 and 1. 3. We may define u D .k/ D 0 for k 2 Z but as we will see this is not essential.
A simple calculation shows that we have We define similarly the extension

Extension corresponding to homogeneous Neumann boundary condition
Similarly to the previous case, the definition is given for functions u W .0; 1/ ! R.
the Neumann type extension of u and we denote it with O u N .

Remarks.
1. This is sometimes called the even extension of u and can be obtained first with reflecting the graph of u to the vertical line given with x D 0 to extend it to . 1; 0/ then to the vertical line given with x D 1 to extend it to .1; 2/ 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 0 and 1.
A simple calculation shows that Similar notations are used for the "extended" vector O v N 2 R Z of v D .v 0 ; v 1 ; : : : ; v N / 2 R N C1 which is defined as follows: v j j D 0; 1; : : : ; N v j 1 j D N 1; N; : : : ; 1 v j C2N C2 j 6 2 f N 1; : : : ; 0; : : : ; N g : (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 3.1. Proof. According to Lemma 2.3 we can express the solution of x y/ˆ˛; t .y/ dy : Sinceˆ˛; t 2 C 1 .R/, the same holds for u.t; /. Accordingly, the right and left limit of @ x u.t; / in 1 coincide. Using this, (15) and the fact thatˆ˛; t is even we obtain which gives that the homogeneous Neumann boundary condition is satisfied in 1. With an obvious modification, using (15), we can also verify the homogeneous Neumann boundary condition @ x u.t; 0/ D 0.
Similarly, we can express the solution of x y/ˆ˛; t .y/ dy : Sinceˆ˛; t 2 C 1 .R/, the same holds for u.t; /. Accordingly, the right and left limit of u.t; / in 1 coincide. Using this and (13) we obtain which gives that the homogeneous Dirichlet boundary condition is satisfied in 1. With an obvious modification, using (13), we can verify homogeneous Dirichlet boundary condition also in 0.
In both cases, the equality @ t u D @j xj u is satisfied in .0; 1/ 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. Proof. We first note that the problem is well-posed and corresponding to Lemma 3.4 its solution can be given by This implies that Therefore, the restriction of the solution u of (18) solves the problem in (17). Here we have also used throughout that for the Neumann extension: u N 2 C I .R/ such that the Riesz derivative @j xj u N makes sense.
To prove the uniqueness, we assume that v solves (17). According to (15) we obviously have that the Neumann extension satisfies To compute the fractional order integral we introduce the function Consequently, (19) and (20) imply and the periodicity obviously gives such that the Neumann extension u N solves (18). This, however, has a unique solution, which gives the uniqueness of the solution of (17).

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 OE0; 1 will be transformed to OE0; , where we use the following grid points: u n j denotes the numerical approximation at time n in the grid point x j and u.nı; / the values of the analytic solution at time nı in the grid points.

Analysis of a finite difference scheme
For the spatial discretization we use the Grünwald-Letnikov approximations in (9) introducing A˛; h 2 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
In the general case, using (22), (9) and (16) we obtain that j D 1; 2; : : : ; N 1 For K D 2mC1 the corresponding matrix can be given with a slight modification. Observe that all of the coefficients g i arise once on the right hand side of (24), (25) and (26).
Proposition 3.6. The matrix A˛; h has negative diagonal and non-negative off-diagonal elements.
Proof. Observe that in (24), (25) and (26) the coefficient of v j ; v 0 and v N , respectively, is g 1 , and g 1 appears only here. Therefore, using also (8) we obtain that A has negative diagonals and positive off-diagonals.
We analyze the properties of the matrix A˛; h and the corresponding differential operator. The technical proofs of these statements are postponed to the Appendix. Proof. It is sufficient to prove that the right hand side of (23) provides a first order approximation for the Riesz derivative of order j˛j. According to the proof of Theorem 3.5 the analytic solution u N of (18) is periodic, smooth and it satisfies the homogeneous Neumann boundary conditions in 0 and 1. Therefore, its cosine Fourier series is pointwise convergent: where we do not denote the time dependence of the cosine Fourier coefficients F k . Using again that u.t; / 2 C 1 .R/ we also have -using the regularity theory of Fourier series -that for any r 2 N there exists a constant C r such that Using (27) componentwise for t D nı we have that u.hı; / D OEu.nı; x 0 / u.nı; x 1 / : : : u.nı; F k cos kx 1 : : : Using Lemma 3.8 for the expansion in (27) and the matrix A˛; h for (29) according to (42) we obtain the following equality: To prove the proposition we first verify that for a mesh-independent constant C the following inequality is valid: which will be applied with s D k h 2 . We first verify that where lim s!0C h.s/ D 0. Therefore, it is sufficient to prove that h 0 is bounded on OE0; 2 . Obviously, Hence, all components in the expansion of h 0 .s/ are bounded, which really verifies (32). Using (32) in (30) and applying (28) we obtain the following estimation: which completes the proof. Proof. Using Proposition 3.9 we only have to verify the stability of (23). For this, we rewrite it into a linear system Using Proposition 3.6 we obtain that the diagonal of I A˛; h is strictly positive and has non-positive off-diagonals. Moreover, using (24), a simple calculation shows that for the indices j D 1; 2; : : : ; N 1 we have h .I A˛; h / .1; 1; : : : Observe that in the brackets in (34) each coefficient g i ; i 6 D 1 appears once (see the Example and the remark after (26)). According to (8), we obtain˛ and therefore, the sum in both brackets on the right hand side of (34) is zero such that the entire right hand side is one. Using (25) and (26), an obvious modification of (34) gives its positivity both for the indices j D 0 and j D N . Therefore, h .I A˛; h / .1; 1; : : : is valid for all j D 0; 1; : : : ; N . In this way, .I A˛; h / 1 elementwise positive such that k.I A˛; h / 1 k 1 D 1, and consequently, the scheme in (23) is unconditionally stable.

Construction of the matrix A˛; h
Whenever the coefficient in the matrix A˛; h are based on an infinite number of grid points, it can be computed in concrete terms. For where on the right hand side all terms can be computed.
Remark. The statement in Lemma 2.3 offers an alternative to our approach: the numerical solution of the extended problem can be performed with an accurate approximation of the fundamental solution and taking numerical convolution. For a two-dimensional approximation of (4) we refer to [36].

Complexity and extension to higher-order methods
Since we have explored the eigenvectors of A˛; h we do not have to compute its components in practice as a series.
We simply obtain the matrix using (35) such that the extension does not result in extra computational costs.
Higher order methods can be obtained in the same fashion. In the semidiscretization, we should then choose a higher order spatial approximation, e.g., the one in (12) and the accuracy of the time stepping can also be increased, e.g., using a Crank-Nicolson scheme. For homogeneous Dirichlet boundary conditions, such a study is performed (even for the multidimensional case) in [8].

A homogeneous model problem
We first investigate the model problem The analytic solution of (37) on .0; 1/ .0; 1 2 / is which has been computed in the grid points with a high accuracy to verify the convergence of the implicit Euler method. The results of the computations are summarized in Table 1. We computed the error ke ;h k 1 of the approximation in maximum-norm for various time steps and discretization parameters at the final time t D 1. One can clearly see the first order convergence which was predicted by the theory, see Theorem 3.10. The convergence rate was estimated in the consecutive refinement steps using the formula log 2 ke 2 ;2h k 1 ke ;h k 1 Á . The model of (fractional order) diffusion predicts that in case of homogeneous Neumann boundary conditions the total mass should be preserved. Accordingly, in the numerical simulations above, the l 1 norm should be constant, which is an easy consequence of the fact, that the sum of the elements in the columns of A˛; h;1 is zero. Therefore, we examine the boundary conditions in course of the simulations. For this, we use the second order accurate approximation and the results are shown in Figure 3. For the simplicity, we applied the same number of grid points in the spatial and the time coordinates. This accurate approximation can be recognized as the numerical equivalent of Lemma 3.4.

An inhomogeneous model problem
Secondly, we investigate the model problem These have been computed in the grid points with a high accuracy to verify the convergence of the implicit Euler method. Then, according to the Duhamel's principle (see [37], p. 49) we have Accordingly, for the numerical solution at t D 1 we compute first the approximation of u 1 .t; / at the grid points. Then with the same time steps and spatial accuracy we approximate  Table 2, where in all cases D 1=1024 such that the numerical integration does not harm the predicted order of convergence.

Conclusion and future work
In this paper, we have developed and analyzed a finite difference scheme for the numerical solution of onedimensional fractional diffusion problems with Neumann type boundary conditions. The main idea is to extend first the problem to the real line using the boundary data and analyze its numerical solution. Whenever the scheme is given in R, owing to the periodicity of the extension, one can give the numerical solution of the original problem (on a bounded domain) in a conventional matrix-vector form. The corresponding method exhibits optimal convergence rate with respect to the maximum norm. In the analysis we did not need any smoothness assumption. The proposed method was also implemented: the convergence results have been confirmed in the numerical experiments. The most exciting questions for the continuation of this work are the multidimensional generalization of this approach and the treatment of Robin boundary conditions. Proof of Lemma 3.7. We first observe that the Neumann extension is the natural one for v k h in the sense that h c v k h N i j D cos kx j j 2 Z: Then according to (22), (41), (21), (6) and (7) we obtain  The definition of C gives then the statement in the lemma.
Proof of Lemma 3.8. In the proof, we use the identities (43) which can be found in [32], 3.761/9. Observe that the even extension of the cos.k /j .0;1/ function to the real axis is the cos.k / function itself. On the other hand, as it was pointed out in [29], we can differentiate the integrals in the Riemann-Liouville formula to obtain @j xj cos.k x/ D .k / 2 2 cos ˛  On the other hand, the system fcos k xg 1 kD1 is complete in L 2 .0; 1/ such that no further eigenfunctions can exist.