A Fast Implicit Finite Difference Method for Fractional Advection-Dispersion Equations with Fractional Derivative Boundary Conditions

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 porousmedium. 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 O(K) and computational cost of O(K logK). Traditionally, the Gaussian elimination method requires storage of O(K2) and computational cost of O(K3). Finally, the accuracy and efficiency of the method are checked with a numerical example.

For a function () over the interval  <  < , the left-sided Riemann-Liouville fractional derivative of order  is defined by where  is an integer, such that  − 1 <  ≤ .See [12] for details.If  is an integer, then (4) gives the standard integer derivative. = 0 corresponds to a fractional Neumann boundary condition and  > 0 corresponds to a fractional Robin boundary condition.In this paper, we only discuss  > 0.
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 ( log ).In Section 5, a numerical experiment is presented to verify the accuracy and efficiency of the proposed scheme.

Lemma 1.
Let  be a positive real number and the integer  ≥ 1.We have 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  +1  for 1 ≤  ≤ .Then, from (7), we have From ( 8), this implicit finite difference scheme is consistent.Equations ( 7) may be rearranged as follows: where   =   Δ/ℎ  ,   = V  Δ/ℎ.Equations ( 9) can be written in matrix form: where where the coefficient matrix A = [ , ] × and its entries are

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.
Matrix  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 p  (1 ≤  ≤ ) the approximate solution of the difference scheme with the initial condition p0 (1 ≤  ≤ ) 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 which implies that ‖  ‖ ∞ < ‖ 0 ‖ ∞ .Thus, the implicit finite difference method defined by ( 7) is unconditionally stable.
|, 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 By induction, we can finally obtain This completes the proof.

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 ( 2 ) memories.
Proof.We express coefficient matrix  defined by (12) in block form as follows: where  −1,−1 is the ( − 1) × ( − 1) matrix that consists of the first  − 1 rows and the first  − 1 columns.  ,−1 is an ( − 1)-dimensional row vector that consists of the first  − 1 entries in the th row of matrix , and  −1, is an ( − 1)-dimensional column vector that consists of the first  − 1 entries in the th column of matrix .Matrix  −1,−1 has the similar form as that in the finite difference method for (1)-( 3) with the fractional derivative boundary conditions [15]. −1,−1 can be shown to have the following decomposition: where  is the identity matrix of the order −1.And  () and  () are ( − 1)th-order Toeplitz matrix given as follows: According to the definition of ( 12), we can write  −1, and   ,−1 as follows: where  () = ( ,−1 ,  ,−2 , . . .,  ,1 ), and  −1 is the ( − 1)th column of .
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 (FFT).First, we embed Toeplitz matrix   into a 2 × 2 circulant matrix  2 as follows [18,19]: It is known that a circulant matrix  2 can be decomposed as follows [19,20]: where  2 is the first column vector of  2 . 2 is the 2 × 2 discrete Fourier transform matrix in which (, )-entry ( 2 ) , of the matrix  2 is given by ) ,  = √ −1, 0 ≤ ,  ≤ 2 − 2. (34) For any vector   ∈   , we extend the vector by 0 to a 2-dimensional vector  2 = (   0 ).It is well known that the matrix-vector multiplication  2  2 can be carried out in (2 log 2) = ( log ) operations via the fast Fourier transform.Equation (33) shows that  2  2 can be evaluated in ( log ) operations.We observe At the same time,  requires storage of ( 2 ) and computational cost of ( 3 ).In this paper, the system matrix  is supertriangular; the traditional Gaussian elimination method requires storage of ( 2 ) and computational cost of ( 2 ).We present a fast iterative method for the implicit finite difference scheme, which only requires storage of () and computational cost of ( log ).

A Numerical Example
The following fractional differential equation was considered: Note that the exact solution to this problem is (, ) = 5 − (1 − ).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  = 1 and  = 1.0.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  = 1.9,
Let   be the sum of absolute values of elements along the th row excluding the diagonal elements  , .According to Lemma 1, we can get Theorem 2 (solvability).If  > 0, (10) is uniquely solvable.Proof.

Table 1 :
Error behavior for the implicit finite difference scheme (7), for example, at time  = 1 and  = 1.0.
36)Thus, the upper half of the vector multiplication  2  2 gives the matrix-vector multiplication     .So we can evaluate     in ( log ) operations.Then, by (28) we evaluate  −1,−1  −1 in ( log ) operations.Then, from (27), for any vector   ∈   , the matrix-vector multiplication is evaluated in ( log ) operations without any lossy compression.Remark 7. If the system matrix  from (10) is a dense and full coefficient matrix, the Gaussian elimination method