Abstract

In this paper, we concentrate on the efficient solvers for the time-space fractional advection-diffusion equations. Firstly, the implicit finite difference schemes with the shifted Grünwald-Letnikov approximations for spatial fractional derivative and unshifted Grünwald-Letnikov approximations for time fractional derivative are employed to discretize time-space fractional advection-diffusion equations. The discretization results in a series of large dense linear systems. Then, a banded preconditioner is proposed and some theoretical properties for the preconditioning matrix are studied. Numerical implementations show that the banded preconditioner may lead to satisfactory experimental results when we choose appropriate bandwidth in the preconditioner.

1. Introduction

Consider the following initial-boundary value problem of time-space fractional advection-diffusion equations (TS-FADE):where is the source and sink term, the coefficient of anomalous diffusion satisfying and , the real constant denotes the drift of the process, that is, the mean advection velocity, and and are the differentiation parameters.

The Caputo time derivative of real order with in (1) is defined as [1, 2]The right- and left-hand spatial fractional derivatives are defined by [3]respectively, where is the gamma function.

The TS-FADE of form (1) arises in variety of research areas such as modeling chaotic dynamics of classical conservative systems [4], turbulent flow [5], groundwater contaminant transport [6], and applications in finance [7], image processing [8], physics [9], biological systems [10], and so on. Though many analytic approaches, such as the Fourier transform method, the Mellin transform method, and the Laplace transform methods, have been used to seek the closed-form solutions [11], there are few available analytical closed-form solutions for fractional differential equations (FDEs). Hence, it is important to find efficient methods to solve fractional differential equations.

Traditional methods for solving FDEs generate the computational cost of and storage of as a result of full coefficient matrices, where is the grid point number. However, when we use the shifted Grünwald-Letnikov discretization scheme proposed by Meerschaert and Tadjeran [12] to approximate the FDE, we will obtain a special Toeplitz-like discretized coefficient matrix [13, 14]. Therefore, according to [13], the storage requirement is instead of , and the complexity of the matrix-vector multiplication only requires operations by fast Fourier transform (FFT). By making use of this advantage, K. Wang and H. Wang [14] showed that the conjugate gradient normal residual (CGNR) method which has computational cost of is very efficient when the diffusion coefficients are very small; that is, the discretized systems are well-conditioned. However, when the diffusion coefficients are not small, the discretized diffusion coefficient matrices become ill-conditioned. Therefore, preconditioning technique must be used to improve the efficiency of the CGNR method. Many references proposed efficient methods. For example, Pang and Sun [15] proposed the multigrid (MG) method to solve the system. Lei and Sun [16] used the circulant preconditioned CGNR method, which is proved to be more efficient than the MG method. After that, Qu et al. [17] introduced a circulant and skew-circulant splitting method to solve the fractional diffusion equations. For more accuracy analysis of the discretization from the fractional difference equations, we refer to the documents [18–20]. Some efficient algorithms for the discretized linear systems from the space fractional difference equations can be found in [16, 21–24] and so on.

In this paper, we will first discretize TS-FADE (1) by using the unshifted Grünwald-Letnikov approximation to the Caputo time fractional derivative and the shifted Grünwald-Letnikov approximation to the spatial fractional derivative. Then this discretization will lead to a system with Toeplitz-like coefficient matrix. Following the strategy proposed in [21], we will propose a banded preconditioning technique for solving the discretized system. Some theoretical results about the preconditioning matrix will be presented. Numerical experiments will also be given to testify the effectiveness of the new preconditioning technique.

The remainder of this paper is organized as follows. In Section 2, we give a detailed description for the discretization of TS-FADE (1). In Section 3, a banded preconditioner is constructed. Then some properties of the discretized system (16) are given and some spectral properties for the preconditioning matrix are studied. Numerical experiments are presented in Section 4 to illustrate the efficiency of the banded preconditioning technique for solving the TS-FADE. Finally, in Section 5, we draw some conclusions to end this paper.

2. Discretization of the TS-FADE

Let and be positive integers and and be the sizes of spatial grid and time step, respectively. We first perform discretization in time using the unshifted Grünwald-Letnikov approximation to the Caputo time fractional derivative . By making use of the abbreviations and , where , ,   ,  and , we can obtain the following semidiscretized scheme:where is identity matrix andrepresents the discrete Caputo derivative. Here is defined asBecause of the zero initial condition, we can rewrite (4) asand here the first line of (7) has the formWe use the central difference approximation to and approximate the spatial derivative of order by the shifted Grünwald-Letnikov approximations [12] for the right- and left-handed spatial fractional derivative; that is,

Denote ,and .According to the boundary condition , the space-time discretization of (7) will lead to the following algebraic linear system in Kronecker form:where using the MATLAB notation, and each is a vector of dimension associated with the point in time. Meanwhile, the right-hand side vector in (13) is defined by .

LetThen we can rewrite the linear system (13) into a matrix equation formwith defined in (14) being a sum of some dense diagonal-times-Toeplitz matrices and being upper triangular matrix.

Note that is a sum of some diagonal-times-Toeplitz matrices and is upper triangular matrix, the right-hand side matrix in (16) is not low rank, and the size of usually has much smaller size than that of . Then the first column of can be obtained by solving the shifted system . All subsequent columns of may be obtained with some substitutions aswhere ; that is, denotes the th-entry of .

Since . Then we have to solve these systems with the same coefficient matrix , where and are positive diagonal matrices, is an identity matrix, and and are described previously.

3. Banded Preconditioning Iteration Method

Definewhere and the matrices , , , and are defined in Section 2. It can be easily found that is a skew-symmetric matrix.

Obviously, the main work for solving the discretized system (17) concentrates on solving the system with coefficient matrix for . Since is a nonsymmetric matrix, then we can use the Krylov subspace methods [25], such as the generalized minimal residual (GMRES) method, to solve the linear systems (17). The preconditioning techniques are often employed to improve the performance and reliability of the Krylov subspace methods.

Denoteand define a banded preconditioner asIt can be easily seen that is a banded Toeplitz matrix with bandwidth containing the central diagonals of . Therefore, is a banded matrix with bandwidth .

To study the property of the preconditioner , we will give the following lemmas first.

Lemma 1 (see [12]). Let and be defined in (6). Then we have

Lemma 2. For , the matrix defined in (19) is a strongly diagonally dominant matrix and the real part of the spectrum of is contained in the open interval .

Proof. According to Lemma 1, it holds that . Therefore, for each row of the matrix , we have or, equivalently, for . That is, (). Hence, is a strongly diagonally dominant matrix.
Moreover, by simple computation, we can obtain the first column of the matrix asHence, the Gershgorin circles of are centered at . Moreover, according to Lemma 1, the largest radius is at mostIt then follows that the real part of the spectrum of is contained in the open interval .

Next, we will concentrate on the uniqueness of the solution of the matrix equation (16).

Theorem 3. Let , where is defined as in (11). Then is a symmetric positive definite matrix.

Proof. Suppose that are the eigenvalues set of the matrix . According to the Gershgorin theorem, all eigenvalues of are contained within a complex disk centered at , with the maximum radius being in line or in line (dependent on whether is odd or even). Therefore, the maximum radius satisfiesThat is, .
Because is a real symmetric matrix, then the eigenvalues are all real. So is symmetric positive definite.

Theorem 4. Suppose that defined in (16) satisfies , and then there exists only one solution for the Sylvester equation (16).

Proof. From [26, 27], we know that when , the continuous Sylvester equation (16) has a unique solution if there is no common eigenvalue between and . It is easy to know that is upper triangular matrix and its eigenvalues are all positive. Note that +. According to Theorem 3, if the real parts of the eigenvalues of are positive, then the real parts of all eigenvalues of are also positive. Furthermore, the matrices and are both symmetric positive definite. Therefore, is symmetric positive definite, and the real parts of the eigenvalues of are positive. Therefore, there is no common eigenvalues between the matrices and .

Theorem 5. For , the preconditioner defined in (20) is nonsingular.

Proof. Since is a positive diagonal matrix and from Lemma 2, the real part of the spectrum of is positive. Then the real part of the spectrum of is also positive. Using the same strategy, we can obtain that the real part of the spectrum of is also positive.
Assume that is a singular matrix and there exists a nonzero eigenvector such that . Then it follows thatMultiplying by the left of the above equation, we getAs is a skew-symmetric matrix and is a real number, then is a pure imaginary. Because the real part of is a positive number, then the real part of is a positive number. According to (26), we can easily find that there is a contradiction because the real part of the right hand side of (26) is zero. That is, if , it must hold that ; that is, the matrix is nonsingular.

Theorem 6. For and , the relative difference between and can be described aswhich implies can be an efficient preconditioner for the coefficient matrix as increases.

Proof. The proof process is similar to [13].

At the end of this section, we will draw the main steps for solving the discretized system (13).

Step 1. Compute the matrices , , and in (16).

Step 2. Solve the system by using the banded preconditioner .

Step 3. Compute .

Step 4. For , solve the systems by using the banded preconditioner .

Step 5. Obtain the numerical solution .

4. Numerical Results

In this section, we will test the feasibility and effectiveness of the banded preconditioning iteration method with the preconditioner for TS-FADE (1). All runs at each time step are started from zero initial guess and terminated if the current iteration satisfies or the elapsed CPU time is exceeding 2000 in seconds, whereAll the experiments are performed in MATLAB R2012a on Intel® Core™ i7-3770 CPU 3.40 GHz and 8.00 GB of RAM, with machine precision . The average iteration counts (denoted by “Iter”), elapsed CPU time in seconds (denoted by “CPU”), and total iteration counts for solving the discretized system (17) (denoted by “IT”) are reported. Here “Iter” is defined as