Rational approximation to the fractional Laplacian operator in reaction-diffusion problemsThis work was partially supported by GNCS-INdAM and FRA-University of Trieste.

Rational approximation to the fractional Laplacian operator in reaction-diffusion problemsthanks: This work was partially supported by GNCS-INdAM and FRA-University of Trieste.

Lidia Aceto Department of Mathematics, University of Pisa, Italy Paolo Novati Department of Mathematics and Geosciencies, University of Trieste, Italy

This paper provides a new numerical strategy to solve fractional in space reaction-diffusion equations on bounded domains under homogeneous Dirichlet boundary conditions. Using the matrix transform method the fractional Laplacian operator is replaced by a matrix which, in general, is dense. The approach here presented is based on the approximation of this matrix by the product of two suitable banded matrices. This leads to a semi-linear initial value problem in which the matrices involved are sparse. Numerical results are presented to verify the effectiveness of the proposed solution strategy.

Keyword: fractional Laplacian operator, matrix functions, Guass-Jacobi rule.

MSC: 65F60, 35R11, 65D32.

1 Introduction

Fractional-order in space mathematical models, in which an integer-order differential operator is replaced by a corresponding fractional one, are becoming increasingly used since they provide an adequate description of many processes that exhibit anomalous diffusion. This is due to the fact that the non-local nature of the fractional operators enable to capture the spatial heterogeneity that characterizes these processes.

There are however some challenges when facing fractional models. First of all, there is no unique way to define fractional in space derivatives and, in general, these definitions are not equivalent especially when more than one spatial dimension is considered [16]. In addition, considering that the value of the solution at a given point depends on the solution behavior on the entire domain, it is intuitive to understand that the boundary conditions deserve a particular attention and should be appropriately chosen in order to model the phenomenon properly.

In this paper we consider the following fractional in space reaction-diffusion differential equation


subject to homogeneous Dirichlet boundary conditions


and the initial condition


where represents the diffusion coefficient and the forcing term and are sufficiently smooth functions. The symmetric space fractional derivative of order () is defined through the spectral decomposition of the homogeneous Dirichlet Laplace operator [9, Definition 2]. Assuming that is a Lipschitz domain, the spectrum of is discrete and positive, and accumulate at infinity. Thus,


where are the Fourier coefficients of and are the eigenvalues and the eigenvectors of respectively.

We remark that the fractional power of the Laplace operator is alternatively defined in the literature using the Fourier transform on an infinite domain [12], with a natural extension to finite domain when the function vanishes on and outside the boundary of the spatial domain. In this case, in fact, it is possible to consider non-local problems on bounded domain by simply assuming that the solution of fractional problem is equal to zero everywhere outside the domain of interest. Using such definition and assuming to work with homogeneous Dirichlet boundary conditions, in [14, Lemma 1] it has been proved that the one-dimensional fractional Laplacian operator as defined in (4) is equivalent to the Riesz fractional derivative. Hence, it can be approximated by the ‘fractional centered derivative’ introduced by Ortigueira in [11]. Çelik and Duman in [5] have used such a method for solving a fractional diffusion equation with the Riesz fractional derivative in a finite domain. Moreover, by exploiting the decay of the coefficients characterizing the method, in [10] a ‘short memory’ version of the scheme has been implemented. However, both the original and the modified methods only work for one-dimensional space cases.

A mainstay in the numerical treatment of partial differential problems of type (1)–(3) is to apply the method of lines. Discretizing in space with a uniform mesh of stepsize in each domain direction and using the matrix transfer technique proposed in [8, 9] by Ilić et al., we obtain

where is the approximate matrix representation of the standard Laplacian obtained by using finite difference methods. Consequently, (1) is transformed into a system of ordinary differential equations


where and denote the vectors of node values of and respectively. The matrix raised to the fractional power is, in general, a dense matrix which could be also very large depending on the numbers of mesh points used for the spatial discretization. Therefore, the computational effort for solving (5) could be really heavy, independently of the integrator used. Recently, some authors have developed techniques for reducing this cost. In particular, an approach which can be equally applicable to fractional-in-space problems in two or three spatial dimensions has been considered in [4]. The computational heart of this approach is the efficient computation of the fractional power of a matrix times a vector. However, its effectiveness depends on the mesh discretization.

In this paper, we propose a solution strategy based on a suitable approximation of In particular, we look for a decomposition of the type

where and are both banded matrices arising from a rational approximation of the function based on the Gauss-Jacobi rule applied to the integral representation of cf. [6]. The poles of the formula depends on a continuous parameter whose choice is crucial for a fast and accurate approximation. The above factorization allow to approximate the solution of (5) by solving


By virtue of the structure of the matrices and the numerical solution of (6) may be computed in a more efficient way with respect to the one of (5). We remark that the approach is independent of the Laplacian working dimension.

The paper is organized as follows. In Section 2, the main results about the matrix transfer technique are recalled. Section 3 is devoted to the construction of the rational approximation together with the analysis of the asymptotically optimal choices of the poles. In Section 4 a theoretical error analysis is presented. Numerical experiments are carried out in Section 5, and the conclusions follow in Section 6.

2 Background on the matrix transfer technique

For an independent reading, in this section the basic facts concerning the matrix transfer technique proposed by Ilić et al. in [8, 9] to discretize the one-dimensional fractional Laplacian operator are recalled. In addition, since in this work we also lead with problems in two spatial dimensions, we refer to the results given in [15] as well.

Working with the basic assumption that the fractional Laplacian operator with Dirichlet boundary conditions can be defined as the fractional power of the standard Laplacian, the matrix transfer technique simply consists in approximating the operator through the matrix , where is any finite-difference approximation of on a uniform mesh of size . The only important requirement is that the matrix is positive definite so that its fractional power is well defined. This requirement is fulfilled by the existing standard central difference schemes. Working like that, the original problem (1)–(3) is then transformed into the semi-linear initial value problem


where denotes the vector of node values of

It is important to remark that while is typically sparse, when the matrix loses its sparsity and becomes dense. Observe moreover that the stiffness property of (7) for is essentially inherited by the fractional counterpart so that an implicit scheme or an exponential integrator is generally needed for solving this initial value problem. In both cases the density of may lead to a computational demanding integrator when the discretization is sharp. In order to overcome the limitations in terms of computational efficiency, we propose a strategy based on a suitable approximate factorization of In the next section we focus on the construction of such approximation.

3 Approximation to the matrix fractional power

From the theory of matrix functions (see [7] for a survey) we know that the fractional power of a generic matrix can be written as a contour integral

where is a suitable closed contour enclosing the spectrum of , , in its interior. The following known result (see, e.g., [2]) expresses in terms of a real integral. The proof is based on a particular choice of and a subsequent change of variable.

Proposition 1

Let be such that For the following representation holds

In order to confine the dependence of to a weight function, we consider the change of variable



and hence


that naturally leads to the use of the -point Gauss-Jacobi rule. Such a formula yields a rational approximation of the type


where the coefficients and are given by

in which and are, respectively, the weights and nodes of the Gauss-Jacobi quadrature rule with weight function . Of course, the above approximation can be used in our case with whenever represents the discrete Laplacian operator with Dirichlet boundary conditions, whose spectrum is contained in In the field of fractional calculus, we need to mention that this technique has already been used in [1] for the approximation of the Caputo’s fractional derivative. At this point, denoting by and the polynomials of degree such that we can approximate the solution of (7) by solving (6) with and We remark that the use of the Gauss-Jacobi rule ensures that and for each , and hence it is immediate to verify that the spectrum of is strictly contained in the positive real axis. This condition is fundamental to preserve the stability properties of (7) whenever is replaced by

3.1 Choice of

The choice of the parameter in the change of variable (8) is crucial for the quality of the approximation attainable by (10). Assuming that the generic matrix is symmetric and positive definite, let and be its smallest and largest eigenvalue, respectively. Let moreover It is well known that


In this view, looking at (9), a good choice of is the one that minimizes, uniformly with respect to the error of the Gauss-Jacobi formula when applied to the computation of

From the theory of best uniform polynomial approximation and its application to the analysis of the Gauss quadrature rules (see e.g. [13] for a recent study) it is known that the position of the poles of the integrand function with respect to the interval of integration defines the quality of the approximation. In our case, we observe that for each the poles of the integrand function

are functions of defined by

and we clearly have for , and for . Our aim is to define in order to maximize the distance of the set

from the interval of integration . We observe that for the worst case is given by or since we have respectively


As consequence, the idea is to set such that

that leads directly to the equation

whose solution is


Formally, is given by

In this way, the set is symmetric with respect to the origin, that is where


in which denotes the spectral condition number of This situation is summarized in an example reported in Figure 1.

Figure 1: Example of function for , The choice of as in (12) ensures the symmetry of the set The minimum distance of the curve from the set is given by and is attained in either or

4 Error analysis

In this section we analyze the error of the rational approximation (10) with the choice of in (8). We start with the following result, whose proof is given in [13, Theorems 4.3 and 4.4].

Theorem 1

Let be a function analytic in an open subset of the complex plane containing the ellipse

Let moreover be the polynomial of degree of best uniform approximation of in and




Theorem 2

Let be a symmetric positive definite matrix and Then for sufficiently large, the error of the rational approximation (10), generated by the Gauss-Jacobi rule applied to the integral (9) for is given by

where is a constant independent of and

Proof. For let


Let moreover be the corresponding -point Gauss-Jacobi approximation with weights , By standard arguments we have that


where, since

Now, independently of the choice of makes possible to use the bound (14) for each where solves

since Thus by (15), (14) and using

we obtain


Now, neglecting the factor and then minimizing with respect to yields

Hence, for large enough (we need ), we can use in (16), obtaining


Indeed, defining such that

we have

Moreover, in (17) we have used the inequalities

Finally, since by (9)

using (17) we obtain the result.

Corollary 1

The asymptotic convergence factor fulfils

Proof. By (13)

Remark 1

From the above analysis it is easy to observe that for the Laplacian operator , discretized with standard central differences (3-points or 5-points in one or two dimensions, respectively), we have

where represents the number of discretization points in one dimension.

In Figure 2 we plot the relative error for the one- and two-dimensional Laplacian discretized as in the previous remark for some values of The geometric convergence theoretically proved in this section is clear in the pictures, together with the substantial independence of which is absorbed by the weight function. It is also quite evident that the method is particularly effective for the two-dimensional case; this represents an important feature since most of the standard techniques for the discretization of the fractional Laplacian only work in one dimension.

Figure 2: Relative error of the rational approximation versus the number of points of the Gauss-Jacobi rule, for some values of The one- and the two-dimensional cases are on the left and on the right, respectively. In the first case the dimension of the problem is and in the second one it is

Independently of the above analysis, we observe that the error with respect to can be easily monitored by using (11) with or and hence working scalarly. This consideration suggests a simple strategy for the choice of when the rational approximation is employed to solve (7). Indeed, neglecting for a moment the forcing term, we have that the most important component of the solution is the one corresponding to the smallest eigenvalue, since it leads to the slowest decay rate, given by (here for simplicity). Using our approximation, this decay will be replaced by where

In this way, working scalarly we can easily define in order to keep under control the error factor and hence the error with respect to each component. While this represents a general indication, in the experiments of the next section we have taken very small in order to show the computational efficiency of the method.

5 Solving fractional in space reaction-diffusion problems

As already said in Section 2, if we discretize on a uniform mesh the fractional Laplacian operator occurring in (1), we obtain the initial value problem


Therefore, the application of the rational approximation (10) of based on the -point Gauss-Jacobi rule and given by leads to the following initial value problem


The sparse structure of the matrices and represents the main advantage of this approach in terms of computational work and memory saving.

Following the examples reported in [14], we first focus on two fractional in space diffusion equations with different initial conditions. Then, we consider a reaction-diffusion equation with forcing term independent of the solution. All of these examples are in one spatial dimension. In each case, discretizing the spatial domain with a uniform mesh having stepsize we consider the standard 3-point central difference discretization of the Laplacian Finally, we also report the results obtained by applying our approach for the numerical solution of a fractional reaction-diffusion example in two space dimensions. In this case, we discretize in space the problem via the 5-point finite difference stencil. The matrix is therefore a block tridiagonal matrix of size having the following form with denoting the identity matrix of size and

In all examples, we solve (18) and (19) by the MATLAB routine ode15s. Moreover, we indicate by ‘exact’ the analytical solution, by ‘MT’ the solution of the problem (18), obtained by applying the matrix transfer approach, and by ‘rational’ the solution arising from (19).

Example 1

Consider the problem (1) on the spatial domain with and According to [14, Section 3.1], the analytic solution corresponding to the initial condition is given by

Setting at time in the left-hand side of Figure 3 the exact solution is compared with the numerical solutions of the semi-discrete problems (18) and (19) with (that is ) and On the right picture, the step-by-step maximum norm of the difference between the analytic solution and the numerical ones is reported. As one can see, the numerical solution provided by the rational approximation is in good agreement with the one obtained by the matrix transfer technique.

Figure 3: Comparison of the analytic solution of the problem in Example 1 with the numerical solutions provided by the rational approach and the matrix transfer technique at (left) and corresponding errors (right).

To illustrate the impact of the fractional order in space we consider the following example which differs from the previous one only for the choice of the initial condition.

Example 2

Consider now the problem (1) on the spatial domain with and

In this example, we use and we compute the numerical solutions provided by the matrix transfer technique and the rational approximation approach with In particular, the solutions profiles corresponding to and are shown in Figure 4 at time It is interesting to see that the diffusion depends on the value of the fractional order

Figure 4: Comparison of the numerical solutions of the problem in Example 2 provided by the MT and rational approaches at (left) and corresponding step-by-step maximum norm of their difference (right) for (top) and (bottom), respectively.
Example 3

Consider the problem (1) on the spatial domain with and

The exact solution is given by

In our experiments, we select the model parameters and discretize the spacial domain using In Figure 5 we report the step-by-step error provided by the numerical solutions obtained by applying the Gauss-Jacobi rule with at compared with the one obtained by solving directly (18). As expected, the approximation of the solution provided by the rational approach improves as increases.

Figure 5: Comparison of the errors provided by solving the problem of the Example 3 using both rational with (blue dashed-dot-line), (red dashed-line) and (black dot-line) and MT (green solid-line).
Example 4

We solve the fractional reaction-diffusion equation in two space dimensions



subject to and homogeneous Dirichlet boundary conditions [3].

The exact solution to this problem is

The numerical solution provided by the rational approach based on the Gauss-Jacobi rule with and the matrix transfer technique are drawn at in Figure 6 using and points in each domain direction. It is worth noting that in order to obtain the same accuracy, the matrix transfer technique costs three times the rational approach.

Figure 6: Comparison of the analytical solution of the problem in Example 4 with the numerical solution provided at by rational (top) and MT (bottom) and corresponding relative errors (right) for and

6 Conclusions

In this paper we have proposed a rational approximation to the discrete fractional Laplacian. When applied for solving the reaction-diffusion equations this leads to a semi-discrete problem which can be solved in an efficient way due to the band structure of the matrices occurring in the definition of the approximation. Another advantage of this approach is that it can be generalized to high dimensions without modifying the overall solution methodology.



  • [1] L. ACETO, C. MAGHERINI, P. NOVATI, On the construction and properties of -step methods for FDEs, SIAM J. Sci. Comput., 37 (2015), pp. A653–A675.
  • [2] D. A. BINI, N. J. HIGHAM, B. MEINI, Algorithms for the matrix pth root, Numer. Algorithms, 39 (2005), pp. 349–378.
  • [3] A. BUENO-OROVIO, D. KAY, K. BURRAGE, Fourier spectral methods for fractional-in-space reaction-diffusion equations, BIT, 54 (2014), pp. 937–954.
  • [4] K. BURRAGE, N. HALE, D. KAY, An efficient implicit FEM scheme for fractional-in-space reaction-diffusion equations, SIAM J. Sci. Comput., 34 (2012), pp. A2145–A2172.
  • [5] C. ÇELIK, M. DUMAN, Crank-Nicolson method for the fractional diffusion equation with the Riesz fractional derivative, J. Comput. Phys., 231 (2012), pp. 1743–1750.
  • [6] A. FROMMER, S. GÜTTEL, M. SCHWEITZER, Efficient and stable Arnoldi restarts for matrix functions based on quadrature, Preprint BUW-IMACM 13/17 (2013), http://www.math.uni-wuppertal.de.
  • [7] N. J. HIGHAM, Functions of matrices. Theory and computation, SIAM, Philadelphia, PA, 2008.
  • [8] M. ILIĆ, F. LIU, I. TURNER, V. ANH, Numerical approximation of a fractional-in-space diffusion equation I, Fract. Calc. Appl. Anal., 8 (2005), pp. 323–341.
  • [9] M. ILIĆ, F. LIU, I. TURNER, V. ANH, Numerical approximation of a fractional-in-space diffusion equation (II)-with nonhomogeneous boundary conditions, Fract. Calc. Appl. Anal., 9 (2006), pp. 333–349.
  • [10] M. POPOLIZIO, A matrix approach for partial differential equations with Riesz space fractional derivatives, Eur. Phys. J. Special Topics, 222 (2013), pp. 1975–1985.
  • [11] M. D. ORTIGUEIRA, Riesz potential operators and inverses via fractional centred derivatives, Int. J. Math. Math. Sci., (2006) Article ID 48391, pp. 1–12.
  • [12] S. G. SAMKO, A. A. KILBAS, O. I. MARICHEV, Fractional integral and derivatives: theory and applications, Gordon and Breach Science Publishers, New York, 1987.
  • [13] L. N. TREFETHEN, Is Gauss quadrature better than Clenshaw-Curtis?, SIAM Rev., 50 (2008), pp. 67–87.
  • [14] Q. YANG, F. LIU, I. TURNER, Numerical methods for the fractional partial differential equations with Riesz space fractional derivatives, Appl. Math. Model., 34 (2010), pp. 200–218.
  • [15] Q. YANG, I. TURNER, F. LIU, M. ILIĆ, Novel numerical methods for solving the time-space fractional diffusion equation in two dimensions, SIAM J. Sci. Comput., 33 (2011), pp. 1159–1180.
  • [16] Q. YU, F. LIU, I. TURNER, K. BURRAGE, Numerical investigation of the three types of space and time fractional Bloch-Torrey equations in 2D, Cent. Eur. J. Phys., 11 (2013), pp. 646–665.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description