- Email: [email protected]
A circulant preconditioner for fractional diffusion equations
Journal of Computational Physics 242 (2013) 715–725
Contents lists available at SciVerse ScienceDirect
Journal of Computational Physics journal homepage: www.elsevier.com/locate/jcp
A circulant preconditioner for fractional diffusion equations q Siu-Long Lei, Hai-Wei Sun ⇑ Department of Mathematics, University of Macau, Macao, China
a r t i c l e
i n f o
Article history: Received 10 September 2012 Received in revised form 4 February 2013 Accepted 12 February 2013 Available online 26 February 2013 Keywords: Fractional diffusion equations Shifted Grünwald discretization Toeplitz Circulant preconditioner Fast Fourier transform CGNR method
a b s t r a c t The implicit finite difference scheme with the shifted Grünwald formula, which is unconditionally stable, is employed to discretize fractional diffusion equations. The resulting systems are Toeplitz-like and then the fast Fourier transform can be used to reduce the computational cost of the matrix–vector multiplication. The preconditioned conjugate gradient normal residual method with a circulant preconditioner is proposed to solve the discretized linear systems. The spectrum of the preconditioned matrix is proven to be clustered around 1 if diffusion coefficients are constant; hence the convergence rate of the proposed iterative algorithm is superlinear. Numerical experiments are carried out to demonstrate that our circulant preconditioner works very well, even though for cases of variable diffusion coefficients. Ó 2013 Elsevier Inc. All rights reserved.
1. Introduction In the last decade or so, the catholicity of anomalous diffusion phenomena in the real world has led to the fractional diffusion equation (FDE). FDEs arise in research topics including modeling chaotic dynamics of classical conservative systems [33], groundwater contaminant transport [3,4], turbulent flow [6,25], and applications in biology [16], finance [24], image processing [1], and physics [26]. As there are very few cases of FDEs in which the closed-form analytical solutions are available, numerical solutions for FDEs become major ways and then have been developed intensively [5,11,12,14,15,17–20,27–29]. Nevertheless, since the fractional differential operator is nonlocal, it was shown that a naive discretization of the FDE, even though implicit, leads to unconditionally unstable [18,19]. Moreover, most numerical methods for FDEs tend to generate full coefficient matrices, which require computational cost of OðN 3 Þ and storage of OðN 2 Þ, where N is the number of grid points [31]. It is quite different from second-order diffusion equations which usually yield sparse coefficient matrices with OðNÞ nonzero entries and can be solved very efficiently by fast iterative methods with OðNÞ complexity. To overcome the difficulty of the stability, Meerschaet and Tadjeran [18,19] proposed a shifted Grünwald discretization to approximate FDEs. Their method has been proven to be unconditionally stable. Later, Wang, Wang, and Sircar [31] discovered that the full coefficient matrix by the Meerschaet–Tadjeran’s method holds a Toeplitz-like structure. More precisely, such a full matrix can be written as the sum of diagonal-multiply-Toeplitz matrices. Thus the storage requirement is significantly reduced from OðN 2 Þ to OðNÞ. It is well known that the matrix–vector multiplication for the Toeplitz matrix can be computed by the fast Fourier transform (FFT) with OðN log NÞ operations [8,21]. With this advantage, Wang and Wang [32] employed the conjugate gradient normal residual (CGNR) method to solve the discretized system of the FDE by the Meerschaet–Tadjeran’s method. Thanks to the Toeplitz-like structure, the cost per iteration by the CGNR method is of
q This work was partially supported by the research grant MYRG038(Y2-L1)-FST12-LSL and MYRG140(Y3-L2)-FST11-SHW from University of Macau, and the research grant 105/2012/A3 from FDCT of Macao. ⇑ Corresponding author. E-mail addresses: [email protected] (S.-L. Lei), [email protected] (H.-W. Sun).
0021-9991/$ - see front matter Ó 2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.jcp.2013.02.025
716
S.-L. Lei, H.-W. Sun / Journal of Computational Physics 242 (2013) 715–725
OðN log NÞ. The convergence of the CGNR method is fast with smaller diffusion coefficients [32] (in that case the discretized system is well-conditioned). Nevertheless, if the diffusion coefficient functions are not small, the resulting system will become ill-conditioned and hence the CGNR method converges very slowly. To overtake this shortcoming, Pang and Sun [22] proposed a multigrid method to solve the discretized system of the FDE by the Meerschaet–Tadjeran’s method. With the damped-Jacobi method as the smoother, the multigrid algorithm can preserve the computational cost per iteration as OðN log NÞ operations. Numerical results showed that their multigrid method converges very fast, even for the ill-conditioned systems. However, from the theoretical point of view, the linear convergence of their multigrid method, despite a very simple case (both diffusion coefficients are equal and constant), has not been proven; see [22] for details. In the literature on Toeplitz or Toeplitz-like systems, circulant preconditioners always played important roles. In fact, circulant preconditioners have been theoretically and numerically studied with numerous applications for over twenty years; see for instances [7,8,21]. Nevertheless, to our knowledge, circulant preconditioners were never particularly applied to the Toeplitz-like matrix with structure as the sum of diagonal-multiply-Toeplitz matrices before. In this paper, we propose the preconditioned CGNR (PCGNR) method with a circulant preconditioner to solve such Toeplitz-like systems. The proposed preconditioner is invertible and its spectrum is theoretically proven to be clustered around 1 under some conditions. Thus the superlinear convergence rate of the PCGNR method is obtained. As the computational cost per iteration is of OðN log NÞ, the total complexity of the PCGNR method at each time step retains OðN log NÞ operations. The paper is organized as follows. In Section 2, the background of the discretization for the FDE is reviewed. In Section 3, the PCGNR method with a circulant preconditioner is proposed to solve the discretized linear system. The convergence rate of the PCGNR method is studied in Section 4. In Section 5, numerical results are reported to demonstrate the efficiency of the proposed method. Concluding remarks are given in Section 6. 2. FDE and finite difference discretization In this paper, we study an initial-boundary value problem of the FDE:
8 @uðx; tÞ @ a uðx; tÞ @ a uðx; tÞ > > > ¼ dþ ðx; tÞ þ d ðx; tÞ þ f ðx; tÞ; > a > @t @þx @ xa > > > < t 2 ð0; T; x 2 ðxL ; xR Þ; > > > > uðxL ; tÞ ¼ uðxR ; tÞ ¼ 0; 0 6 t 6 T; > > > > : x 2 ½xL ; xR ; uðx; 0Þ ¼ u0 ðxÞ;
ð1Þ
where a 2 ð1; 2Þ is the order of the fractional derivative, f ðx; tÞ is the source term, and diffusion coefficient functions d ðx; tÞ a
a
@ uðx;tÞ are nonnegative; i.e., d ðx; tÞ P 0. The left-sided and the right-sided fractional derivatives @ @uðx;tÞ a and @ xa are defined in the þx
Grünwald–Letnikov form [23]
@ a uðx; tÞ 1 ¼ limþ a @ þ xa Dx!0 Dx
bðxx L Þ=Dxc X
@ a uðx; tÞ 1 ¼ limþ a @ xa Dx!0 Dx
bðxRX xÞ=Dxc
ðaÞ
g k uðx kDx; tÞ;
k¼0
ðaÞ
g k uðx þ kDx; tÞ;
k¼0 ðaÞ
where bc denotes the floor function, and g k is the alternating fractional binomial coefficient given as
(
ðaÞ
g 0 ¼ 1; ðaÞ
k
aða 1Þ ða k þ 1Þ; g k ¼ ð1Þ k!
k ¼ 1; 2; 3; . . . ;
ð2Þ
which can be evaluated by the recurrence relation ðaÞ
g kþ1 ¼
a þ 1 ðaÞ g ; 1 kþ1 k
k ¼ 0; 1; 2; . . . :
R xL Let N and M be positive integers, and Dx ¼ xNþ1 and Dt ¼ T=M be the sizes of spatial grid and time step, respectively. We define a spatial and temporal partition xi ¼ xL þ iDx for i ¼ 0; 1; . . . ; N þ 1 and tm ¼ mDt for m ¼ 0; 1; . . . ; M. Let
ðmÞ
ui
ðmÞ
ðmÞ
¼ uðxi ; t m Þ; d;i ¼ d ðxi ; tm Þ, and fi
¼ f ðxi ; t m Þ. The shifted Grünwald approximation in [18,19] is as follows,
iþ1 @ a uðxi ; tm Þ 1 X ðaÞ ðmÞ ¼ a g u þ OðDxÞ; a @þx Dx k¼0 k ikþ1
X ðaÞ ðmÞ @ a uðxi ; tm Þ 1 Niþ2 ¼ g u þ OðDxÞ; @ xa Dxa k¼0 k iþk1
S.-L. Lei, H.-W. Sun / Journal of Computational Physics 242 (2013) 715–725
717
ðaÞ
where g k is defined in (2), and the corresponding implicit finite difference scheme ðmÞ
ui
ðmÞ
ðm1Þ
ui Dt
¼
ðmÞ
iþ1 X ðaÞ ðmÞ dþ;i X d;i Niþ2 ðaÞ ðmÞ ðmÞ g k uikþ1 þ a g u þ fi a Dx k¼0 Dx k¼0 k iþk1
ð3Þ ðaÞ
is unconditionally stable. The alternating fractional binomial coefficient g k have some useful properties, that are observed in [18,19,31], and are summarized in the following proposition. ða Þ
Proposition 1. Let 1 < a < 2 and g k
be defined in (2). We have
8 ðaÞ ðaÞ ðaÞ ðaÞ > < g 0 ¼ 1; g 1 ¼ a < 0; g 2 > g 3 > > 0; 1 n X X ðaÞ ðaÞ > g k ¼ 0; g k < 0; 8n P 1: :
ð4Þ
k¼0 h k¼0 iT h iT ðmÞ ðmÞ ðmÞ ðmÞ ðmÞ ðmÞ Let uðmÞ ¼ u1 ; u2 ; ; uN ; f ðmÞ ¼ f1 ; f2 ; ; fN , and I be the identity matrix with an appropriate size. Then the
numerical scheme (3) can be written in the following matrix form
a Dx Dxa ðm1Þ þ Dxa f ðmÞ ; I þ AðmÞ uðmÞ ¼ u Dt Dt
ð5Þ
with ðmÞ T AðmÞ ¼ DðmÞ þ Ga þ D Ga ;
where
DðmÞ
¼
ð6Þ
ðmÞ ðmÞ diagðd;1 ; ; d;N Þ
2
ðaÞ g1
6 6 ðaÞ 6 g2 6 6 . 6 . 6 . 6 Ga ¼ 6 6 .. 6 . 6 6 6 ðaÞ 6 g N1 4 ðaÞ
gN
ðaÞ g0
0
ðaÞ
g0
g2
ðaÞ
g1
..
.
..
..
.
..
g1
ðaÞ
g N1
and
0
0
..
.
..
.
.
..
.
..
.
.
..
.
g1
ðaÞ
ðaÞ
ðaÞ ðaÞ
g2
0
3
7 7 0 7 7 .. 7 7 . 7 7 7 7 0 7 7 7 ðaÞ 7 g0 7 5
:
ð7Þ
ðaÞ
g1
NN
It is obvious that Ga is a Toeplitz matrix (see [8,9]). Therefore, it can be stored with N þ 1 entries [31]. Furthermore, the matrix–vector multiplication for the Toeplitz-like matrix AðmÞ in (6) can be obtained in OðN log NÞ operations by the FFT; see [8,22]. Denote
mN;M ¼
Dxa M ¼ ðxR xL Þa T 1 Dt ðN þ 1Þa
ð8Þ
which is related to the numbers of time steps and grid points. The linear system (5) can be written as
MðmÞ uðmÞ ¼ b
ðm1Þ
ð9Þ
;
where
MðmÞ ¼ mN;M I þ AðmÞ ¼ mN;M I þ DþðmÞ Ga þ DðmÞ GTa
ð10Þ
and ðm1Þ
b
ðmÞ ¼ mN;M uðm1Þ þ Dtf :
The coefficient matrix MðmÞ in (10) is a strictly diagonally dominant M-matrix (see [31]) and therefore it is nonsingular. 3. PCGNR method with circulant preconditioner The conjugate gradient (CG) method is a popular and effective iterative method for solving symmetric positive definite linear systems. Furthermore, for the nonsymmetric linear system (9), the CG method can also be applied to its equivalent normalized linear system
ðMðmÞ ÞT MðmÞ uðmÞ ¼ ðMðmÞ ÞT b
ðm1Þ
;
which is called the CGNR method [2]. The algorithm of the CGNR method is given in Algorithm 1 below.
718
S.-L. Lei, H.-W. Sun / Journal of Computational Physics 242 (2013) 715–725
ðm1Þ
Algorithm 1. CGNR for MðmÞ uðmÞ ¼ b ðmÞ
1.
Given initial guess u0
2.
r0 ¼ b
3. 4.
ðm1Þ
ðmÞ
MðmÞ u0
ðmÞ T
z0 ¼ ðM Þ r 0 ; p0 ¼ z0 For i ¼ 0; 1; . . ., until convergence, Do
5.
wi ¼ MðmÞ pi
6. 7. 8.
ai ¼ jjzi jj22 =jjwi jj22 ðmÞ ðmÞ uiþ1 ¼ ui þ ai pi r iþ1 ¼ ri ai wi
9.
ziþ1 ¼ ðMðmÞ ÞT riþ1
10. 11. 12.
bi ¼ jjziþ1 jj22 =jjzi jj22 piþ1 ¼ ziþ1 þ bi pi EndDo
In [32], the CGNR method is employed to solve (5). However, the drawback of the CGNR method is at its slow convergence rate, due to the fact that the condition number of the coefficient matrix ðMðmÞ ÞT MðmÞ is usually large. In fact, the numerical results given in Section 5 show that the CGNR method really has a poor convergence rate for ill-conditioned systems. Therefore, preconditioner technique could be exploited to accelerate the convergence rate of the CGNR method. When the CG method, with a matrix P ðmÞ , is employed to solve the following normalized preconditioned system
h iT h i h iT ðm1Þ ðP ðmÞ Þ1 MðmÞ ðP ðmÞ Þ1 MðmÞ uðmÞ ¼ ðP ðmÞ Þ1 MðmÞ ðP ðmÞ Þ1 b ; ðm1Þ
it is equivalent to the original system MðmÞ uðmÞ ¼ b mathematically. This method is called the PCGNR method with a preconditioner P ðmÞ . The algorithm of the PCGNR method is given in Algorithm 2 below, which is the direct application of the ðm1Þ Algorithm 2 to the system ðP ðmÞ Þ1 MðmÞ uðmÞ ¼ ðP ðmÞ Þ1 b . ðm1Þ
Algorithm 2. PCGNR for MðmÞ uðmÞ ¼ b
with preconditioner P ðmÞ .
ðmÞ
1. Given initial guess u0 ðmÞ 1
ðm1Þ
ðmÞ
2. r0 ¼ ðP Þ ðb MðmÞ u0 Þ h iT 3. z0 ¼ ðP ðmÞ Þ1 MðmÞ r0 , p0 ¼ z0 4. For i ¼ 0; 1; . . ., until convergence, Do 5.
wi ¼ MðmÞ pi
6.
wi ¼ ðP ðmÞ Þ1 wi
7.
ai ¼ kzi k22 =kwi k22 ðmÞ ðmÞ uiþ1 ¼ ui þ ai pi r iþ1 ¼ ri ai wi
8. 9. 10.
ziþ1 ¼ ðP ðmÞ ÞT riþ1
11.
ziþ1 ¼ ðMðmÞ ÞT ziþ1
12. bi ¼ kziþ1 k22 =kzi k22 13. piþ1 ¼ ziþ1 þ bi pi 14. EndDo
Now we propose a circulant preconditioner, which is generated from the Strang’s circulant preconditioner [9], in the PCGNR method to solve (9). The Strang’s circulant matrix sðBÞ ¼ ½sjk 06j;k
8 bj ; > > > < 0; sj ¼ > bjN ; > > : sjþN ;
0 6 j < N=2; j ¼ N=2 if N is even; N=2 < j < N; 0 < j < N:
Recall the formula (10) that
719
S.-L. Lei, H.-W. Sun / Journal of Computational Physics 242 (2013) 715–725
MðmÞ ¼ mN;M I þ AðmÞ ¼ mN;M I þ DþðmÞ Ga þ DðmÞ GTa ; where ðmÞ
ðmÞ
DðmÞ ¼ diagðd;1 ; ; d;N Þ; and Ga is the Toeplitz matrix given in (7). Then our circulant preconditioner is defined as
ðmÞ sðGa Þ þ d ðmÞ sðGT Þ; S ðmÞ ¼ mN;M I þ d þ a
ð11Þ
where N X ðmÞ ðmÞ ¼ 1 d d : N i¼1 ;i
More precisely, the first columns of both sðGa Þ and sðGTa Þ are given by
2
ða Þ g1
2
3
ðaÞ
g1
3
6 ðaÞ 7 6 g0 7 6 . 7 7 6 6 . 7 6 0 7 6 . 7 7 6 6 ðaÞ 7 6 .. 7 6 g Nþ1 7 6 . 7 6 b 2 c7 7 6 7 6 7; 6 0 7 and 6 0 7 6 6 . 7 6 ðaÞ 7 6 . 7 6 g Nþ1 7 6 . 7 c b 6 2 7 7 6 7 6 4 0 5 6 .. 7 4 . 5 ða Þ g0 ðaÞ g2 respectively. We know that a circulant matrix C can be diagonalized by the Fourier matrix F [9]; i.e.,
C ¼ F KF; where the entries of F are given by
1 F j;k ¼ pffiffiffiffi e2pijk=N ; N
0 6 j; k 6 N 1;
with the imaginary unit i, and K is a diagonal matrix holding the eigenvalues of C. a is the complex conjugate of Ka . The following lemma is the key to prove a F, where K Let sðGa Þ ¼ F Ka F. Then sðGTa Þ ¼ F K the invertibility of S ðmÞ in (11). Lemma 1. All eigenvalues of sðGa Þ and sðGTa Þ fall inside the open disc
fz 2 C : jz aj < ag: ðaÞ
Proof. All the Gershgorin disc of the circulant matrices sðGa Þ and sðGTa Þ are centered at g 1 ¼ a with radius ðaÞ
rN ¼ g 0 þ
bNþ1 c 2
X
ðaÞ
gk <
1 X
ðaÞ
ðaÞ
g k ¼ g 1 ¼ a;
k¼0;k–1
k¼2
ðaÞ
by the properties of the sequence fg k g; see Proposition 1. h Remark 1. It is worth to note that: 1. The real parts of all eigenvalues of sðGa Þ and sðGTa Þ are strictly positive for all N. 2. The absolute values of all eigenvalues of sðGa Þ and sðGTa Þ are bounded above by 2a < 4 for all N. Ka þ d ðmÞ K a . Then S ðmÞ is invertible Decompose the circulant matrix S ðmÞ ¼ F Ks F with the diagonal matrix Ks ¼ mN;M I þ d þ if all diagonal entries of Ks are nonzero. ðmÞ
Lemma 2. Let 1 < a < 2. The preconditioner S ðmÞ in (11) is invertible and
kðS ðmÞ Þ1 k2 6
1
mN;M
:
ð12Þ
720
S.-L. Lei, H.-W. Sun / Journal of Computational Physics 242 (2013) 715–725
Proof. By Part 1 of Remark 1, we have Reð½Ka k;k Þ > 0. Noting that
mN;M > 0 and dðmÞ P 0, we have
ðmÞ Reð½Ka Þ þ d ðmÞ Reð½K a Þ P mN;M > 0 j½Ks k;k j P Reð½Ks k;k Þ ¼ mN;M þ d þ k;k k;k for each k ¼ 1; . . . ; N. Therefore, S ðmÞ is invertible. Furthermore, we have
kðS ðmÞ Þ1 k2 ¼
1 1 : 6 min j½Ks k;k j mN;M
16k6N
h The computational cost per iteration of the PCGNR method with circulant preconditioners is OðN log NÞ complexity, as the main workload is the matrix–vector multiplication. Therefore, if the iterative method is of linear or superlinear convergence, the total cost for solving (9) at each time step by the proposed method will keep OðN log NÞ operations. Before to analyze the convergence of the proposed method, we prove the following theorem which gives a bound of kS ðmÞ k2 . This theorem will be exploited to prove the superlinear convergence rate of the PCGNR method in next section. ~ We have Lemma 3. Let 1 < a < 2 and 0 6 d ðx; tÞ 6 d.
~ kS ðmÞ k2 6 mN;M þ 8d:
ð13Þ
Proof. For any k ¼ 1; . . . ; N, by Part 2 of Remark 1,
½ Ka þ d ðmÞ ½K j½Ka j þ d ðmÞ j½K þ 4d ðmÞ 6 mN;M þ 8d: ~ a j 6 mN;M þ d a j < mN;M þ 4d j½Ks k;k j ¼ jmN;M þ d þ þ þ k;k k;k k;k k;k ðmÞ
ðmÞ
ðmÞ
Then
~ kS ðmÞ k2 ¼ max j½Ks k;k j < mN;M þ 8d: 16k6N
h
4. Spectrum of the preconditioned matrix In this section, we study the convergence rate of the PCGNR method with the proposed circulant preconditioner S ðmÞ in (11) for the linear system (9). For convenience of our investigation, we assume that the left and the right diffusion coefficient functions d ðx; tÞ of the FDE (1) are two nonnegative constants and M is properly chosen, depending on N, such that mN;M in and m ^ such that (8) is bounded away from 0; i.e., there exist two real numbers m
6 mN;M 6 m ^; 0
8N:
ð14Þ
By this assumption, for all i ¼ 1; . . . ; N, ðmÞ
dþ;i ¼ dþ P 0;
ðmÞ
d;i ¼ d P 0;
and dþ þ d – 0:
ð15Þ
Throughout this section, we add a subscript N to each matrix to denote the matrix size. Under the assumption in (15), the matrix AðmÞ defined in (6), MðmÞ defined in (10), and S ðmÞ defined in (11) are independent of m, and we therefore simply denote them as AN ; MN , and S N , respectively. Now the coefficient matrix MðmÞ in (10) becomes
MN ¼ mN;M I þ AN ;
ð16Þ
which is a nonsymmetric Toeplitz matrix, where
2
ðaÞ
ðaÞ
dþ g 1 þ d g 1
6 6 6 dþ g ð2aÞ þ d g ð0aÞ 6 6 ðaÞ AN ¼ dþ Ga þ d GTa ¼ 6 6 dþ g 3 6 6 .. 6 . 4 ðaÞ
dþ g N
ðaÞ
ðaÞ
ðaÞ
ðaÞ
dþ g 0 þ d g 2
ðaÞ
dþ g 1 þ d g 1 .. . .. .
.. . .. . .. .
...
dþ g 2 þ d g 0
ða Þ
3
ðaÞ
d g 3 .. . .. . .. .
d g N .. . .. .
ðaÞ
7 7 7 7 7 7 ½ajk : NN 7 7 7 ðaÞ ðaÞ dþ g 0 þ d g 2 7 5 ðaÞ
ðaÞ
dþ g 1 þ d g 1
To study the spectrum of the preconditioned matrix
1 T 1 S N MN S N MN ; we first introduce the generating function of the sequence of Toeplitz matrices fAN g1 N¼1 [8]:
ð17Þ
S.-L. Lei, H.-W. Sun / Journal of Computational Physics 242 (2013) 715–725
pðhÞ ¼
1 X
ak eikh ;
721
ð18Þ
k¼1
where ak is the k-th diagonal of AN . The generating function pðhÞ is in the Wiener class [8,9] if and only if 1 X
jak j < 1:
k¼1
For AN defined in (17), we have
pðhÞ ¼
1 X
ak eikh ¼
k¼1
We remark that
1 X
ðaÞ
g kþ1 ðdþ eikh þ d eikh Þ:
ð19Þ
k¼1
mN;M þ pðhÞ cannot be a generating function of MN in (16) since mN;M is dependent on N.
T 1 Lemma 4. Let p and q be the generating functions of fAN g1 N¼1 and fAN þ AN gN¼1 , respectively, we have
1. p is in the Wiener class; 2. q is real-valued and nonnegative.
Proof ðaÞ
1. By the properties of the sequence fg k g1 k¼0 given in (4), we have 1 X
jak j ¼ ðdþ þ
k¼1
1 X
ðaÞ d Þ jg kþ1 j k¼1
¼ ðdþ þ d Þ
ðaÞ 2g 1
1 X ðaÞ þ gk
!
¼ 2aðdþ þ d Þ < 1:
k¼0
Thus p is in the Wiener class. 2. From (4) and the definition of generating function in (18), we have
qðhÞ ¼
1 X
" # 1 X ðaÞ ða Þ ðaÞ ih ðaÞ ðak þ ak Þeikh ¼ ðdþ þ d Þ 2g 1 þ g 0 þ g 2 g kþ1 eikh þ eikh e þ eih þ
k¼1
k¼2
"
# 1 1 X X ðaÞ ðaÞ ðaÞ ðaÞ ðaÞ ¼ ðdþ þ d Þ 2g 1 þ 2 g 0 þ g 2 cos h þ 2 g kþ1 cos kh P 2ðdþ þ d Þ g kþ1 ¼ 0: k¼2
k¼1
h
From Part 2 of Lemma 4, and according to the Grenander–Szegö Theorem (see Theorem 1.12 in [7]), the matrix AN þ ATN is symmetric positive semidefinite. Furthermore, by noting that
MN MTN ¼ m2N;M I þ mN;M AN þ ATN þ AN ATN ; we conclude that the smallest singular value of the matrix MN is bounded below by
rmin ðMN Þ P mN;M :
mN;M ; i.e., ð20Þ
By Part 1 of Lemma 4, we have the following lemma. Lemma 5 (see [8]). If pðhÞ, the generating function of AN , is in the Wiener class, then for any > 0, there exist N0 and M0 > 0, such that for all N > N0,
AN sðAN Þ ¼ U N þ V N where
rank ðU N Þ 6 M0 and
kV N k2 < : Now we consider the spectrum of S 1 N MN I. By Lemma 5, we have 1 1 1 S 1 N MN I ¼ S N ðMN S N Þ ¼ S N U N S N V N :
722
S.-L. Lei, H.-W. Sun / Journal of Computational Physics 242 (2013) 715–725
0 Note that rankðS 1 N U N Þ 6 rankðU N Þ 6 M and
1 kS 1 N V N k2 6 kS N k2 kV N k2 <
6 ; mN;M m
by the assumption (14) and the inequality (12) in Lemma 2. Therefore, the following corollary is resulted. Corollary 1. If
mN;M satisfies the assumption (14), for any > 0, there exists N00 and M00 > 0 such that, for any N > N00 ,
e e S 1 N MN I ¼ U N þ V N ;
ð21Þ
e N Þ 6 M 00 and k V e N k < . where rankð U 2 Now we show that the spectrum of the normalized preconditioned matrix
1 T 1 S N MN S N MN is clustered around 1. Theorem 1. If mN;M satisfies the assumption (14), for any 0 < < 1, there exists N 00 and M 00 > 0 such that, for all N > N 00 , at most 2M00 eigenvalues of the matrix
1 T 1 S N MN S N MN I have absolute values larger than 3. Proof. By (21), we have
1 T 1 e N ÞT ðI þ U eN þ V e NÞ ¼ I þ U bN þ V b N; eN þ V S N MN S N MN ¼ ðI þ U where
bN ¼ U e T ðI þ U eN þ V e N Þ þ ðI þ V e T ÞU eN U N N and
eN þ V eT þ V eT V bN ¼ V e V N N N: b N Þ 6 2M00 and k V b N k < 3. Thus, we have Then by Corollary 1, we see that rankð U 2
1 T 1 b N; bN þ V S N MN S N MN I ¼ U b N and V b N are real and symmetric. By the well-known Weyl’s theorem, we conclude that at most 2M00 eigenwhere both U values of the matrix
1 T 1 S N MN S N MN I have absolute values larger than 3.
h
For the superlinear convergence of the PCGNR method, according to the Corollary 1.11 in [7], it is further required that the smallest singular value of the matrix S 1 N MN is uniformly bounded away from zero. In fact, using (13), (20), and the assumption (14), we have
rmin ðS 1 N MN Þ P
rmin ðMN Þ kS N k2
P
mN;M m^ P > 0; mN;M þ 8d~ m þ 8d~
~ ¼ maxfdþ ; d g. Therefore, the PCGNR method for solving (9) converges superlinearly. where d 5. Numerical results In this section, we solve the FDE problem (1) numerically by the implicit finite difference method given in Section 2. After ðm1Þ the finite difference discretization, at each time step, the nonsymmetric linear system MðmÞ uðmÞ ¼ b is solved by the CGNR method (Algorithm 1), the PCGNR method (Algorithm 2), and the multigrid method [22], respectively. Number of iterations required for convergence and CPU time of those methods are reported. The stopping criterion of those methods is
krðkÞ k2 < 107 ; kr ð0Þ k2 where rðkÞ is the residual vector of the linear system after k iterations, and the initial guess at each time step is chosen as the zero vector. In the following tables, ‘‘N’’ denotes the number of spatial grid points, ‘‘M’’ denotes the number of time steps,
723
S.-L. Lei, H.-W. Sun / Journal of Computational Physics 242 (2013) 715–725
‘‘Error’’ denotes the error between the true solution and the approximation under the infinity norm, ‘‘CPU (s)’’ denotes the total CPU time in seconds for solving the whole FDE problem, and ‘‘Iter’’ denotes the average number of iterations required for solving the FDE problem; i.e., M 1X IterðmÞ; M m¼1
Iter ¼
where IterðmÞ denotes the number of iterations required for solving (9). For the PCGNR method, besides the proposed circulant preconditioner S ðmÞ in (11), we also test the Chan’s preconditioner [10] which can be written as
cðGa Þ þ d ðmÞ cðGT Þ; cðMðmÞ Þ ¼ mN;M I þ d þ a ðmÞ
where cðBÞ denotes the Chan’s preconditioner for an arbitrary matrix B. More precisely, the first columns of the circulant matrices cðGa Þ and cðGTa Þ are given as
2
6 6 6 16 6 N6 6 6 4
3
ðaÞ
Ng 1
ða Þ
ðN 1Þg 2 .. . ðaÞ
2g N1 ðaÞ gN
ðaÞ
2
3
ðaÞ
Ng 1
7 7 6 7 6 ðN 1Þg ðaÞ þ g ðaÞ 7 N 7 0 7 6 7 7 6 ðaÞ 7 and 1 6 7; 2g N1 7 7 6 N 7 7 6 . 7 7 6 . 5 5 4 . ðaÞ
þ ðN 1Þg 0
ðN 1Þg 2
respectively. We remark that the superlinear convergence of the PCGNR method with the preconditioner cðMðmÞ Þ, as that with S ðmÞ in Section 4, can be analogously proven if the assumptions (14) and (15) hold. In the following tables, we denote ‘‘I ’’ as the CGNR method, ‘‘C’’ as the PCGNR method with the T. Chan’s circulant preconditioner cðMðmÞ Þ, ‘‘S’’ as the PCGNR method with the proposed circulant preconditioner S ðmÞ , and ‘‘MGM’’ as the multigrid method in [22]. In Table 2, ‘‘Error’’ denotes the error between the exact solution and the numerical solution obtained by the PCGNR method with the proposed preconditioner SðmÞ . All numerical experiments are run in MATLAB 7.12 (R2011a) on a Dell workstation Precision T7500 with the configuration: Intel (R) Xeon (R) CPU X5690 3.47 GHz and 12 GB RAM. Example 1. In this example, we solve the initial-boundary value problem of FDE (1) with source term f ðx; tÞ 0, for order of fractional derivatives a ¼ 1:2; 1:5 and 1:8. The spatial domain is ½xL ; xR ¼ ½0; 2 and the time interval is ½0; T ¼ ½0; 1. The initial condition uðx; 0Þ is the following Gaussian pulse
uðx; 0Þ ¼ exp
! ðx xc Þ2 ; 2r 2
xc ¼ 1:2;
r ¼ 0:08;
and the diffusion coefficients are
dþ ðx; tÞ 0:6;
and d ðx; tÞ 0:5:
Table 1 Comparisons for solving Example 1 by the CGNR method, the PCGNR method with preconditioners cðMðmÞ Þ and S ðmÞ , and the multigrid method for a ¼ 1:2; 1:5, and 1:8. Here Dt 2Dxa .
a 1.2
1.5
1.8
Nþ1
M
I
C
MGM
S
Iter
CPU (s)
Iter
CPU (s)
Iter
CPU (s)
Iter
CPU (s)
6
2
32
37.6
0.10
6.0
0.05
5.8
0.04
10.1
0.30
27
74
34.4
0.36
6.0
0.12
5.3
0.11
8.8
0.65
28
169
31.4
0.96
5.0
0.28
5.0
0.28
7.3
1.56
29
388
28.5
4.08
5.0
1.29
5.0
1.27
6.3
3.93
210
891
25.7
15.22
5.0
4.49
5.0
4.44
5.5
10.54
2
91
40.9
0.30
6.0
0.11
5.6
0.10
7.0
0.52
27
256
39.2
1.40
6.0
0.43
5.2
0.37
6.1
1.59
28
724
35.8
4.70
5.4
1.28
5.0
1.20
5.6
5.25
29
2048
32.3
24.21
5.0
6.74
5.0
6.73
5.1
16.91
210
5793
29.0
110.65
5.0
28.92
5.0
28.93
4.8
60.20
6
2
256
42.6
0.85
7.0
0.31
5.8
0.26
6.1
1.25
27
891
41.0
5.01
6.0
1.45
5.5
1.36
6.0
5.50
28
3104
36.3
20.59
6.0
5.92
5.3
5.38
5.0
19.93
29
10809
31.8
126.14
5.2
36.55
5.1
35.88
5.0
88.35
210
37641
27.5
685.80
5.0
187.06
5.0
187.07
4.1
339.35
6
724
S.-L. Lei, H.-W. Sun / Journal of Computational Physics 242 (2013) 715–725
Table 2 Comparisons for solving Example 2 by the CGNR method, the PCGNR method with preconditioners cðMðmÞ Þ and S ðmÞ , and the multigrid method for a ¼ 1:2; 1:5 and 1:8.
a 1:2
1:5
1:8
Nþ1
M
Error
I
C
MGM
S
Iter
CPU (s)
Iter
CPU (s)
Iter
CPU (s)
Iter
CPU (s) 0.19
2
5
2
3.1501e2
33.8
0.08
8.0
0.05
8.0
0.05
8.0
27
26
1.5983e2
65.5
0.53
8.0
0.12
8.0
0.12
7.1
0.41
28
27
8.0488e3
82.0
1.70
8.0
0.29
7.0
0.26
7.1
1.03
29
28
4.0377e3
89.9
7.92
8.0
1.18
8.0
1.18
7.0
2.60
210
29
2.0214e3
96.0
30.90
7.0
3.32
8.0
3.73
7.0
6.91
26
25
2.2529e2
46.6
0.11
10.0
0.06
8.0
0.05
8.0
0.19
27
26
1.1164e2
111.6
0.87
10.4
0.16
9.0
0.14
7.0
0.41
28
27
5.5563e3
264.6
5.34
10.9
0.38
9.3
0.33
7.0
1.03
29
28
2.7721e3
568.3
48.74
9.9
1.44
9.9
1.44
7.0
2.62
210
29
1.3838e3
903.7
288.51
11.0
5.01
10.0
4.52
6.0
5.93
2
5
2
1.7434e2
70.6
0.16
16.0
0.08
13.0
0.07
7.0
0.16
27
26
8.3524e3
202.0
1.58
18.0
0.26
14.0
0.21
7.0
0.40
28
27
4.0838e3
587.2
11.68
18.9
0.63
14.0
0.48
7.0
1.01
29
28
2.0186e3
1703.0
146.05
21.0
2.89
14.0
1.97
7.0
2.60
210
29
1.0035e3
4900.0
1541.66
20.0
8.66
13.0
5.78
6.0
5.92
6
6
For satisfying (14), we choose M ¼ bð2N þ 2Þa =2c such that Dt 2Dxa . Thus, the superlinear convergence rate by the PCGNR method with the proposed circulant preconditioner is guaranteed as both (14) and (15) hold. We remark that the linear convergence rate by the multigrid method in [22] has not been proven theoretically. The numerical results of Example 1 are shown in Table 1. It shows that both the average number of iterations and the CPU time by the PCGNR method with circulant preconditioners are much less than those by the CGNR method. It also shows that the CPU times by the PCGNR methods are less than that by the multigrid method.
Example 2. In this example, we study the case for which the diffusion coefficients are variable. We solve the FDE problem (1) of order a with the following data. The spatial domain is ½xL ; xR ¼ ½0; 2 and the time interval is ½0; T ¼ ½0; 1. The left and right diffusion coefficients are
dþ ðx; tÞ ¼ Cð3 aÞxa
and d ðx; tÞ ¼ Cð3 aÞð2 xÞa ;
respectively. The source term is
1 3 3 f ðx; tÞ ¼ 32et x2 þ ð2 xÞ2 ð8 þ x2 Þ ½x3 þ ð2 xÞ3 þ ½x4 þ ð2 xÞ4 8 3a ð4 aÞð3 aÞ and the initial condition is
uðx; 0Þ ¼ 4x2 ð2 xÞ2 : The exact solution of this problem is
uðx; tÞ ¼ 4et x2 ð2 xÞ2 for any a 2 ð1; 2Þ. With the exact solution, we calculate the exact error of the numerical solution. In the following tests, we consider a ¼ 1:2; 1:5, and 1:8. We remark that the same example in [19,22,31] is only considered for a ¼ 1:8. In the numerical test, we naturally choose Dt ¼ Dx. Table 2 reports the history for numerically solving Example 2. Note that the average number of iterations for the CGNR method is increasing as N is increasing rapidly, while the average numbers of iterations by the PCGNR methods almost keep constant, and therefore their CPU times are much less than that by the CGNR method. We remark that both diffusion coefficients are not constant and mN;M tends to 0. Therefore, both assumptions (14) and (15) are not satisfied and hence the superlinear convergence rate is not theoretically guaranteed. However, from the point of view on CPU time, the performance of the proposed circulant preconditioner still works well and is the best among all the methods.
6. Concluding remarks In this paper, we have employed the PCGNR method with a circulant preconditioner to solve the discretized linear systems of (1). Theoretically, we prove the superlinear convergence rate of the proposed method under the conditions that the
S.-L. Lei, H.-W. Sun / Journal of Computational Physics 242 (2013) 715–725
725
a
diffusion coefficients d ðx; tÞ are constant and the ratio DDxt is bounded away from zero. Numerical tests have shown the efficiency of the PCGNR method. Moreover, for variable diffusion coefficients in the FDE (1), the proposed algorithm still works very well. In our future consideration, more effective preconditioners will be developed to solve the two or more dimensional FDE problem [30] and the steady-state fractional partial differential equations. Acknowledgement The authors thank Dr. Hong-kui Pang for his helpful discussions. We are also grateful to the anonymous referees for the useful suggestions and comments that improved the presentation of this paper. References [1] J. Bai, X. Feng, Fractional-order anisotropic diffusion for image denoising, IEEE Trans. Image Proc. 16 (2007) 2492–2502. [2] R. Barrett, M. Berry, T. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo, C. Romine, H. van der Vorst, Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods, SIAM, Philadelphia, 1994. [3] D. Benson, S.W. Wheatcraft, M.M. Meerschaert, Application of a fractional advection–dispersion equation, Water Resour. Res. 36 (2000) 1403–1413. [4] D. Benson, S.W. Wheatcraft, M.M. Meerschaert, The fractional-order governing equation of Lévy motion, Water Resour. Res. 36 (2000) 1413–1423. [5] B. Beumer, M. Kovács, M.M. Meerschaert, Numerical solutions for fractional reaction–diffusion equations, Comput. Math. Appl. 55 (2008) 2212–2226. [6] B.A. Carreras, V.E. Lynch, G.M. Zaslavsky, Anomalous diffusion and exit time distribution of particle tracers in plasma turbulence models, Phys. Plasma 8 (2001) 5096–5103. [7] R. Chan, X. Jin, An Introduction to Iterative Toeplitz Solvers, SIAM, Philadelphia, 2007. [8] R. Chan, M. Ng, Conjugate gradient methods for Toeplitz systems, SIAM Rev. 38 (1996) 427–482. [9] R. Chan, G. Strang, Toeplitz equations by conjugate gradients with circulant preconditioner, SIAM J. Sci. Stat. Comput. 10 (1989) 104–119. [10] T. Chan, An optimal circulant preconditioner for Toeplitz systems, SIAM J. Sci. Stat. Comput. 9 (1988) 766–771. [11] W. Deng, Finite element method for the space and time fractional Fokker–Planck equation, SIAM J. Numer. Anal. 47 (2008) 204–226. [12] V.J. Ervin, N. Heuer, J.P. Roop, Numerical approximation of a time dependent, nonlinear, space-fractional diffusion equation, SIAM J. Numer. Anal. 45 (2007) 572–591. [13] X. Jin, Preconditioning Techniques for Teoplitz Systems, Higher Education Press, Beijing, 2010. [14] T.A.M. Langlands, B.I. Henry, The accuracy and stability of an implicit solution method for the fractional diffusion equation, J. Comput. Phys. 205 (2005) 719–736. [15] F. Liu, V. Anh, I. Turner, Numerical solution of the space fractional Fokker–Planck equation, J. Comput. Appl. Math. 166 (2004) 209–219. [16] R.L. Magin, Fractional Calculus in Bioengineering, Begell House Publishers, 2006. [17] M.M. Meerschaert, H.P. Scheffler, C. Tadjeran, Finite difference methods for two-dimensional fractional dispersion equation, J. Comput. Phys. 211 (2006) 249–261. [18] M.M. Meerschaert, C. Tadjeran, Finite difference approximations for fractional advection–dispersion flow equations, J. Comput. Appl. Math. 172 (2004) 65–77. [19] M.M. Meerschaert, C. Tadjeran, Finite difference approximations for two-sided space-fractional partial differential equations, Appl. Numer. Math. 56 (2006) 80–90. [20] D.A. Murio, Implicit finite difference approximation for time fractional diffusion equations, Comput. Math. Appl. 56 (2008) 1138–1145. [21] M. Ng, Iterative Methods for Toeplitz Systems, Numerical Mathematics and Scientific Computation, Oxford University Press, New York, 2004. [22] H. Pang, H. Sun, Multigrid method for fractional diffusion equations, J. Comput. Phys. 231 (2012) 693–703. [23] I. Podlubny, Fractional Differential Equations, Academic Press, New York, 1999. [24] M. Raberto, E. Scalas, F. Mainardi, Waiting-times and returns in high-frequency financial data: an empirical study, Physica A 314 (2002) 749–755. [25] M.F. Shlesinger, B.J. West, J. Klafter, Lévy dynamics of enhanced diffusion: application to turbulence, Phys. Rev. Lett. 58 (1987) 1100–1103. [26] I.M. Sokolov, J. Klafter, A. Blumen, Fractional kinetics, Phys. Today Nov. (2002) 28–53. [27] E. Sousa, Finite difference approximates for a fractional advection diffusion problem, J. Comput. Phys. 228 (2009) 4038–4054. [28] L. Su, W. Wang, Z. Yang, Finite difference approximations for the fractional advection–diffusion equation, Phys. Lett. A 373 (2009) 4405–4408. [29] C. Tadjeran, M.M. Meerschaert, H.P. Scheffler, A second-order accurate numerical approximation for the fractional diffusion equation, J. Comput. Phys. 213 (2006) 205–213. 2 [30] H. Wang, K. Wang, An OðNlog NÞ alternating-direction finite difference method for two-dimensional fractional diffusion equations, J. Comput. Phys. 230 (2011) 7830–7839. 2 [31] H. Wang, K. Wang, T. Sircar, A direct OðNlog NÞ finite difference method for fractional diffusion equations, J. Comput. Phys. 229 (2010) 8095–8104. [32] K. Wang, H. Wang, A fast characteristic finite difference method for fractional advection–diffusion equations, Adv. Water Resour. 34 (2011) 810–816. [33] G.M. Zaslavsky, D. Stevens, H. Weitzner, Self-similar transport in incomplete chaos, Phys. Rev. E 48 (1993) 1683–1694.
Contents lists available at SciVerse ScienceDirect
Journal of Computational Physics journal homepage: www.elsevier.com/locate/jcp
A circulant preconditioner for fractional diffusion equations q Siu-Long Lei, Hai-Wei Sun ⇑ Department of Mathematics, University of Macau, Macao, China
a r t i c l e
i n f o
Article history: Received 10 September 2012 Received in revised form 4 February 2013 Accepted 12 February 2013 Available online 26 February 2013 Keywords: Fractional diffusion equations Shifted Grünwald discretization Toeplitz Circulant preconditioner Fast Fourier transform CGNR method
a b s t r a c t The implicit finite difference scheme with the shifted Grünwald formula, which is unconditionally stable, is employed to discretize fractional diffusion equations. The resulting systems are Toeplitz-like and then the fast Fourier transform can be used to reduce the computational cost of the matrix–vector multiplication. The preconditioned conjugate gradient normal residual method with a circulant preconditioner is proposed to solve the discretized linear systems. The spectrum of the preconditioned matrix is proven to be clustered around 1 if diffusion coefficients are constant; hence the convergence rate of the proposed iterative algorithm is superlinear. Numerical experiments are carried out to demonstrate that our circulant preconditioner works very well, even though for cases of variable diffusion coefficients. Ó 2013 Elsevier Inc. All rights reserved.
1. Introduction In the last decade or so, the catholicity of anomalous diffusion phenomena in the real world has led to the fractional diffusion equation (FDE). FDEs arise in research topics including modeling chaotic dynamics of classical conservative systems [33], groundwater contaminant transport [3,4], turbulent flow [6,25], and applications in biology [16], finance [24], image processing [1], and physics [26]. As there are very few cases of FDEs in which the closed-form analytical solutions are available, numerical solutions for FDEs become major ways and then have been developed intensively [5,11,12,14,15,17–20,27–29]. Nevertheless, since the fractional differential operator is nonlocal, it was shown that a naive discretization of the FDE, even though implicit, leads to unconditionally unstable [18,19]. Moreover, most numerical methods for FDEs tend to generate full coefficient matrices, which require computational cost of OðN 3 Þ and storage of OðN 2 Þ, where N is the number of grid points [31]. It is quite different from second-order diffusion equations which usually yield sparse coefficient matrices with OðNÞ nonzero entries and can be solved very efficiently by fast iterative methods with OðNÞ complexity. To overcome the difficulty of the stability, Meerschaet and Tadjeran [18,19] proposed a shifted Grünwald discretization to approximate FDEs. Their method has been proven to be unconditionally stable. Later, Wang, Wang, and Sircar [31] discovered that the full coefficient matrix by the Meerschaet–Tadjeran’s method holds a Toeplitz-like structure. More precisely, such a full matrix can be written as the sum of diagonal-multiply-Toeplitz matrices. Thus the storage requirement is significantly reduced from OðN 2 Þ to OðNÞ. It is well known that the matrix–vector multiplication for the Toeplitz matrix can be computed by the fast Fourier transform (FFT) with OðN log NÞ operations [8,21]. With this advantage, Wang and Wang [32] employed the conjugate gradient normal residual (CGNR) method to solve the discretized system of the FDE by the Meerschaet–Tadjeran’s method. Thanks to the Toeplitz-like structure, the cost per iteration by the CGNR method is of
q This work was partially supported by the research grant MYRG038(Y2-L1)-FST12-LSL and MYRG140(Y3-L2)-FST11-SHW from University of Macau, and the research grant 105/2012/A3 from FDCT of Macao. ⇑ Corresponding author. E-mail addresses: [email protected] (S.-L. Lei), [email protected] (H.-W. Sun).
0021-9991/$ - see front matter Ó 2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.jcp.2013.02.025
716
S.-L. Lei, H.-W. Sun / Journal of Computational Physics 242 (2013) 715–725
OðN log NÞ. The convergence of the CGNR method is fast with smaller diffusion coefficients [32] (in that case the discretized system is well-conditioned). Nevertheless, if the diffusion coefficient functions are not small, the resulting system will become ill-conditioned and hence the CGNR method converges very slowly. To overtake this shortcoming, Pang and Sun [22] proposed a multigrid method to solve the discretized system of the FDE by the Meerschaet–Tadjeran’s method. With the damped-Jacobi method as the smoother, the multigrid algorithm can preserve the computational cost per iteration as OðN log NÞ operations. Numerical results showed that their multigrid method converges very fast, even for the ill-conditioned systems. However, from the theoretical point of view, the linear convergence of their multigrid method, despite a very simple case (both diffusion coefficients are equal and constant), has not been proven; see [22] for details. In the literature on Toeplitz or Toeplitz-like systems, circulant preconditioners always played important roles. In fact, circulant preconditioners have been theoretically and numerically studied with numerous applications for over twenty years; see for instances [7,8,21]. Nevertheless, to our knowledge, circulant preconditioners were never particularly applied to the Toeplitz-like matrix with structure as the sum of diagonal-multiply-Toeplitz matrices before. In this paper, we propose the preconditioned CGNR (PCGNR) method with a circulant preconditioner to solve such Toeplitz-like systems. The proposed preconditioner is invertible and its spectrum is theoretically proven to be clustered around 1 under some conditions. Thus the superlinear convergence rate of the PCGNR method is obtained. As the computational cost per iteration is of OðN log NÞ, the total complexity of the PCGNR method at each time step retains OðN log NÞ operations. The paper is organized as follows. In Section 2, the background of the discretization for the FDE is reviewed. In Section 3, the PCGNR method with a circulant preconditioner is proposed to solve the discretized linear system. The convergence rate of the PCGNR method is studied in Section 4. In Section 5, numerical results are reported to demonstrate the efficiency of the proposed method. Concluding remarks are given in Section 6. 2. FDE and finite difference discretization In this paper, we study an initial-boundary value problem of the FDE:
8 @uðx; tÞ @ a uðx; tÞ @ a uðx; tÞ > > > ¼ dþ ðx; tÞ þ d ðx; tÞ þ f ðx; tÞ; > a > @t @þx @ xa > > > < t 2 ð0; T; x 2 ðxL ; xR Þ; > > > > uðxL ; tÞ ¼ uðxR ; tÞ ¼ 0; 0 6 t 6 T; > > > > : x 2 ½xL ; xR ; uðx; 0Þ ¼ u0 ðxÞ;
ð1Þ
where a 2 ð1; 2Þ is the order of the fractional derivative, f ðx; tÞ is the source term, and diffusion coefficient functions d ðx; tÞ a
a
@ uðx;tÞ are nonnegative; i.e., d ðx; tÞ P 0. The left-sided and the right-sided fractional derivatives @ @uðx;tÞ a and @ xa are defined in the þx
Grünwald–Letnikov form [23]
@ a uðx; tÞ 1 ¼ limþ a @ þ xa Dx!0 Dx
bðxx L Þ=Dxc X
@ a uðx; tÞ 1 ¼ limþ a @ xa Dx!0 Dx
bðxRX xÞ=Dxc
ðaÞ
g k uðx kDx; tÞ;
k¼0
ðaÞ
g k uðx þ kDx; tÞ;
k¼0 ðaÞ
where bc denotes the floor function, and g k is the alternating fractional binomial coefficient given as
(
ðaÞ
g 0 ¼ 1; ðaÞ
k
aða 1Þ ða k þ 1Þ; g k ¼ ð1Þ k!
k ¼ 1; 2; 3; . . . ;
ð2Þ
which can be evaluated by the recurrence relation ðaÞ
g kþ1 ¼
a þ 1 ðaÞ g ; 1 kþ1 k
k ¼ 0; 1; 2; . . . :
R xL Let N and M be positive integers, and Dx ¼ xNþ1 and Dt ¼ T=M be the sizes of spatial grid and time step, respectively. We define a spatial and temporal partition xi ¼ xL þ iDx for i ¼ 0; 1; . . . ; N þ 1 and tm ¼ mDt for m ¼ 0; 1; . . . ; M. Let
ðmÞ
ui
ðmÞ
ðmÞ
¼ uðxi ; t m Þ; d;i ¼ d ðxi ; tm Þ, and fi
¼ f ðxi ; t m Þ. The shifted Grünwald approximation in [18,19] is as follows,
iþ1 @ a uðxi ; tm Þ 1 X ðaÞ ðmÞ ¼ a g u þ OðDxÞ; a @þx Dx k¼0 k ikþ1
X ðaÞ ðmÞ @ a uðxi ; tm Þ 1 Niþ2 ¼ g u þ OðDxÞ; @ xa Dxa k¼0 k iþk1
S.-L. Lei, H.-W. Sun / Journal of Computational Physics 242 (2013) 715–725
717
ðaÞ
where g k is defined in (2), and the corresponding implicit finite difference scheme ðmÞ
ui
ðmÞ
ðm1Þ
ui Dt
¼
ðmÞ
iþ1 X ðaÞ ðmÞ dþ;i X d;i Niþ2 ðaÞ ðmÞ ðmÞ g k uikþ1 þ a g u þ fi a Dx k¼0 Dx k¼0 k iþk1
ð3Þ ðaÞ
is unconditionally stable. The alternating fractional binomial coefficient g k have some useful properties, that are observed in [18,19,31], and are summarized in the following proposition. ða Þ
Proposition 1. Let 1 < a < 2 and g k
be defined in (2). We have
8 ðaÞ ðaÞ ðaÞ ðaÞ > < g 0 ¼ 1; g 1 ¼ a < 0; g 2 > g 3 > > 0; 1 n X X ðaÞ ðaÞ > g k ¼ 0; g k < 0; 8n P 1: :
ð4Þ
k¼0 h k¼0 iT h iT ðmÞ ðmÞ ðmÞ ðmÞ ðmÞ ðmÞ Let uðmÞ ¼ u1 ; u2 ; ; uN ; f ðmÞ ¼ f1 ; f2 ; ; fN , and I be the identity matrix with an appropriate size. Then the
numerical scheme (3) can be written in the following matrix form
a Dx Dxa ðm1Þ þ Dxa f ðmÞ ; I þ AðmÞ uðmÞ ¼ u Dt Dt
ð5Þ
with ðmÞ T AðmÞ ¼ DðmÞ þ Ga þ D Ga ;
where
DðmÞ
¼
ð6Þ
ðmÞ ðmÞ diagðd;1 ; ; d;N Þ
2
ðaÞ g1
6 6 ðaÞ 6 g2 6 6 . 6 . 6 . 6 Ga ¼ 6 6 .. 6 . 6 6 6 ðaÞ 6 g N1 4 ðaÞ
gN
ðaÞ g0
0
ðaÞ
g0
g2
ðaÞ
g1
..
.
..
..
.
..
g1
ðaÞ
g N1
and
0
0
..
.
..
.
.
..
.
..
.
.
..
.
g1
ðaÞ
ðaÞ
ðaÞ ðaÞ
g2
0
3
7 7 0 7 7 .. 7 7 . 7 7 7 7 0 7 7 7 ðaÞ 7 g0 7 5
:
ð7Þ
ðaÞ
g1
NN
It is obvious that Ga is a Toeplitz matrix (see [8,9]). Therefore, it can be stored with N þ 1 entries [31]. Furthermore, the matrix–vector multiplication for the Toeplitz-like matrix AðmÞ in (6) can be obtained in OðN log NÞ operations by the FFT; see [8,22]. Denote
mN;M ¼
Dxa M ¼ ðxR xL Þa T 1 Dt ðN þ 1Þa
ð8Þ
which is related to the numbers of time steps and grid points. The linear system (5) can be written as
MðmÞ uðmÞ ¼ b
ðm1Þ
ð9Þ
;
where
MðmÞ ¼ mN;M I þ AðmÞ ¼ mN;M I þ DþðmÞ Ga þ DðmÞ GTa
ð10Þ
and ðm1Þ
b
ðmÞ ¼ mN;M uðm1Þ þ Dtf :
The coefficient matrix MðmÞ in (10) is a strictly diagonally dominant M-matrix (see [31]) and therefore it is nonsingular. 3. PCGNR method with circulant preconditioner The conjugate gradient (CG) method is a popular and effective iterative method for solving symmetric positive definite linear systems. Furthermore, for the nonsymmetric linear system (9), the CG method can also be applied to its equivalent normalized linear system
ðMðmÞ ÞT MðmÞ uðmÞ ¼ ðMðmÞ ÞT b
ðm1Þ
;
which is called the CGNR method [2]. The algorithm of the CGNR method is given in Algorithm 1 below.
718
S.-L. Lei, H.-W. Sun / Journal of Computational Physics 242 (2013) 715–725
ðm1Þ
Algorithm 1. CGNR for MðmÞ uðmÞ ¼ b ðmÞ
1.
Given initial guess u0
2.
r0 ¼ b
3. 4.
ðm1Þ
ðmÞ
MðmÞ u0
ðmÞ T
z0 ¼ ðM Þ r 0 ; p0 ¼ z0 For i ¼ 0; 1; . . ., until convergence, Do
5.
wi ¼ MðmÞ pi
6. 7. 8.
ai ¼ jjzi jj22 =jjwi jj22 ðmÞ ðmÞ uiþ1 ¼ ui þ ai pi r iþ1 ¼ ri ai wi
9.
ziþ1 ¼ ðMðmÞ ÞT riþ1
10. 11. 12.
bi ¼ jjziþ1 jj22 =jjzi jj22 piþ1 ¼ ziþ1 þ bi pi EndDo
In [32], the CGNR method is employed to solve (5). However, the drawback of the CGNR method is at its slow convergence rate, due to the fact that the condition number of the coefficient matrix ðMðmÞ ÞT MðmÞ is usually large. In fact, the numerical results given in Section 5 show that the CGNR method really has a poor convergence rate for ill-conditioned systems. Therefore, preconditioner technique could be exploited to accelerate the convergence rate of the CGNR method. When the CG method, with a matrix P ðmÞ , is employed to solve the following normalized preconditioned system
h iT h i h iT ðm1Þ ðP ðmÞ Þ1 MðmÞ ðP ðmÞ Þ1 MðmÞ uðmÞ ¼ ðP ðmÞ Þ1 MðmÞ ðP ðmÞ Þ1 b ; ðm1Þ
it is equivalent to the original system MðmÞ uðmÞ ¼ b mathematically. This method is called the PCGNR method with a preconditioner P ðmÞ . The algorithm of the PCGNR method is given in Algorithm 2 below, which is the direct application of the ðm1Þ Algorithm 2 to the system ðP ðmÞ Þ1 MðmÞ uðmÞ ¼ ðP ðmÞ Þ1 b . ðm1Þ
Algorithm 2. PCGNR for MðmÞ uðmÞ ¼ b
with preconditioner P ðmÞ .
ðmÞ
1. Given initial guess u0 ðmÞ 1
ðm1Þ
ðmÞ
2. r0 ¼ ðP Þ ðb MðmÞ u0 Þ h iT 3. z0 ¼ ðP ðmÞ Þ1 MðmÞ r0 , p0 ¼ z0 4. For i ¼ 0; 1; . . ., until convergence, Do 5.
wi ¼ MðmÞ pi
6.
wi ¼ ðP ðmÞ Þ1 wi
7.
ai ¼ kzi k22 =kwi k22 ðmÞ ðmÞ uiþ1 ¼ ui þ ai pi r iþ1 ¼ ri ai wi
8. 9. 10.
ziþ1 ¼ ðP ðmÞ ÞT riþ1
11.
ziþ1 ¼ ðMðmÞ ÞT ziþ1
12. bi ¼ kziþ1 k22 =kzi k22 13. piþ1 ¼ ziþ1 þ bi pi 14. EndDo
Now we propose a circulant preconditioner, which is generated from the Strang’s circulant preconditioner [9], in the PCGNR method to solve (9). The Strang’s circulant matrix sðBÞ ¼ ½sjk 06j;k
8 bj ; > > > < 0; sj ¼ > bjN ; > > : sjþN ;
0 6 j < N=2; j ¼ N=2 if N is even; N=2 < j < N; 0 < j < N:
Recall the formula (10) that
719
S.-L. Lei, H.-W. Sun / Journal of Computational Physics 242 (2013) 715–725
MðmÞ ¼ mN;M I þ AðmÞ ¼ mN;M I þ DþðmÞ Ga þ DðmÞ GTa ; where ðmÞ
ðmÞ
DðmÞ ¼ diagðd;1 ; ; d;N Þ; and Ga is the Toeplitz matrix given in (7). Then our circulant preconditioner is defined as
ðmÞ sðGa Þ þ d ðmÞ sðGT Þ; S ðmÞ ¼ mN;M I þ d þ a
ð11Þ
where N X ðmÞ ðmÞ ¼ 1 d d : N i¼1 ;i
More precisely, the first columns of both sðGa Þ and sðGTa Þ are given by
2
ða Þ g1
2
3
ðaÞ
g1
3
6 ðaÞ 7 6 g0 7 6 . 7 7 6 6 . 7 6 0 7 6 . 7 7 6 6 ðaÞ 7 6 .. 7 6 g Nþ1 7 6 . 7 6 b 2 c7 7 6 7 6 7; 6 0 7 and 6 0 7 6 6 . 7 6 ðaÞ 7 6 . 7 6 g Nþ1 7 6 . 7 c b 6 2 7 7 6 7 6 4 0 5 6 .. 7 4 . 5 ða Þ g0 ðaÞ g2 respectively. We know that a circulant matrix C can be diagonalized by the Fourier matrix F [9]; i.e.,
C ¼ F KF; where the entries of F are given by
1 F j;k ¼ pffiffiffiffi e2pijk=N ; N
0 6 j; k 6 N 1;
with the imaginary unit i, and K is a diagonal matrix holding the eigenvalues of C. a is the complex conjugate of Ka . The following lemma is the key to prove a F, where K Let sðGa Þ ¼ F Ka F. Then sðGTa Þ ¼ F K the invertibility of S ðmÞ in (11). Lemma 1. All eigenvalues of sðGa Þ and sðGTa Þ fall inside the open disc
fz 2 C : jz aj < ag: ðaÞ
Proof. All the Gershgorin disc of the circulant matrices sðGa Þ and sðGTa Þ are centered at g 1 ¼ a with radius ðaÞ
rN ¼ g 0 þ
bNþ1 c 2
X
ðaÞ
gk <
1 X
ðaÞ
ðaÞ
g k ¼ g 1 ¼ a;
k¼0;k–1
k¼2
ðaÞ
by the properties of the sequence fg k g; see Proposition 1. h Remark 1. It is worth to note that: 1. The real parts of all eigenvalues of sðGa Þ and sðGTa Þ are strictly positive for all N. 2. The absolute values of all eigenvalues of sðGa Þ and sðGTa Þ are bounded above by 2a < 4 for all N. Ka þ d ðmÞ K a . Then S ðmÞ is invertible Decompose the circulant matrix S ðmÞ ¼ F Ks F with the diagonal matrix Ks ¼ mN;M I þ d þ if all diagonal entries of Ks are nonzero. ðmÞ
Lemma 2. Let 1 < a < 2. The preconditioner S ðmÞ in (11) is invertible and
kðS ðmÞ Þ1 k2 6
1
mN;M
:
ð12Þ
720
S.-L. Lei, H.-W. Sun / Journal of Computational Physics 242 (2013) 715–725
Proof. By Part 1 of Remark 1, we have Reð½Ka k;k Þ > 0. Noting that
mN;M > 0 and dðmÞ P 0, we have
ðmÞ Reð½Ka Þ þ d ðmÞ Reð½K a Þ P mN;M > 0 j½Ks k;k j P Reð½Ks k;k Þ ¼ mN;M þ d þ k;k k;k for each k ¼ 1; . . . ; N. Therefore, S ðmÞ is invertible. Furthermore, we have
kðS ðmÞ Þ1 k2 ¼
1 1 : 6 min j½Ks k;k j mN;M
16k6N
h The computational cost per iteration of the PCGNR method with circulant preconditioners is OðN log NÞ complexity, as the main workload is the matrix–vector multiplication. Therefore, if the iterative method is of linear or superlinear convergence, the total cost for solving (9) at each time step by the proposed method will keep OðN log NÞ operations. Before to analyze the convergence of the proposed method, we prove the following theorem which gives a bound of kS ðmÞ k2 . This theorem will be exploited to prove the superlinear convergence rate of the PCGNR method in next section. ~ We have Lemma 3. Let 1 < a < 2 and 0 6 d ðx; tÞ 6 d.
~ kS ðmÞ k2 6 mN;M þ 8d:
ð13Þ
Proof. For any k ¼ 1; . . . ; N, by Part 2 of Remark 1,
½ Ka þ d ðmÞ ½K j½Ka j þ d ðmÞ j½K þ 4d ðmÞ 6 mN;M þ 8d: ~ a j 6 mN;M þ d a j < mN;M þ 4d j½Ks k;k j ¼ jmN;M þ d þ þ þ k;k k;k k;k k;k ðmÞ
ðmÞ
ðmÞ
Then
~ kS ðmÞ k2 ¼ max j½Ks k;k j < mN;M þ 8d: 16k6N
h
4. Spectrum of the preconditioned matrix In this section, we study the convergence rate of the PCGNR method with the proposed circulant preconditioner S ðmÞ in (11) for the linear system (9). For convenience of our investigation, we assume that the left and the right diffusion coefficient functions d ðx; tÞ of the FDE (1) are two nonnegative constants and M is properly chosen, depending on N, such that mN;M in and m ^ such that (8) is bounded away from 0; i.e., there exist two real numbers m
6 mN;M 6 m ^; 0
8N:
ð14Þ
By this assumption, for all i ¼ 1; . . . ; N, ðmÞ
dþ;i ¼ dþ P 0;
ðmÞ
d;i ¼ d P 0;
and dþ þ d – 0:
ð15Þ
Throughout this section, we add a subscript N to each matrix to denote the matrix size. Under the assumption in (15), the matrix AðmÞ defined in (6), MðmÞ defined in (10), and S ðmÞ defined in (11) are independent of m, and we therefore simply denote them as AN ; MN , and S N , respectively. Now the coefficient matrix MðmÞ in (10) becomes
MN ¼ mN;M I þ AN ;
ð16Þ
which is a nonsymmetric Toeplitz matrix, where
2
ðaÞ
ðaÞ
dþ g 1 þ d g 1
6 6 6 dþ g ð2aÞ þ d g ð0aÞ 6 6 ðaÞ AN ¼ dþ Ga þ d GTa ¼ 6 6 dþ g 3 6 6 .. 6 . 4 ðaÞ
dþ g N
ðaÞ
ðaÞ
ðaÞ
ðaÞ
dþ g 0 þ d g 2
ðaÞ
dþ g 1 þ d g 1 .. . .. .
.. . .. . .. .
...
dþ g 2 þ d g 0
ða Þ
3
ðaÞ
d g 3 .. . .. . .. .
d g N .. . .. .
ðaÞ
7 7 7 7 7 7 ½ajk : NN 7 7 7 ðaÞ ðaÞ dþ g 0 þ d g 2 7 5 ðaÞ
ðaÞ
dþ g 1 þ d g 1
To study the spectrum of the preconditioned matrix
1 T 1 S N MN S N MN ; we first introduce the generating function of the sequence of Toeplitz matrices fAN g1 N¼1 [8]:
ð17Þ
S.-L. Lei, H.-W. Sun / Journal of Computational Physics 242 (2013) 715–725
pðhÞ ¼
1 X
ak eikh ;
721
ð18Þ
k¼1
where ak is the k-th diagonal of AN . The generating function pðhÞ is in the Wiener class [8,9] if and only if 1 X
jak j < 1:
k¼1
For AN defined in (17), we have
pðhÞ ¼
1 X
ak eikh ¼
k¼1
We remark that
1 X
ðaÞ
g kþ1 ðdþ eikh þ d eikh Þ:
ð19Þ
k¼1
mN;M þ pðhÞ cannot be a generating function of MN in (16) since mN;M is dependent on N.
T 1 Lemma 4. Let p and q be the generating functions of fAN g1 N¼1 and fAN þ AN gN¼1 , respectively, we have
1. p is in the Wiener class; 2. q is real-valued and nonnegative.
Proof ðaÞ
1. By the properties of the sequence fg k g1 k¼0 given in (4), we have 1 X
jak j ¼ ðdþ þ
k¼1
1 X
ðaÞ d Þ jg kþ1 j k¼1
¼ ðdþ þ d Þ
ðaÞ 2g 1
1 X ðaÞ þ gk
!
¼ 2aðdþ þ d Þ < 1:
k¼0
Thus p is in the Wiener class. 2. From (4) and the definition of generating function in (18), we have
qðhÞ ¼
1 X
" # 1 X ðaÞ ða Þ ðaÞ ih ðaÞ ðak þ ak Þeikh ¼ ðdþ þ d Þ 2g 1 þ g 0 þ g 2 g kþ1 eikh þ eikh e þ eih þ
k¼1
k¼2
"
# 1 1 X X ðaÞ ðaÞ ðaÞ ðaÞ ðaÞ ¼ ðdþ þ d Þ 2g 1 þ 2 g 0 þ g 2 cos h þ 2 g kþ1 cos kh P 2ðdþ þ d Þ g kþ1 ¼ 0: k¼2
k¼1
h
From Part 2 of Lemma 4, and according to the Grenander–Szegö Theorem (see Theorem 1.12 in [7]), the matrix AN þ ATN is symmetric positive semidefinite. Furthermore, by noting that
MN MTN ¼ m2N;M I þ mN;M AN þ ATN þ AN ATN ; we conclude that the smallest singular value of the matrix MN is bounded below by
rmin ðMN Þ P mN;M :
mN;M ; i.e., ð20Þ
By Part 1 of Lemma 4, we have the following lemma. Lemma 5 (see [8]). If pðhÞ, the generating function of AN , is in the Wiener class, then for any > 0, there exist N0 and M0 > 0, such that for all N > N0,
AN sðAN Þ ¼ U N þ V N where
rank ðU N Þ 6 M0 and
kV N k2 < : Now we consider the spectrum of S 1 N MN I. By Lemma 5, we have 1 1 1 S 1 N MN I ¼ S N ðMN S N Þ ¼ S N U N S N V N :
722
S.-L. Lei, H.-W. Sun / Journal of Computational Physics 242 (2013) 715–725
0 Note that rankðS 1 N U N Þ 6 rankðU N Þ 6 M and
1 kS 1 N V N k2 6 kS N k2 kV N k2 <
6 ; mN;M m
by the assumption (14) and the inequality (12) in Lemma 2. Therefore, the following corollary is resulted. Corollary 1. If
mN;M satisfies the assumption (14), for any > 0, there exists N00 and M00 > 0 such that, for any N > N00 ,
e e S 1 N MN I ¼ U N þ V N ;
ð21Þ
e N Þ 6 M 00 and k V e N k < . where rankð U 2 Now we show that the spectrum of the normalized preconditioned matrix
1 T 1 S N MN S N MN is clustered around 1. Theorem 1. If mN;M satisfies the assumption (14), for any 0 < < 1, there exists N 00 and M 00 > 0 such that, for all N > N 00 , at most 2M00 eigenvalues of the matrix
1 T 1 S N MN S N MN I have absolute values larger than 3. Proof. By (21), we have
1 T 1 e N ÞT ðI þ U eN þ V e NÞ ¼ I þ U bN þ V b N; eN þ V S N MN S N MN ¼ ðI þ U where
bN ¼ U e T ðI þ U eN þ V e N Þ þ ðI þ V e T ÞU eN U N N and
eN þ V eT þ V eT V bN ¼ V e V N N N: b N Þ 6 2M00 and k V b N k < 3. Thus, we have Then by Corollary 1, we see that rankð U 2
1 T 1 b N; bN þ V S N MN S N MN I ¼ U b N and V b N are real and symmetric. By the well-known Weyl’s theorem, we conclude that at most 2M00 eigenwhere both U values of the matrix
1 T 1 S N MN S N MN I have absolute values larger than 3.
h
For the superlinear convergence of the PCGNR method, according to the Corollary 1.11 in [7], it is further required that the smallest singular value of the matrix S 1 N MN is uniformly bounded away from zero. In fact, using (13), (20), and the assumption (14), we have
rmin ðS 1 N MN Þ P
rmin ðMN Þ kS N k2
P
mN;M m^ P > 0; mN;M þ 8d~ m þ 8d~
~ ¼ maxfdþ ; d g. Therefore, the PCGNR method for solving (9) converges superlinearly. where d 5. Numerical results In this section, we solve the FDE problem (1) numerically by the implicit finite difference method given in Section 2. After ðm1Þ the finite difference discretization, at each time step, the nonsymmetric linear system MðmÞ uðmÞ ¼ b is solved by the CGNR method (Algorithm 1), the PCGNR method (Algorithm 2), and the multigrid method [22], respectively. Number of iterations required for convergence and CPU time of those methods are reported. The stopping criterion of those methods is
krðkÞ k2 < 107 ; kr ð0Þ k2 where rðkÞ is the residual vector of the linear system after k iterations, and the initial guess at each time step is chosen as the zero vector. In the following tables, ‘‘N’’ denotes the number of spatial grid points, ‘‘M’’ denotes the number of time steps,
723
S.-L. Lei, H.-W. Sun / Journal of Computational Physics 242 (2013) 715–725
‘‘Error’’ denotes the error between the true solution and the approximation under the infinity norm, ‘‘CPU (s)’’ denotes the total CPU time in seconds for solving the whole FDE problem, and ‘‘Iter’’ denotes the average number of iterations required for solving the FDE problem; i.e., M 1X IterðmÞ; M m¼1
Iter ¼
where IterðmÞ denotes the number of iterations required for solving (9). For the PCGNR method, besides the proposed circulant preconditioner S ðmÞ in (11), we also test the Chan’s preconditioner [10] which can be written as
cðGa Þ þ d ðmÞ cðGT Þ; cðMðmÞ Þ ¼ mN;M I þ d þ a ðmÞ
where cðBÞ denotes the Chan’s preconditioner for an arbitrary matrix B. More precisely, the first columns of the circulant matrices cðGa Þ and cðGTa Þ are given as
2
6 6 6 16 6 N6 6 6 4
3
ðaÞ
Ng 1
ða Þ
ðN 1Þg 2 .. . ðaÞ
2g N1 ðaÞ gN
ðaÞ
2
3
ðaÞ
Ng 1
7 7 6 7 6 ðN 1Þg ðaÞ þ g ðaÞ 7 N 7 0 7 6 7 7 6 ðaÞ 7 and 1 6 7; 2g N1 7 7 6 N 7 7 6 . 7 7 6 . 5 5 4 . ðaÞ
þ ðN 1Þg 0
ðN 1Þg 2
respectively. We remark that the superlinear convergence of the PCGNR method with the preconditioner cðMðmÞ Þ, as that with S ðmÞ in Section 4, can be analogously proven if the assumptions (14) and (15) hold. In the following tables, we denote ‘‘I ’’ as the CGNR method, ‘‘C’’ as the PCGNR method with the T. Chan’s circulant preconditioner cðMðmÞ Þ, ‘‘S’’ as the PCGNR method with the proposed circulant preconditioner S ðmÞ , and ‘‘MGM’’ as the multigrid method in [22]. In Table 2, ‘‘Error’’ denotes the error between the exact solution and the numerical solution obtained by the PCGNR method with the proposed preconditioner SðmÞ . All numerical experiments are run in MATLAB 7.12 (R2011a) on a Dell workstation Precision T7500 with the configuration: Intel (R) Xeon (R) CPU X5690 3.47 GHz and 12 GB RAM. Example 1. In this example, we solve the initial-boundary value problem of FDE (1) with source term f ðx; tÞ 0, for order of fractional derivatives a ¼ 1:2; 1:5 and 1:8. The spatial domain is ½xL ; xR ¼ ½0; 2 and the time interval is ½0; T ¼ ½0; 1. The initial condition uðx; 0Þ is the following Gaussian pulse
uðx; 0Þ ¼ exp
! ðx xc Þ2 ; 2r 2
xc ¼ 1:2;
r ¼ 0:08;
and the diffusion coefficients are
dþ ðx; tÞ 0:6;
and d ðx; tÞ 0:5:
Table 1 Comparisons for solving Example 1 by the CGNR method, the PCGNR method with preconditioners cðMðmÞ Þ and S ðmÞ , and the multigrid method for a ¼ 1:2; 1:5, and 1:8. Here Dt 2Dxa .
a 1.2
1.5
1.8
Nþ1
M
I
C
MGM
S
Iter
CPU (s)
Iter
CPU (s)
Iter
CPU (s)
Iter
CPU (s)
6
2
32
37.6
0.10
6.0
0.05
5.8
0.04
10.1
0.30
27
74
34.4
0.36
6.0
0.12
5.3
0.11
8.8
0.65
28
169
31.4
0.96
5.0
0.28
5.0
0.28
7.3
1.56
29
388
28.5
4.08
5.0
1.29
5.0
1.27
6.3
3.93
210
891
25.7
15.22
5.0
4.49
5.0
4.44
5.5
10.54
2
91
40.9
0.30
6.0
0.11
5.6
0.10
7.0
0.52
27
256
39.2
1.40
6.0
0.43
5.2
0.37
6.1
1.59
28
724
35.8
4.70
5.4
1.28
5.0
1.20
5.6
5.25
29
2048
32.3
24.21
5.0
6.74
5.0
6.73
5.1
16.91
210
5793
29.0
110.65
5.0
28.92
5.0
28.93
4.8
60.20
6
2
256
42.6
0.85
7.0
0.31
5.8
0.26
6.1
1.25
27
891
41.0
5.01
6.0
1.45
5.5
1.36
6.0
5.50
28
3104
36.3
20.59
6.0
5.92
5.3
5.38
5.0
19.93
29
10809
31.8
126.14
5.2
36.55
5.1
35.88
5.0
88.35
210
37641
27.5
685.80
5.0
187.06
5.0
187.07
4.1
339.35
6
724
S.-L. Lei, H.-W. Sun / Journal of Computational Physics 242 (2013) 715–725
Table 2 Comparisons for solving Example 2 by the CGNR method, the PCGNR method with preconditioners cðMðmÞ Þ and S ðmÞ , and the multigrid method for a ¼ 1:2; 1:5 and 1:8.
a 1:2
1:5
1:8
Nþ1
M
Error
I
C
MGM
S
Iter
CPU (s)
Iter
CPU (s)
Iter
CPU (s)
Iter
CPU (s) 0.19
2
5
2
3.1501e2
33.8
0.08
8.0
0.05
8.0
0.05
8.0
27
26
1.5983e2
65.5
0.53
8.0
0.12
8.0
0.12
7.1
0.41
28
27
8.0488e3
82.0
1.70
8.0
0.29
7.0
0.26
7.1
1.03
29
28
4.0377e3
89.9
7.92
8.0
1.18
8.0
1.18
7.0
2.60
210
29
2.0214e3
96.0
30.90
7.0
3.32
8.0
3.73
7.0
6.91
26
25
2.2529e2
46.6
0.11
10.0
0.06
8.0
0.05
8.0
0.19
27
26
1.1164e2
111.6
0.87
10.4
0.16
9.0
0.14
7.0
0.41
28
27
5.5563e3
264.6
5.34
10.9
0.38
9.3
0.33
7.0
1.03
29
28
2.7721e3
568.3
48.74
9.9
1.44
9.9
1.44
7.0
2.62
210
29
1.3838e3
903.7
288.51
11.0
5.01
10.0
4.52
6.0
5.93
2
5
2
1.7434e2
70.6
0.16
16.0
0.08
13.0
0.07
7.0
0.16
27
26
8.3524e3
202.0
1.58
18.0
0.26
14.0
0.21
7.0
0.40
28
27
4.0838e3
587.2
11.68
18.9
0.63
14.0
0.48
7.0
1.01
29
28
2.0186e3
1703.0
146.05
21.0
2.89
14.0
1.97
7.0
2.60
210
29
1.0035e3
4900.0
1541.66
20.0
8.66
13.0
5.78
6.0
5.92
6
6
For satisfying (14), we choose M ¼ bð2N þ 2Þa =2c such that Dt 2Dxa . Thus, the superlinear convergence rate by the PCGNR method with the proposed circulant preconditioner is guaranteed as both (14) and (15) hold. We remark that the linear convergence rate by the multigrid method in [22] has not been proven theoretically. The numerical results of Example 1 are shown in Table 1. It shows that both the average number of iterations and the CPU time by the PCGNR method with circulant preconditioners are much less than those by the CGNR method. It also shows that the CPU times by the PCGNR methods are less than that by the multigrid method.
Example 2. In this example, we study the case for which the diffusion coefficients are variable. We solve the FDE problem (1) of order a with the following data. The spatial domain is ½xL ; xR ¼ ½0; 2 and the time interval is ½0; T ¼ ½0; 1. The left and right diffusion coefficients are
dþ ðx; tÞ ¼ Cð3 aÞxa
and d ðx; tÞ ¼ Cð3 aÞð2 xÞa ;
respectively. The source term is
1 3 3 f ðx; tÞ ¼ 32et x2 þ ð2 xÞ2 ð8 þ x2 Þ ½x3 þ ð2 xÞ3 þ ½x4 þ ð2 xÞ4 8 3a ð4 aÞð3 aÞ and the initial condition is
uðx; 0Þ ¼ 4x2 ð2 xÞ2 : The exact solution of this problem is
uðx; tÞ ¼ 4et x2 ð2 xÞ2 for any a 2 ð1; 2Þ. With the exact solution, we calculate the exact error of the numerical solution. In the following tests, we consider a ¼ 1:2; 1:5, and 1:8. We remark that the same example in [19,22,31] is only considered for a ¼ 1:8. In the numerical test, we naturally choose Dt ¼ Dx. Table 2 reports the history for numerically solving Example 2. Note that the average number of iterations for the CGNR method is increasing as N is increasing rapidly, while the average numbers of iterations by the PCGNR methods almost keep constant, and therefore their CPU times are much less than that by the CGNR method. We remark that both diffusion coefficients are not constant and mN;M tends to 0. Therefore, both assumptions (14) and (15) are not satisfied and hence the superlinear convergence rate is not theoretically guaranteed. However, from the point of view on CPU time, the performance of the proposed circulant preconditioner still works well and is the best among all the methods.
6. Concluding remarks In this paper, we have employed the PCGNR method with a circulant preconditioner to solve the discretized linear systems of (1). Theoretically, we prove the superlinear convergence rate of the proposed method under the conditions that the
S.-L. Lei, H.-W. Sun / Journal of Computational Physics 242 (2013) 715–725
725
a
diffusion coefficients d ðx; tÞ are constant and the ratio DDxt is bounded away from zero. Numerical tests have shown the efficiency of the PCGNR method. Moreover, for variable diffusion coefficients in the FDE (1), the proposed algorithm still works very well. In our future consideration, more effective preconditioners will be developed to solve the two or more dimensional FDE problem [30] and the steady-state fractional partial differential equations. Acknowledgement The authors thank Dr. Hong-kui Pang for his helpful discussions. We are also grateful to the anonymous referees for the useful suggestions and comments that improved the presentation of this paper. References [1] J. Bai, X. Feng, Fractional-order anisotropic diffusion for image denoising, IEEE Trans. Image Proc. 16 (2007) 2492–2502. [2] R. Barrett, M. Berry, T. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo, C. Romine, H. van der Vorst, Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods, SIAM, Philadelphia, 1994. [3] D. Benson, S.W. Wheatcraft, M.M. Meerschaert, Application of a fractional advection–dispersion equation, Water Resour. Res. 36 (2000) 1403–1413. [4] D. Benson, S.W. Wheatcraft, M.M. Meerschaert, The fractional-order governing equation of Lévy motion, Water Resour. Res. 36 (2000) 1413–1423. [5] B. Beumer, M. Kovács, M.M. Meerschaert, Numerical solutions for fractional reaction–diffusion equations, Comput. Math. Appl. 55 (2008) 2212–2226. [6] B.A. Carreras, V.E. Lynch, G.M. Zaslavsky, Anomalous diffusion and exit time distribution of particle tracers in plasma turbulence models, Phys. Plasma 8 (2001) 5096–5103. [7] R. Chan, X. Jin, An Introduction to Iterative Toeplitz Solvers, SIAM, Philadelphia, 2007. [8] R. Chan, M. Ng, Conjugate gradient methods for Toeplitz systems, SIAM Rev. 38 (1996) 427–482. [9] R. Chan, G. Strang, Toeplitz equations by conjugate gradients with circulant preconditioner, SIAM J. Sci. Stat. Comput. 10 (1989) 104–119. [10] T. Chan, An optimal circulant preconditioner for Toeplitz systems, SIAM J. Sci. Stat. Comput. 9 (1988) 766–771. [11] W. Deng, Finite element method for the space and time fractional Fokker–Planck equation, SIAM J. Numer. Anal. 47 (2008) 204–226. [12] V.J. Ervin, N. Heuer, J.P. Roop, Numerical approximation of a time dependent, nonlinear, space-fractional diffusion equation, SIAM J. Numer. Anal. 45 (2007) 572–591. [13] X. Jin, Preconditioning Techniques for Teoplitz Systems, Higher Education Press, Beijing, 2010. [14] T.A.M. Langlands, B.I. Henry, The accuracy and stability of an implicit solution method for the fractional diffusion equation, J. Comput. Phys. 205 (2005) 719–736. [15] F. Liu, V. Anh, I. Turner, Numerical solution of the space fractional Fokker–Planck equation, J. Comput. Appl. Math. 166 (2004) 209–219. [16] R.L. Magin, Fractional Calculus in Bioengineering, Begell House Publishers, 2006. [17] M.M. Meerschaert, H.P. Scheffler, C. Tadjeran, Finite difference methods for two-dimensional fractional dispersion equation, J. Comput. Phys. 211 (2006) 249–261. [18] M.M. Meerschaert, C. Tadjeran, Finite difference approximations for fractional advection–dispersion flow equations, J. Comput. Appl. Math. 172 (2004) 65–77. [19] M.M. Meerschaert, C. Tadjeran, Finite difference approximations for two-sided space-fractional partial differential equations, Appl. Numer. Math. 56 (2006) 80–90. [20] D.A. Murio, Implicit finite difference approximation for time fractional diffusion equations, Comput. Math. Appl. 56 (2008) 1138–1145. [21] M. Ng, Iterative Methods for Toeplitz Systems, Numerical Mathematics and Scientific Computation, Oxford University Press, New York, 2004. [22] H. Pang, H. Sun, Multigrid method for fractional diffusion equations, J. Comput. Phys. 231 (2012) 693–703. [23] I. Podlubny, Fractional Differential Equations, Academic Press, New York, 1999. [24] M. Raberto, E. Scalas, F. Mainardi, Waiting-times and returns in high-frequency financial data: an empirical study, Physica A 314 (2002) 749–755. [25] M.F. Shlesinger, B.J. West, J. Klafter, Lévy dynamics of enhanced diffusion: application to turbulence, Phys. Rev. Lett. 58 (1987) 1100–1103. [26] I.M. Sokolov, J. Klafter, A. Blumen, Fractional kinetics, Phys. Today Nov. (2002) 28–53. [27] E. Sousa, Finite difference approximates for a fractional advection diffusion problem, J. Comput. Phys. 228 (2009) 4038–4054. [28] L. Su, W. Wang, Z. Yang, Finite difference approximations for the fractional advection–diffusion equation, Phys. Lett. A 373 (2009) 4405–4408. [29] C. Tadjeran, M.M. Meerschaert, H.P. Scheffler, A second-order accurate numerical approximation for the fractional diffusion equation, J. Comput. Phys. 213 (2006) 205–213. 2 [30] H. Wang, K. Wang, An OðNlog NÞ alternating-direction finite difference method for two-dimensional fractional diffusion equations, J. Comput. Phys. 230 (2011) 7830–7839. 2 [31] H. Wang, K. Wang, T. Sircar, A direct OðNlog NÞ finite difference method for fractional diffusion equations, J. Comput. Phys. 229 (2010) 8095–8104. [32] K. Wang, H. Wang, A fast characteristic finite difference method for fractional advection–diffusion equations, Adv. Water Resour. 34 (2011) 810–816. [33] G.M. Zaslavsky, D. Stevens, H. Weitzner, Self-similar transport in incomplete chaos, Phys. Rev. E 48 (1993) 1683–1694.