A spectrally accurate step-by-step method for the numerical solution of fractional differential equations

In this paper we consider the numerical solution of fractional differential equations. In particular, we study a step-by-step graded mesh procedure based on an expansion of the vector field using orthonormal Jacobi polynomials. Under mild hypotheses, the proposed procedure is capable of getting spectral accuracy. A few numerical examples are reported to confirm the theoretical findings.


Introduction
In recent years fractional differential equation modelling has become more and more frequent -see, for example, the two monographs [37] and [19].In fact the use of fractional models has quite a long history and early applications considered their behaviour in describing anomalous diffusion in a randomly hetergeneous porous medium [31] and in water resources modelling [6].In mathematical biology, a number of researchers have considered fractional reaction diffusion equations -see, for example, [16] who have studied via a spectral approach the interfacial properties of the Allen-Cahn equation with a quartic double well potential, a model often used for alloy kinetics and the dynamics of phase boundaries.Henry, Langlands and Wearne [27] analysed Turing pattern formation in fractional activator-inhibitor systems.In a computational medicine setting, Henry and Langlands [26] developed fractional cable models for spiny neuronal dendrites.Orsingher and Behgin [36] have considered the dynamics of fractional chemical kinetics for unimolecular systems using time change.The authors in [17] have made a series of arguments based on Riesz potential theory and the wide range of heterogeneous cardiac tissues for the use of fractional models in cardiac electrophysiology, while [18] have explored the effects of fractional diffusion on the cardiac action potential shape and electrical propagation.In addition, [34], [28] and [15] have considered the calibration of fractional cardiac models through a fractional Bloch-Torrey equation.
Due to these application advances there is clearly a need to develop advanced numerical methods for fractional differential equations.In a purely time fractional setting the standard approach has been to use a first order approximation to the Caputo derivative.Lubich [32] has extended this and developed higher order approximations based on an underlying linear multistep method that are convergent of the order of this multistep method.
There are two important issues when designing effective numerical methods for Fractional Differential equations: memory and non-locality.In the case of memory, in principle one needs to evaluate the numerical approximations all the way back to the initial time point at each time step and this becomes very computationally intensive over long time intervals.Techniques such as fixed time windowing have been proposed [37] and more recently a number of approaches have been considered to compute the discrete convolutions to approximate the fractional operator [39], [43].Zeng et al. [42] developed fast memory-saving time stepping methods in which the fractional operator is decoupled into a local part with fixed memory length and the history part is computed by Laguerre-Gauss quadrature, whose error is independent of the stepsize.
In the second situation the solutions to FDEs are often non-smooth and may have a strong singularity at t = 0.In these cases the vector field may inherit this singularity behaviour and so a constant timestep would be very inefficient.Thus graded meshes have been proposed by a number of authors [30,44,40].Lubich [33] on the other hand deals with singularities by introducing certain correction terms.Other approaches have been considered such as Adams methods [21], trapezoidal methods [23] and Bernstein splines [38] while Garrappa [24] has given a survey and software tutorial on numerical methods for FDEs.
In this paper, we consider a major improvement of the recent solution approach described in [2,10], based on previous work on Hamiltonian Boundary value Methods (HBVMs) [12,8,9,7] (also used as spectral methods in time [14,13,3,5]), for solving fractional initial value problems in the form: where, without loss of generality, we omit t as a formal argument, at the r.h.s.Further, for sake of brevity, the dimension m of the state vector will be also omitted, when unnecessary.
Here, for α ∈ (0, 1], y (α) (t) ≡ D α y(t) is the Caputo fractional derivative: The Riemann-Liouville integral associated to (2) is given by: It is known that [37]: Consequently, the solution of ( 1) is formally given by: We shall generalize the approach given in [2], by defining a step-by-step procedure, whereas the procedure given in [2] was of a global nature.This strategy, when combined with a graded mesh selection, enables better handling of singularities in the derivative of the solution to (1) at the origin, which occur, for example, when f is a suitably regular function.Additionally, our approach partially mitigates issues stemming from the non-local nature of the problem.Indeed, at a given time-step of the integration procedure, it requires the solution of a differential problem within the current time interval along with a finite sequence of pure quadrature problems that account for the contribution brought by the history term.
With this premise, the structure of the paper is as follows: in Section 2 we recall some basic facts about Jacobi polynomials, for later use; in Section 3 we describe a piecewise quasi-polynomial approximation to the solution of (1); in Section 4 the analysis of the new method is carried out; in Section 5 we report a few numerical tests for assessing the theoretical achievements; at last, in Section 6 a few concluding remarks are given.

Orthonormal Jacobi polynomials
Jacobi polynomials form an orthogonal set of polynomials w.r.t. a given weighting function.In more detail, for r, ν > −1: where, as is usual, Π i is the vector space of polynomials of degree at most i and δ ij is the Kronecker symbol.Consequently, the shifted and scaled Jacobi polynomials In particular, hereafter we shall consider the polynomials with α ∈ (0, 1] the same parameter in (1), such that: where we have slightly changed the weighting function, in order that it has a unit integral: We refer to [25] and the accompanying software, for their computation, and for computing the nodes and weights of the Gauss-Jacobi quadrature of order 2k.
As is well known, the polynomials (6) form an orthonormal basis for the Hilbert space L 2 [0, 1], equipped with the scalar product and the associated norm Consequently, from the Parseval identity, it follows that any square summable function can be expanded along such a basis, and the corresponding expansion is convergent to the function itself.That said, with reference to the vector field in (1), let us consider the interval [0, h], and the expansion with the Fourier coefficients given by: Consequently, on the interval [0, h], (1) can be rewritten as and, by virtue of (3)-(4), one obtains: In particular, due to ( 7) and ( 8), one has: Remark 2 By considering that P 0 (c) ≡ 1, from (12) one derives that Consequently, because of (3), (15) becomes i.e., (5) with t = h.This clearly explains the use of the Jacobi basis for the expansion (11).Further, we observe that, when α = 1, one retrieves the approach described in [12] for ODE-IVPs.
We also need the following preliminary results.
Lemma 1 Let G : [0, h] → V , with V a vector space, admit a Taylor expansion at t = 0.Then, Proof By the hypothesis and ( 7), one has: Concerning the above result, we observe that, by considering ( 9)- (10), from the Cauchy-Schwartz theorem one derives that From Lemma 1, one has the following: Corollary 1 Assume that the r.h.s. of problem (1) admits a Taylor expansion at t = 0, then the coefficients defined at (12) satisfy: The result of the previous lemma can be weakened as follows (the proof is similar and, therefore, omitted).
Lemma 2 Let G : [0, h] → V , with V a vector space, admit a Taylor polynomial expansion of degree k with remainder at t = 0.Then, In such a case, the result of Corollary 1 is weakened accordingly.
Corollary 2 Assume that the r.h.s. of problem (1) admits a Taylor polynomial expansion of degree k with remainder at t = 0.Then, the coefficients defined at (12) satisfy: However, in general, at t = 0 the solution may be not regular, since the derivative may be singular (this is, indeed, the case, when f is a regular function).In such a case, we shall resort to the following result.
Lemma 3 Assume that the r.h.s. in ( 1) is continuous for all t ∈ [0, h].Then, for a convenient ξ t ∈ (0, t), one has: Proof In fact, from (4), and by using the weighted mean-value theorem for integrals, one has: As a consequence, we derive the following weaker results concerning the Fourier coefficients (12).
Corollary 3 In the hypotheses of Lemma 3, and assuming f is continuously differentiable in a neighborhood of y 0 , one has: Proof In fact, from Lemma 3, one has: From this relation, and from (12), the statement then follows.□ For later use, we also need to discuss the quadrature error in approximating the first s Fourier coefficients (12) by means of the Gauss-Jacobi formula of order 2k, for a convenient k ≥ s, whose abscissae are the zeros of P k , as defined in (7), and with weights: That is, Concerning the quadrature errors the following result holds true.
Theorem 1 Assume the function G j,h (c) := P j (c)f (y(ch)) be of class C (2k) [0, h].Then, with reference to (18), one has, for a suitable ξ = (0, 1) : with the constants K j independent of both G j,h and h.
Proof The statement follows by considering that the quadrature is exact for polynomials in Π 2k−1 .□ However, as recalled above, when f is a smooth function, the derivative of the solution of problem (1) may have a singularity at t = 0.In such a case, estimates can be also derived (see, e.g., [20,41] and [35,Theorem 5.1.8]).However, for our purposes, because of ( 7), for j = 0, . . ., s − 1, from Corollary 3 we can easily derive the following one, assuming that k ≥ s and ( 13) hold true: In order to derive a polynomial approximation to (13), it is enough to truncate the infinite series in (11) to a finite sum, thus obtaining: with γ j (σ) formally given by ( 12) upon replacing y by σ.In so doing, ( 14) and ( 15) respectively become: and It can be shown that Corollary 1, Corollary 2, Lemma 3, and Corollary 3 continue formally to hold for σ.Further, by considering the approximation of the Fourier coefficients obtained by using the Gauss-Jacobi quadrature of order 2k, the result of Theorem 1 continues to hold.Moreover, we shall assume that the expansion holds true, similarly as (11), from which also (20) follows for the quadrature error Σ j (σ, h), when σ has a singular derivative at 0. In the next sections, we shall better detail, generalize, and analyze this approach.

Piecewise quasi-polynomial approximation
To begin with, in order to obtain a piecewise quasi-polynomial approximation to the solution of (1), we consider a partition of the integration interval in the form: where Then, according to ( 21)-( 23), on the first subinterval [0, h 1 ] we can derive a polynomial approximation of degree s − 1 to (11), thus getting the fractional initial value problem where γ j (σ 1 ) is formally given by ( 12), upon replacing y by σ 1 at the right-hand side: The solution of ( 28) is a quasi-polynomial of degree s − 1 + α, formally given by: However, in order to actually compute the Fourier coefficients {γ j (σ 1 )}, one needs to approximate them by using the Gauss-Jacobi quadrature (16): In so doing, (30) now formally becomes Moreover, one solves the system of equations, having (block) dimension s independently of k: This kind of procedure is typical of HBVM(k, s) methods, in the case of ODE-IVPs [8]: the main difference, here, derives from the non-locality of the operator.As is clear (compare with ( 23)), the approximation at t = h 1 will be now given by:3 Remark 4 For an efficient and stable evaluation of the integrals we refer to the procedure described in [4].

The general procedure
We now generalize the previous procedure for the subsequent integration steps.For later use, we shall denote: the restriction of the solution on the interval [t n−1 , t n ].Similarly, we shall denote by σ(t) ≈ y(t) the piecewise quasi-polynomial approximation such that: Then, by using an induction argument, let us suppose one knows the quasi-polynomial approximations where ϕ α i−1 (c, σ) denotes a history term, to be defined later, such that: • in the first subinterval, according to (31), The corresponding approximations at the grid-points are given by ȳi : and assume we want to compute Hereafter, we shall assume that the time-steps in (26) define a graded mesh.In more detail, for a suitable r > 1: Remark 5 As is clear, in the limit case where r = 1, (40) reduces to a uniform mesh with a constant time-step h = T /N .
In order to derive the approximation (39), we start computing the solution of the problem in the subinterval [t n , t n+1 ].Then, for t ≡ t n + ch n+1 , c ∈ [0, 1], one has: having set the history term which is a known quantity, since it only depends on y 1 , . . ., y n (see (35)).Consequently, we have obtained which reduces to (5), when n = 0 and t ∈ [0, h].Further, by considering the expansion with the Fourier coefficients given by: one obtains that (42) can be rewritten as: with the value at t = h given by: The corresponding approximation is obtained by truncating the series in (42) after s terms, and approximating the corresponding Fourier coefficients via the Gauss-Jacobi quadrature of order 2k: with ϕ α n (c, σ) an approximation of G α n (c, y) in ( 41), defined as follows: having set which, as we shall see in Section 3.3, can be accurately and efficiently computed.
As is clear, in (49) we have used the approximation having (block)-size s, independently of the considered value of k.

Solving the discrete problems
The discrete problem (55), to be solved at each integration step, can be better cast in vector form by introducing the matrices I α P 0 (c 1 ) . . .I α P s−1 (c 1 ) . . . . . .
and the (block) vectors: In fact, in so doing (55) can be rewritten as: This formulation is very similar to that used for HBVMs in the case α = 1 [11].The formulation (56) has a twofold use: • it shows that, assuming, for example, f Lipschitz continuous, then the solution exists and is unique, for all sufficiently small timesteps h n+1 ; • it induces a straightforward fixed-point iteration: which we shall use for the numerical tests.
In fact, the following result holds true.
Theorem 2 Assume f be Lipchitz with constant L in in the interval [t n , t n+1 ].Then, the iteration (57) is convergent for all timesteps h n+1 such that Proof In fact, one has: hence the iteration function defined at (57) is a contraction, when (58) holds true.□ A simplified Newton-type iteration, akin to that defined in [11] for HBVMs, will be the subject of future investigations.

Remark 6
We observe that the discrete problem (56) can be cast in a Runge-Kutta type form.In fact, the vector in argument to the function f at the r.h.s. in (56), can be regarded as the stage vector of a Runge-Kutta method, tailored for the problem at hand.Substituting (56) into (59), and using Y n+1 as argument of f , then gives It is worth noticing that (60) has (block) dimension k, instead of s.However, by considering that usually k ≫ s (see Section 5), it follows that solving (60) is less efficient than solving (56).This fact is akin to what happens for a HBVM(k, s) method [12,8,9], which is the k-stage Runge-Kutta method obtained when α = 1.In fact, for such a method, the discrete problem can also be cast in the form (56), then having (block) dimension s, independently of k (which is usually much larger than s).

Computing the integrals J α j (x)
From ( 55) and (49), it follows that one needs to compute the integrals with {c 1 , . . ., c k } the abscissae of the Gauss-Jacobi quadrature (16).As an example, in Figure 1 we plot the Gauss-Jacobi abscissae for α = 0.5 and k = 5, 10, 15, 20, 25, 30.Further, there is numerical evidence that a sufficiently high-order Gauss-Legendre formula can compute the integrals (50) up to round-off, when x ≥ 1.5.Namely, From the last two formulae, one derives that Jjalfa(a,b,d,alfa,1+ci*r) computes all the integrals in (64) corresponding to the abscissa ci, with Jjalfa the Matlab © function listed in Table 1, and the vectors a,b,d containing the scalars (67).An implementation of the previous function, using the standard variable precision arithmetic (i.e., using vpa in Matlab © ), allows handling values of s up to 20, at least.

Analysis of the method
From (45) and (47), and with reference to (51)-(54), one derives: Assuming f Lipschitz with constant L in a suitable neighborhood of the solution, and setting one then obtains that, for h α n+1 sufficiently small, and a suitable constant K 1 > 0: Considering that 1 Γ(α) and (see (54)) so that (see (40)), for a suitable constant K 2 > 0, one eventually obtains: with (see (40) and ( 68)) At last, by considering that, for n = 1, . . ., N − 1, setting and considering that e 1 ≤ O(h 2α 1 ), from the Gronwall lemma (see, e.g., [29]) one eventually obtains, for a suitable K > 0: Considering that (see (27) and ( 40)) N = log r (1 + T (r − 1)/h 1 ), the error is optimal when h 2α 1 ∼ h s+α N .Remark 8 In the case where the vector field of problem (1) is suitably smooth at t = 0, so that a constant timestep h = T /N can be used, the estimate can be derived for the error, by using similar arguments as above.

Numerical Tests
We here report a few numerical tests to illustrate the theoretical findings.For all tests, when not otherwise specified, we use k = 30 for the Gauss-Jacobi quadrature (16), and p = 30 for the Gauss-Legendre quadrature formula (65).Also, the following values of s will be considered: The first problem, taken from [24], is given by: whose solution is given by the following Mittag-Leffler function:6 E 0.6 (t) = j≥0 (−10t 0.6 ) j Γ(0.6j + 1) .
In Table 2, we list the obtained results, in terms of maximum error, when using different values of s and h 1 .We use the value r = 1.01 in (40), and the number of timesteps N is chosen in order that t N is the closest point to the end of the integration interval.The ***, in the line corresponding to s = 1, means that the solution is not properly evaluated: this is due to the failure of the fixedpoint iteration (57) in the last integration steps.As is expected from (71), the error decreases as s increases and the initial timestep h 1 decreases.Further, having a large value of s is not effective, if h 1 is not suitably small, and vice versa (again from the bound (71)).Remarkably, by choosing s large enough and h 1 suitably small, full machine accuracy is gained (cf. the last 4 entries in the last column of the table, having the same error).
In Table 3 we list the results obtained by using k = s, which is the minimum value allowed.In such a case, the accuracy is generally slightly worse, due to the fact that the quadrature error is of the same order as the truncation error.It is, however, enough choosing k only slightly larger than s, to achieve a comparable accuracy, as is shown in Table 4 for k = s + 5. Nevertheless, it must be stressed that choosing larger values of k is not an issue, since the discrete problem (56), to be solved at each integration step, has (block) size s, independently of k.The next problem that we consider is from [21] (see also [24]): y (0.5) = −y According to [24], "this problem is surely of interest because, unlike several other problems often proposed in the literature, it does not present an artificial smooth solution, which is indeed not realistic in most of the fractional-order applications."Despite this, a constant timestep h = 1/N turns out to be appropriate since, unlike the solution (75), the vector field (74) is very smooth at the origin, as one may see in Figure 2. In Table 5 we list the maximum error by using different values of N : as is clear, by using s ≥ 8, only 32 steps are needed to gain full machine accuracy for the computed solution.
Next, we consider the following problems, taken from [38]: whose solution is y(t) = t 2/3 + 1, and whose solution is y(t) = t 4/3 .We solve Problem (76) by using h 1 = 10 −11 and r = 1.2.In so doing 130 timesteps are needed to cover approximately the integration interval.The obtained results are listed in Table 6.Also in this case, we obtain full accuracy starting from s = 8.
In addition, we solve Problem (77) by using a constant timestep h = 1/N : in fact, the vector field can be seen to be a polynomial of degree 1.The obtained results are listed in Table 7: an order 1 convergence can be observed for s = 1 (which is consistent with (72)), whereas full machine accuracy is obtained for s ≥ 2, due to the fact that, as anticipated above, the vector field of problem (77) is a polynomial of degree 1 in t and, consequently, ( 13) and ( 21) coincide, for all s ≥ 2.
Finally, we reformulate the two problems (76) and (77) as a system of two equations, having  We solve Problem (78) by using the same parameters used for Problem (76): h 1 = 10 −11 , r = 1.2, and 130 timesteps.The obtained results are again listed in Table 6: it turns out that they are similar to those obtained for Problem (76) and, also in this case, we obtain full accuracy starting from s = 8.

Conclusions
In this paper we have devised a novel step-by-step procedure for solving fractional differential equations.The procedure, which generalizes that given in [2], relies on the expansion of the vector field along a suitable orthonormal basis, here chosen as the shifted and orthonormal Jacobi polynomial basis.The analysis of the method has been given, along with its implementation details.These latter details show that the method can be very efficiently implemented.A few numerical tests confirm the theoretical findings.
It is worth mentioning that systems with FDEs of different orders can be also solved by slightly adapting the previous framework: as matter of fact, it suffices using different Jacobi polynomials, corresponding to the different orders of the FDEs.This, in turn, allows easily managing fractional differential equations of order α > 1, by casting them as a system of ⌊α⌋ ODEs, coupled with a fractional equation of order β := α − ⌊α⌋ ∈ (0, 1).
As anticipated in Section 3.2, a future direction of investigation will concern the efficient numerical solution of the generated discrete problems, which is crucial when using larger timesteps.Also the optimal choice of the parameters for the methods will be further investigated, as well as the possibility of adaptively defining the computational mesh.Further, the parallel implementation of the methods could be also considered, following an approach similar to [1].

Table 2 :
Maximum error for Problem (73), r = 1.01 and k = 30.All the numerical tests have been performed in Matlab © .As anticipated, we use the fixed-point iteration (57) to solve the generated discrete problems (56).The iteration is carried out until full machine accuracy is gained.