Fast finite volume methods for space-fractional diffusion equations

Fractional diffusion equations model phenomena exhibiting anomalous diffusion that is characterized by a heavy tail or an inverse power law decay, which cannot be modeled accurately by second-order diffusion equations that is well known to model Brownian motions that are characterized by an exponential decay. However, fractional differential equations introduce new mathematical and numerical difficulties that have not been encountered in the context of traditional second-order differential equations. For instance, because of the nonlocal property of fractional differential operators, the corresponding numerical methods have full coefficient matrices which require storage of $O(N^2)$ and computational cost of $O(N^3)$ where $N$ is the number of grid points. 
   
We develop a fast locally conservative finite volume method for a time-dependent variable-coefficient conservative space-fractional diffusion equation. This method requires only a computational cost of $O(N \log N)$ at each iteration and a storage of $O(N)$. Numerical experiments are presented to investigate the performance of the method and to show the strong potential of these methods.


1.
Introduction. Classical diffusion equations describe the spreading of a cloud of tracer particles, which can be modeled by Fick's law or equivalently the random process of Brownian motion. The underlying assumption is that a particle's motion has little or no spatial correlation, so the variance of a particle excursion distance is finite. Then the central limit theorem concludes that the probability of finding a particle somewhere in space satisfies a Gaussian distribution, which gives a probabilistic description of a diffusion precess. In last few decades more and more diffusion processes, which range from the signalling of biological cells to the solute transport in groundwater, were found not Fickian. In these processes, a particle's long walk is not necessarily independent of each other and may have long correlation length so the processes can have large deviations from the stochastic process of Brownian motion or even have an infinite variance. Fractional diffusion equations provide a more adequate and accurate description of these processes than classical second-order diffusion equations do [4,24,25].
In recent years extensive research has been carried out on the development and analysis of numerical methods for fractional partial differential equations [3,6,8,9,14,16,17,18,19,22,23,29,30]. However, fractional differential equations introduce new mathematical and numerical difficulties that have not been encountered in the context of traditional second-order differential equations. Because of the relative simplicity of the Grünwald-Letnikov fractional derivatives, finite difference methods are probably the most developed and widely used methods for fractional diffusion equations [22,23]. However, finite difference methods also have serious limitations: (i) Since they are derived based upon the Grünwald-Letnikov fractional derivatives, finite difference methods require a uniform partition. (ii) A finite difference method derived via a naive truncation of the infinite series of the Grünwald-Letnikov fractional derivatives is unconditionally unstable [22]. A shifted Grünwald-Letnikov discretization was used in finite differenmce methods, which ensures the stability and convergence of the numerical methods but leads to only first-order spatial accuracy. A Richardson extrapolation has been proposed to retain second-order spatial accuracy, which requires a solution at a finer spatial mesh [32]. (iii) Finite difference methods apply to nonconservative space-fractional diffusion equations.
The governing equations arising in applications such as solute transport in groundwater hydrology are conservative space-fractional diffusion equations [4], which are not equivalent to nonconservative space-fractional diffusion equations for variable diffusivity coefficients. While finite difference methods do not apply to conservative space-fractional diffusion equations, the Galerkin finite element methods and finite volume methods directly apply. In [10] Ervin and Roop proved the wellposedness of a Galerkin weak formulation to a conservative space-fractional diffusion equation with a constant diffusivity coefficient and a homogeneous Dirichlet boundary condition. In [40] we presented a counterexample to show that the the Galerkin weak formulation loses its coercivity for a variable-coefficient steady-state conservative space-fractional diffusion equation. A Petrov-Galerkin weak formulation was proposed for the problem and was proved to be weakly coercive. Hence the Petrov-Galerkin formulation is well posed. Numerical experiments show that Galerkin finite element methods do not necessarily converge for a steady-state conservative space-fractional diffusion equation with a piecewise-contant diffusivity coefficient, but do converge for the same equation with a smoothly varying diffusivity coefficient [41].
In this paper we develop a fast finite volume method for the efficient numerical solution of time-dependent variable-coefficient conservative space-fractional diffusion equations. The movitation is as follows: Finite volume methods (i) conserve mass locally, which is a crucial property in many applications; (ii) are derived based upon Riemann-Liouville or Caputo fractional derivatives, and so can naturally be used on nonuniform meshes, which is particularly suited to handle the singularity of space-fractional diffusion equations; and (iii) retain second-order space accuracy without an extrapolation needed [2,13,20,41]. However, because of the nonlocal property of fractional differential operators, finite volume methods, like the finite difference methods and finite element methods, have full stiffness matrices that require O(N 2 ) of memory to store and O(N 3 ) of computational cost to invert per time step where N is the number of grid points. This imposes computational challenges in the application of fractional diffusion equations, especially in the context of three dimensional problems.
In [39] we first broke this bottleneck in the context of finite difference methods by decomposing the stiffness matrix as a combination of Toeplitz-like matrices [27]. Consequently, we significantly reduced the storage requirement from O(N 2 ) to O(N ) and a computational cost from O(N 2 ) to O(N log N ) per iteration. As the derivation was conducted without any compression, but by fully exploring the structure of the stiffness matrix, the resulting fast method is not lossy. In this paper we utilize a similar idea to develop a fast and faithful finite volume method for the efficient numerical solution of time-dependent variable-coefficient conservative space-fractional diffusion equations. The rest of the paper is organized as follows. In section 2 we present the fractional diffusion equation and derive its finite volume approximation. In section 3 we study the matrix structure and efficient storage of these finite volume methods. In section 4 we derive a fast solution method based on the conjugate gradient squared method (CGS). In section 5 we carry out numerical experiments to study the spatial and temporal convergence behavior of the finite volume methods. We also compare their performance with that of the well established finite difference methods. In section 6 we draw some concluding remarks.

2.
Model problem and its finite volume discretization. We consider the initial-boundary value problem of a fractional diffusion equation with an anomalous diffusion of order 2 − β with 0 < β < 1 [4,21,26,28]  are the left and right fractional integral operators defined by [26,28] where Γ(·) is the Gamma function.
We define a temporal partition on [0, T ] by t m = m∆t for m = 0, 1, . . . , M with ∆t = T /M . We then approximate the time derivative in Eq. (2.1) by a time difference quotient and drop the local truncation error term to write the governing equation in Eq. (2.1) at time instant t m,θ = t m−1 + θ∆t for as follows: +f (x, t m,θ ).

(2.3)
Here θ ∈ [0, 1] is a weighting parameter. In particular, (2.3) with θ = 1, 1/2, and 0 corresponds to the backward Euler, Crank-Nicolson, and forward Euler scheme, respectively. We focus on the backward Euler and Crank-Nicolson schemes in this paper, although the method developed in this paper also works for any value of θ ∈ [0, 1] in particular the forward Euler scheme. We define a uniform spatial partition on [a, b] by x j = jh for j = 0, 1, . . . , N + 1 with h = (b − a)/(N + 1). We let x j− 1 2 = (x j−1 + x j )/2 for j = 1, 2, . . . , N + 1 to be the mid point of the interval [x j−1 , x j ]. Then we integrate the governing equation on (x j− 1 2 , x j+ 1 2 ) for j = 1, 2, . . . , N to obtain the following weak formulation for Eq. (2.3): We let S h (a, b) denote the space of continuous and piecewise-linear functions with respect to the spatial partition, which vanish at the boundary x = a and The corresponding finite volume scheme can be formulated as (2.7) The finite volume schemes (2.6) can be written in a matrix form 3. Properties of the coefficient matrix and its efficient storage. The finite volume schemes (2.6) formally look similar to their counterpart for the classical diffusion equation. In fact, the mass matrix has exactly the same form The fundamental difference comes from the stiffness matrix B m . Even though the support of Thus, the stiffness matrix B m is a full matrix. It is well known that inverting a general full matrix requires O(N 3 ) operations while storing such a matrix requires O(N 2 ) storage. The significant increase of computational cost and storage requirement of numerical methods for fractional diffusion equations has been one major reason which seriously limits their applications, even though fractional diffusion equations provide an adequate and accurate description of anomalous diffusion processes.
To develop an efficient algorithm for the finite volume schemes (2.6), we have to explore the structure of the stiffness matrix B m carefully. We go through algebraic manipulations to find out that the entries of the stiffness matrix B m are of the form: are defined as follows: In other words, the stiffness matrix B m has the following decomposition . . , N . B L and B R are Toeplitz matrices of order N given by In order to store the stiffness matrix B m we need only to store K(x i− 1 2 , t m ) for i = 1, . . . , N + 1, g (β) = [g  (5N −1) parameters. This is in contrast the storage of the full matrix B m which requires storing N 2 parameters. In particular, the vector g (β) and q (β) are independent of space or time and can be stored a priori for a given β, γ, and N .

4.
A fast solution method for the finite volume schemes. In this section we base on the features of the stiffness matrix B m to develop a fast iterative algorithm for the efficient solution of the finite volume scheme. A fast O(N log N ) algorithm for the evaluation of B m v. Note that the Toeplitz matrices B L and B R can be embedded into 2N × 2N circulant matrices C 2N,L and C 2N,R as follows [5,12]

4.1.
It is known that a circulant matrix C 2N can be decomposed as follows [7,12] where c is the first column vector of C 2N and F 2N is the 2N × 2N discrete Fourier transform matrix in which the (j, l)-entries F 2N (j, l) of the matrix F 2N are given by The fact that the stiffness matrix-vector multiplication B m v N can be evaluated in O(N log N ) operations motivates that we should use a conjugate gradient type method to solve the system. Since the stiffness matrix B m is not symmetric, we need to use a nonsymmetric version of conjugate gradient method. To explore the idea, we use the conjugate gradient squared method (CGS) as an example to solve finite volume scheme Eq. (2.8) [1]: At each time step t m , we choose u (0) = u m−1 Compute r (0) = f m − A m u (0) Chooser (for example,r = r (0) ) for i = 1, 2, · · ·

Numerical experiments.
In this section we carry out numerical experiments to investigate the performance of the backward Euler and Crank-Nicolson finite volume methods developed in this paper. These experiments serve for two purposes: (i) We observe the spatial and temporal convergence behavior of these finite volume methods. (ii) We investigate the reduction of the computational cost and storage of the fast conjugate gradient squared method (FCGS) over the standard method for fractional diffusion equations. To gain a better understanding of these newly developed methods, we compare them with the well developed backward Euler and Crank-Nicolson finite difference methods [22,23,32]. In the numerical experiments we use the following widely tested examples.
Example 5.1. The fundamental solution of the fractional diffusion equation. We solve a constant-coefficient homogeneous fractional diffusion equation (2.1) which is subject to an initial condition of a Dirac δ function located at x = 0. The analytical solution u(x, t) of this problem is expressed in terms of the inverse Fourier transform [4] u(x, t) = 1 π Example 5.2. Gaussian pulse subject to an anomalous diffusion. We simulate an anomalous diffusive process of a Gaussian pulse subject to the homogeneous fractional diffusion equation (2.1) with the initial condition In the numerical experiments, the data is chosen as follows: β = 0.5, γ = 0.5, (x L , x R ) = (0, 2), (0, T ) = (0, 1), the mean x c = 1.0, the standard deviation σ = 0.08 and f (x, t) ≡ 0. The diffusion coefficient is chosen to be: K(x, t) = 0.002(1 + x(2 − x) + t 2 ) .

5.1.
Spatial and temporal convergence of the finite volume methods. We apply the backward Euler and Crank-Nicolson finite volume methods (2.6) and the corresponding finite difference methods in [22,23,32] to solve Examples 5.1 and 5.2 to investigate their spatial and temporal convergence behavior in the context of fractional diffusion equations. We use a linear regression to fit the convergence rates and the associated constants in the error estimates where u h refers to the numerical solution. We use u F V , u CF V , u F D , and u CF D to denote the numerical solutions by the backward Euler finite volume method, the Crank-Nicolson finite volume method, the backward Euler finite difference method, and the Crank-Nicolson finite difference method, respectively. We start with Example 5.1. In the context of the backward Euler finite volume method, we expect α = 2 and ν = 1 in (5.21). We accordingly perform two types of computations in the numerical experiments. The first tests the spatial convergence rate of the finite volume scheme, where we fix a small time step ∆t and compute the constant M α and the rate α with respect to the spatial mesh size h; the other tests the temporal convergence rate of the finite volume scheme, where we choose a small Table 3. The L 2 errors of the numerical solutions by the backward Euler finite volume method (FV), Crank-Nicolson finite volume method (CFV), the backward Euler finite difference method (FD), and Crank-Nicolson finite difference method (CFD) at time t = T for Example 5.1   Tables 1 and 2. We clearly observe a second-order convergence rate in space and a first-order convergence rate in time in both the L 2 norm and the L ∞ norm. As for the Crank-Nicolson finite volume method, we anticipate that α = ν = 2 in (5.21). This suggests a choice of h = ∆t. In this case (5.21) reduces to We present the L 2 and the L ∞ errors of the numerical solutions by the backward Euler finite volume method, the Crank-Nicolson finite volume method, the backward Euler finite difference method, and the Crank-Nicolson finite difference method in Tables 3 and 4, respectively. We observe the second-order spatial and temporal convergence rates of the Crank-Nicolson finite volume method. We observe only the first-order accuracy for the remaining methods. The reason that the Crank-Nicolson finite difference method is only first-order accurate lies in the fact that it is second-order accurate in time but only first-order accurate in space. We now turn to Example 5.2. In this case no closed-form analytical solution is available. Hence we have to use the Crank-Nicolson finite volume method with extremely fine time step ∆t = 2 −13 and fine spatial grid h = 2 −13 to generate a     Table 7. L 2 and L ∞ errors of the finite volume method (FV) and Crank-Nicolson finite volume method (CFV) for Example 5.2 fine-scale numerical solution as a reference solution u. In addition, we note that the mass-conservative finite volume methods were developed for a fractional diffusion equation in a conservative form [9,10], while the finite difference methods were developed for a fractional diffusion equation in a nonconservative form [22,23,32]. These two equations are not equivalent in general for a variable diffusion coefficient. So we can run the numerical experiments for both the backward Euler and Crank-Nicolson finite volume methods in a similar manner as in Example 5.1. We present the corresponding L 2 and L ∞ errors of the numerical solutions in Tables 5-7. These results show that the backward Euler finite volume method possesses second-order accuracy in space and first-order accuracy in time, while the Crank-Nicolson finite volume method has second-order accuracy in both space and time. This is in contrast to the backward Euler and Crank-Nicolson finite difference methods for fractional diffusion equations, both of which have only first-order accuracy in space [22,23,32].  [22,23,32]. In Table 8 we present the CPU time consumed by these methods for the same spatial mesh sizes and time steps as used in Tables 3 and 4. We notice that with a spatial mesh size of h = 2 −9 , i.e., 4096 nodes and 512 time steps, the CPU times consumed by the backward Euler and Crank-Nicolson finite difference methods with Gaussian elimination are about 8.7 × 10 3 seconds. In contrast, the backward Euler and the Crank-Nicolson finite volume methods with the fast CGS acceleration take 22 and 28.44 seconds, respectively, which represent a saving of CPU of more than 300 times. If we take the corresponding errors in Tables 3  and 4 into account, we note that the Crank-Nicolson finite volume method with h = ∆t = 2 −7 is already less than those by the backward Euler and Crank-Nicolson finite difference methods with h = ∆t = 2 −9 . Thus, to reach the same accuracy, the Crank-Nicolson finite volume method with the fast CGS acceleration has a CPU saving of more than 10,000 times! We now look at the storage savings. With a spatial mesh size h = 2 −9 or N = 4096 nodes, the backward Euler and Crank-Nicolson finite difference methods require a storage of O(N 2 ) or 1.68×10 7 parameters. The Crank-Nicolson finite volume method with the fast acceleration requires no more than 20N storage and generates a more accurate numerical solution with h = 2 −7 . That is, the method requires the storage around 10,000 parameters. This represents a saving of the storage about 1,680 times! More importantly, as the problem size N increases the savings on both the computational cost and the storage increases significantly too.
We now turn the issue of CPU and storage for Example 5.2. As we explained in section 5.1, we cannot compare the finite volume methods with the finite difference methods for Example 5.2. In this case we compare the performance of the finite volume methods with or without the fast CGS acceleration. With the same number of nodes N = 2096, the savings of CPU are between 70-110 times and the saving of storage is more than 100 times. The savings of computational cost and storage will increase significantly as the problem size N increases.
6. Concluding remarks. In this paper we developed a fast finite volume method for the efficient numerical solution of time-dependent variable-coefficient conservative space-fractional diffusion equations. Finite volume methods naturally apply to conservative fractional diffusion eqautions, conserve mass locally, and retain secondorder space accuracy without an extrapolation needed [13,20,41]. Furthermore, finite volume methods allow the use of nonuniform meshes and so can be used on a locally refined mesh to treat the singularity of fractional diffusion equations [13,41]. We utilized our idea in [39] to develop a fast and faithful solution algorithm for the finite volume method, which does not require any compression but by carefully exploring the structure of the stiffness matrix and hence is not lossy. This algorithm significantly reduced the storage requirement from O(N 2 ) to O(N ) and a computational cost from O(N 2 ) to O(N log N ) per iteration.
In the context of steady-state space-fractional diffusion equations, the stiffness matrix of corresponding numerical methods is ill-posed. Even though the computational cost of each iteration is reduced significantly, the number of iterations could still be large. In [35] we developed a full symmetric and positive-definite Toeplitz matrix preconditioner for the variable-coefficient conservative Riesz spacefractional diffusion equation. The development was relied heavily on the symmetry of the fractional derivatives. Numerical experiments show that the preconditioner is optimal. For a general time-dependent two-sided varaible-coefficient conservative space-fractional diffusion equations, the coefficient matrix is not Toeplitz. The development of an efficient preconditioner is currently under investigation.
We now turn to time-dependent variable-coefficient space-fractional diffusion equations in multiple space dimensions. In [37,38] we developed fast alternatingdirection implicit finite difference methods for time-dependent nonconservative space-fractional diffusion equations in two and three space dimensions, respectively. The fundamental idea is that via the alternating-direction approach, the multidimensional finite difference methods are reduced to families of independent one-dimensional systems, so the fast solultion algorithms for one-dimensional problems apply. Numerical results show the strength of these methods. Nevertheless, alternating-direction implicit methods have proved stability and convergence only if the diffusivity coefficients have certain special forms [31,33,37]. To overcome these limitations, in [34,36] we developed fast solution algorithms for the truly multidimensional finite difference methods for time-dependent nonconservative spacefractional diffusion equations in two and three space dimensions, respectively. Numerical results show the utility of these methods. We plan to follow the same ideas to develop fast solution methods for finite volume methods for conservative space-fractional diffusion equations in the future.