# Computational Challenge of Fractional Differential Equations and the Potential Solutions: A Survey

Authors: Chunye Gong; Weimin Bao; Guojian Tang; Yuewen Jiang; Jie Liu
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

We present a survey of fractional differential equations and in particular of the computational cost for their numerical solutions from the view of computer science. The computational complexities of time fractional, space fractional, and space-time fractional equations areO(N2M), O(NM2), and O(NM(M + N)) compared with O(MN) for the classical partial differential equations with finite difference methods, where M, N are the number of space grid points and time steps. The potential solutions for this challenge include, but are not limited to, parallel computing, memory access optimization (fractional precomputing operator), short memory principle, fast Fourier transform (FFT) based solutions, alternating direction implicit method, multigrid method, and preconditioner technology. The relationships of these solutions for both space fractional derivative and time fractional derivative are discussed. The authors pointed out that the technologies of parallel computing should be regarded as a basic method to overcome this challenge, and some attention should be paid to the fractional killer applications, high performance iteration methods, high order schemes, and Monte Carlo methods. Since the computation of fractional equations with high dimension and variable order is even heavier, the researchers from the area of mathematics and computer science have opportunity to invent cornerstones in the area of fractional calculus.


---

## Body

## 1. Introduction

The idea of fractional is natural. If∂u/∂x and ∂2u/∂x2 exist, ∂1.5u/∂x1.5 maybe exists too. Fractional equations can be used to describe some physical phenomena more accurately than the classical integer order differential equation [1]. Fractional differential equations (FDEs) provide a powerful instrument for the description of memory and hereditary properties of different substances. The fractional diffusion equations play an important role in dynamical systems of semiconductor research, hydrogeology, bioinformatics, finance [2], and other research areas [3–6]. Rajeev and Kushwaha [7] presented a mathematical model describing the time fractional anomalous diffusion process of a generalized Stefan problem which is a limiting case of a shoreline problem. Space fractional advection-diffusion equations arise when velocity variations are heavy tailed and describe particle motion that account for variation in the flow field over the entire system [8]. FDEs may be divided into two fundamental types: time fractional differential equations and space fractional differential equations. For the fractional ordinary equations and fractional order control systems are also studied [9, 10]. The stability of fractional order control systems attracts many attentions [11, 12]. For example, Laguerre continued fraction expansion of the Tustin fractional discrete-time operator was investigated by Maione [13].Some analytical methods were proposed for fractional differential equation [14, 15]. Saha Ray [16] presented the analytical solutions of the space fractional diffusion equations by two-step Adomian decomposition method. Momani and Odibat [17] gave a comparison between the homotopy perturbation method and the variational iteration method for linear fractional partial differential equations. By using initial conditions, the explicit solutions of the equations have been presented in the closed form and then their solutions have been represented graphically. Because most of fractional problems cannot be solved analytically, more and more works focus on their numerical solutions. There are many numerical solutions proposed for fractional equations [18], such as finite difference method (FDM) [18]. FDM is intuitive to understand and easy to learn for inexperienced researcher from the areas rather than mathematics. So this survey focuses on FDM for fractional equations.For the numerical solutions of different differential equations, the area of mathematics pays much attention to approximating the equation more accurately and faster (accuracy and speed). The area of computer science mainly focuses on the runtime (speed) and code reuse (software). The numerical methods (mathematic area) have eternal value and may exist along with the existence of human culture [19–21]. The implementations (computer area) are closer to the human society and the real applications but vary quickly along with the computer architecture [22–28]. The computational cost of the numerical solutions for fractional equations is much heavier than that for the traditional integer order equations. In the near future, the fractional problems with high dimension, long time iterations, and huge grid points will need to be solved. These problems are real challenge for today’s computer technologies and algorithms.


## 2. Fractional Differential Equations

### 2.1. Origins

In 1695, L’Hopital wrote to Leibniz asking him about a particular notation he had used in his publications for thenth-derivative of the linear function f(x)=x, Dnx/Dxn [29, 30]. L’Hopital’s posed the question to Leibniz: what would the result be if n=1/2? Leibniz’s response was “An apparent paradox, from which one day useful consequences will be drawn.” In these words fractional calculus was born [31]. Later, Fourier, Euler, and Laplace dabbled with fractional calculus [32].


### 2.2. A Short Summary of FDE

There are mainly three kinds of FDEs:(1)
time fractional [33, 34],(2)
space fractional [35],(3)
space-time fractional [36].One of the main challenges in fractional differential equation is the nonlocality of the fractional operator. In the case of a time fractional derivative, one needs to store all the history, whereas in the case of a space fractional derivative, one needs to deal with almost-dense matrices.The partial differential equations (PDEs) mainly have three categories: parabolic, hyperbolic, and elliptic. There are corresponding FDEs to deal with the fractional conditions. There are various FDEs listed below:(1)
the parabolic, widely studied fractional diffusion equation [37, 38],(2)
the hyperbolic Telegraph equation [39],(3)
the elliptic fractional Laplace equation [40],(4)
fractional Black-Scholes equation in computational finance [41],(5)
fractional wave equation for sound, light, and water waves [42],(6)
fractional Fokker-Planck equations describing the time evolution of the probability density function of the velocity of a particle [43],(7)
fractional Euler equations and fractional Navier-Stokes equations [44] for fluid dynamics [45],(8)
fractional kinetic equations for motion of objects [46],(9)
the perfect, fractional Maxwell’s equations for electrodynamics [47],(10)
fractional Boltzmann equations for particle transport [48].We can suppose thatif there is a PDE, there will be a FDE. The reason is that replacing the integer derivative with fractional derivative and solving the FDE with different numerical methods is not a hard job.The numerical schemes [49] for these FDEs are listed below:(1)
finite difference method (FDM) [18, 50, 51],(2)
finite element method [52–54],(3)
finite volume method [55–57],(4)
Adomian decomposition method [58],(5)
Fourier transform method [59],(6)
spectral method [60],(7)
meshless method [61–63],(8)
exponential difference method [64, 65],(9)
Monte Carlo method [66].The Monte Carlo (MC) method, which uses repeated random sampling to obtain numerical results, belongs to undetermined computational methods.


### 2.3. A Hot Topic in Recent Years

Machado et al. [29] collected 20 special issues on fractional calculus. In 2011 and 2012, there are about 7 special issues on this topic. There are more special issues in the nearby four years (2011–2014) than those of 1999 to 2010 [29]. And there are additional 20 special issues in the years of 2013 and 2014. So the researches on fractional calculus can be regarded as a collective revelry.


## 2.1. Origins

In 1695, L’Hopital wrote to Leibniz asking him about a particular notation he had used in his publications for thenth-derivative of the linear function f(x)=x, Dnx/Dxn [29, 30]. L’Hopital’s posed the question to Leibniz: what would the result be if n=1/2? Leibniz’s response was “An apparent paradox, from which one day useful consequences will be drawn.” In these words fractional calculus was born [31]. Later, Fourier, Euler, and Laplace dabbled with fractional calculus [32].


## 2.2. A Short Summary of FDE

There are mainly three kinds of FDEs:(1)
time fractional [33, 34],(2)
space fractional [35],(3)
space-time fractional [36].One of the main challenges in fractional differential equation is the nonlocality of the fractional operator. In the case of a time fractional derivative, one needs to store all the history, whereas in the case of a space fractional derivative, one needs to deal with almost-dense matrices.The partial differential equations (PDEs) mainly have three categories: parabolic, hyperbolic, and elliptic. There are corresponding FDEs to deal with the fractional conditions. There are various FDEs listed below:(1)
the parabolic, widely studied fractional diffusion equation [37, 38],(2)
the hyperbolic Telegraph equation [39],(3)
the elliptic fractional Laplace equation [40],(4)
fractional Black-Scholes equation in computational finance [41],(5)
fractional wave equation for sound, light, and water waves [42],(6)
fractional Fokker-Planck equations describing the time evolution of the probability density function of the velocity of a particle [43],(7)
fractional Euler equations and fractional Navier-Stokes equations [44] for fluid dynamics [45],(8)
fractional kinetic equations for motion of objects [46],(9)
the perfect, fractional Maxwell’s equations for electrodynamics [47],(10)
fractional Boltzmann equations for particle transport [48].We can suppose thatif there is a PDE, there will be a FDE. The reason is that replacing the integer derivative with fractional derivative and solving the FDE with different numerical methods is not a hard job.The numerical schemes [49] for these FDEs are listed below:(1)
finite difference method (FDM) [18, 50, 51],(2)
finite element method [52–54],(3)
finite volume method [55–57],(4)
Adomian decomposition method [58],(5)
Fourier transform method [59],(6)
spectral method [60],(7)
meshless method [61–63],(8)
exponential difference method [64, 65],(9)
Monte Carlo method [66].The Monte Carlo (MC) method, which uses repeated random sampling to obtain numerical results, belongs to undetermined computational methods.


## 2.3. A Hot Topic in Recent Years

Machado et al. [29] collected 20 special issues on fractional calculus. In 2011 and 2012, there are about 7 special issues on this topic. There are more special issues in the nearby four years (2011–2014) than those of 1999 to 2010 [29]. And there are additional 20 special issues in the years of 2013 and 2014. So the researches on fractional calculus can be regarded as a collective revelry.


## 3. Computational Challenge

### 3.1. A Classical Partial Differential Equation (PDE)

#### 3.1.1. Heat Equation/Diffusion Equation

The heat equation is a parabolic PDE which describes the distribution of heat (or variation in temperature) in a given region over time, shown in (1). This equation is also known as the diffusion equation:(1)∂ux,t∂t=dx∂2ux,t∂x2+fx,t,ux,0=ϕx,x∈xL,xR,uxL,t=φt,t∈0,T,uxR,t=ψt,t∈0,T,where the d(x) is the diffusion constant.


#### 3.1.2. Numerical Method

Definetn=nτ, xi=ih for 0≤n≤N, 0≤i≤M, where M and N are positive integer and τ=T/N, h=(xR-xL)/M are time step size and space step size, respectively. Define uin, fin, di as the numerical approximation to u(xi,tn), f(xi,tn), and d(xi). Using a forward difference at time tn and a second-order central difference for the space derivative at position xi, we get the explicit finite difference approximation for the one-dimensional heat equation:(2)uin+1-uinτ=diui+1n-2uin+ui-1nh2+fin.Theuin+1 can be obtained by this way:(3)uin+1=a1ui+1n+a2uin+a3ui-1n+τfin,where a1=a3=-diτ/h2 and  a2=1-2diτ/h2.If we use the backward difference at timetn+1 and a second-order central difference for the space derivative at position xi we get the explicit finite difference approximation for (1):(4)uin+1-uinτ=diui+1n+1-2uin+1+ui-1n+1h2+fin+1.Theuin+1 can be obtained by this way:(5)a1ui+1n+1+a2uin+1+a3ui-1n+1=uin+τfin,where a1=a3=diτ/h2 and a2=1+2diτ/h2.The classical PDE is used to compare the fractional equations. For convenience, there are several hypotheses for this paper.(1)
We just focus on the explicit FDM of (3), because these kinds of comparisons between fractional equations and classical equations are straightforward. The comparison between implicit FDE of (5) and implicit schemes of fractional equations is more complex [67].(2)
The computational cost of diffusion coefficientdi and source fin (or fin+1) is ignored, because they only need compute once. And they are different with different equations and make the comparison complicated.


#### 3.1.3. Computational Cost

Each grid point of time steptn+1 needs 4 multiplications and 3 additions. There are M-1 grid points in each time step. So each time step needs 7(M-1) arithmetic logic operations. There are about N time steps. So the total computational cost is about 7N(M-1). The computational cost will vary linearly along the number of time steps and grid points.


### 3.2. Time Fractional Diffusion Equation

#### 3.2.1. Numerical Method

Liu et al. [63] developed an implicit radial basis function (RBF) meshless approach for time fractional diffusion equations and found that the presented meshless formulation is very effective for modeling and simulation of fractional differential equations. Murillo and Yuste developed an explicit difference method for solving fractional diffusion with Caputo form [68]:(6)∂αux,t∂αt=∂2ux,t∂x20<α<1,ux,0=gx,x∈0,xR,u0,t=uxR,t=0on a finite domain 0<x<xR and 0≤t≤T.The explicit finite difference approximation for (6) is described as follows [69]:(7)τ-α∑i=0n+1‍wiujn+1-i-uj0=1h2uj+1n-2ujn+uj-1n,where wi is the normalized Grünwald weight defined with(8)wi=-1iαi.Defineμ=τα/h2, bn=∑i=0n-1wi, and Un=(u1n,u2n,…,uM-1n)T; we can get(9)Un+1=bn+1U0-∑i=1n+1‍wiUn+1-i+μAUn,where matrix A is a tridiagonal matrix, defined by(10)AM-1×M-1=-211-21••••11-2.


#### 3.2.2. Computational Cost

In order to getUn+1, the right-sided computation of (9) should be performed. There are mainly one tridiagonal matrix-vector multiplication, many constant-vector multiplications, and many vector-vector additions in the right-sided computation.(1)
The tridiagonal matrix-vector multiplication isAUn and a new vector Un+1=AUn is got.(2)
The constant-vector multiplications areV′=bn+1U0, V′′=μV1, Vi=wiUn+1-i.(3)
The vector-vector additions areUn+1=Un+1+V′+V′′-∑i=1n+1Vi.There are about5(M-1) operations for tridiagonal matrix-vector multiplication. For time step tn+1 needs (n+8)(M-1) arithmetic logic operations with n=1→N. So the total computational cost for (9) is about(11)8NM-1+1+2+⋯+N-1M-1=8NM-1+NN-1M-12=NN2+7.5M-1.The computational cost varies linearly along the number of grid points but squares with the number of time steps.


### 3.3. Space Fractional Diffusion Equation

#### 3.3.1. Numerical Method

A classical numerical scheme for the space fractional diffusion equation is the second-order fractional Crank-Nicolson method proposed by Tadjeran et al. [38, 70], where the Richardson extrapolation technique is used to the first order shift Grünwald formula for space fractional derivative. Tadjeran et al. [70] presented a practical numerical method in time and space to solve a class of initial-boundary value with variable coefficients on a finite domain for (12)∂ux,t∂t=dx∂αux,t∂xα+qx,ton a finite domain xL<x<xR with 1<α<2 and 0≤t≤T. The case of 1<α<2 models a super-diffusive flow in which a cloud of diffusing particles spreads at a faster rate than the classical diffusion model predictions [71, 72].The explicit finite difference approximation for (12) is(13)uin+1-uinτ=dihα∑k=0i+1‍gα,kui-k+1n+qin,where uin=u(xi,tn), di=d(xi), and qin=q(xi,tn). The gα,k is the normalized Grünwald weight [70, 73]:(14)gα,k=Γk-αΓ-αΓk+1.These normalized weights only depend on the orderα and the index k. The resulting equation can be explicitly solved for uin+1 to give(15)uin+1=uin+diτhα∑k=0i+1‍gα,kui-k+1n+τqin.


#### 3.3.2. Computational Cost

Theith grid point of time step tn+1 needs i+6 multiplication, addition, and division. There are M-1 grid points in each time step. So each time step needs (M-1)6+(1+2+⋯+M-1)=(M-1)((M/2)+6) arithmetic logic operations. There are about N time steps. So the total computational cost is about N(M-1)((M/2)+6).


### 3.4. Riesz Space Fractional Diffusion Equation

#### 3.4.1. Numerical Method

Shen et al. [74] investigated the Green function and a discrete random walk model for Riesz fractional advection-dispersion equation on infinite domain with an initial condition and also presented implicit and explicit finite differences to this problem on a finite domain. Çelik and Duman [75] applied the Crank-Nicolson method to a fractional diffusion equation which has the Riesz fractional derivative and obtained that the method is unconditionally stable and convergent. The Riesz space fractional reaction-diffusion equation [76, 77] is(16)∂ux,t∂t=-ux,t+xD0αux,t+qx,t,ux,0=gx,x∈xL,xR,uxL,t=uxR,t=0with 1<α≤2 and 0≤t≤T. Both u(x,t) and g(x) are real valued and sufficiently well-behaved function. Dx0αu(x,t) is the Riesz space fractional derivative.With adopting an Euler approximation in time, the explicit difference approximation can be got:(17)uin+1-uinτ=-uin-1chα×∑k=0i+1‍gkui+1-kn+∑k=0M-i+1‍gkui-1+kn+fin,where g0=1, gk=(-1)kαα-1⋯α-k+1/k!, k=1,2,3,… is the normalized Grünwald weight, c=1/2cos⁡(απ/2), and fin=q(xi,tn).Equation (17) results in a linear system of equations(18)Un+1=AUn+Qn,where Un=(u1n,u2n,…,uM-1n)T,  Qn=(τf1n,τf2n,…,τfM-1n)T with qin=τfin (1≤i<M-1), and A=(aij)(M-1)×(M-1) is a matrix of coefficient. A is defined by(19)ai,j=-τchαgi+1-j,for1≤j<i-1,-τchαg0+g2,forj=i-1,1-τ-2τchαg1,forj=i,-τchαg0+g2,forj=i+1,-τchαgj-i+1,fori+1<j<M-1.


#### 3.4.2. Computational Cost

From (18), the M-1 results are produced by AUn for each time step. The inner product of vectors v1 and v2 with size m has m multiplications and m-1 additions. AUn for each time step needs (M-1)(2M-3) operations. Assuming τfin is pre-performed, each time step needs (M-1)(2M-2) arithmetic logic operations. There are about N time steps. So the total computational cost is about N(M-1)(2M-2). The computational cost will vary linearly along the number of time steps but square with the number of grid points. From analytical view, the computation of (18) is about four times heavier than that of (15).


### 3.5. Space-Time Riesz-Caputo Fractional Convection-Diffusion Equation

The fractional advection-diffusion (dispersion) equation has been applied to many problems. Fractional advection-dispersion equations are used in groundwater hydrology to model the transport of passive tracers carried out by fluid flow in a porous medium [78]. Shen et al. [79] presented an explicit difference approximation and an implicit difference approximation for the space-time Riesz-Caputo fractional advection-diffusion equation with initial and boundary conditions in a finite domain. They proved that the implicit difference approximation is unconditionally stable and convergent, but the explicit difference approximation is conditionally stable and convergent. Liu et al. [80] proposed an implicit difference method and an explicit difference method to solve the space-time fractional advection-diffusion equation and discussed the stability and convergence of the method.


#### 3.5.1. Numerical Method

We consider the following space-time Riesz-Caputo fractional advection-diffusion equation (STRCFADE) studied by Shen et al. [79]. This equation is obtained by replacing the space-derivative in the advection-diffusion equation with a generalized derivative of order β1, β2 with 0<β1<1, 1<β2<2 and time-derivative with a generalized derivative of order α with 0<α<1. We consider(20)D0tαux,t=B1∂β1ux,t∂xβ1+B2∂β2ux,t∂xβ2+fx,t,ux,0=ϕx,x∈0,xR,u0,t=uxR,t=0,t∈0,Ton a finite domain 0≤x≤xR and 0≤t≤T. The coefficients B1 and B2 are both positive constants and represent the diffusion (dispersion) coefficient and the average fluid velocity. The ∂β1u(x,t)/∂|x|β1 and ∂β2u(x,t)/∂|x|β2 are Riesz space fractional derivatives of order β1 and β2, respectively. Then we can get the explicit finite differential approximation for (20) [79]:(21)uin+1=∑j=0n-1‍bj-bj+1uin-j+bnui0+B1μ1∑k=-iN-i‍ωkβ1ui+knuin+1=+B2μ2∑k=-iN-i‍ωkβ2ui+kn+μ0φin,i=1,…,N-1,u0n+1=uNn+1=0,ui0=ϕi,where μ0=ταΓ(2-α), μ1=μ0/hβ1, and μ2=μ0/hβ2. More information can be referred to in [79].


#### 3.5.2. Computational Cost

From Section3.4.2, we know that the computational cost of B1μ1(∑k=-iN-iωkβ1ui+kn)+B2μ2(∑k=-iN-iωkβ2ui+kn)+μ0φin of N time steps is about N(M-1)(2M-2) operations. From Section 3.2.2, we know that the computational cost of ∑j=0n-1(bj-bj+1)uin-j+bnui0 of N time steps with n ranging from 1 to N is about N(N/2)(M-1). So the total computational cost is about N(M-1)(2M-2+(N/2)). The computational cost varies quadratically with the number of time steps or the number of grid points.For fixedM=4097, the comparison between the computational costs of the numerical solution among classical PDE, time fractional, space fractional, Riesz space fractional, and space-time Riesz-Caputo fractional equations is shown in Figure 1. For fixed N=2048, the computational cost is shown in Figure 2.Figure 1
Comparison of the computational cost with different grid points for fixedM=4097 among classical, time fractional, space fractional, Riesz space fractional, and space-time Riesz-Caputo fractional equations.Figure 2
Comparison of the computational cost with different time step for fixedN=2048 among classical, time fractional, space fractional, Riesz space fractional, and space-time Riesz-Caputo fractional equations.


### 3.6. More Challenges

There are more computational challenges listed below:(1)
high dimensional problems [55, 81],(2)
implicit schemes [67, 82],(3)
high order schemes [33, 83],(4)
variable order problems [84],(5)
huge memory space requirement.The high dimensional problems are more computation expensive [55, 81, 85, 86]. For example, the two-dimensional time fractional diffusion equation (2D-TFDE) [67, 87] is(22)Dtα0ux,y,t=ax,y,t∂2ux,y,t∂x2Dtα0ux,y,t=+ax,y,t∂2ux,y,t∂y2+fx,y,t,ux,y,0=ϕx,y,x,y∈Ω,ux,y,t∂Ω=0,t∈0,T,where Ω={x,y∣0≤x≤L1, 0≤y≤L2, a(x,y,t)>0, and b(x,y,t)>0}.The approximating scheme is [87](23)ui,jn+1-ui,jn+∑s=1n‍bsui,jn+1-s-ui,jn-s=μ1Γ2-αai,jnui+1,jn-2ui,jn+ui-1,jn+μ2Γ2-αbi,jnui,j+1n-2ui,jn+ui,j-1n+ταΓ2-αfi,jn,where μ1=τα/hx2 and μ2=τα/hy2. The hx and hy are the step size along X and Y directions. The computational complexity is about O(N2M2), which is much bigger than the computational complexity of one-dimensional problems when the number of grid points is bigger enough.The direct Gauss elimination for implicit scheme of FDE is not convenient. Solving unsteady FDE with implicit schemes relies on the iteration method at each time step. The variable order problems [84, 88] have complex coefficients, which need more arithmetic logic operations. The high order schemes [33, 38, 70] need more arithmetic logic operations too.The memory usage belongs to the computing resources. So the huge memory requirement of FDE is a kind of computational challenge in the broader sense. This is especially true for time fractional problems. Ignoring the memory usage of the coefficients and source terms, the one-dimensional equation (9) needs 8(M-1)N bytes of memory space and the two-dimensional equation (23) needs 8M2N bytes of memory space. For three-dimensional problems, the memory usage is 8M3N at least. It needs 15.625 PB (1 PB = 10245 bytes) memory space with M=10240, N=2048 for three-dimensional problems. Maybe only the most powerful supercomputer [89] can satisfy the huge memory requirement.


## 3.1. A Classical Partial Differential Equation (PDE)

### 3.1.1. Heat Equation/Diffusion Equation

The heat equation is a parabolic PDE which describes the distribution of heat (or variation in temperature) in a given region over time, shown in (1). This equation is also known as the diffusion equation:(1)∂ux,t∂t=dx∂2ux,t∂x2+fx,t,ux,0=ϕx,x∈xL,xR,uxL,t=φt,t∈0,T,uxR,t=ψt,t∈0,T,where the d(x) is the diffusion constant.


### 3.1.2. Numerical Method

Definetn=nτ, xi=ih for 0≤n≤N, 0≤i≤M, where M and N are positive integer and τ=T/N, h=(xR-xL)/M are time step size and space step size, respectively. Define uin, fin, di as the numerical approximation to u(xi,tn), f(xi,tn), and d(xi). Using a forward difference at time tn and a second-order central difference for the space derivative at position xi, we get the explicit finite difference approximation for the one-dimensional heat equation:(2)uin+1-uinτ=diui+1n-2uin+ui-1nh2+fin.Theuin+1 can be obtained by this way:(3)uin+1=a1ui+1n+a2uin+a3ui-1n+τfin,where a1=a3=-diτ/h2 and  a2=1-2diτ/h2.If we use the backward difference at timetn+1 and a second-order central difference for the space derivative at position xi we get the explicit finite difference approximation for (1):(4)uin+1-uinτ=diui+1n+1-2uin+1+ui-1n+1h2+fin+1.Theuin+1 can be obtained by this way:(5)a1ui+1n+1+a2uin+1+a3ui-1n+1=uin+τfin,where a1=a3=diτ/h2 and a2=1+2diτ/h2.The classical PDE is used to compare the fractional equations. For convenience, there are several hypotheses for this paper.(1)
We just focus on the explicit FDM of (3), because these kinds of comparisons between fractional equations and classical equations are straightforward. The comparison between implicit FDE of (5) and implicit schemes of fractional equations is more complex [67].(2)
The computational cost of diffusion coefficientdi and source fin (or fin+1) is ignored, because they only need compute once. And they are different with different equations and make the comparison complicated.


### 3.1.3. Computational Cost

Each grid point of time steptn+1 needs 4 multiplications and 3 additions. There are M-1 grid points in each time step. So each time step needs 7(M-1) arithmetic logic operations. There are about N time steps. So the total computational cost is about 7N(M-1). The computational cost will vary linearly along the number of time steps and grid points.


## 3.1.1. Heat Equation/Diffusion Equation

The heat equation is a parabolic PDE which describes the distribution of heat (or variation in temperature) in a given region over time, shown in (1). This equation is also known as the diffusion equation:(1)∂ux,t∂t=dx∂2ux,t∂x2+fx,t,ux,0=ϕx,x∈xL,xR,uxL,t=φt,t∈0,T,uxR,t=ψt,t∈0,T,where the d(x) is the diffusion constant.


## 3.1.2. Numerical Method

Definetn=nτ, xi=ih for 0≤n≤N, 0≤i≤M, where M and N are positive integer and τ=T/N, h=(xR-xL)/M are time step size and space step size, respectively. Define uin, fin, di as the numerical approximation to u(xi,tn), f(xi,tn), and d(xi). Using a forward difference at time tn and a second-order central difference for the space derivative at position xi, we get the explicit finite difference approximation for the one-dimensional heat equation:(2)uin+1-uinτ=diui+1n-2uin+ui-1nh2+fin.Theuin+1 can be obtained by this way:(3)uin+1=a1ui+1n+a2uin+a3ui-1n+τfin,where a1=a3=-diτ/h2 and  a2=1-2diτ/h2.If we use the backward difference at timetn+1 and a second-order central difference for the space derivative at position xi we get the explicit finite difference approximation for (1):(4)uin+1-uinτ=diui+1n+1-2uin+1+ui-1n+1h2+fin+1.Theuin+1 can be obtained by this way:(5)a1ui+1n+1+a2uin+1+a3ui-1n+1=uin+τfin,where a1=a3=diτ/h2 and a2=1+2diτ/h2.The classical PDE is used to compare the fractional equations. For convenience, there are several hypotheses for this paper.(1)
We just focus on the explicit FDM of (3), because these kinds of comparisons between fractional equations and classical equations are straightforward. The comparison between implicit FDE of (5) and implicit schemes of fractional equations is more complex [67].(2)
The computational cost of diffusion coefficientdi and source fin (or fin+1) is ignored, because they only need compute once. And they are different with different equations and make the comparison complicated.


## 3.1.3. Computational Cost

Each grid point of time steptn+1 needs 4 multiplications and 3 additions. There are M-1 grid points in each time step. So each time step needs 7(M-1) arithmetic logic operations. There are about N time steps. So the total computational cost is about 7N(M-1). The computational cost will vary linearly along the number of time steps and grid points.


## 3.2. Time Fractional Diffusion Equation

### 3.2.1. Numerical Method

Liu et al. [63] developed an implicit radial basis function (RBF) meshless approach for time fractional diffusion equations and found that the presented meshless formulation is very effective for modeling and simulation of fractional differential equations. Murillo and Yuste developed an explicit difference method for solving fractional diffusion with Caputo form [68]:(6)∂αux,t∂αt=∂2ux,t∂x20<α<1,ux,0=gx,x∈0,xR,u0,t=uxR,t=0on a finite domain 0<x<xR and 0≤t≤T.The explicit finite difference approximation for (6) is described as follows [69]:(7)τ-α∑i=0n+1‍wiujn+1-i-uj0=1h2uj+1n-2ujn+uj-1n,where wi is the normalized Grünwald weight defined with(8)wi=-1iαi.Defineμ=τα/h2, bn=∑i=0n-1wi, and Un=(u1n,u2n,…,uM-1n)T; we can get(9)Un+1=bn+1U0-∑i=1n+1‍wiUn+1-i+μAUn,where matrix A is a tridiagonal matrix, defined by(10)AM-1×M-1=-211-21••••11-2.


### 3.2.2. Computational Cost

In order to getUn+1, the right-sided computation of (9) should be performed. There are mainly one tridiagonal matrix-vector multiplication, many constant-vector multiplications, and many vector-vector additions in the right-sided computation.(1)
The tridiagonal matrix-vector multiplication isAUn and a new vector Un+1=AUn is got.(2)
The constant-vector multiplications areV′=bn+1U0, V′′=μV1, Vi=wiUn+1-i.(3)
The vector-vector additions areUn+1=Un+1+V′+V′′-∑i=1n+1Vi.There are about5(M-1) operations for tridiagonal matrix-vector multiplication. For time step tn+1 needs (n+8)(M-1) arithmetic logic operations with n=1→N. So the total computational cost for (9) is about(11)8NM-1+1+2+⋯+N-1M-1=8NM-1+NN-1M-12=NN2+7.5M-1.The computational cost varies linearly along the number of grid points but squares with the number of time steps.


## 3.2.1. Numerical Method

Liu et al. [63] developed an implicit radial basis function (RBF) meshless approach for time fractional diffusion equations and found that the presented meshless formulation is very effective for modeling and simulation of fractional differential equations. Murillo and Yuste developed an explicit difference method for solving fractional diffusion with Caputo form [68]:(6)∂αux,t∂αt=∂2ux,t∂x20<α<1,ux,0=gx,x∈0,xR,u0,t=uxR,t=0on a finite domain 0<x<xR and 0≤t≤T.The explicit finite difference approximation for (6) is described as follows [69]:(7)τ-α∑i=0n+1‍wiujn+1-i-uj0=1h2uj+1n-2ujn+uj-1n,where wi is the normalized Grünwald weight defined with(8)wi=-1iαi.Defineμ=τα/h2, bn=∑i=0n-1wi, and Un=(u1n,u2n,…,uM-1n)T; we can get(9)Un+1=bn+1U0-∑i=1n+1‍wiUn+1-i+μAUn,where matrix A is a tridiagonal matrix, defined by(10)AM-1×M-1=-211-21••••11-2.


## 3.2.2. Computational Cost

In order to getUn+1, the right-sided computation of (9) should be performed. There are mainly one tridiagonal matrix-vector multiplication, many constant-vector multiplications, and many vector-vector additions in the right-sided computation.(1)
The tridiagonal matrix-vector multiplication isAUn and a new vector Un+1=AUn is got.(2)
The constant-vector multiplications areV′=bn+1U0, V′′=μV1, Vi=wiUn+1-i.(3)
The vector-vector additions areUn+1=Un+1+V′+V′′-∑i=1n+1Vi.There are about5(M-1) operations for tridiagonal matrix-vector multiplication. For time step tn+1 needs (n+8)(M-1) arithmetic logic operations with n=1→N. So the total computational cost for (9) is about(11)8NM-1+1+2+⋯+N-1M-1=8NM-1+NN-1M-12=NN2+7.5M-1.The computational cost varies linearly along the number of grid points but squares with the number of time steps.


## 3.3. Space Fractional Diffusion Equation

### 3.3.1. Numerical Method

A classical numerical scheme for the space fractional diffusion equation is the second-order fractional Crank-Nicolson method proposed by Tadjeran et al. [38, 70], where the Richardson extrapolation technique is used to the first order shift Grünwald formula for space fractional derivative. Tadjeran et al. [70] presented a practical numerical method in time and space to solve a class of initial-boundary value with variable coefficients on a finite domain for (12)∂ux,t∂t=dx∂αux,t∂xα+qx,ton a finite domain xL<x<xR with 1<α<2 and 0≤t≤T. The case of 1<α<2 models a super-diffusive flow in which a cloud of diffusing particles spreads at a faster rate than the classical diffusion model predictions [71, 72].The explicit finite difference approximation for (12) is(13)uin+1-uinτ=dihα∑k=0i+1‍gα,kui-k+1n+qin,where uin=u(xi,tn), di=d(xi), and qin=q(xi,tn). The gα,k is the normalized Grünwald weight [70, 73]:(14)gα,k=Γk-αΓ-αΓk+1.These normalized weights only depend on the orderα and the index k. The resulting equation can be explicitly solved for uin+1 to give(15)uin+1=uin+diτhα∑k=0i+1‍gα,kui-k+1n+τqin.


### 3.3.2. Computational Cost

Theith grid point of time step tn+1 needs i+6 multiplication, addition, and division. There are M-1 grid points in each time step. So each time step needs (M-1)6+(1+2+⋯+M-1)=(M-1)((M/2)+6) arithmetic logic operations. There are about N time steps. So the total computational cost is about N(M-1)((M/2)+6).


## 3.3.1. Numerical Method

A classical numerical scheme for the space fractional diffusion equation is the second-order fractional Crank-Nicolson method proposed by Tadjeran et al. [38, 70], where the Richardson extrapolation technique is used to the first order shift Grünwald formula for space fractional derivative. Tadjeran et al. [70] presented a practical numerical method in time and space to solve a class of initial-boundary value with variable coefficients on a finite domain for (12)∂ux,t∂t=dx∂αux,t∂xα+qx,ton a finite domain xL<x<xR with 1<α<2 and 0≤t≤T. The case of 1<α<2 models a super-diffusive flow in which a cloud of diffusing particles spreads at a faster rate than the classical diffusion model predictions [71, 72].The explicit finite difference approximation for (12) is(13)uin+1-uinτ=dihα∑k=0i+1‍gα,kui-k+1n+qin,where uin=u(xi,tn), di=d(xi), and qin=q(xi,tn). The gα,k is the normalized Grünwald weight [70, 73]:(14)gα,k=Γk-αΓ-αΓk+1.These normalized weights only depend on the orderα and the index k. The resulting equation can be explicitly solved for uin+1 to give(15)uin+1=uin+diτhα∑k=0i+1‍gα,kui-k+1n+τqin.


## 3.3.2. Computational Cost

Theith grid point of time step tn+1 needs i+6 multiplication, addition, and division. There are M-1 grid points in each time step. So each time step needs (M-1)6+(1+2+⋯+M-1)=(M-1)((M/2)+6) arithmetic logic operations. There are about N time steps. So the total computational cost is about N(M-1)((M/2)+6).


## 3.4. Riesz Space Fractional Diffusion Equation

### 3.4.1. Numerical Method

Shen et al. [74] investigated the Green function and a discrete random walk model for Riesz fractional advection-dispersion equation on infinite domain with an initial condition and also presented implicit and explicit finite differences to this problem on a finite domain. Çelik and Duman [75] applied the Crank-Nicolson method to a fractional diffusion equation which has the Riesz fractional derivative and obtained that the method is unconditionally stable and convergent. The Riesz space fractional reaction-diffusion equation [76, 77] is(16)∂ux,t∂t=-ux,t+xD0αux,t+qx,t,ux,0=gx,x∈xL,xR,uxL,t=uxR,t=0with 1<α≤2 and 0≤t≤T. Both u(x,t) and g(x) are real valued and sufficiently well-behaved function. Dx0αu(x,t) is the Riesz space fractional derivative.With adopting an Euler approximation in time, the explicit difference approximation can be got:(17)uin+1-uinτ=-uin-1chα×∑k=0i+1‍gkui+1-kn+∑k=0M-i+1‍gkui-1+kn+fin,where g0=1, gk=(-1)kαα-1⋯α-k+1/k!, k=1,2,3,… is the normalized Grünwald weight, c=1/2cos⁡(απ/2), and fin=q(xi,tn).Equation (17) results in a linear system of equations(18)Un+1=AUn+Qn,where Un=(u1n,u2n,…,uM-1n)T,  Qn=(τf1n,τf2n,…,τfM-1n)T with qin=τfin (1≤i<M-1), and A=(aij)(M-1)×(M-1) is a matrix of coefficient. A is defined by(19)ai,j=-τchαgi+1-j,for1≤j<i-1,-τchαg0+g2,forj=i-1,1-τ-2τchαg1,forj=i,-τchαg0+g2,forj=i+1,-τchαgj-i+1,fori+1<j<M-1.


### 3.4.2. Computational Cost

From (18), the M-1 results are produced by AUn for each time step. The inner product of vectors v1 and v2 with size m has m multiplications and m-1 additions. AUn for each time step needs (M-1)(2M-3) operations. Assuming τfin is pre-performed, each time step needs (M-1)(2M-2) arithmetic logic operations. There are about N time steps. So the total computational cost is about N(M-1)(2M-2). The computational cost will vary linearly along the number of time steps but square with the number of grid points. From analytical view, the computation of (18) is about four times heavier than that of (15).


## 3.4.1. Numerical Method

Shen et al. [74] investigated the Green function and a discrete random walk model for Riesz fractional advection-dispersion equation on infinite domain with an initial condition and also presented implicit and explicit finite differences to this problem on a finite domain. Çelik and Duman [75] applied the Crank-Nicolson method to a fractional diffusion equation which has the Riesz fractional derivative and obtained that the method is unconditionally stable and convergent. The Riesz space fractional reaction-diffusion equation [76, 77] is(16)∂ux,t∂t=-ux,t+xD0αux,t+qx,t,ux,0=gx,x∈xL,xR,uxL,t=uxR,t=0with 1<α≤2 and 0≤t≤T. Both u(x,t) and g(x) are real valued and sufficiently well-behaved function. Dx0αu(x,t) is the Riesz space fractional derivative.With adopting an Euler approximation in time, the explicit difference approximation can be got:(17)uin+1-uinτ=-uin-1chα×∑k=0i+1‍gkui+1-kn+∑k=0M-i+1‍gkui-1+kn+fin,where g0=1, gk=(-1)kαα-1⋯α-k+1/k!, k=1,2,3,… is the normalized Grünwald weight, c=1/2cos⁡(απ/2), and fin=q(xi,tn).Equation (17) results in a linear system of equations(18)Un+1=AUn+Qn,where Un=(u1n,u2n,…,uM-1n)T,  Qn=(τf1n,τf2n,…,τfM-1n)T with qin=τfin (1≤i<M-1), and A=(aij)(M-1)×(M-1) is a matrix of coefficient. A is defined by(19)ai,j=-τchαgi+1-j,for1≤j<i-1,-τchαg0+g2,forj=i-1,1-τ-2τchαg1,forj=i,-τchαg0+g2,forj=i+1,-τchαgj-i+1,fori+1<j<M-1.


## 3.4.2. Computational Cost

From (18), the M-1 results are produced by AUn for each time step. The inner product of vectors v1 and v2 with size m has m multiplications and m-1 additions. AUn for each time step needs (M-1)(2M-3) operations. Assuming τfin is pre-performed, each time step needs (M-1)(2M-2) arithmetic logic operations. There are about N time steps. So the total computational cost is about N(M-1)(2M-2). The computational cost will vary linearly along the number of time steps but square with the number of grid points. From analytical view, the computation of (18) is about four times heavier than that of (15).


## 3.5. Space-Time Riesz-Caputo Fractional Convection-Diffusion Equation

The fractional advection-diffusion (dispersion) equation has been applied to many problems. Fractional advection-dispersion equations are used in groundwater hydrology to model the transport of passive tracers carried out by fluid flow in a porous medium [78]. Shen et al. [79] presented an explicit difference approximation and an implicit difference approximation for the space-time Riesz-Caputo fractional advection-diffusion equation with initial and boundary conditions in a finite domain. They proved that the implicit difference approximation is unconditionally stable and convergent, but the explicit difference approximation is conditionally stable and convergent. Liu et al. [80] proposed an implicit difference method and an explicit difference method to solve the space-time fractional advection-diffusion equation and discussed the stability and convergence of the method.


### 3.5.1. Numerical Method

We consider the following space-time Riesz-Caputo fractional advection-diffusion equation (STRCFADE) studied by Shen et al. [79]. This equation is obtained by replacing the space-derivative in the advection-diffusion equation with a generalized derivative of order β1, β2 with 0<β1<1, 1<β2<2 and time-derivative with a generalized derivative of order α with 0<α<1. We consider(20)D0tαux,t=B1∂β1ux,t∂xβ1+B2∂β2ux,t∂xβ2+fx,t,ux,0=ϕx,x∈0,xR,u0,t=uxR,t=0,t∈0,Ton a finite domain 0≤x≤xR and 0≤t≤T. The coefficients B1 and B2 are both positive constants and represent the diffusion (dispersion) coefficient and the average fluid velocity. The ∂β1u(x,t)/∂|x|β1 and ∂β2u(x,t)/∂|x|β2 are Riesz space fractional derivatives of order β1 and β2, respectively. Then we can get the explicit finite differential approximation for (20) [79]:(21)uin+1=∑j=0n-1‍bj-bj+1uin-j+bnui0+B1μ1∑k=-iN-i‍ωkβ1ui+knuin+1=+B2μ2∑k=-iN-i‍ωkβ2ui+kn+μ0φin,i=1,…,N-1,u0n+1=uNn+1=0,ui0=ϕi,where μ0=ταΓ(2-α), μ1=μ0/hβ1, and μ2=μ0/hβ2. More information can be referred to in [79].


### 3.5.2. Computational Cost

From Section3.4.2, we know that the computational cost of B1μ1(∑k=-iN-iωkβ1ui+kn)+B2μ2(∑k=-iN-iωkβ2ui+kn)+μ0φin of N time steps is about N(M-1)(2M-2) operations. From Section 3.2.2, we know that the computational cost of ∑j=0n-1(bj-bj+1)uin-j+bnui0 of N time steps with n ranging from 1 to N is about N(N/2)(M-1). So the total computational cost is about N(M-1)(2M-2+(N/2)). The computational cost varies quadratically with the number of time steps or the number of grid points.For fixedM=4097, the comparison between the computational costs of the numerical solution among classical PDE, time fractional, space fractional, Riesz space fractional, and space-time Riesz-Caputo fractional equations is shown in Figure 1. For fixed N=2048, the computational cost is shown in Figure 2.Figure 1
Comparison of the computational cost with different grid points for fixedM=4097 among classical, time fractional, space fractional, Riesz space fractional, and space-time Riesz-Caputo fractional equations.Figure 2
Comparison of the computational cost with different time step for fixedN=2048 among classical, time fractional, space fractional, Riesz space fractional, and space-time Riesz-Caputo fractional equations.


## 3.5.1. Numerical Method

We consider the following space-time Riesz-Caputo fractional advection-diffusion equation (STRCFADE) studied by Shen et al. [79]. This equation is obtained by replacing the space-derivative in the advection-diffusion equation with a generalized derivative of order β1, β2 with 0<β1<1, 1<β2<2 and time-derivative with a generalized derivative of order α with 0<α<1. We consider(20)D0tαux,t=B1∂β1ux,t∂xβ1+B2∂β2ux,t∂xβ2+fx,t,ux,0=ϕx,x∈0,xR,u0,t=uxR,t=0,t∈0,Ton a finite domain 0≤x≤xR and 0≤t≤T. The coefficients B1 and B2 are both positive constants and represent the diffusion (dispersion) coefficient and the average fluid velocity. The ∂β1u(x,t)/∂|x|β1 and ∂β2u(x,t)/∂|x|β2 are Riesz space fractional derivatives of order β1 and β2, respectively. Then we can get the explicit finite differential approximation for (20) [79]:(21)uin+1=∑j=0n-1‍bj-bj+1uin-j+bnui0+B1μ1∑k=-iN-i‍ωkβ1ui+knuin+1=+B2μ2∑k=-iN-i‍ωkβ2ui+kn+μ0φin,i=1,…,N-1,u0n+1=uNn+1=0,ui0=ϕi,where μ0=ταΓ(2-α), μ1=μ0/hβ1, and μ2=μ0/hβ2. More information can be referred to in [79].


## 3.5.2. Computational Cost

From Section3.4.2, we know that the computational cost of B1μ1(∑k=-iN-iωkβ1ui+kn)+B2μ2(∑k=-iN-iωkβ2ui+kn)+μ0φin of N time steps is about N(M-1)(2M-2) operations. From Section 3.2.2, we know that the computational cost of ∑j=0n-1(bj-bj+1)uin-j+bnui0 of N time steps with n ranging from 1 to N is about N(N/2)(M-1). So the total computational cost is about N(M-1)(2M-2+(N/2)). The computational cost varies quadratically with the number of time steps or the number of grid points.For fixedM=4097, the comparison between the computational costs of the numerical solution among classical PDE, time fractional, space fractional, Riesz space fractional, and space-time Riesz-Caputo fractional equations is shown in Figure 1. For fixed N=2048, the computational cost is shown in Figure 2.Figure 1
Comparison of the computational cost with different grid points for fixedM=4097 among classical, time fractional, space fractional, Riesz space fractional, and space-time Riesz-Caputo fractional equations.Figure 2
Comparison of the computational cost with different time step for fixedN=2048 among classical, time fractional, space fractional, Riesz space fractional, and space-time Riesz-Caputo fractional equations.


## 3.6. More Challenges

There are more computational challenges listed below:(1)
high dimensional problems [55, 81],(2)
implicit schemes [67, 82],(3)
high order schemes [33, 83],(4)
variable order problems [84],(5)
huge memory space requirement.The high dimensional problems are more computation expensive [55, 81, 85, 86]. For example, the two-dimensional time fractional diffusion equation (2D-TFDE) [67, 87] is(22)Dtα0ux,y,t=ax,y,t∂2ux,y,t∂x2Dtα0ux,y,t=+ax,y,t∂2ux,y,t∂y2+fx,y,t,ux,y,0=ϕx,y,x,y∈Ω,ux,y,t∂Ω=0,t∈0,T,where Ω={x,y∣0≤x≤L1, 0≤y≤L2, a(x,y,t)>0, and b(x,y,t)>0}.The approximating scheme is [87](23)ui,jn+1-ui,jn+∑s=1n‍bsui,jn+1-s-ui,jn-s=μ1Γ2-αai,jnui+1,jn-2ui,jn+ui-1,jn+μ2Γ2-αbi,jnui,j+1n-2ui,jn+ui,j-1n+ταΓ2-αfi,jn,where μ1=τα/hx2 and μ2=τα/hy2. The hx and hy are the step size along X and Y directions. The computational complexity is about O(N2M2), which is much bigger than the computational complexity of one-dimensional problems when the number of grid points is bigger enough.The direct Gauss elimination for implicit scheme of FDE is not convenient. Solving unsteady FDE with implicit schemes relies on the iteration method at each time step. The variable order problems [84, 88] have complex coefficients, which need more arithmetic logic operations. The high order schemes [33, 38, 70] need more arithmetic logic operations too.The memory usage belongs to the computing resources. So the huge memory requirement of FDE is a kind of computational challenge in the broader sense. This is especially true for time fractional problems. Ignoring the memory usage of the coefficients and source terms, the one-dimensional equation (9) needs 8(M-1)N bytes of memory space and the two-dimensional equation (23) needs 8M2N bytes of memory space. For three-dimensional problems, the memory usage is 8M3N at least. It needs 15.625 PB (1 PB = 10245 bytes) memory space with M=10240, N=2048 for three-dimensional problems. Maybe only the most powerful supercomputer [89] can satisfy the huge memory requirement.


## 4. Potential Solutions

There are several ways which can be used to overcome the computational challenge of FDEs.


### 4.1. Parallel Computing

Large scale applications and algorithms in science and engineering such as neutron transport [90–92], light transport [93], computational fluid dynamics [94, 95], molecular dynamics [96], and computational finance and different numerical methods [97] rely on parallel computing [98, 99].Gong et al. [77] present a parallel algorithm for Riesz space fractional diffusion equation based on MPI parallel programming model at the first time. The parallel algorithm is as accurate as the serial algorithm and achieves 79.39% parallel efficiency compared with 8 processes on a distributed memory cluster system. The parallel implicit iterative algorithm was studied for two-dimensional time fractional problem and a task distribution model is shown in Figure 3 [67]. Mx, My, Px, Py stand for the discrete grid points, parallel processes along X, Y directions.Figure 3
The two-dimensional task distribution model for 2D-TFDE.Domain decomposition method is regarded as the basic mathematical background for many parallel applications [100–102]. A domain decomposition algorithm for time fractional reaction-diffusion equation with implicit finite difference method was proposed [103]. The domain decomposition algorithm keeps the same parallelism as Jacobi iteration but needs much fewer iterations.Diethelm [104] implemented the fractional version of the second-order Adams-Bashforth-Moulton method on a parallel computer and discussed the precise nature of the parallelization concept. Parallel computing has already appeared in some studies on FDEs, but until today their power for approximating fractional derivatives and solving FDEs has not been fully recognized.


### 4.2. Memory Access Optimization (Fractional Precomputing Operator)

Memory access optimization is generally regarded as a technology of computer architecture. It is also useful from the view of applications. One feature of the modern computer architecture is multilayer memory. For example, the traditional CPU has fast cache and slow main memory. The new graphics processing unit (GPU) has fast shared memory and slow memory. Reusing the data in shared memory is a key point to improve the performance of GPU applications [105–109]. A very useful optimization is presented for fractional derivative [69, 110]. We can name this fractional precomputing operator, which will be a very basic and useful optimizing technology for the implementation of fractional algorithms and applications. The example of fractional precomputing operator can refer to Algorithms 4 and 5 in reference [69].On CPU platform, an optimization of the sum of constant vector multiplication is presented and 2-time speedup can be got for both serial and parallel algorithm for the time fractional equation [69]. The key technology is reusing the data in cache through loop unrolling. Zhang et al. [110] presented code optimization for fractional Adams-Bashforth-Moulton method by loop fusion and loop unrolling.On GPU platform, Liu et al. [111] presented an optimized CUDA kernel for the numerical solution of time fractional equation and 1.2–2.7-time performance improvement is achieved.


### 4.3. Short Memory Principle

The short memory principle [3] means the unknown grid points only rely on the recent past (in time) or near neighbors (in space). This principle has been proved to be an easy and powerful way for various kinds of fractional differential equations [112, 113]. The short memory principle is also called fixed memory principle or logarithmic memory principle. The idea of short memory principle is simple. Taking time fractional equation as an example in Section 3.2.1, the value of wi becomes smaller while i is increasing. In (9), the accumulation of ∑i=1n+1wiUn+1-i is replaced with ∑i=min⁡{1,n+1-K}n+1wiUn+1-i. K is a positive integer, which means if n is very big, the computation of the accumulation is fixed.Based on short memory principle, some principles with better balance of computation speed and computation accuracy are presented. The logarithmic memory principle is developed with a good approximation to the true solution at reasonable computational cost [114]. Another principle is equal-weight memory principle, in which an equal-weight is applied to all past data in history, and the result is reserved instead of being discarded [115]. The equal-weight memory principle is an interesting and useful approximation method to fractional derivative.


### 4.4. FFT Based Solution

Because of the nonlocal property of fractional differential operators, the numerical methods have full coefficient matrices which require computational cost ofO(M3) for implicit scheme per time step, where M is the number of grid points. Ford and Simpson [114] developed a faster scheme for the calculation of fractional integrals. A reduction in the amount of computational work can be achieved by using a graded mesh, thereby making the O(N2) method to a O(Nlog⁡N) method. The underlying idea is based on the fact that the fractional integral possesses a fundamental scaling property that can be exploited in a natural way [114].Wang et al. [116, 117] developed a fast finite difference method for fractional diffusion equations, which only requires computational cost of O(Nlog2N) while retaining the same accuracy and approximation property as the regular finite difference method. Numerical experiments show that the fast method has a significant reduction of CPU time [86]. The fast method should have a banded coefficient matrix instead of the full matrix. The properties of Toeplitz and circulant matrices, fast Fourier transform (FFT), and inverse FFT are used to reduce the computational cost. The method is also called superfast-preconditioned iterative method for steady-state space fractional diffusion equations [118].An efficient iteration method for Toeplitz-plus-band triangular systems, which may produced by fractional ordinary differential equations, was developed. Some methods such as matrix splitting, FFT, compress memory storage, and adjustable matrix bandwidth are used in the presented solution. The interesting technologies are the adjustable matrix bandwidth and solving fractional ordinary differential equations with iteration method. The experimental results show that the presented efficient iteration method is 4.25 times faster than the regular solution [10].


### 4.5. Alternating Direction Implicit Method

The alternating direction implicit (ADI) method is a finite difference method for solving traditional PDEs. The approximation methods for fractional equations result in a very complicated set of equations in multiple dimensions, which are hard to solve [67, 81, 85]. So the ADI method is developed for high dimensional problems [33, 117]. The advantage of the fractional ADI method is that the equations that have to be solved in each step have a simpler structure. The time fractional problems can be solved efficiently with the tridiagonal matrix algorithm [33]. The space fractional problems can use the FFT to accelerate the computation [117].


### 4.6. Multigrid Method

The multigrid method is usually exploited for solving ill-conditioned systems. The main idea of multigrid is to accelerate the convergence of a basic iterative method by global correction form, accomplished by solving a coarse problem. Pang and Sun [119] proposed a multigrid method to solve the initial-boundary value problem of a fractional diffusion equation. The experimental results show that the multigrid + FFT method runs hundred times faster than Gaussian elimination method and the conjugate gradient normal residual (CGNR) method. A full V-cycle multigrid method is proposed for the stationary fractional advection dispersion equation [120] and ten-time performance improvement is achieved.


### 4.7. Preconditioning Technology

Preconditioning is typically related to reducing a condition number of the problem with iterative methods. It shows that both the average number of iterations and the CPU time by the PCGNR (preconditioner CGNR) method with circulant preconditioners are much less than those by the CGNR method and less than that by the multigrid method [121]. The circulant preconditioner [121], banded preconditioner [122], fast Poisson preconditioners [123], and preconditioned conjugate gradient squared method plus FFT [118] are developed for different FDEs.


### 4.8. Relationships among These Potential Solutions

The potential solutions for the computational challenge of FDE are investigated above. Many people will be curious about these relationships: can we combine these methods to develop a fastest solution for FDEs? The answer is still unknown. The performance of these solutions varies from different FDE applications. Their relationships are shown in Table1. The PC, MAO, SMP, FFT, ADI, MGM, and PT stand for parallel computing, memory access optimization, short memory principle, FFT based solution, alternating direction implicit method, multigrid method, and preconditioning technology. The score means the degrees of the two solutions are harmonious. Higher score means using two solutions can achieve better performance.Table 1
Relationships among the potential solutions for space fractional derivative.
PC
MAO
SMP
FFT
ADI
MGM
FT
PC
—
3
3
2a
3
2a
2a
MAO
—
3
1b
3
3
3
SMP
—
1c
3
3
3
FFT
—
3
3
3
ADI
—
3
3
MGM
—
2d
PT
—
aThe parallel efficiency of FFT, multigrid, preconditioner is limiting.
bIt is hard to use MAO within FFT.
cBecause of the Toeplitz structure, FFT cannot use short memory principle.
dThe philosophy of multigrid and preconditioning technology are different. But there are some multigrid preconditioners.For time fractional derivative [124], only memory access optimization [69, 111] and short memory principle [3] are useful. Here, the time fractional derivative is different from time fractional problems. For example, the one- and two-dimensional FDEs are parallelized by Gong et al. [67, 69]. These parallel algorithms are based on the partition of space not time.


## 4.1. Parallel Computing

Large scale applications and algorithms in science and engineering such as neutron transport [90–92], light transport [93], computational fluid dynamics [94, 95], molecular dynamics [96], and computational finance and different numerical methods [97] rely on parallel computing [98, 99].Gong et al. [77] present a parallel algorithm for Riesz space fractional diffusion equation based on MPI parallel programming model at the first time. The parallel algorithm is as accurate as the serial algorithm and achieves 79.39% parallel efficiency compared with 8 processes on a distributed memory cluster system. The parallel implicit iterative algorithm was studied for two-dimensional time fractional problem and a task distribution model is shown in Figure 3 [67]. Mx, My, Px, Py stand for the discrete grid points, parallel processes along X, Y directions.Figure 3
The two-dimensional task distribution model for 2D-TFDE.Domain decomposition method is regarded as the basic mathematical background for many parallel applications [100–102]. A domain decomposition algorithm for time fractional reaction-diffusion equation with implicit finite difference method was proposed [103]. The domain decomposition algorithm keeps the same parallelism as Jacobi iteration but needs much fewer iterations.Diethelm [104] implemented the fractional version of the second-order Adams-Bashforth-Moulton method on a parallel computer and discussed the precise nature of the parallelization concept. Parallel computing has already appeared in some studies on FDEs, but until today their power for approximating fractional derivatives and solving FDEs has not been fully recognized.


## 4.2. Memory Access Optimization (Fractional Precomputing Operator)

Memory access optimization is generally regarded as a technology of computer architecture. It is also useful from the view of applications. One feature of the modern computer architecture is multilayer memory. For example, the traditional CPU has fast cache and slow main memory. The new graphics processing unit (GPU) has fast shared memory and slow memory. Reusing the data in shared memory is a key point to improve the performance of GPU applications [105–109]. A very useful optimization is presented for fractional derivative [69, 110]. We can name this fractional precomputing operator, which will be a very basic and useful optimizing technology for the implementation of fractional algorithms and applications. The example of fractional precomputing operator can refer to Algorithms 4 and 5 in reference [69].On CPU platform, an optimization of the sum of constant vector multiplication is presented and 2-time speedup can be got for both serial and parallel algorithm for the time fractional equation [69]. The key technology is reusing the data in cache through loop unrolling. Zhang et al. [110] presented code optimization for fractional Adams-Bashforth-Moulton method by loop fusion and loop unrolling.On GPU platform, Liu et al. [111] presented an optimized CUDA kernel for the numerical solution of time fractional equation and 1.2–2.7-time performance improvement is achieved.


## 4.3. Short Memory Principle

The short memory principle [3] means the unknown grid points only rely on the recent past (in time) or near neighbors (in space). This principle has been proved to be an easy and powerful way for various kinds of fractional differential equations [112, 113]. The short memory principle is also called fixed memory principle or logarithmic memory principle. The idea of short memory principle is simple. Taking time fractional equation as an example in Section 3.2.1, the value of wi becomes smaller while i is increasing. In (9), the accumulation of ∑i=1n+1wiUn+1-i is replaced with ∑i=min⁡{1,n+1-K}n+1wiUn+1-i. K is a positive integer, which means if n is very big, the computation of the accumulation is fixed.Based on short memory principle, some principles with better balance of computation speed and computation accuracy are presented. The logarithmic memory principle is developed with a good approximation to the true solution at reasonable computational cost [114]. Another principle is equal-weight memory principle, in which an equal-weight is applied to all past data in history, and the result is reserved instead of being discarded [115]. The equal-weight memory principle is an interesting and useful approximation method to fractional derivative.


## 4.4. FFT Based Solution

Because of the nonlocal property of fractional differential operators, the numerical methods have full coefficient matrices which require computational cost ofO(M3) for implicit scheme per time step, where M is the number of grid points. Ford and Simpson [114] developed a faster scheme for the calculation of fractional integrals. A reduction in the amount of computational work can be achieved by using a graded mesh, thereby making the O(N2) method to a O(Nlog⁡N) method. The underlying idea is based on the fact that the fractional integral possesses a fundamental scaling property that can be exploited in a natural way [114].Wang et al. [116, 117] developed a fast finite difference method for fractional diffusion equations, which only requires computational cost of O(Nlog2N) while retaining the same accuracy and approximation property as the regular finite difference method. Numerical experiments show that the fast method has a significant reduction of CPU time [86]. The fast method should have a banded coefficient matrix instead of the full matrix. The properties of Toeplitz and circulant matrices, fast Fourier transform (FFT), and inverse FFT are used to reduce the computational cost. The method is also called superfast-preconditioned iterative method for steady-state space fractional diffusion equations [118].An efficient iteration method for Toeplitz-plus-band triangular systems, which may produced by fractional ordinary differential equations, was developed. Some methods such as matrix splitting, FFT, compress memory storage, and adjustable matrix bandwidth are used in the presented solution. The interesting technologies are the adjustable matrix bandwidth and solving fractional ordinary differential equations with iteration method. The experimental results show that the presented efficient iteration method is 4.25 times faster than the regular solution [10].


## 4.5. Alternating Direction Implicit Method

The alternating direction implicit (ADI) method is a finite difference method for solving traditional PDEs. The approximation methods for fractional equations result in a very complicated set of equations in multiple dimensions, which are hard to solve [67, 81, 85]. So the ADI method is developed for high dimensional problems [33, 117]. The advantage of the fractional ADI method is that the equations that have to be solved in each step have a simpler structure. The time fractional problems can be solved efficiently with the tridiagonal matrix algorithm [33]. The space fractional problems can use the FFT to accelerate the computation [117].


## 4.6. Multigrid Method

The multigrid method is usually exploited for solving ill-conditioned systems. The main idea of multigrid is to accelerate the convergence of a basic iterative method by global correction form, accomplished by solving a coarse problem. Pang and Sun [119] proposed a multigrid method to solve the initial-boundary value problem of a fractional diffusion equation. The experimental results show that the multigrid + FFT method runs hundred times faster than Gaussian elimination method and the conjugate gradient normal residual (CGNR) method. A full V-cycle multigrid method is proposed for the stationary fractional advection dispersion equation [120] and ten-time performance improvement is achieved.


## 4.7. Preconditioning Technology

Preconditioning is typically related to reducing a condition number of the problem with iterative methods. It shows that both the average number of iterations and the CPU time by the PCGNR (preconditioner CGNR) method with circulant preconditioners are much less than those by the CGNR method and less than that by the multigrid method [121]. The circulant preconditioner [121], banded preconditioner [122], fast Poisson preconditioners [123], and preconditioned conjugate gradient squared method plus FFT [118] are developed for different FDEs.


## 4.8. Relationships among These Potential Solutions

The potential solutions for the computational challenge of FDE are investigated above. Many people will be curious about these relationships: can we combine these methods to develop a fastest solution for FDEs? The answer is still unknown. The performance of these solutions varies from different FDE applications. Their relationships are shown in Table1. The PC, MAO, SMP, FFT, ADI, MGM, and PT stand for parallel computing, memory access optimization, short memory principle, FFT based solution, alternating direction implicit method, multigrid method, and preconditioning technology. The score means the degrees of the two solutions are harmonious. Higher score means using two solutions can achieve better performance.Table 1
Relationships among the potential solutions for space fractional derivative.
PC
MAO
SMP
FFT
ADI
MGM
FT
PC
—
3
3
2a
3
2a
2a
MAO
—
3
1b
3
3
3
SMP
—
1c
3
3
3
FFT
—
3
3
3
ADI
—
3
3
MGM
—
2d
PT
—
aThe parallel efficiency of FFT, multigrid, preconditioner is limiting.
bIt is hard to use MAO within FFT.
cBecause of the Toeplitz structure, FFT cannot use short memory principle.
dThe philosophy of multigrid and preconditioning technology are different. But there are some multigrid preconditioners.For time fractional derivative [124], only memory access optimization [69, 111] and short memory principle [3] are useful. Here, the time fractional derivative is different from time fractional problems. For example, the one- and two-dimensional FDEs are parallelized by Gong et al. [67, 69]. These parallel algorithms are based on the partition of space not time.


## 5. Future Directions

### 5.1. Fractional Killer Applications

Killer application is a kind of application that is so necessary or desirable that it proves the core value of some larger technology. A killer application is something like Project Apollo in space technology, the IBM 360 in personal computer industry, and the IPhone in the smart phone industry. The fractional research still lacks these kinds of killer applications. It needs fractional applications to solve scientific or engineering problems in physical world, such as the fractional flow/control of hypersonic vehicle, not only the academical problems. The fractional killer application should be proved that it is more useful than the traditional classical application. Solving real problems in physical world will build an economic foundation for fractional researches.


### 5.2. Parallel Computing

The technologies of parallel computing should be regarded as a basic method to overcome the computational challenge for FDEs. There are three potential solutions for the computational challenge of FDEs. The short memory principle is an experimental method, which is useful in real fractional applications. TheO(Nlog⁡N) methods [117] used the property of Toeplitz matrices and FFT technology. Parallel computing is a foundational technology for scientific and engineering computation. Fortunately, the short memory principle and O(Nlog⁡N) methods are compatible with parallel computing. So it is interesting to develop algorithms which is faster than O(Nlog⁡N) methods. The numerical methods of space fractional equations with global dependence are much harder to be parallelized than that of the time fractional equations. So the task distribution model, load balance of the parallel algorithm for space fractional equations should be paid attention to.


### 5.3. High Performance Iteration Methods

Different kind of numerical methods is very easily found for FDEs. The direct methods such as Gauss elimination are not suitable for large scale fractional applications. The iterative methods for these numerical methods are not fully studied and very few works can be found [118, 125]. We think that different iterative methods, such as Jacobi method and Gauss-Seidel and Krylov subspace method, are effective for fractional linear systems which are produced from FDEs. Does the coefficient matrix of FDEs have some other special properties or not? The answer is still unknown. The convergence and stability of these iterative methods should be proved as well.


### 5.4. High Order Schemes for Fractional Derivatives

The traditional partial derivatives have many high order schemes. For time fractional equation, the high order schemes for traditional integer derivative are not hard to build. But for fractional derivatives, the high order schemes are under developing [70, 126–128]. High order schemes will be used in the numerical solutions of FDEs where high accuracy is required in the presence of shocks or discontinuities.


### 5.5. Monte Carlo Method

The Monte Carlo method has advantage in solving nonlinear, high dimensional, complex geometry problems. In order to get the approximation of a small domain, the determined methods, such as FDM, must resolve the total definition domain with boundary conditions. The Monte Carlo method only focuses on the small domain. This property is very useful if we are only interested in this small domain. The Monte Carlo method is easy to be parallelized and needs much less memory space than determined methods. The fractional equations are a kind of nonlinear problem and their high dimensional problems are very computation intensive. So Monte Carlo method for FDEs needs to be studied in the future. Because of the high nonlinear and nonlocal property of FDEs, the sampling efficiency will be the key point of Monte Carlo method for FDEs.


## 5.1. Fractional Killer Applications

Killer application is a kind of application that is so necessary or desirable that it proves the core value of some larger technology. A killer application is something like Project Apollo in space technology, the IBM 360 in personal computer industry, and the IPhone in the smart phone industry. The fractional research still lacks these kinds of killer applications. It needs fractional applications to solve scientific or engineering problems in physical world, such as the fractional flow/control of hypersonic vehicle, not only the academical problems. The fractional killer application should be proved that it is more useful than the traditional classical application. Solving real problems in physical world will build an economic foundation for fractional researches.


## 5.2. Parallel Computing

The technologies of parallel computing should be regarded as a basic method to overcome the computational challenge for FDEs. There are three potential solutions for the computational challenge of FDEs. The short memory principle is an experimental method, which is useful in real fractional applications. TheO(Nlog⁡N) methods [117] used the property of Toeplitz matrices and FFT technology. Parallel computing is a foundational technology for scientific and engineering computation. Fortunately, the short memory principle and O(Nlog⁡N) methods are compatible with parallel computing. So it is interesting to develop algorithms which is faster than O(Nlog⁡N) methods. The numerical methods of space fractional equations with global dependence are much harder to be parallelized than that of the time fractional equations. So the task distribution model, load balance of the parallel algorithm for space fractional equations should be paid attention to.


## 5.3. High Performance Iteration Methods

Different kind of numerical methods is very easily found for FDEs. The direct methods such as Gauss elimination are not suitable for large scale fractional applications. The iterative methods for these numerical methods are not fully studied and very few works can be found [118, 125]. We think that different iterative methods, such as Jacobi method and Gauss-Seidel and Krylov subspace method, are effective for fractional linear systems which are produced from FDEs. Does the coefficient matrix of FDEs have some other special properties or not? The answer is still unknown. The convergence and stability of these iterative methods should be proved as well.


## 5.4. High Order Schemes for Fractional Derivatives

The traditional partial derivatives have many high order schemes. For time fractional equation, the high order schemes for traditional integer derivative are not hard to build. But for fractional derivatives, the high order schemes are under developing [70, 126–128]. High order schemes will be used in the numerical solutions of FDEs where high accuracy is required in the presence of shocks or discontinuities.


## 5.5. Monte Carlo Method

The Monte Carlo method has advantage in solving nonlinear, high dimensional, complex geometry problems. In order to get the approximation of a small domain, the determined methods, such as FDM, must resolve the total definition domain with boundary conditions. The Monte Carlo method only focuses on the small domain. This property is very useful if we are only interested in this small domain. The Monte Carlo method is easy to be parallelized and needs much less memory space than determined methods. The fractional equations are a kind of nonlinear problem and their high dimensional problems are very computation intensive. So Monte Carlo method for FDEs needs to be studied in the future. Because of the high nonlinear and nonlocal property of FDEs, the sampling efficiency will be the key point of Monte Carlo method for FDEs.


## 6. Conclusions

In this paper, we give a comprehensive review of FDEs and its computational challenge. This kind of challenge will become an incoming problem for the computer industry if the real fractional problems need to be approximated. We reviewed a wide range of computational costs that come from different kinds of fractional equations. While we have collected several potential solutions on this challenge, we believe that the long-term legacy of solutions will allow the real world scientific and engineering applications come true.



---

*Source: 258265-2015-03-26.xml*