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

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 are O(NM), O(NM), 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), shortmemory principle, fast Fourier transform (FFT) based solutions, alternating direction implicitmethod,multigridmethod, 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.


Introduction
The idea of fractional is natural.If / and  2 / 2 exist,  1.5 / 1.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][4][5][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 2 Mathematical Problems in Engineering 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][20][21].The implementations (computer area) are closer to the human society and the real applications but vary quickly along with the computer architecture [22][23][24][25][26][27][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.

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 the nth-derivative of the linear function () = ,   /  [29,30].L'Hopital's posed the question to Leibniz: what would the result be if  = 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].
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.
We can suppose that if 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 Monte Carlo (MC) method, which uses repeated random sampling to obtain numerical results, belongs to undetermined computational methods.

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.

Computational Challenge
where the () is the diffusion constant.
The  +1  can be obtained by this way: where  1 =  3 = −  /ℎ 2 and  2 = 1 − 2  /ℎ 2 .If we use the backward difference at time  +1 and a second-order central difference for the space derivative at position   we get the explicit finite difference approximation for (1): The  +1  can be obtained by this way: where 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 coefficient   and source    (or  +1

𝑖
) is ignored, because they only need compute once.And they are different with different equations and make the comparison complicated.

Computational Cost.
Each grid point of time step  +1 needs 4 multiplications and 3 additions.There are  − 1 grid points in each time step.So each time step needs 7( − 1) arithmetic logic operations.There are about  time steps.So the total computational cost is about 7( − 1).The computational cost will vary linearly along the number of time steps and grid points.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]:

Time Fractional Diffusion Equation
on a finite domain 0 <  <   and 0 ≤  ≤ .

Computational Cost.
In order to get  +1 , the rightsided 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 is   and a new vector  +1 =   is got.
(3) The vector-vector additions are There are about 5( − 1) operations for tridiagonal matrix-vector multiplication.For time step  +1 needs ( + 8)( − 1) arithmetic logic operations with  = 1 → .So the total computational cost for ( 9) is about The computational cost varies linearly along the number of grid points but squares with the number of time steps.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 initialboundary value with variable coefficients on a finite domain for

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.

Numerical Method.
We consider the following spacetime Riesz-Caputo fractional advection-diffusion equation (STRCFADE) studied by Shen et al. [79].This equation is obtained by replacing the space-derivative in the advectiondiffusion 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 on a finite domain 0 ≤  ≤   and 0 ≤  ≤ .The coefficients  1 and  2 are both positive constants and represent the diffusion (dispersion) coefficient and the average fluid velocity.The   1 (, )/||  1 and   2 (, )/||  2 are Riesz space fractional derivatives of order  1 and  2 , respectively.Then we can get the explicit finite differential approximation for (20) [79]: where  0 =   Γ(2 − ),  1 =  0 /ℎ  1 , and  2 =  0 /ℎ  2 .More information can be referred to in [79].For fixed  = 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  = 2048, the computational cost is shown in Figure 2.
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(−1) bytes of memory space and the twodimensional equation ( 23) needs 8 2  bytes of memory space.For three-dimensional problems, the memory usage is 8 3  at least.It needs 15.625 PB (1 PB = 1024 5 bytes) memory space with  = 10240,  = 2048 for three-dimensional problems.Maybe only the most powerful supercomputer [89] can satisfy the huge memory requirement.

Potential Solutions
There are several ways which can be used to overcome the computational challenge of FDEs.
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].  ,   ,   ,   stand for the discrete grid points, parallel processes along ,  directions.Domain decomposition method is regarded as the basic mathematical background for many parallel applications [100][101][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.

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][106][107][108][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-timeperformance improvement is achieved.

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   becomes smaller while  is increasing.In ( 9), the accumulation of is a positive integer, which means if  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.

FFT Based Solution.
Because of the nonlocal property of fractional differential operators, the numerical methods have full coefficient matrices which require computational cost of ( 3 ) for implicit scheme per time step, where  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 ( 2 ) method to a ( log ) 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 (log 2 ) 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.

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.

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 Table 1.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.
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 twodimensional FDEs are parallelized by Gong et al. [67,69].These parallel algorithms are based on the partition of space not time.

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.

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.The ( log ) 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 ( log ) methods are compatible with parallel computing.So it is interesting to develop algorithms which is faster than ( log ) 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.

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.

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][127][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.

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 longterm legacy of solutions will allow the real world scientific and engineering applications come true.

Table 1 :
Relationships among the potential solutions for space fractional derivative.The parallel efficiency of FFT, multigrid, preconditioner is limiting.b It is hard to use MAO within FFT. c Because of the Toeplitz structure, FFT cannot use short memory principle.d The philosophy of multigrid and preconditioning technology are different.But there are some multigrid preconditioners. a