Fast permutation preconditioning for fractional diffusion equations.

In this paper, an implicit finite difference scheme with the shifted Grünwald formula, which is unconditionally stable, is used to discretize the fractional diffusion equations with constant diffusion coefficients. The coefficient matrix possesses the Toeplitz structure and the fast Toeplitz matrix-vector product can be utilized to reduce the computational complexity from [Formula: see text] to [Formula: see text], where N is the number of grid points. Two preconditioned iterative methods, named bi-conjugate gradient method for Toeplitz matrix and bi-conjugate residual method for Toeplitz matrix, are proposed to solve the relevant discretized systems. Finally, numerical experiments are reported to show the effectiveness of our preconditioners.

Background

In the last few decades, many anomalous diffusion phenomena have been found in the real world, which lead to the generation of the fractional diffusion equations (FDEs). The FDEs emerge from numerous research fields such as modeling chaotic dynamics of classical conservative systems (Zaslavsky et al. 1993), groundwater contaminant transport (Benson et al. 2000a, b), turbulent flow (Carreras et al. 2001; Shlesinger et al. 1987), and applications in physics (Sokolov et al. 2002), finance (Raberto et al. 2002), biology (Magin 2006), hydrology (Baeumer et al. 2001) and image processing (Blackledge 2009; Bai and Feng 2007). Usually, it is unavailable to access the closed-form analytical solutions of the FDEs (Alquran et al. 2015; Allan and Al-Khaled 2006; Sababheh et al. 2003). Therefore, many numerical approaches for the FDEs have been proposed and developed intensively in the last decade, for instance Zhang et al. (2010), Ervin et al. (2007), Langlands and Henry (2005), Liu et al. (2004), Meerschaert and Tadjeran (2004, 2006), Tian et al. (2015), Gu et al. (2015). However, even if the discretized approach of the FDEs is implicit, it still can result in unconditionally unstable (Meerschaert and Tadjeran 2004, 2006) because of the nonlocality of the fractional differential operators.

In order to overcome the difficulty of the stability, (Meerschaert and Tadjeran 2004, 2006) put forward a shifted Grünwald discretization to approximate FDEs with a left-sided fractional derivative and the FDEs with two-sided fractional derivatives, respectively, and their method has been proved to be unconditionally stable. However, it is worth noting that most of the numerical methods for FDEs tend to generate full coefficient matrices, which require computational cost of and storage of , there is no doubt that it will certainly increase the computational work; see, e.g., Wang et al. (2010) for a discussion of these issues.

Recently, there is some progress on fast numerical solutions of FDEs. Wang et al. (2010) discovered that the full coefficient matrix generated by Meerschaert and Tadjeran (2006) method has a good feature, i.e., it can be written as a sum of diagonal-multiply-Toeplitz matrices. Thus the storage requirement is reduced from to . As we know, the Toeplitz matrix-vector product (MVP) can be computed in operations by the fast Fourier transforms (FFTs) (Chan and Ng 1996). Fast methods Lei and Sun (2013), Popolizio (2015), Gu et al. (2015, 2015) have been developed to solve FDEs with the shifted Grünwald formula. Wang and Wang (2011) proposed a conjugate gradient normal residual (CGNR) to solve the discretized system by Meerschact and Tadjeran’s method with the computational cost of . The preconditioned CGNR with a circulant preconditioner is proposed by Lei and Sun (2013), to solve FDEs by Meerschact and Tadjeran’s method with constant diffusion coefficients.

In this paper, we employ the implicit finite difference method to discretize the FDEs and the problem is transformed to solve a linear nonsymmetric Toeplitz system in each time step. Since, the Bi-Conjugate Gradient (BiCG) (Saad 2003, pp. 234–236) and Bi-Conjugate Residual (BiCR) (Sogabe et al. 2009) can be regarded as two effective methods for solving the nonsymmetric system. There is no doubt that the two methods with Toeplitz fast MVP can be used to solve such discretized Toeplitz linear systems. However, from Sogabe et al. (2005), Pestana and Wathen (2015), if the two iterative methods are employed directly, then we indeed fail to make full use of Toeplitz structure of the discretized system, it also means that their computational cost fail to attain optimality. Hence, it is still worth finding more effective methods to reduce the computational complexity. Recently, in Sogabe et al. (2005), Pestana and Wathen (2015), a permutation matrix P was introduced to transform the nonsymmetric matrix into a symmetric one so as to improve the performance of iterative methods. In view of this point, we re-explain the ideas in Sogabe et al. (2005), Pestana and Wathen (2015) as a kind of preconditioning techniques for solving the discretized system of the FDEs by the method of Meerschaet and Tadjeran. More precisely, we do equivalent transformation for the original discretized system, left multiplying by a permutation matrix (Pestana and Wathen 2015) at the same time, then we obtain a new symmetric linear system with the coefficient matrix being a Hankle matrix, which has the same solution with the original discretized system. As we know, the symmetric linear systems are usually simpler to be solved than the nonsymmetric cases. Conjugate Gradient (CG) and Conjugate Residual (CR) are two effective methods for solving symmetric linear system. In this paper, we extend CG and CR to BiCGT and BiCRT, respectively, which are proposed to solve the equivalent equation. The numerical results show that both BiCGT and BiCRT are more competitive than CGNR.

The paper is organized as follows. in Sect. 2, we briefly introduce the discretization of FDEs by finite difference method. In Sect. 3, we construct the permutation preconditioner and propose BiCGT and BiCRT to solve the equivalent system of linear equations. In Sect. 4, numerical results are reported to illustrate the efficiency of the proposed methods. Concluding remarks are given in Sect. 5.

Discretization of FDEs by finite difference method

In this section, we are interested in solving an initial-boundary value problem of the following FDEs,where is the order of the space fractional derivative, f(x, t) is the source term, and the diffusion coefficients satisfying , , and . In this paper, we use the Grünwald-Letnikov form Podlubny (1999) to define the left-sided and the right-sided fractional derivatives and :where denotes the floor function, and the Grünwald coefficients are defined as followswhich can be evaluated by the recurrence relationIn order to derive the proposed scheme, let and be the sizes of spatial grid and time step, respectively (N, M are positive integers). Define () and the temporal partition (). Let , . We employ the shifted Grünwald approximation (Meerschaert and Tadjeran 2004, 2006):where is defined in Eq. (2). Then the corresponding finite difference schemeis unconditionally stable, see Meerschaert and Tadjeran (2004, 2006) for details.

Let Then we can rewrite (3) into the matrix formwithwhereWe can note that is a nonsymmetric Toeplitz matrix, thus it can be stored with entries (Wang et al. 2010). The Toeplitz matrix-vector product (MVP) can be computed in complexities with the aid of FFTs (Pang and Sun 2012).

Define which is related to the number of time steps and grid points. The above linear system (4) can be rewritten in the following matrix formwhereand the right hand vectorIn order to illustrate the convergence and stability of the implicit difference scheme (3), we note that satisfy the following proposition.

Proposition 1

(Meerschaert and Tadjeran 2004; 2006; Wang et al. 2010) Letandbe defined in (2). Then we have

Since is a strongly diagonally dominant and nonsingular Toeplitz matrix, and thus the scheme (3) is monotone; refer to Wang et al. (2010).

The BiCGT method and the BiCRT method

We will show how to construct the permutation preconditioners for accelerating the iterative solver and describe the derivation process of BiCGT and BiCRT. Furthermore, an analysis of computational cost for each iteration step is also proposed. In the linear system (5), is a nonsymmetric Toeplitz matrix. As previously mentioned, we cannot exploit BiCG and BiCR for resulting linear systems (5) directly, otherwise it is not possible to take advantage of the Toeplitz structure of coefficient matrix. So we need to modify and improve the BiCG and BiCR methods particularly for solving (5). Recently, Sogabe et al. (2005), proposed a preconditioner of permutation matrix for improving the performance of the Krylov subspace method, which is used to solve a nonsymmetric Toeplitz linear system. Later, Pestana and Wathen rigorously establish a circulant preconditioned MINRES method (Paige and Saunders 1975) for nonsymmetric Toeplitz systems. Inspired by their pioneer work, we construct a preconditioner, which is a permutation matrix P (Sogabe et al. 2005) with the form of

We would like to solve the Toeplitz system (5) by the CG-like method. This goal can be achieved with little additional computing cost, since we can get the equivalent system:which can be regarded as a left preconditioning technique (Saad 2003) and also has the same solution with (5). Define then (6) can be rewritten intoan equivalent statement is that is self-adjoint with respect to the bilinear form defined by P (Paige and Saunders 1975). P is symmetric positive definite, a nonsymmetric Toeplitz matrix is exactly changed into a symmetric matrix , so that (7) can be solved by the modified BiCG, where the additional operations for dual systems have been eliminated.

Firstly, we consider CG and BiCG (Saad 2003). Then we followed the philosophy behind the derivation of iterative method in Sogabe et al. (2005 Algorithm1), i.e., in the CG Algorithm, we replace A and b with and , respectively. Then we get the following new algorithm:

In BiCGT, we only need one MVP, i.e. and two inner products, three vector additions/subtractions per iteration. The rewritten algorithm is more effective than CGNR, because P multiply an arbitrary vector is to reorder the vector in its reversed order (Sogabe et al. 2005). Therefore, it can greatly reduce the required number of MVPs.

In a similar way, we can get the algorithm of BiCRT for symmetric linear system (7), BiCRT reduces to CR if we get rid of the permutation preconditioner P, and the algorithm of BiCRT is presented as blew.

If we employ BiCG or BiCR directly, it is impossible to minimize the residual vector in some special conditions. However, BiCRT could realize this goal to some extent. The approximation is generated from by moving from in a certain direction to a minimum point of the residual function , . In other words, for , is chosen to minimize E(u).

It is useful to consider the computational cost. We give a table to illustrate the computational cost of BiCG, BiCR, BiCGT, BiCRT and CGNR. “AXPY” denotes addition of scaled vectors, and “” denotes 1 product with the matrix and 1 with its transpose. From Table 1, it is remarkable that BiCG and BiCR, BiCGT and BiCRT require almost the same memory and computational cost in each iteration step. More precisely, BiCGT is the best method to solve the above system in terms of computational cost (i.e, AXPYs and MVPs). For BiCRT, the number of AXPYs is one more than BiCGT, and for CGNR, the number of MVPs is one more than BiCGT. As we know, the computational complexity of one Toeplitz MVP is by FFTs, but one AXPY can be computed in complexity. So BiCRT is more efficient than CGNR from this perspective.

Table 1

Summary of algorithmic cost per iteration step

Method	Dot product	AXPY	MVP
BiCG	2	5	1+1
BiCR	2	6	1+1
BiCGT	2	3	1
BiCRT	2	4	1
CGNR	2	3	2

Numerical results

We solve the FDEs (1) numerically by the implicit finite difference method given in Sect. 2. After the finite difference discretization and the equivalent transformation, the symmetric linear system (7) is solved by BiCGT (Algorithm 1), BiCRT (Algorithm 2), and CGNR (Wang and Wang 2011), respectively. All the MVPs are done by FFTs in operations (Lei and Sun 2013) and the initial guess is chosen to be zero vector at each time step. The stopping criterion of BiCGT, BiCRT and CGNR is set to bewhere is the residual vector of the linear system after k iterations and is the initial residual vector.

In the following tables, “N” denotes the number of spatial grid points, “M” denotes the number of time steps, CPU(s) denotes the total CPU time (in seconds) for solving the whole discretized system. “Error” denotes the infinity norm of the difference between the true solution and the approximation at the last time step. “Iter” denotes the average number of iterations over all time discretized level for solving the FDEs, i.e., Iter Iter (m), where Iter(m) is the number of iterations required for solving the linear system (7) in the mth time discretized level. All experiments are run in MATLAB R2010a on a PC with the following configuration: Windows 7 (32 bit), Iter(R) Core(TM) i3-2130 CPU 3.40 GHz and 4 GB RAM.

Example 1

We consider FDEs (1) on space interval and time interval with diffusion coefficients initial condition and source termThe exact solution is For the finite difference discretization, the space step and time step are taken to be , , i.e., .

The numerical results are listed in Table 2, as for comparisons, we also carry out CGNR without preconditioner. From Table 2, it is remarkable that the error is improved for CGNR, BiCGT and BiCRT as the increasing of . However, BiCGT and BiCRT are more effective than CGNR in terms of CPU time. More precisely, the performance of BiCGT is the best in terms of the CPU time except the cases of , , , and . In addition, the average number of iterations of BiCGT is less than that by CGNR and BiCRT sometimes. For instance, look at these cases in the numerical results at discretized size , , and for . BiCGT and BiCRT have faster convergence speed, less computational time expenditure than CGNR. Meanwhile, we can see that the CPU time increases as the order of the time derivative increases. To explain this phenomenon, we list the spectra of the original matrix with different in Fig. 1. As we can see from the figure, most of eigenvalues is close to zero with the increasing of the valve of , it means that the coefficient matrix become increasingly ill-conditioned, it also implies that the linear systems will be difficult to solve.

Table 2

Comparisons for solving Example 1 by different methods with and 1.8 at

α	N+1	CGNR			BiCGT			BiCRT
		CPU(s)	Error	Iter	CPU(s)	Error	Iter	CPU(s)	Error	Iter
1.4	64	0.129	3.7873e−4	69.6	0.101	3.7873e–4	63.6	0.103	3.7873e–4	63.7
	128	1.038	1.9163e–4	163.1	0.632	1.9163e−4	144.0	0.646	1.9163e−4	143.7
	256	4.828	9.6389e−5	272.0	3.963	9.6389e−5	354.0	3.898	9.6390e−5	338.7
	512	28.391	4.8338e−5	375.8	25.382	4.8338e−5	590.5	23.993	4.8342e−5	546.4
	1024	136.318	2.4205e−5	476.0	129.990	2.4206e−5	838.5	114.42	2.4212e−5	726.3
1.5	64	0.140	2.7756e−4	75.7	0.107	2.7756e−4	65.4	0.150	2.7756e−4	65.4
	128	1.206	1.4046e−4	188.2	0.760	1.4046e−4	162.6	0.844	1.4046e−4	162.2
	256	7.339	7.0658e−5	410.5	4.888	7.0654e−5	421.0	5.662	7.0654e−5	418.9
	512	51.154	3.5436e−5	681.4	42.872	3.5435e−5	955.4	45.933	3.5439e−5	904.0
	1024	289.79	1.7747e−5	1015.0	267.921	1.7744e−5	1657.2	277.638	1.7751e−5	1464.7
1.8	64	0.186	8.0708e−5	103.8	0.131	8.0708e−5	85.2	0.140	8.0708e−5	85.0
	128	1.852	4.0979e−5	292.8	1.042	4.0979e−5	234.3	1.100	4.0978e−5	234.0
	256	14.583	2.0659e−5	841.5	8.570	2.0659e−5	746.9	8.629	2.0659e−5	729.4
	512	182.390	1.0375e−5	2472.5	96.153	1.0374e−5	2205.7	96.526	1.0374e−5	2164.3
	1024	1708.245	5.2057e−6	6053.9	1075.756	5.1984e−6	6778.4	1058.220	5.1985e−6	6604.0

Fig. 1

The spectra of the matrix with for Example 1

Example 2

Consider FDEs (1) on space interval and time interval with diffusion coefficients initial condition source termThe exact solution of this example is For the finite difference discretization, the space step and time step are taken to be , and , respectively.

Table 3 shows the numerical results for solving Example 2 by different methods. The error is decreased for those methods as the increasing of , and the accuracy is almost the same as Example 1. Similar to Example 1, BiCGT and BiCRT are more effective than CGNR in terms of CPU time elapsed. Besides, the CPU time increases as the order of the time derivative increases is similar to Example 1, and the reason is the same as Example 1. We also list the spectra of the matrix with different in Fig. 2.

Table 3

Comparisons for solving Example 2 by different methods with and 1.8 at

α	N+1	CGNR			BiCGT			BiCRT
		CPU(s)	Error	Iter	CPU(s)	Error	Iter	CPU(s)	Error	Iter
1.4	64	0.131	4.0801e−4	70	0.101	4.0801e−4	64.1	0.103	4.0801e−4	64.2
	128	1.007	2.0542e−4	157.3	0.671	2.0542e−4	151.9	0.682	2.0542e−4	151.7
	256	4.764	1.0283e−4	274.0	3.985	1.0283e−4	353.7	3.977	1.0283e−4	349.0
	512	30.947	5.1397e−5	408	27.848	5.1395e−5	643.2	26.268	5.1398e−5	601.4
	1024	148.506	2.5696e−5	519	147.36	2.5677e−5	948.6	130.875	2.5682e−5	837.8
1.5	64	0.147	2.7980e−4	76.3	0.100	2.7980e−4	66.9	0.111	2.7980e−4	66.6
	128	1.286	1.4222e−4	190.0	0.745	1.4222e−4	171.5	0.791	1.4222e−4	171.0
	256	7.410	7.1577e−5	397.5	4.806	7.1582e−5	427.2	4.962	7.1582e−5	426.0
	512	56.933	3.5883e−5	713.8	43.722	3.5877e−5	1012.6	42.825	3.5883e−5	972.4
	1024	331.186	1.7968e−5	1106	286.337	1.7953e−5	1841.4	261.763	1.7948e−5	1659.9
1.8	64	0.182	9.8722e−5	105.5	0.144	9.8722e−5	91.3	0.237	9.8722e−5	90.2
	128	1.890	4.6159e−5	298.2	1.088	4.6159e−5	241.3	1.592	4.6159e−5	240.1
	256	15.027	2.2331e−5	858.4	8.798	2.2331e−5	758.0	10.687	2.2331e−5	746.2
	512	189.009	1.0989e−5	2503.8	105.062	1.0989e−5	2332.5	133.211	2.2331e−5	2283.1
	1024	1733.04	5.4529e−6	6104.3	1158.655	5.4529e−6	7074.2	1007.543	5.4529e−6	6900.6

Fig. 2

The spectra of the matrix with for Example 2

Concluding remarks

Two new iterative methods, named BiCGT and BiCRT, are presented to solve the resulting linear system of the FDEs (1), which are discretized by the implicit finite difference method. Namely, with the help of the permutation matrix P, we transform the difficult nonsymmetric linear systems into the symmetric cases, which are often simpler to be solved. The computational complexity can be reduced from to by utilizing FFTs. Numerical experiments illustrate the effectiveness of the proposed methods. In our future work, we will apply BiCGT and BiCRT to solve other (two dimensional) fractional differential equations (Wang and Basu 2012), such as fractional advection-diffusion equations; and we will investigate and develop some suitable preconditioning, see e.g. Lei and Sun (2013), Gu et al. (2015) to further accelerate our proposed methods.