It is nowadays well established that several real–life phenomena are better described by fractional differential equations (FDEs), where the term fractional, used for historical reasons, refers to derivative operators of any real positive order. Applications of FDEs are commonly found in bioengineering, chemistry, control theory, electronic circuit theory, mechanics, physics, seismology, signal processing and so on (e.g., [3, 7, 5, 33, 22, 41, 50]). We refer to  for an historical perspective on fractional calculus.
As a consequence of the growing interest for fractional order models, the analysis and the development of numerical methods for FDEs and partial differential equations with time or space fractional derivatives have become an active area of research (see, for instance,[4, 18, 17, 13, 25, 26, 32, 35, 38, 44, 45, 47, 48, 46, 53, 54]).
Some of the methods for FDEs are directly derived from methods for integral equations; this is the case, for instance, of product integration (PI) rules. A completely different approach, named as fractional linear multistep methods (FLMMs) , has instead been specifically tailored to provide a solid theoretical background for the numerical treatment of FDEs.
Interestingly, these two approaches (PI and FLMM) lead to the same methods when the fractional order
converges to the nearest integer; indeed, different methods for FDEs are often generalizations of the same corresponding method for ordinary differential equations (ODEs).
Some important differences are however observed in the pure fractional case. For instance, the PI and FLMM generalizations of the implicit Euler method behave in a different way when applied to problems with discontinuous right–hand side [10, 11, 12] and only the method devised in the framework of FLMMs succeeds in stabilizing the solution on the switching surface, as expected from theoretical considerations, without generating unwanted spurious oscillations .
The presence of distinguishing characteristics in methods for FDEs developed from the same method for ODEs deserves to be more thoroughly investigated. Some differences are of theoretical interest (convergence, stability, etc.); nevertheless, the two approaches also present distinctive features of computational nature that affect the efficiency of the solution process.
The main aim of this paper is to describe and compare some methods for FDEs, highlight strengths and weaknesses and investigate the circumstances under which a certain method is preferable to some others.
The analysis will be confined to those PI rules and FLMMs that generalize, in some respects, the implicit second–order trapezoidal rule. This –stable one–step method has the smallest possible error constant  and is appreciated since provides a satisfactory balance between accuracy, stability properties and computational complexity. Our interest is in exploring when these strengths are inherited by the corresponding methods for FDEs.
There are also some more specific motivations for restricting the attention to second–order methods: a) PI rules very seldom converge with an order equal or greater than 2; b) the implementation of higher order FLMMs is usually computationally demanding; c) as for multistep methods for ODEs, order 2 represents a barrier for stability.
Despite this apparently limited extent, we will however be able to include in our analysis a sufficiently large number of methods.
The choice of focusing on implicit methods is also motivated by the fact that they have been so far very seldom investigated (very often, implicit methods are taken into account in the framework of predictor–corrector algorithms in which the implicit scheme is actually implemented as an explicit one) and explicit methods are unduly preferred in most of the works on FDEs, regardless of stability issues which discourage their use.
The paper is organized as follows. In Section 2 we present some basic definitions and properties concerning FDEs. Section 3 is devoted to describe PI rules and discuss different implementations on uniform and graded grids. In Section 4 we present FLMMs and we describe three different generalizations of the implicit trapezoidal rule; the use of an effective formula for the evaluation of the weights is also proposed. In Section 5 linear stability is investigated and some implementation details are discussed in Section 6. Finally, in Section 7 the results of some numerical tests are presented to compare the methods under investigation and some concluding remarks are presented at the end of the paper.
2 Fractional differential equations
Let us consider the initial value problem for a system of FDEs in the form
where is the smallest integer such that , is assumed to be sufficiently smooth, is the unknown solution and denotes classical derivative of integer order .
with the Riemann–Liouville (RL) integral of order on the interval
A classical result in fractional calculus  allows to compute, for any , the RL integral of the power function as
It is a well known result that (1) can be equivalently written in terms of the integral equation
where is the Taylor expansion of centered at
One of the major difficulties in the treatment of FDEs derives from the lack of smoothness of the true solution on the whole interval of integration. Indeed, as proved in , the true solution of (1) possesses an expansion in mixed (integer and real) powers of the form
with , and
with the consequence that the derivatives of are unbounded at .
3 Product integration rules
) by approximating the vector fieldwith suitable polynomials. For the integer–order case this is the way in which Adams multistep methods are devised.
To generalize the –step trapezoidal rule consider a (not necessarily uniform) grid on the whole interval of integration and, in the piecewise decomposition of (4)
in each subinterval by the first–degree polynomial interpolant
where , and is the numerical approximation to the solution . After exactly solving the integrals
where, for shortness, we denoted
it is elementary, after a simple change of the summation index, to obtain the numerical approximation of the solution at as
According to the analysis presented in , this scheme is convergent of second order whenever the true solution is assumed sufficiently smooth. Unfortunately, this is a rather optimistic result compared to the expansion (5) of : due to the presence of unbounded derivatives at the convergence rate of (6) is actually lower than the theoretical value , as it can be formally stated by a direct application of the main result by Dixon in .
Let be Lipschitz continuous with respect to the second variable and the numerical approximation obtained by applying the PI trapezoidal rule (6) on the interval . There exists a constant , which does not depend on , such that
Unless very restrictive (and usually unrealistic) assumptions are made to ensure the smoothness of also at the left endpoint of , the rate by which the error reduces to is therefore, away from , proportional to when . The same drop in the convergence order occurs also when higher degree polynomials are employed; for this reason PI rules of higher order are never taken into account for FDEs when . When order of convergence is instead obtained.
3.1 Uniform meshes
Using (6) on uniform grids , with a constant step–size, offers some advantages of computational type. The evaluation of the weights is indeed less expensive and just involves the computation of real powers of integer numbers. Furthermore, (6) presents a convolution structure, in which each weight depends just on , and it can be simplified as
3.2 Graded meshes
Unlike FLMMs, PI rules can be implemented on non uniform grids, with the aim of overcoming the order reduction stated in Theorem 1. This issue has been investigated, in the context of weakly singular integral equations, by several authors [1, 9] and graded grids are highly praised to this purpose.
By clustering the grid–points at the left boundary of the integration interval, graded grids are able to cope with the lack of regularity and counteract the inefficiency of polynomials to satisfactorily approximate real powers.
A graded mesh on the interval is defined according to
where is a suitable grading exponent. After denoting and observing that , a simple manipulation allows to reformulate (6) according to
Let , . For it is
The use of PI rules on nonuniform grids is more demanding; not only the weights must be recomputed at each step but, since the absence of a convolution structure, it is not possible to employ fast algorithms. By means of some numerical experiments we will verify, at the end of the paper, whether or not this further amount of computation is, in some way, balanced by the enhancement in the accuracy.
4 Fractional linear multistep methods
In the pioneering work on discretized fractional calculus , Lubich proposed an elegant and effective strategy for classical linear multistep methods (LMMs) originally devised for integer–order ODEs.
Although the potentials of FLMMs have been recognized in several papers (see, for instance, the review article by Gorenflo ), their use for practical computation has not captured so far widespread attention, perhaps due to some difficulties in explicitly formulating FLMMs.
The key aspect in FLMMs is the approximation of the RL integral (2) by means of the convolution quadrature
on uniform grids , , and where convolution and starting quadrature weights and do not depend on .
Starting weights play an important role, especially in the first part of the integration interval, in order to cope with the possibly singular character of the integrand function at ; we will discuss the problem of their evaluation later on.
Convolution quadrature weights are the main components of the quadrature rule and characterize the specific FLMM. They are obtained starting from any LMM for ODEs
and the function goes therefore under the name of generating function of the LMM. The idea proposed in  is to generate a quadrature rule for fractional problems (1) by evaluating the convolution weights as the coefficients in the FPS of the fractional–order power of the generating function
Methods of this kind, named as FLMMs, when applied to (1) read as
Let be a stable and consistent of order implicit LMM with the zeros of having absolute value . The FLMM (10) is convergent of order .
One of the major difficulties in FLMMs (10) is in evaluating the weights as the coefficients in the FPS (9). Although some sophisticated algorithms are available for manipulating FPS , for most of the methods an efficient tool is the J.C.P. Miller formula stated by the following theorem [30, Theorem 1.6c].
Let be a FPS. Then for any ,
where coefficients can be recursively evaluated as
The Miller formula allows to evaluate, in the most general case, the first coefficients of with a number of operations proportional to . In most cases of practical interest, this computational effort is actually reduced by an order of magnitude. This is the case of the evaluation of the coefficients in the FPS of ; indeed, since and , it is elementary to see that the application of Theorem 4 leads to
involving just multiplications and additions.
4.1 Fractional trapezoidal rule
In order to generalize to FDEs the classical trapezoidal rule, let us first consider its classical formulation for ODEs
with characteristic polynomials and and generating function
To practically compute the first weights in the FPS of we advise against the direct application of the Miller formula to the infinite series in since it would involve a number of operations proportional to . It is instead preferable to first apply twice Theorem 4 to and , by means of (11), and hence evaluate the coefficient of their product by an FFT algorithm, with a number of operations proportional to when is a power of . This last task can be efficiently performed by means of the very few lines of Matlab code
x = fft([omega1,zeros(size(omega1))]) ; y = fft([omega2,zeros(size(omega2))]) ; omega = ifft(x.*y) ; omega = omega(1:N) ;
where omega1 and omega2 are the row arrays with coefficients of and respectively, evaluated in a fast way thanks to (11).
4.2 Newton–Gregory formula
The trapezoidal rule (12) belongs to the more general family of –step Adams–Moulton methods
with the classical backward differences (i.e., and ); the coefficients are the first terms in the truncation
of the power series expansion, with respect to , of and the resulting generating function is .
It is easy to see that and , thus immediately leading to (12) when . Therefore, a straightforward generalization to fractional–order problems by means of the –power of leads, in the case of the –step method, to the trapezoidal rule discussed in the previous subsection.
It has been shown in  that alternative methods, named as Newton–Gregory formulas, can be devised by first evaluating the power of and successively truncating to the first terms as
The interchange between the truncation series symbol and the –power is important and generates a different family of methods. By operating for simplicity the change of variable ,
Since and , the generating function of the FLLM obtained by the –step Adams method is therefore
and, after a simple manipulation, it is possible to evaluate the first coefficients in , with a cost proportional to , as
with the coefficients in the FPS evaluated by means of (11).
4.3 Fractional BDF formula
The second order BDF formula for ODEs is given by
with characteristic polynomials and and generating function
BDF method (13) can be considered in no way as coming from the trapezoidal rule; it is indeed obtained by differentiating a second–degree interpolant polynomial. However, it is equally of interest in our analysis due to its second order of convergence as the trapezoidal rule and the excellent stability properties that we are interested in verifying also for FDEs.
The coefficients in the FPS of can be profitably evaluated again by the Miller formula applied to , with negative power , to first produce
and hence by computing with an overall number of operations that still remains proportional to .
5 Linear stability analysis
The analysis of stability is of primary importance to understand the qualitative behavior of a numerical method. By assuming the existence of a stable steady state, stability analysis aims to verify whether or not approximated trajectories move towards the steady state.
As usual, stability issues for FDEs can be addressed by studying the scalar linear test equation
Let . The steady–state of (14) is stable if and only if , where .
In Figure 1 we present the sector of stability for two exemplifying cases. It is immediate to verify that coincides with the left complex semi–plane when , thus generalizing the classical result for ODEs.
When PI rules or FLMMs are applied to (14), the numerical solution is evaluated by means of a convolution quadrature
Let . Under the assumption that the sequence is convergent and that the quadrature weights satisfy
the stability region of the convolution quadrature (15) is
In the numerical treatment of ODEs –stability means that the whole complex semi–plane is included in the stability region of the method. –stable methods are therefore particularly attractive since they ensure that the numerical solution of stable linear problems converges to regardless of the size of the step . Because of the result in Theorem 5, for FDEs it is more useful to use the concept of –stability, which means that the whole sector must be included in the stability region.
While is a priori known for FLMMs, the generating function of PI rules does not possess an analytical representation and, therefore, it can be just numerically approximated.
In Figure 2 we present a comparison of the stability regions (the gray areas) of the methods discussed in the previous sections for various values (dotted lines denote the corresponding sectors of stability ). To identify each method, here and throughout the following section we make use of the acronyms listed in Table 1.
|PI U||PI rule on uniform grid||3.1|
|PI G||PI rule on graded grid||3.2|
|FT||fractional trapezoidal rule||4.1|
|FBDF||fractional BDF formula||4.3|
We can distinctly observe that, as increases from to , the stability region of each method decreases but always includes ; we can draw the conclusion that all the methods are –stable for .
Within ODEs, the stability region of the trapezoidal rule is the negative semi–plane and exactly corresponds to the region in which the true solution converges to the steady–state. Among the methods presented in this paper, just FT preserves this feature and indeed its stability regions coincide with for any value of . NG and PI have regions of almost the same size but larger with respect to FT; as in the ODE case, FBDF possesses the largest regions of stability.
Interestingly, the situation is completely different for as shown in Figure 3. Whilst the stability regions of FT still coincide with the stability sectors and the stability properties of FBDF still remain the most remarkable ones, we observe that NG and PI completely lose the –stability. This is an important aspect to be taken into account: despite the implicit nature of these methods, numerical instability can arise when a not sufficiently small step–size is used and, therefore, they seem not suited to solve stiff problems when .
6 Implementation details and Matlab codes
In this section we discuss some implementation issues that are essential for the development of efficient codes.
6.1 Evaluation of the lag term
The evaluation of the (usually long) lag term is one of the most expensive tasks in convolution quadratures like (7) or (10). The integration along a large number of steps is usually a critical bottleneck and, if not well designed, it can be an impracticable task involving a number of operations proportional to .
In our codes we implement the FFT algorithm described in  by which the computational effort downgrades from to . We do not further discuss technical details of this algorithm, for which we just refer to the above mentioned paper; we simply observe that, without efficient algorithms of this kind, the numerical solution of FDEs very quickly degenerates in an unworkable task as the number of grid–points increases.
6.2 Starting weights
Starting weights are introduced in (10) to cope with the singular behavior of the solution close to the left endpoint of the integration interval. They are evaluated after imposing that (8) is exact for with and the order of convergence of the FLMM .
Because of (3), and since the methods investigated in this paper have order , it is immediate to see that weights must satisfy the relation
and, hence, they are evaluated by solving, at each step , a system of linear equations, with the cardinality of .
This problem has been in–depth discussed in , where the authors stressed the attention on the mildly ill–conditioned nature of the coefficient matrix of (16). Nevertheless, for the methods under investigation this is not a serious issue since the size of is small, being (in most cases will be 3 or 4 and the inversion of can be done analytically). The only exception is for very small values of which, however, are extremely uncommon in applications.
6.3 Starting the computation
The initialization of (10) requires the knowledge of the first approximations of the solution. Since just is given by the problem, the remaining values must be evaluated in some other way, for instance by using some different methods.
To avoid mixing methods of different nature, we prefer to evaluate all the approximations at the same time by setting the single implicit system
This nonlinear system has size , where is the size of the system of FDEs we want to solve. No particular problems arise when the size of the problem is small; otherwise some memory allocation problems could occur in presence of very small values of , involving a very large number of initial approximations to compute.
6.4 Solution of the nonlinear systems
where and do not depend on . To preserve the good stability properties, it is convenient to use Newton iterations
where is solution of the linear system
and, usually, . We refer to  for more insights on this issue and for a discussion on the stopping criteria.
We just highlight here that the convergence can depend on , instead of as in ODEs. Whilst this is not an issue when , some unexpected problems could arise whenever is small. For this reason it is necessary to use a quite small step–size with low values of .
7 Numerical experiments
By means of numerical experiments we compare the methods discussed in the paper. We perform all the experiments in Matlab, version 22.214.171.1249, on a computer equipped with the Intel Dual Core E5400 processor running at 2.70 GHz under the Windows XP operating system; all the Matlab codes have been optimized, to get the best performance, according to the discussion reported in the previous sections.
As a first test problem we consider the linear FDE (14) for . For each method, and for an increasing number of grid–points, we report the error at
, with respect to a reference solution, and we provide an estimation of the order of convergence (EOC) as.
The numerical results in Table 2 for clearly confirm theoretical findings on the order of convergence. As expected from Theorem 1, the PI rule on uniform grids fails to obtain order 2 and converges with order ; on the other hand, the PI rule on graded grids converges with full order .
As in the ODE case, FT produces the smallest error. Although also NG and PI are generalizations of the trapezoidal rule, it seems that only FT inherits this feature; by taking into account, from the analysis of Section 5, that FT preserves the other peculiar character of the trapezoidal rule of having the region of stability exactly equal to the region in which the solution converges towards the steady state, it is tempting to conclude that FT is the most natural generalization of the trapezoidal rule.
Unfortunately, the situation is more complicate since when it is NG that provides the lowest error, as shown in Table 3 for .
This phenomenon can be more clearly observed in Figure 4, where we show the errors as a function of (we used the same values for and from the previous experiment and a fixed number of grid points).
We must highlight that this is a not completely unexpected result. Stability and accuracy are indeed conflicting requirements in most of numerical methods; thus, FT shows the smallest regions of stability but the lowest error when , exactly as it happens with NG for . What we observe is a kind of roles swap between FT and NG when moves from the left to the right of 1.
The very bad performance of PI on graded grids for low values of is due to the fact that the grids generated by the graded exponent are very coarse close to right endpoint of the integration interval. Furthermore, graded grids do not add any improvement for , since the maximum expected order is already achieved by using uniform grids. We can draw the conclusion that the utility of graded grids is confined only to .
It is also interesting to analyze the computational effort in relation to the accuracy. To this purpose, the plots in Figure 5, which are related to the experiments presented in Tables 2 and 3, show that FT, NG and FBDF behave in a more efficient way with respect to PI when ; indeed, PI U pays for the reduction in the accuracy whilst PI G pays for the absence of a convolution structure which precludes the application of FFT algorithms. The minor efficiency of FBDF with respect to FT and NG is due to the major complexity in the evaluation of weights. It is worthwhile to note that the preeminence of FLMMs with respect to PI rules holds despite the need of solving linear systems for evaluating the starting weights. For the absence of order reduction increases the efficiency of PI on uniform nodes but however it is below that of NG which instead has lower errors.
To observe the behavior of the methods under investigation for nonlinear problems, we consider the fractional reaction diffusion system of Brusselator type
It has been shown  that there exists a marginal value