Abstract
Fractional advection-dispersion equations, as generalizations of classical integer-order advection-dispersion equations, are used to model the transport of passive tracers carried by fluid flow in a porous medium. In this paper, we develop an implicit finite difference method for fractional advection-dispersion equations with fractional derivative boundary conditions. First-order consistency, solvability, unconditional stability, and first-order convergence of the method are proven. Then, we present a fast iterative method for the implicit finite difference scheme, which only requires storage of and computational cost of . Traditionally, the Gaussian elimination method requires storage of and computational cost of . Finally, the accuracy and efficiency of the method are checked with a numerical example.
1. Introduction
Fractional derivatives are almost as old as their more familiar integer-order counterparts [1, 2]. Fractional advection-dispersion equations are used to model many problems in physics, biology, and finance [3–5]. Fractional differential equations have attracted fantastic attention of many authors in recent years [6–11].
The fractional space derivatives are used to establish anomalous diffusion and dispersion models, where the propagation velocity of the particle beam is different from the classical Brown motion model. When the space second derivative in a diffusion or dispersion model is replaced with a fractional space derivative, it leads to superdiffusion. In this paper, we consider a one-dimensional fractional advection-dispersion equation with fractional derivative boundary conditions:where ; is the distribution function of the diffusion quantity; is flow velocity; is diffusion coefficient; is the source term. We assume and , so that the flow is from left to right. The space fractional partial derivatives in (1) and (2) are left-sided Riemann-Liouville fractional derivatives over the domain. For a function over the interval , the left-sided Riemann-Liouville fractional derivative of order is defined by where is an integer, such that . See [12] for details. If is an integer, then (4) gives the standard integer derivative. corresponds to a fractional Neumann boundary condition and corresponds to a fractional Robin boundary condition. In this paper, we only discuss .
For a one-dimensional advection-dispersion equation with constant coefficients, Benson et al. [13] got available analytical solutions using Fourier transform method. However, many problems require a model with variable coefficients. Meerschaert and Tadjeran [14] developed finite difference approximations for one-dimensional fractional advection-dispersion equations with Dirichlet boundary conditions. And, recently, some authors have discussed numerical method for fractional equations with fractional derivative boundary conditions. Jia and Wang [15] developed fast implicit finite difference methods for two-sided space fractional diffusion equations with fractional derivative boundary conditions. Guo et al. [16] developed implicit finite difference method for fractional percolation equation with Dirichlet and fractional boundary conditions. As far as we know, the study on numerical method for one-dimensional advection-dispersion equations with fractional derivative boundary conditions is still limited. This motivates us to examine a fast implicit finite difference method for fractional advection-dispersion equations with fractional derivative boundary conditions in this paper.
The rest of the paper is organized as follows. In Section 2, we develop an implicit finite difference method and analyze its consistency. In Section 3, the solvability, unconditional stability, and convergence of the method are proven. In Section 4, we present a fast iterative method for the implicit finite difference scheme, which only requires storage of and computational cost of . In Section 5, a numerical experiment is presented to verify the accuracy and efficiency of the proposed scheme.
2. The Implicit Finite Difference Method for the Advection-Dispersion Equations with Fractional Derivative Boundary Conditions and Its Consistency
First, let and be the spatial mesh width and the time step, where is a positive integer. We define a spatial partition by for and a temporal partition by for . Let , , , , . Denote the exact and numerical solutions at the mesh point by and , respectively. and .
Then, we employ the backward Euler scheme to approximate the first-order time derivative; a right-shifted Grünwald formula to approximate left-sided -order Riemann-Liouville fractional derivative and standard Grünwald formula to approximate -order Riemann-Liouville fractional derivative are defined as follows [12]:where is normalized Grünwald weights:Moreover, the normalized Grünwald weights satisfy the properties in the following Lemma 1.
Lemma 1. Let be a positive real number and the integer . We have(i), for ;(ii), for ;(iii), for .
Therefore implicit finite difference method for the advection-dispersion equations with fractional derivative boundary conditions (1)–(3) can be written as follows:
Denote the local truncation error by for . Then, from (7), we have From (8), this implicit finite difference scheme is consistent.
Equations (7) may be rearranged as follows:where , . Equations (9) can be written in matrix form: wherewhere the coefficient matrix and its entries are
3. Solvability, Stability, and Convergence Analysis
Having developed a numerical scheme and shown that it is consistent, in the following theorems, we show this method is solvable, unconditionally stable, and convergent.
Theorem 2 (solvability). If , (10) is uniquely solvable.
Proof. Let be the sum of absolute values of elements along the th row excluding the diagonal elements . According to Lemma 1, we can getMatrix is strictly diagonal, which guarantees the invertibility of matrix . Hence, (10) is uniquely solvable.
To discuss the stability of the numerical method, we denote by the approximate solution of the difference scheme with the initial condition and define , , .
Theorem 3. The implicit finite difference method for the advection-dispersion equations with fractional derivative boundary conditions defined by (7) is unconditionally stable.
Proof. From (9) and the definition of , we have By Lemma 1, we have It follows from (15) that Therefore, we suppose that , where . By Lemmas 1(i) and (iii), then which implies that .
Thus, the implicit finite difference method defined by (7) is unconditionally stable.
For the convergence analysis, we define and , .
Theorem 4. There exists a positive constant , independent of and , such that
Proof. From (9), we have with , and , . If , we have By Lemma 1(ii), then With Lemma 1 and the Stirling formula for the Gamma function [17], we have as . Combining (22) and (23), we get If , , then By induction, we can finally obtain This completes the proof.
4. Storage and Fast Krylov Subspace Method of the Coefficient Matrix
The following theorem shows that the storage of matrix can be stored in memories, instead of memories.
Theorem 5. Coefficient matrix can be stored in memories.
Proof. We express coefficient matrix defined by (12) in block form as follows: where is the matrix that consists of the first rows and the first columns. is an -dimensional row vector that consists of the first entries in the th row of matrix , and is an -dimensional column vector that consists of the first entries in the th column of matrix . Matrix has the similar form as that in the finite difference method for (1)–(3) with the fractional derivative boundary conditions [15]. can be shown to have the following decomposition: where is the identity matrix of the order . And and are th-order Toeplitz matrix given as follows:According to the definition of (12), we can write and as follows: where , and is the th column of .
Thus, in order to store , we need only to store , , , , , and , which totally have parameters, instead of storing the super triangular matrix in a straightforward but naive way that required memories.
Theorem 6. For any vector and to be a Toeplitz matrix of order , the matrix-vector multiplication can be carried out in operations without any lossy compression.
Proof. A Toeplitz matrix is a matrix in which each descending diagonal from left to right is constant. Let be a Toeplitz matrix of order of the following form: The matrix-vector multiplication can be carried out via the fast Fourier transform . First, we embed Toeplitz matrix into a circulant matrix as follows [18, 19]: It is known that a circulant matrix can be decomposed as follows [19, 20]: where is the first column vector of . is the discrete Fourier transform matrix in which -entry of the matrix is given by For any vector , we extend the vector by 0 to a -dimensional vector . It is well known that the matrix-vector multiplication can be carried out in operations via the fast Fourier transform. Equation (33) shows that can be evaluated in operations. We observe At the same time, Thus, the upper half of the vector multiplication gives the matrix-vector multiplication . So we can evaluate in operations. Then, by (28) we evaluate in operations. Then, from (27), for any vector , the matrix-vector multiplication is evaluated in operations without any lossy compression.
Remark 7. If the system matrix from (10) is a dense and full coefficient matrix, the Gaussian elimination method requires storage of and computational cost of . In this paper, the system matrix is supertriangular; the traditional Gaussian elimination method requires storage of and computational cost of . We present a fast iterative method for the implicit finite difference scheme, which only requires storage of and computational cost of .
5. A Numerical Example
The following fractional differential equation was considered:
On a finite domain . The flow velocity is . Diffusion coefficient is . The source term is . The initial function is . And the boundary conditions are .
Note that the exact solution to this problem is . In the numerical experiments, we consider three different in the case, respectively. Table 1 shows the maximum error of the numerical solution of at the time and . The stability and its order of convergence are proved.
Figures 1–3 show the numerical solutions obtained by applying the Euler method (10) discussed above, with , , and respectively, with , and . The numerical solution compares well with the exact analytic solution to the fractional partial differential equation in these test cases.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This work was supported by the National Science Foundation of China under Grants 61375063, 61271355, 11301549, and 11271378.