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 are O(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.

1. Introduction

The idea of fractional is natural. If and exist, 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 [36]. 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 [1921]. The implementations (computer area) are closer to the human society and the real applications but vary quickly along with the computer architecture [2228]. 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 the nth-derivative of the linear function , [29, 30]. L’Hopital’s posed the question to Leibniz: what would the result be if ? 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 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 numerical schemes [49] for these FDEs are listed below:(1)finite difference method (FDM) [18, 50, 51],(2)finite element method [5254],(3)finite volume method [5557],(4)Adomian decomposition method [58],(5)Fourier transform method [59],(6)spectral method [60],(7)meshless method [6163],(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:where the is the diffusion constant.

3.1.2. Numerical Method

Define , for , , where and are positive integer and , are time step size and space step size, respectively. Define , , as the numerical approximation to , , and . Using a forward difference at time and a second-order central difference for the space derivative at position , we get the explicit finite difference approximation for the one-dimensional heat equation:

The can be obtained by this way:where and  .

If we use the backward difference at time and a second-order central difference for the space derivative at position we get the explicit finite difference approximation for (1):

The can be obtained by this way:where and .

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 ) 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 step needs 4 multiplications and 3 additions. There are grid points in each time step. So each time step needs arithmetic logic operations. There are about time steps. So the total computational cost is about . 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]:on a finite domain and .

The explicit finite difference approximation for (6) is described as follows [69]:where is the normalized Grünwald weight defined with

Define , , and ; we can getwhere matrix is a tridiagonal matrix, defined by

3.2.2. Computational Cost

In order to get , 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 is and a new vector is got.(2)The constant-vector multiplications are , , .(3)The vector-vector additions are .

There are about operations for tridiagonal matrix-vector multiplication. For time step needs arithmetic logic operations with . 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. 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 on a finite domain with and . The case of 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) iswhere , , and . The is the normalized Grünwald weight [70, 73]:

These normalized weights only depend on the order and the index . The resulting equation can be explicitly solved for to give

3.3.2. Computational Cost

The th grid point of time step needs multiplication, addition, and division. There are grid points in each time step. So each time step needs arithmetic logic operations. There are about time steps. So the total computational cost is about .

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] iswith and . Both and are real valued and sufficiently well-behaved function. is the Riesz space fractional derivative.

With adopting an Euler approximation in time, the explicit difference approximation can be got:where , , is the normalized Grünwald weight, , and .

Equation (17) results in a linear system of equationswhere ,   with (), and is a matrix of coefficient. is defined by

3.4.2. Computational Cost

From (18), the results are produced by for each time step. The inner product of vectors and with size has multiplications and additions. for each time step needs operations. Assuming is pre-performed, each time step needs arithmetic logic operations. There are about time steps. So the total computational cost is about . 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 , with , and time-derivative with a generalized derivative of order with . We consideron a finite domain and . The coefficients and are both positive constants and represent the diffusion (dispersion) coefficient and the average fluid velocity. The and are Riesz space fractional derivatives of order and , respectively. Then we can get the explicit finite differential approximation for (20) [79]:where , , and . More information can be referred to in [79].

3.5.2. Computational Cost

From Section 3.4.2, we know that the computational cost of of time steps is about operations. From Section 3.2.2, we know that the computational cost of of time steps with ranging from to is about . So the total computational cost is about . The computational cost varies quadratically with the number of time steps or the number of grid points.

For fixed , 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 , the computational cost is shown in Figure 2.

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] iswhere , , , and .

The approximating scheme is [87]where and . The and are the step size along and directions. The computational complexity is about , 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 bytes of memory space and the two-dimensional equation (23) needs bytes of memory space. For three-dimensional problems, the memory usage is at least. It needs 15.625 PB (1 PB = bytes) memory space with , 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 [9092], 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]. , , , stand for the discrete grid points, parallel processes along , directions.

Domain decomposition method is regarded as the basic mathematical background for many parallel applications [100102]. 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 [105109]. 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 becomes smaller while is increasing. In (9), the accumulation of is replaced with . 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.

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 of 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 method to a 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 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 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 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. The 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 methods are compatible with parallel computing. So it is interesting to develop algorithms which is faster than 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, 126128]. 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.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This research work is supported in part by the National Natural Science Foundation of China under Grants no. 61402039, no. 61303265, and no. 60970033, in part by Specialized Research Fund for the Doctoral Program of Higher Education under Grant no. 20114307110001, and in part by 973 Program of China under Grants no. 2014CB430205 and no. 61312701001. The authors would like to thank the anonymous reviewers for their helpful comments as well.