A GPIU method for fractional diffusion equations

*Correspondence: hailong_shen@126.com 1Department of Mathematics, College of Sciences, Northeastern University, Shenyang, P.R. China Abstract The fractional diffusion equations can be discretized by applying the implicit finite difference scheme and the unconditionally stable shifted Grünwald formula. Hence, the generating linear system has a real Toeplitz structure when the two diffusion coefficients are non-negative constants. Through a similarity transformation, the Toeplitz linear system can be converted to a generalized saddle point problem. We use the generalization of a parameterized inexact Uzawa (GPIU) method to solve such a kind of saddle point problem and give a new algorithm based on the GPIU method. Numerical results show the effectiveness and accuracy for the new algorithm.


Introduction
The fractional differential operator is suitable for describing the memory, genetic, mechanical and electrical properties of various materials. Compared with the classical integer differential operator, it can more concisely and accurately describe the biological, mechanical and physical processes with historical memory and spatial global correlation characteristics, such as abnormal diffusion of particles, quantization problem of non-local field theory, the fractional capacitance theory, universal voltage shunt, chaotic circuit analysis, physical semiconductors field, dispersion in porous media, physical and engineering issues related to fractal dimensions, and non-Newtonian fluid mechanics.
The fractional diffusion equations can be abstracted from many practical problems, such as a random walk model describing the competition between sub-diffusion and superdiffusion Three different types of fractional diffusion equations can be derived according to the difference in particle waiting time and jumping steps: when the mean of the waiting time of each step is infinite and the squared mean of the jump steps is limited, the random walk model describes the abnormal sub-diffusion phenomenon, and the time fractional diffusion equation is derived accordingly; when the mean of the waiting time is limited and the squared mean of the jump steps is infinite, the random walk model describes the super-diffusion phenomenon, and the spatial fractional-order diffusion equation is derived accordingly; when the mean of the waiting time and the squared mean of the jump steps are both infinite, the random walk model describes the phenomenon of competition between sub-diffusion and super-diffusion, and the fractional diffusion equation in time and space is derived accordingly; using the fractional convection-diffusion equations to simulate Levy's motion is very effective, it can more accurately simulate the solute motion process with long tail behavior; when studying Eulerian estimation of solute transport in porous media, the fractional-order Fick law can be used to replace the traditional Fick law to obtain the fractional-order convection-diffusion equations; in the application of non-Newtonian fluid mechanics, to solve the analytical solution of a material model problem between the Hook solid model and Newtonian fluid model, a special function related to fractional calculus needs to be applied; in the early days, the problem of dynamic response to broadband crust induced by earthquakes that was hardly affected by velocity can be described by a fractional linear rheological solid model. It is precisely because of the advantages of fractional derivatives and the wide application of fractional diffusion equations (FDEs) (1) in groundwater contaminant transport [1,2], turbulent flow [3], image processing [4], finance [5], physics [6] and other fields, that people have drawn extensive attention to its research.
Considerable numerical methods have been developed as the major ways for solving FDEs [7][8][9][10][11][12]. However, some salient features of fractional differential operators lead to the unconditional instability [13]. In addition, most numerical methods for finite difference equations (FDEs) can generate a full coefficient matrix, which requires the computational cost of O(N 3 ) and the storage of O(N 2 ), where N is the number of grid points [14]. In contrast, the second-order diffusion equations which usually contribute to sparse coefficient matrices can be solved efficiently by fast iteration methods with computational complexity O(N).
To guarantee the stability, Meerschaert and Tadjeran [13] proposed a shifted Grünwald discretization to approximate fractional diffusion equations and it has been proved that it is actually unconditionally stable. Then, based on the Meerschaert-Tadjeran method, a special structure of the full coefficient matrices was presented in [14]. The method maintained a Toeplitz-like structure and reduced the storage requirement from O(N 2 ) to O(N). Furthermore, the matrix-vector multiplication for the Toeplitz matrix can be calculated by the fast Fourier transform (FFT) with O(N log N) operations [15]. As a result, fast numerical methods for solving fractional equations with the shifted Grünwald formula have been developed, including the conjugate gradient normal residual (CGNR) methods [16,17], the preconditioned CGNR (PCGNR) methods which employ some circulant preconditioners [17][18][19] and the multigrid method [20]. If the diffusion coefficients are constant, the full coefficient matrix originating from the Meerschaert-Tadjeran method will be non-Hermitian positive definite. Therefore, the Hermitian and skew-Hermitian splitting (HSS) iteration and the corresponding preconditioned HSS (PHSS) iteration [21] have also been developed to solve fractional equations.
Since the full coefficient matrix is also a real Toeplitz matrix, the FDE can be transformed into a generalized saddle point problem based on some properties of the real Toeplitz matrix [22]. Therefore, the generalized parameterized inexact Uzawa (GPIU) method [23,24] is suitable for solving this new saddle point problem. However, it shows some differences from the traditional GPIU method. This method should store three N/2 order matrices. Although these matrices have symmetric structures (which will require less storage space), the complexity is less than O(N 2 ). In the meantime, the complexity of calculation is also O(N 2 ) since the method needs to compute the multiplication of several N/2 order matrices.
The framework of this paper is as follows. In Sect. 2, we introduce the background of the discretization for the FDEs. Then the real Toeplitz linear system is converted to a generalized saddle point problem by similarity transformation in Sect. 3. In Sect. 4, the GPIU method is considered to solve this kind of generalized saddle point problem and its convergence is analyzed. Besides, a GPIU algorithm is presented. In Sect. 5, numerical examples are performed to illustrate the effectiveness and power of our new method. Finally, conclusions are given in Sect. 6.

The depiction of FDEs and finite difference discretization
Firstly, an initial-boundary value problem of the FDEs is given by where α ∈ (1, 2) is the order of fractional derivatives, f (x, t) is the source term and diffusion coefficient functions d ± ≥ 0 with d + + d -= 0. The left-sided and the right-sided fractional derivatives ∂ α u(x,t) ∂ + x α and ∂ α u(x,t) ∂ -x α of (1) have the same definition as the Grünwald form [25] ∂ α u(x, t) where x denotes the floor of x and the Grünwald weights g (α) k are the alternating fractional binomial coefficient given as which can be evaluated by the following recursive form: Besides, the coefficients g (α) k satisfy the following properties [13,14,17,21] when 1 < α < 2: Let N and M be positive integers, then x = x R -x L N+1 and t = T M are the sizes of spatial grids and time steps, respectively. Additionally, the spatial and temporal partition are defined as x i = x L + i x for i = 0, 1, . . . , N + 1 and t m = m t for m = 0, 1, . . . , M. We also let u (m) . Consider the shifted Grünwald approximations [13,17,21]: where g (α) k is defined in (2) and the corresponding implicit finite difference scheme in (1) can be modified as follows: . . , f (m) N ] T and I ∈ R N×N be an identity matrix. Then the scheme (3) can construct the following equations: and where It is obvious that G α is a Toeplitz matrix [15] and the coefficient matrix A (m) is a full Toeplitz-like matrix, which is non-symmetric when d + = dand the matrix-vector multiplication for the A (m) can be obtained in O(N log N) operations by the fast Fourier transform (FFT). We N+1) α , related to the numbers of time steps and grid points. The linear system (4) is represented as follows: where As the coefficient matrix M (m) in (8)  is a symmetric positive definite matrix. Therefore, M (m) is also positive definite (see [21,26]).

The transformation for the Toeplitz matrix M (m)
As mentioned in Sect. 2, the coefficient matrix M (m) is a real positive definite Toeplitz-like matrix. We assume that the diffusion coefficient functions are two non-negative constants, i.e., d (m) Then M (m) becomes a real non-symmetric positive definite Toeplitz matrix. In addition, M is appropriately chosen according to N and satisfies the limit of v N,M away from 0.
Then we transform the Toeplitz linear system (7) into a generalized saddle point problem by constructing an orthogonal matrix to deal with M (m) . The first row and column of the real N × N Toeplitz matrix M (m) can be indicated as respectively. We use the Hermitian and skew-Hermitian splitting (HSS) method to split M (m) and get where Considering that N is even or odd, we discuss the structures of H (m) and S (m) in detail.
Case I: N is even, let N = 2m. Then H (m) can be rewritten as where E ∈ R m×m is a symmetric Toeplitz matrix whose first column can be given as and F ∈ R m×m is a Toeplitz matrix whose first row and column are as follows: respectively.
Analogously, S (m) can be rewritten as where G ∈ R m×m is a skew-symmetric Toeplitz matrix whose first column can be given as and L ∈ R m×m is a Toeplitz matrix whose first row and column are as follows: respectively. We give an orthogonal matrix P defined as then multiply the left and right sides of (9) by P T from the left and P from the right, respectively, we have the new forms ofH (m) andS (m) : Hence, letM (m) = P T M (m) P,ũ (m) = P T u (m) , andb (m-1) = P T b (m-1) , we can rewrite Eq. (7) as the following linear system: and we havẽ For simplicity, let B = E + J m F, C = E -J m F and W = -G + J m L, then we can find that The coefficient matrixM (m) can be expressed as where E, F and J m are similar to the above, and the column vector z 1 is expressed as follows: In addition, S (m) can be expressed as where G and L are similar to the above, and the column vector z 2 is expressed as follows: Then the orthogonal matrix is defined as Multiplying the left and right sides of (9) by P T from the left and P from the right, respectively, we have the new forms ofH (m) andS (m) as The Toeplitz linear system (7) can also be constructed as a generalized saddle point problem which has a form similar to (10), where All in all, no matter what positive integer N we choose, the scheme for the FDEs is able to be transformed into a generalized saddle point problem.

A GPIU method on the transformed saddle point problem
By splitting the new coefficient matrixM (m) in the linear system (10), we get where here the definitions of B, C and W are given in the above section, Q 1 ∈ R N/2 is a symmetric positive definite matrix, Q 2 ∈ R N/2 is a symmetric positive definite matrix, and Q 3 ∈ R N/2 × N/2 is arbitrary, where N/2 and N/2 represents the ceil and the floor of N/2, respectively.
Considering the GPIU method [23,24], we choose Q 2 = C/δ where δ > 0, so Y (m) can be rewritten as Thus, the linear system (10) is redefined by the splitting form of (12) as The form of the iteration matrix corresponding to (13) is as follows: Denote by ρ(W) the spectral radius of the iteration matrix W, so Eq. (13) converges if and only if ρ(W) < 1.
Let λ be an eigenvalue of W and [ũ H x ,ũ H y ] H be the corresponding eigenvector withũ x ∈ C N/2 andũ y ∈ C N/2 , then we obtain Or, equivalently, Aiming to obtain the convergence condition, the lemmas and theorem as follows are presented. We can rewrite the system as Evidently, the coefficient matrix is non-singular, so [ũ H x ,ũ H y ] H = 0, which is contradictory, so λ = 1.
(b) Ifũ x = 0, then according to (15) we have When W is full of row rank, according to the first equation of (16), we haveũ y = 0. This is a contradiction.
(c) When W is not full of row rank, letũ x = 0.
(i) If λ = 1δ, according to the second equation of (16) and Q 2 is symmetric positive definite, we obtainũ y = 0, which is a contradiction.
(ii) If λ = 1δ, the fundamental set of solution can be obtained by the first equation of (16) onũ y which consists of N/2r(W ) non-zero vector solutions, i.e., λ = 1δ is at least N/2r(W ) eigenvalues, where r(W ) represents the rank of matrix W and r(W ) = r(W T ). Proof (a) Ifũ y = 0, from (15), we obtain When W is full of row rank, since λ = 1 andũ x = 0, from Lemma 4.1, we denote Through multiplying the both sides of the first equation in (17) (b) When W is not full of row rank, if λ = 1δ, thenũ x = 0 from Lemma 4.1. As λ = 1, in the same way as with (a), we get 0 ≤ λ < 1.
In order to ensure the convergence of the iteration method (13), the following theorems are used to give the necessary and sufficient conditions.
x ,ũ H y ] H be a corresponding eigenvector whereũ x ∈ C N/2 andũ y ∈ C N/2 .
Proof Let λ be an eigenvalue of the iteration matrix W and [ũ H x ,ũ H y ] H be its corresponding eigenvector withũ x ∈ C N/2 andũ y ∈ C N/2 . By Lemma 4.1, we know λ = 1.
If λ = 1δ, then (15) results in represents the null space of the corresponding matrix. Then λ = 1δ is an eigenvalue of the iteration matrix W.
When W is not full of row rank,ũ x ∈ null[(1δ)W + δQ 3 ] and there exist N/2r(W ) non-zero eigenvectors ofũ y corresponding to λ = 1δ which are at least N/2r(W ) eigenvalues of the iteration matrix W. Particularly, whenũ x is the zero vector, it has been introduced in Lemma 4.1.
If λ = 1δ, we haveũ x = 0 from the Lemma 4.1, and (15) can be rewritten as Eliminatingũ y in (19), we obtain Then multiplying the left and right sides of (20) byũ where Then (21) can be modified as When γ = 0, we can get Wũ x = 0 and multiply the left and right sides of the first equation in (19) byũ x from the left, we have When Wũ x = 0, we obtain γ > 0. In accordance with Lemma 4.3, the roots of the quadratic Eq. (22) satisfy |λ| < 1 if and only if Solving the inequalities of (23), we have In addition, we consider the case (b) from Lemma 4.2 whenũ y = 0, it actually satisfies the inequalities (24). Therefore, we can get (24) where we have considered that α ≥ 0 and β > 0. By choosing different matrices Q 1 , Q 2 and Q 3 in (12), we get considerable iteration algorithms from the GPIU method for solving the generalized saddle point problem (11).
Then an efficient GPIU algorithm can be presented.

The GPIU Algorithm
the GPIU method gives the following iteration scheme: where the matrices B, C and W rely on N defined in Sect. 4.

Numerical examples
The FDEs problem (1) has been solved by pre-multiplying the coefficient matrix with P T and post-multiplying it with P which is given in Sect. 3, firstly. Then the FDEs problem is changed into a generalized saddle point problem, so in Sect. 4 an algorithm is proposed based on the GPIU method. In this section, we discuss the GPIU algorithm with different parameters and comparison between the GPIU algorithm and the CGNR method (3). The initial value is defined as the zero vector at each time step, and each iteration process terminates when the current residual satisfies r k 2 / r 0 2 < 10 -7 , where r k is the residual vector of the linear system after k iterations. Then we collect data in the table as follows and define "N " as the number of spatial grid points, "M" as the number of time steps, "Error" as the difference between the exact solution and the approximate solution under the infinity norm, "CPU(s)" as the total CPU(s) time (seconds) to solve the whole FDEs problem, and "Iter" as the average number of iterations needed to solve the FDEs problem, i.e.,

Iter
For satisfying v N,M is bounded away from 0, let t ≈ 2 x α . From Table 1, we show the advantage of efficiency according to the CPU time, Iter and Error. When N and M become large, the CPU time and Error will be slow and show inaccuracy, respectively. But the average number of iterations remains stabilized. Table 2 demonstrates that not only the Iter but also CPU(s) by the CGNR method with optimal parameters are much higher than the GPIU method. The iteration steps and CPU  time of the new method is obviously less than HSS. It also indicates that the Error is obviously improved for the FDEs over CGNR and HSS. Figure 1 shows the convergence rate of GPIU and HSS. In order to get the least number of steps which can makes the error decrease to the ε times initial error, we normally define ln(ρ(M)) as the asymptotic convergence rate, where M is the iteration matrix of the iteration method for solving the equations. The minimum number of steps of the iteration is k ≥ -ln ε ln(ρ(M)) . In Fig. 1, the X-axis stands for the scale of the matrix, and the Y-axis stands for the theoretical minimum steps for iteration. It is shown that, for GPIU, the minimum number of steps increases slowly and almost tends to a stable value as N increases. But for HSS, the minimum number of steps increases sharply.    Table 3 illustrates the numerical results of Example 2. Here we specially select the source term f (x, t) as non-zero and nonlinear. Obviously, the Iter for the CGNR method is larger and the Iter of the new method is smaller than HSS. Moreover, the CPU time and the error of the GPIU method are much smaller than the CGNR method and HSS method.
It indicates that the GPIU algorithm is still highly efficient and accurate for the case of nonlinear source term, so the GPIU algorithm can be extendable to the case of nonlinear source term.
Note. The memory usage for all tests in Tables 1, 2 and 3 is O(n), where n denotes the order of the coefficient matrix of equations [7].

Concluding remarks
In this paper, a GPIU method is put forward to solve the generalized saddle point problem generated by the fractional diffusion equations with constant coefficients by premultiplying the orthogonal matrix we constructed and post-multiplying the transposition of it, respectively. Then the convergence condition of the GPIU method is analyzed and we give a new GPIU algorithm.
The numerical results show that the proposed method is more effective than the CGNR method and HSS, and its advantages become more obvious as N increases. Compared with the CGNR method and the HSS method, this GPIU method can achieve a faster convergence rate in practice and theory, and can obtain a more accurate convergence solution in a shorter time. Unlike the CGNR method and the HSS method, with the increase of the matrix order N, the number of iteration steps of the GPIU method grows slowly, it is almost stable near an extremely low value, and the error is significantly improved. The better convergence results of GPIU are due to the fact that the storage of the new algorithm is greatly reduced to O(n), to be compared with the CGNR and HSS method, where n denotes the order of the coefficient matrix of equations.
For further consideration, more effective stationary methods like the GPIU method will be promoted to solve the FDE problems with constant or even variable diffusion coefficients.