INTERNATIONAL JOURNAL OF

c 2012 Institute for Scientific

NUMERICAL ANALYSIS AND MODELING Computing and Information
Volume 9, Number 3, Pages 658–666

A FAST SECOND-ORDER FINITE DIFFERENCE METHOD FOR

SPACE-FRACTIONAL DIFFUSION EQUATIONS

TREENA S. BASU AND HONG WANG

Abstract. Fractional diffusion equations provide an adequate and accurate

description of transport processes that exhibit anomalous diffusion that cannot
be modeled accurately by classical second-order diffusion equations. However,
numerical discretizations of fractional diffusion equations yield full coefficient
matrices, which require a computational operation of O(N 3 ) per time step and
a memory of O(N 2 ) for a problem of size N . In this paper we develop a fast
second-order finite difference method for space-fractional diffusion equations,
which only requires memory of O(N ) and computational work of O(N log2 N ).
Numerical experiments show the utility of the method.

Key Words. circulant and Toeplitz matrix, fast direct solver, fast finite dif-
ference methods, fractional diffusion equations

1. Introduction
Fractional diffusion equations model phenomena exhibiting anomalous diffusion
that cannot be modeled accurately by classical second-order diffusion equations. For
instance, in contaminant transport in groundwater flow the solutes moving through
aquifers do not generally follow a Fickian, second-order partial differential equation
because of large deviations from the stochastic process of Brownian motion. Instead,
a governing equation with a fractional-order anomalous diffusion provides a more
adequate and accurate description of the movement of the solutes [4].
Compared to the classical second-order diffusion equations, the fractional dif-
fusion equations have salient features which introduce new difficulties. From a
computational point of view, fractional differential operators are nonlocal and so
raise subtle stability issues on the corresponding numerical approximations. Numer-
ical methods for space-fractional diffusion equations yield full coefficient matrices,
which require a computational operation of O(N 3 ) per time step and a memory
of O(N 2 ) for a problem of size N . This is in contrast to numerical methods for
second-order diffusion equations which usually generate banded coefficient matri-
ces of O(N ) nonzero entries and can be solved by fast solution methods such as
multigrid methods, domain decomposition methods, and wavelet methods in O(N )
(or O(N log N )) operations per time step with O(N ) memory requirement.
Meerschaert and Tadjeran [7, 8] showed that a direct truncation of the Grünwald-
Letnikov form of fractional derivative, even though discretized implicitly in time,
leads to unstable discretizations. They proposed a shifted Grünwald discretiza-
tion to approximate the fractional diffusion equation and proved the unconditional
stability and convergence of the corresponding finite difference scheme. Numeri-
cal experiments showed that these methods generate satisfactory numerical results.

Received by the editors February 27, 2011 and, in revised form, August 1, 2011.
2000 Mathematics Subject Classification. 35R35, 49J40, 60G40.
658
 A FAST FINITE DIFFERENCE METHOD FOR FRACTIONAL DIFFUSION EQUATION 659

However, the shifted Grünwald discretization is only first-order accurate in space.

Tadjeran et al [11] developed a Crank-Nicolson scheme which is second-order ac-
curate in time. They recovered second-order spatial accuracy by a Richardson
extrapolation. However, these methods still generate full coefficient matrices and
so require storage of O(N 2 ) and computational work of O(N 3 ) per time step.
In this paper we develop a fast second-order finite difference method for two-sided
space-fractional diffusion equations. The method has a significantly reduced mem-
ory requirement of O(N ) and computational work of O(N log2 N ) per time step.
The method is an extension of the fast solution method developed in [13] and can
also viewed as an extension of the superfast method [1, 2, 3], which was a direct so-
lution method of O(N log2 N ) operations for a symmetric positive-definite Toeplitz
system. The rest of the paper is organized as follows. In Section 2 we present the
fractional diffusion equation and its Crank-Nicolson finite difference approximation.
In Section 3 we develop the fast second-order finite difference method. In Section 4
we carry out numerical experiments to compare the performance of the fast finite
difference method with the Crank-Nicolson finite difference method developed and
analyzed in [11].

2. Fractional diffusion equations and its finite difference approximation

We consider the following initial-boundary value problem of a two-sided space-
fractional diffusion equation with an anomalous diffusion of order 1 < α < 2
∂u(x, t) ∂ α u(x, t) ∂ α u(x, t)
− d+ (x, t) − d − (x, t) = f (x, t),
∂t ∂+ xα ∂− xα
(1) xL < x < xR , 0 < t ≤ T,
u(xL , t) = 0, u(xR , t) = 0, 0 ≤ t ≤ T,
u(x, 0) = u0 (x), xL ≤ x ≤ xR .
∂ α u(x,t) ∂ α u(x,t)
The left-sided (+) and the right-sided (−) fractional derivatives ∂+ x α and ∂− xα
of equation (1) are defined in the Grünwald-Letnikov form
b(x−xL )/hc
∂ α u(x, t) 1 X (α)
= lim+ α gk u(x − kh, t),
∂+ xα h→0 h
k=0
(2) b(xR −x)/hc
α
∂ u(x, t) 1 X (α)
= lim α gk u(x + kh, t)
∂− xα h→0+ h
k=0

(α)
where bxc represents the floor of x and the Grünwald weights gk are defined as
gk = (−1)k αk where α
(α)


k represents fractional binomial coefficients. We note
(α)
that the Grünwald weights gk have the recursive relation
 
(α) (α) α+1 (α)
(3) g0 = 1, gk = 1 − gk−1 for k ≥ 1.
k
(α)
Moreover, for 1 < α < 2 the coefficients gk satisfy the following properties:

(α) (α) (α) (α)
 g0 = 1,

 g1 = −α < 0, 1 ≥ g2 ≥ g3 ≥ · · · ≥ 0,
∞ m
(4) X (α)
X (α)


 gk = 0, gk < 0 (m ≥ 1).
k=0 k=0
660 TREENA S. BASU AND HONG WANG

We now focus on the development of a fast numerical method for problem (1).
Let N and M be positive integers and h = (xR −xL )/N and ∆t = T /M be the sizes
of spatial grid and time step, respectively. The spatial and temporal partitions are
defined as xi = xL + ih for i = 0, 1, . . . , N and tm = m∆t for m = 0, 1, . . . , M .
Let um m m m m m
i = u(xi , t ), d+,i = d+ (xi , t ), d−,i = d− (xi , t ), and fi
m
= f (xi , tm ).
We discretize the first-order time derivative in (1) by a standard first-order time
difference quotient, but the for the discretization of the fractional spatial derivative
we use the shifted Grünwald approximations
i+1
∂ α u(xi , tm ) 1 X (α) m
= α gk ui−k+1 + a1 h + O(h2 ),
∂+ xα h
k=0
(5)
α m −i+1
NX
∂ u(xi , t ) 1 (α)
= α gk um 2
i+k−1 + b1 h + O(h )
∂− xα h
k=0

where a1 and b1 do not depend on the grid size h. The Crank-Nicolson finite
difference scheme is formulated as follows
m+1 i+1 i+1
um+1
i − um
i 1  d+,i X (α) m+1 dm
+,i
X (α) m

− gk ui−k+1 − gk ui−k+1
∆t 2 hα hα
k=0 k=0
−i+1
m+1 NX −i+1
NX
(6) 1  d−,i (α) dm
−,i (α) m

− gk um+1
i+k−1 − gk ui+k−1
2 hα hα
k=0 k=0
1
= fim+1 + fim ,


2
which was proved to be unconditionally stable and convergent [11]. Numerical
experiments show that this scheme generates very satisfactory numerical approx-
T m T
, Am =
 
imations. Let um = um m m
1 , u2 , . . . , uN −1 , f = f1m , f2m , . . . , fN
m
−1
 m N −1
ai,j i,j=1 , and I be the identity matrix of order N − 1. Then the numerical scheme
(6) can be expressed in the following matrix form
   
∆t ∆t ∆t m
(7) I + α Am+1 um+1 = I − α Am um + (f + f m+1 ).
2h 2h 2

Here the entries of matrix Am+1 are given by

 m+1 m+1

(α)
 − d+,i + d−,i g1 , 


j = i,
 − dm+1 g (α) + dm+1 g (α) ,

j = i − 1,

  +,i 2

 −,i 0

m+1 m+1 (α) m+1 (α)
(8) ai,j = − d+,i g0 + d −,i g2 , j = i + 1,


 −dm+1 g
 (α)


 +,i i−j+1 , j < i − 1,

 m+1 (α)
−d−,i gj−i+1 , j > i + 1,

It is clear that am+1

i,j ≤ 0 for all i 6= j and that the coefficient matrix I+(∆t/(2hα ))Am+1
is a nonsingular, strictly diagonally dominant M-matrix. (8) implies that the Crank-
Nicolson scheme has a full coefficient matrix, which has a memory requirement of
O(N 2 ) and and computational work of O(N 3 ) per time step.
To develop a fast solution method, we carefully explore the structure of the
coefficient matrices. We conclude from (8) that the stiffness matrices Am+1 and
 A FAST FINITE DIFFERENCE METHOD FOR FRACTIONAL DIFFUSION EQUATION 661

Am can be decomposed as follows

Am+1 = −diag dm+1 AL − diag dm+1



+ − AR ,
(9)
Am = −diag dm m



+ AL − diag d− AR .

Here diag dm+1 , diag dm+1 , diag dm m

+ − + , and diag d− are diagonal matrices of
order N − 1 with their ith entries dm+1 m+1 m m
+,i , d−,i , d+,i , and d−,i for i = 1, 2, . . . , N − 1.
The matrices AL and AR are matrices of order N − 1 and are defined by

(α) (α)
 
g1 g0 0 ... 0 0

 g (α) (α) (α) .. .. 
 2 g1 g0 . . 0 

 . .. .. .. 
 . (α) (α)
. .
 . g2 g1 .

AL =  . ,

 . .. .. .. ..
 . . . . . 0


.. .. ..
 
 (α) (α) (α)
 gN −2 . . . g1 g0


(α) (α) (α) (α)
gN −1 gN −2 ... ... g2 g1
 (α) (α) (α) (α)

g1 g2 ... ... gN −2 gN −1

 g (α) (α) (α) .. (α)

 0 g1 g2 ... . gN −2 


(α) (α) .. .. .. 
 0

g0 g1 . . .

AR =  . .

 . .. .. .. .. ..
 . . . . . .


..
 
(α) (α)
.
 
 0 ... 0 g1 g2 
(α) (α)
0 0 ... 0 g0 g1

Instead of storing two full matrices Am+1 and Am which have 2(N − 1)2 parame-
T
ters we need only store the 5N −4 parameters, dm+1 = dm+1 m+1 m+1

+ +,1 , d+,2 , . . . , d+,N −1 ,
m+1
 m+1 m+1 m+1
T m  m m m
T m  m
d− = d−,1 , d−,2 , . . . , d−,N −1 , d+ = d+,1 , d+,2 , . . . , d+,N −1 , d− = d−,1 ,
T  (α) (α) (α) T
dm m
−,2 , . . . , d−,N −1 , and g(α) = g0 , g1 , . . . , gN −1 . In particular, the fractional
binomial coefficient vector g(α) depends only on the size of the spatial partition and
the order of the anomalous diffusion but is independent of time or space. So it can
be preprocessed and stored in advance.

3. A fast second-order finite difference method

In this section we extend the idea of our previous work [13] to develop a fast
second-order finite difference method for problem (1). The development essentially
consists of three steps: (i) We apply an operator-splitting technique to split the
stiffness matrix Am+1 as the sum of a banded matrix Am+1k and a remaining matrix
Am+1
O = Am+1 − Am+1k . (ii) We then move Am+1
O to the right-hand side of the
numerical scheme and approximate the unknown solution by an extrapolation in
time to retain a second-order accuracy. (iii) We carefully explore the structure of
the coefficient matrices on the right-hand side of the numerical scheme to develop
a fast algorithm to evaluate the right-hand side.
Let Am+1
k contains the 2k + 1 diagonals of Am+1 and zero entries elsewhere, and
m+1 m+1
AO = A − Am+1
k contains the remaining nonzero entries of Am+1 . Then we
662 TREENA S. BASU AND HONG WANG

split the stiffness matrix Am+1 as

(10) Am+1 = Am+1
k + Am+1
O .
It was shown in [13] that if the bandwidth k = log N , then
m+1
− Am+1

A
k
(11) m+1 ∞
= O(log−α N ) → 0 as N → ∞.
A
∞
In other words, as the number of unknowns N increases, the relative weight of
the banded matrix Am+1 k over the full matrix Am+1 increases too, even if the
bandwidth k increases only as log N in contrast to the linear increase of the width
of the full matrix Am+1 . Moreover, such choice of the bandwidth k guarantees the
computational work of inverting I + (∆t/(2hα ))Am+1
k is O(N log2 N ).
We hence split the scheme (7) as follows
   
∆t m+1 ∆t m ∆t
I + α Ak um+1
= I − αA um − α Am+1 um+1
2h 2h 2h O
(12)
∆t m
f + f m+1 .


+
2
The issue that remains is how to approximate the um+1 in the second term on the
right-hand side of the scheme (12). To enhance the accuracy of the approximation,
we evaluate the um+1 on the right-hand side of (12) by a quadratic extrapolation
in time which yields the following approximation ûm+1 of um+1
(13) ûm+1 ≈ 3um − 3um−1 + um−2 , m≥2
Now it remains to approximate the solution u at the first and second time step
i.e., u1 and u2 . We use the fast implicit Euler method in [13] to compute u1 and
u2 . Even though these approximations have only second-order local truncation er-
ror, the numerical analysis theory of time-dependent problems tells us that this
treatment still retains the same second-order global truncation error since these
approximations are used only for the first two steps. Thus after substituting ûm+1
of (13) for um+1 in the second term on the right-hand side of the scheme (12) and
evaluating u1 and u2 via the fast implicit Euler method, we obtain the following
fast second-order finite difference scheme (F2FD)
   
∆t m+1 ∆t m ∆t m+1
I + α Ak um+1 = I − α A − 3 α AO um
2h 2h 2h
∆t ∆t
+3 α Am+1
O um−1 − α Am+1 um−2
2h 2h O
∆t m
f + fm+1 , m ≥ 2,


(14) +
2
   
∆t 2 ∆t ∆t
I+ A u2 = I − 2 α A2O u1 + α A2O u0 + ∆tf 2
hα k h h
   
∆t ∆t
I + α A1k u1 = I − α A1O u0 + ∆tf 1 .
h h
Finally, we turn to issue (iii), i.e., the fast evaluation of the right-hand side of
the scheme. We recall from the matrix decomposition (9) that a fast evaluation
of the right-hand side of the scheme (14) boils down to the fast evaluation of the
matrix-vector multiplication of AL u and AR u. Note that both AL and AR are
 A FAST FINITE DIFFERENCE METHOD FOR FRACTIONAL DIFFUSION EQUATION 663

Toeplitz matrices, which can be embedded into (2N − 2)-by-(2N − 2) circulant

matrix C2N −2,L and C2N −2,R , respectively. A Toeplitz matrix is a matrix in which
each descending diagonal from left to right is constant, while a circulant matrix is
a matrix in which each row vector is rotated one element to the right relative to
the preceding row vecor. In general, an n × n Toeplitz matrix Tn is completely
determined by a sequence of 2n − 1 numbers {ti }n−1 i=1−n such that the (i, j)-entry of
the matrix Tn (i, j) = tj−i for i, j = 1, . . . , n, i.e.,

t0 t1 t2 . . . tn−2 tn−1
 
 t−1 t0 t1 . . . tn−3 tn−2 
 
 .. .. 
 t−2 t−1 t0 . . tn−3 
Tn =  . .
 
 .. .. .. .. .. .. 
 . . . . . 

 .. .. 
 t
2−n t3−n . . t0 t1 
t1−n t2−n t3−n . . . t−1 t0

while an n × n circulant matrix Cn is completely determined by a sequence of n

numbers {ci }n−1 i=0 such that the (i, j)-entry of the matrix Cn (i, j) = c(j−i) mod n for
i, j = 1, . . . , n.

c0 c1 c2 . . . cn−2 cn−1
 
 cn−1 c0 c1 . . . cn−3 cn−2 
 
 .. .. 
 cn−2 cn−1 c0 . . cn−3 
Cn =  . .
 
 .. .. .. .. .. .. 
 . . . . . 

 .. .. 
 c
2 c3 . . c0 c1 
c1 c2 c3 ... cn−1 c0

It is known that a circulant matrix C2N −2 can be decomposed as follows [5, 6]

(15) C2N −2 = F−1

2N −2 diag (F2N−2 c) F2N−2

where c is the first column vector of C2N −2 and F2N −2 is the (2N − 2)-by-(2N − 2)
discrete Fourier transform matrix. For any given N − 1 dimensional vector u,
 T
u2N −2 = u, 0 . (15) shows that C2N −2 u2N −2 can be computed efficiently via the
fast Fourier transform in O((2N ) log(2N )) = O(N log N ) operations. Thus, the
matrix-vector products AL u and AR u can be computed in O(N log N ) operations
and so are Am+1 u and Am u.

4. Numerical experiment
In this section we carry out numerical experiments to investigate the performance
of the fast second-order finite difference method. We consider the fractional diffu-
sion equation (1) with an anomalous diffusion of order α = 1.8 and the left-sided
and right-sided diffusion coefficients

(16) d+ (x, t) = 1.32Γ(1.2)x1.8 , d− (x, t) = 1.32Γ(1.2)(2 − x)1.8 .

664 TREENA S. BASU AND HONG WANG

The spatial domain is [xL , xR ] = [0, 1], the time interval is [0, T ] = [0, 1]. The
source term and the initial condition are given by
h
f (x, t) = −16e1−t x2 (1 − x)2 + 2.64(x2 + (1 − x)2 )

i
−13.2 x3 + (1 − x)3 + 12 x4 + (1 − x)4 ,


(17)

u0 (x) = 16ex2 (1 − x)2 .

The true solution to the fractional diffusion equation (1) is given by [8]
(18) u(x, t) = 16e1−t x2 (1 − x)2 .
Tadjeran et al. proved that Crank-Nicolson scheme (7) is unconditionally stable
and is convergent with the accuracy of O((∆t)2 + h) [11]. They further used the
Richardson extrapolation to recover the second-order spatial convergence, which
involves finding the numerical solution uh on a coarse grid h and then finding the
numerical solution uh/2 on a fine grid h/2 and then computing the extrapolated
solution on the coarse spatial gird h by uh = 2uh/2 − uh . It was shown that the
extrapolated solution has a second-order accuracy in space and time O((∆t)2 +h)2 ).

N =M kuM M
CN − u kL∞ kuM M
CN − u kL2 CPU(seconds)
−3
26 7.98656 × 10 6.80845 × 10−3 1.95312 × 101
−3 −3
27 3.79054 × 10 3.23428 × 10 1.50718 × 102
−3 −3
28 1.84519 × 10 1.57456 × 10 1.23123 × 103
−4 −4
29 9.10152 × 10 7.76630 × 10 9.86920 × 103
−4 −4
210 4.51975 × 10 3.85650 × 10 7.65370 × 104
c=0.58, r = 1.03 c = 0.49, r = 1.03
kuM M
F 2F D − u kL∞ kuM M
F 2F D − u kL2
−3
26
8.00233 × 10 6.82157 × 10−3 5.04687
−3
27
3.79848 × 10 3.24100 × 10−3 2.28593 × 101
28 1.84957 × 10−3 1.57831 × 10−3 1.02171 × 102
9 −4 −4
2 9.13258 × 10 7.79294 × 10 4.13281 × 102
10 −4 −4
2 4.53797 × 10 3.87218 × 10 1.83257 × 103
c=0.58, r = 1.03 c = 0.49, r = 1.03
Table 1. Comparison of the fast second order finite difference
(F2FD) method with the Crank-Nicolson method (CN) with
Gaussian elimination

In the numerical experiment, we solve the problem by the Crank-Nicolson method

(6) (CN), the fast second-order finite difference method (F2FD), the Crank-Nicolson
method with extrapolation (CNE), and the fast second-order finite difference method
with extrapolation (F2FDE) and denote their respective solutions by um m
CN , uCN E ,
m m m m m m
uF 2F D , and uF 2F DE . Let uh be the numerical solution uCN , uCN E , uF 2F D , or
um
F 2F DE at time step t
m
and um = u(x, tm ) be the true solution to problem (1). In
Table 1 we choose h = ∆t and present the errors kuM M M M
CN −u kL∞ , kuF 2F D −u kL∞ ,
M M M M
and kuF 2F D − u kL2 , kuCN − u kL2 for different mesh sizes. We then use a linear
regression to fit the convergence rate r and the associated constant c in the L2 and
L∞ norm
(19) kuM M r
h − u kLp ≤ ch , p = 2, ∞.
 A FAST FINITE DIFFERENCE METHOD FOR FRACTIONAL DIFFUSION EQUATION 665

N =M kuM M
CN E − u kL∞ kuM M
CN E − u kL2 CPU(seconds)
−4 −4
26 4.827901 × 10 3.46162 × 10 1.57656 × 102
−4 −5
27 1.21258 × 10 8.68653 × 10 1.37660 × 103
−5 −5
28 3.03853 × 10 2.17572 × 10 9.91226 × 103
−6 −6
29 7.60523 × 10 5.44443 × 10 7.992601 × 104
c=1.94, r = 1.99 c= 1.40, r = 1.99
kuM M
F 2F DE − u kL∞ kuM M
F 2F DE − u kL2
−4
26
4.82081 × 10 3.45777 × 10−4 2.58125 × 101
7 −4 −5
2 1.20027 × 10 8.60525 × 10 1.14453 × 102
8 −5 −5
2 2.79488 × 10 2.01151 × 10 4.69015 × 102
9 −6 −6
2 6.88048 × 10 4.95583 × 10 2.47515 × 103
c=2.44, r = 2.04 c =1.73, r = 2.04
Table 2. Comparison of the fast second order finite difference
with extrapolation (F2FDE) with the Crank-Nicolson method with
extrapolation (CNE)

In Table 1 we also present the corresponding CPU times of the Crank-Nicolson

method (6) (CN) and the fast second-order finite difference method (F2FD) from
N = M = 64 to N = M = 1024. Each time we reduce the size of the spatial meshes
by half, the total number of unknowns per time step is doubled. Consequently, the
required memory increases 4 times and computational cost increases 8 times. If
the time step size is reduced by half too, then the overall consumed CPU time for
solving the finite difference method is expected to increase 2 × 23 = 16 times as
predicted by the leading order behavior. We can in fact observe that the CPU time
increases around 10 times or so each time we reduce the spatial mesh and time step
size by half. The CPU time increases not exactly 16 times due to the overhead effect
from the computations of other terms. In contrast, the computational work of the
fast second-order finite difference method is O(N log2 N ). Each time we refine the
size of the spatial meshes by half the required memory increases twice and when
reducing the time step by half,the CPU time increases by 4 times. Our numerical
experiments seem to coincide with this analysis.
Finally, in Table 2 we present the numerical results of the Crank-Nicolson method
with Richardson extrapolation (CNE) side by side with fast second-order finite
difference method with Richardson extrapolation (F2FDE), which indeed shows a
second-order convergence rate in space and time.
In summary, the numerical experiments in this section show significant reduction
of computational time, which coincides with the analysis. For example, with 1024
computational nodes, the new scheme developed in this paper has about 40 times
of CPU reduction than the standard scheme. This is in addition to the significan-
t reduction in the storage. In short, these results indeed show the utility of the
method. Finally, even though we did not present a theoretical proof of the stability
of the proposed numerical scheme, the numerical results presented in Table 1 indi-
cate that the new scheme has the same stability constraint as the standard finite
difference scheme which was proven to be unconditionally stable [7].

Acknowledgments
This work was supported in part by the National Science Foundation under
Grant No. EAR-0934747.
666 TREENA S. BASU AND HONG WANG

References
1. G. Ammar and P. Gader, A variant of the Gohberg-Semencul formula involving circulant
matrices, SIAM Journal on Matrix Analysis and Applications, 12 (1991), 534–540
2. G.S. Ammar and W.B. Gragg, The generalized Schur algorithm for the superfast solution of
Toeplitz systems, Rational Approximation and its Applications in Mathematics and Physics,
Lecture Notes in Mathematics, 1237 (1987), 315–330
3. G.S. Ammar and W.B. Gragg, Superfast solution of real positive definite Toeplitz systems,
SIAM Journal on Matrix Anal. Appl., 9 (1988), 61–76
4. D. Benson, S.W. Wheatcraft, and M.M. Meerschaert, The fractional-order governing equation
of Lévy motion, Water Resour. Res. 36 (2000), 1413-1423.
5. P.J. Davis, Circulant Matrices, Wiley-Intersciences, New York, 1979.
6. R.M. Gray, Toeplitz and Circulant Matrices: A Review, Foundations and Trends in Commu-
nications and Information Theory, Vol 2, Issue 3, pp 155-239, 2006.
7. M.M. Meerschaert and C. Tadjeran, Finite difference approximations for fractional advection-
dispersion flow equations, J. Comput. Appl. Math., 172 (2004) 65-77.
8. M.M. Meerschaert and C. Tadjeran, Finite difference approximations for two-sided space-
fractional partial differential equations, Appl. Numer. Math., 56 (2006) 80-90.
9. K.B. Oldham and J. Spanier, The Fractional Calculus, Academic Press, New York, 1974.
10. I. Podlubny, Fractional Differential Equations, Academic Press, New York, 1999.
11. C. Tadjeran, M.M. Meerschaert, and H.P. Scheffler, A second-order accurate numerical ap-
proximation for the fractional diffusion equation, J. Comput. Phys., 213 (2006) 205-213.
12. R.S. Varga, Matrix iterative Analysis, Second Edition, Springer-Verlag, Berlin Heideberg,
2000.
13. H. Wang, K. Wang, and T. Sircar, A direct O(N log 2 N ) finite difference method for fractional
diffusion equations, Journal of Computational Physics, 229 (2010), 8095-8104

Department of Mathematics, University of South Carolina, Columbia South Carolina 29208,

USA
E-mail: sircar@[Link] hwang@[Link]