Numerical solution of FDE-IVPs by using Fractional HBVMs: the fhbvm code

In this paper we describe the efficient numerical implementation of Fractional HBVMs, a class of methods recently introduced for solving systems of fractional differential equations. The reported arguments are implemented in the Matlab code fhbvm, which is made available on the web. An extensive experimentation of the code is reported, to give evidence of its effectiveness.


Introduction
Fractional differential equations have become a common description tool across a variety of applications (see, e.g., the classical references [23,36] for an introduction).For this reason, their numerical solution has been the subject of many researches (see, e.g.[1,24,25,34,35,37,38]), with the development of corresponding software (see, e.g., [21,27,28]).In this context, the present contribution is addressed for solving initial value problems for fractional differential equations (FDE-IVPs) in the form where, for the sake of brevity, we have omitted the argument t for f .Here, for α ∈ (0, 1), y (α) (t) ≡ D α y(t) is the Caputo fractional derivative: The Riemann-Liouville integral associated to (2) is given by: Consequently, the solution of (1) can be formally written as: The numerical method we shall consider, relies on Fractional HBVMs (FH-BVMs), a class of methods recently introduced in [8], as an extension of Hamiltonian Boundary Value Methods (HBVMs), special low-rank Runge-Kutta methods originally devised for Hamiltonian problems (see, e.g., [10,11]), and later extended along several directions (see, e.g., [2,4,6,8,9,12]), including the numerical solution of FDEs.A main feature of HBVMs is the fact that they can gain spectrally accuracy, when approximating ODE-IVPs [3,19,20], and such a feature has been recently extended to the FDE case [8].
With this premise, the structure of the paper is as follows: in Section 2 we recall the main facts about the numerical solution of FDE-IVPs proposed in [8]; in Section 3 we provide full implementation details for the Matlab © code fhbvm used in the numerical tests; in Section 4 we report an extensive experimentation of the code, providing some comparisons with another existing one; at last, a few conclusions are given in Section 5.

Fractional HBVMs
To begin with, in order to obtain a piecewise approximation to the solution of the problem, we consider a partition of the integration interval in the form: where h n > 0, n = 1, . . ., N , and such that Further, by setting the restriction of the solution of (1) on the interval [t n−1 , t n ], and taking into account ( 4) and ( 5)- (6), one obtains, for t ≡ t n−1 + ch n , c ∈ [0, 1], To make manageable the handling of the ratios h i /h ν , i = ν, . . ., n, ν = 1, . . ., n − 1, we shall hereafter consider the following choices for the mesh ( 5)-( 6): -Graded mesh.In order to cope with possible singularities in the derivative of the vector field at the origin, we consider the graded mesh where r > 1 and h 1 > 0 satisfy, by virtue of ( 5)-( 6), As a result, one obtains that (8) becomes: -Uniform mesh.This case is equivalent to allowing r = 1 in (9).Consequently, from ( 5)-( 6) we derive As a result, (8) reads: Remark 1 As is clear, in order to obtain an accurate approximation of the solution, it is important to establish which kind of mesh (graded or uniform) is appropriate.Besides this, also a proper choice of the parameters h 1 , r, and N in ( 9)-( 10) is crucial.Both aspects will be studied in Section 3.1.

Quasi-polynomial approximation
We now discuss a piecewise quasi-polynomial approximation to the solution of (1), such that is the approximation to y n (ch n ), defined in (7).According to [8], such approximation will be derived through the following steps: 1. expansion of the vector field, in each sub-interval [t n−1 , t n ], n = 1, . . ., N , (recall ( 5)-( 6)) along a suitable orthonormal polynomial basis; 2. truncation of the infinite expansion, in order to obtain a local polynomial approximation.
Let us first consider the expansion of the vector field along the orthonormal polynomial basis, w.r.t. the weight function resulting into a scaled and shifted family of Jacobi polynomials: 2 In so doing, for n = 1, . . ., N , one obtains: with (see (15)) Consequently, the FDE (1), can be rewritten as: The approximation is derived by truncating the infinite series in (17) to a finite sum with s terms, thus replacing (19) with a series of local problems, whose vector field is a polynomial of degree s: with γ j (σ n ) defined similarly as in (18): As a consequence, (11) will be approximated as: s−1 j=0 s−1 j=0 where (see (3)) is the Riemann-Liouville integral of P j (c), and, for x ≥ 1: The efficient numerical evaluation of the integrals ( 23) and ( 24) will be explained in detail in Section 3.2.
Remark 2 We observe that, for c = x = 1, by virtue of ( 15)-( 16) one has: Similarly, when using a uniform mesh, by means of similar steps as above, (13) turns out to be approximated by: In both cases, the approximation at t n is obtained, by setting c = 1 and taking into account (25), as: which formally holds also in the case of a uniform mesh (see (12)) by setting r = 1.

The fully discrete method
Quoting Dahlquist and Björk [22], "as is well known, even many relatively simple integrals cannot be expressed in finite terms of elementary functions, and thus must be evaluated by numerical methods".In our framework, this obvious statement means that, in order to obtain a numerical method, at step n the Fourier coefficients γ j (σ n ) in (21) need to be approximated by means of a suitable quadrature formula.Fortunately enough, this can be done up to machine precision by using a Gauss-Jacobi formula of order 2k based at the zeros of P k (c), c 1 , . . ., c k , with corresponding weights (see (15)) by choosing a value of k, k ≥ s, large enough.In other words, where .= means equal within machine precision.Because of this, and for sake of brevity, we shall continue using σ n to denote the fully discrete approximation (compare with (22), or with (26), in case r = 1): As is clear, the coefficientes γ ν j , j = 0, . . ., s − 1, ν = 1, . . ., n − 1, needed to evaluate φ α,s n−1 (c; h 1 , r, σ), have been already computed at the previous timesteps.
We observe that, from (28), in order to compute the (discrete) Fourier coefficients, it is enough evaluating (29) only at the quadrature abscissae c 1 , . . ., c k .In so doing, by combining ( 28) and ( 29), one obtains a discrete problem in the form: Once it has been solved, according to (27), the approximation of y(t n ) is given by: It is worth noticing that the discrete problem ( 30) can be cast in vector form, by introducing the (block) vectors and the matrices I α P 0 (c 1 ) . . .I α P s−1 (c 1 ) . . . . . . as: with the obvious notation for the function f , evaluated in a vector of (block) dimension k, of denoting the (block) vector of dimension k containing the function f evaluated in each (block) entry.As observed in [8], the (block) vector appearing at the r.h.s. in (32) as argument of f , satisfies the equation obtained combining (32) and (33).Consequently, it can be regarded as the stage vector of a Runge-Kutta type method with Butcher tableau Remark 3 Though the two formulations ( 32) and ( 34) are equivalent each other, nevertheless, the former has (block) dimension s independently of the considered value of k (k ≥ s), which is the (block) dimension of the latter.Consequently, in the practical implementation of the method, the discrete problem (32) is the one to be preferred, since it is independent of the number of stages.
The efficient numerical solution of the discrete problem (32) will be considered in Section 3.4.

Implementation issues
In this section we report all the implementation details used in the Matlab © code fhbvm, which will be used in the numerical tests reported in Section 4.

Graded or uniform mesh?
The first relevant problem to face is the choice between a graded or a uniform mesh and, in the former case, also choosing appropriate values for the parameters r, N , and h 1 in ( 9)- (10).We start considering a proper choice of this latter parameter, i.e., h 1 which, in turn, will allow us to choose the type of mesh, too.For this purpose, the user is required to provide, in input, a convenient integer value M > 1, such that, if a uniform mesh is appropriate, then the stepsize is given by: It is to be noted that, since the code is using a method with spectral accuracy in time, the value of M should be as small as possible.That said, we set h 0 1 = h and apply the code on the interval [0, h 0 1 ] both in a single step, and by using a graded mesh of 2 sub-intervals defined by a value of r := r ≡ 3.As a result, the two sub-intervals, obtained by solving (10) for h 1 , with In so doing, we obtain two approximations to y(h 0 1 ), say y 1 and y 2 .According to the analysis in [8], if5 with tol a suitably small tolerance, 6 then this means the stepsize h 1 := h 0 1 is appropriate.Moreover, in such a case, a uniform mesh with N ≡ M turns out to be appropriate, too.
Conversely, the procedure is repeated on the sub-interval [0, h 1 1 ] ≡ [0, h 0 1 /4], and so forth, until a convenient initial stepsize h 1 is obtained.This process can be repeated up to a suitable maximum number of times that, considering that at each iteration the current guess of h 1 is divided by 4, allows for a substantial reduction of the initial value (36).Clearly, in such a case, a graded mesh is needed.At the end of this procedure, we have then chosen the initial stepsize h 1 which, if the procedure ends at the ℓ-th iteration, for a convenient ℓ > 1, is given by: We need now to appropriately choose the parameters r and N in ( 9)- (10).To simplify this choice, we require that the last stepsize in the mesh be equal to (36).Consequently, we would ideally satisfy the following requirements: which, by virtue of (38), become Combining the two equations then gives: As is clear, this value of N is not an integer, in general, so that we shall choose, instead: Remark 4 From (40), and considering that N h 1 < T (conversely, a uniform mesh could have been used) and, by virtue of (38), one has: As is clear, using the value of N in (40) in place of that in (39) implies that we need to recompute r, so that the requirement (now h 1 , N , and T are given) is again fulfilled.Equation (42) can be rewritten as thus inducing the iterative procedure As a convenient choice for r 0 , one can use the guess for r defined in (39).The following result can be proved.
From the two latter properties we conclude that, for any r 0 > 1, the sequence generated by (44) converges monotonically to r.
Remark 5 In the actual implementation of the code, we allow the use of a uniform mesh also when ℓ = 2 steps of the previous procedure are required, provided that M has a moderate value (say, M ≤ 5).Consequently, for the final mesh, h 1 = h/4, r = 1, and N = 4M .

Approximating the fractional integrals
The practical implementation of the method requires the evaluation of the following integrals (recall (30) and the definitions ( 23)-( 24)): and in case a graded mesh is used, or in case a uniform mesh is used.
It must be emphasized that all such integrals ((47) and (48), or (47) and ( 49)), can be pre-computed once for all, for later use.For their computation, in the first version of the software, we adapted an algorithm based on [5] which, however, required a quadruple precision, for the considered values of k and s.In this respect, the use of the standard vpa of Matlab © , which is based on a symbolic computation, turned out to be too slow.For this reason, hereafter we describe two new algorithms for computing the above integrals, which result to be quite satisfactory, when using the standard double precision IEEE.Needless to say that, since vpa is no more required, they turn out to be much faster than the previous ones.
Let us describe, at first, the computation of (47).One has, by considering that k ≥ s and c i ∈ (0, 1): where the last equality holds because the Jacobi quadrature formula has order 2k.Clearly, for each i = 1, . . ., k, all the above integrals can be computed in vector form in "one shot", by using the usual three-term recurrence to compute the Jacobi polynomials.
Concerning the integrals (48)-( 49), let us now consider the evaluation of a generic J α j (x), j = 0, . . ., s − 1, for x > 1.In this respect, there is numerical evidence that, for x ≥ 1.1, a high-order Gauss-Legendre formula is able to approximate the required integral to full machine precision.Since we will use a value of s = 20, we consider, for this purpose, a Gauss-Legendre formula of order 60, which turns out to be fully accurate.Instead, for x ∈ (1, 1.1), one has: again, due to the fact that the quadrature is exact for polynomials of degree at most s − 1.Also in this case, for each fixed x > 1, all the integrals can be computed in "one shot" by using the three-term recurrence of the Jacobi polynomials.

The nonlinear iteration
At last, we describe the efficient numerical solution of the discrete problem (32), which has to be solved at the n-th integration step.As is clear, the very formulation of the problem induces a straightforward fixed-point iteration: which can be conveniently started from γ n,0 = 0.The following straightforward result holds true.
Theorem 2 Assume f be Lipchitz with constant L in in the interval Then, the iteration (50) is convergent for all timesteps h n such that Nevertheless, as is easily seen, also the simple equation whose solution is almost everywhere close to 0, after an initial transient, suffers from stepsize limitations, if the fixed-point iteration (50) is used, since it has to be everywhere proportional to λ −1/α .In order to overcome this drawback, a Newton-type iteration is therefore needed.Hereafter, we consider the so-called blended iteration which has been at first studied in a series of papers [7,14,16,17].It has been implemented in the Fortran codes BIM [15], for ODE-IVPs, and BIMD [18], for ODE-IVPs and linearly implicit DAEs, and in the Matlab code hbvm [10,13], for solving Hamiltonian problems.We here consider its adaption for solving (32).By neglecting, for sake of brevity, the time-step index n, we then want to solve the equation: By setting f ′ 0 the Jacobian of f evaluated at the first entry of φ α,s , I = I s ⊗I m , and the application of the simplified Newton method then reads: set : Even though this iteration has the advantage of using a coefficient matrix which is constant at each time-step, nevertheless, its dimension may be large, when either s or m are large.To study a different iteration, able to get rid of this problem, let us decouple the linear system into the various eigenspaces of f ′ 0 , thus studying the simpler problem solve : with all involved vectors of dimension s, and µ ∈ σ(f ′ 0 ) a generic eigenvalue of f ′ 0 , and an obvious meaning of g(γ ℓ ).By setting q = h α µ, the iteration then reads: solve : set : Hereafter, we consider the iterative solution of the linear system in (54).A linear analysis of convergence (in the case the r.h.s. is constant) is then made, as at first suggested in [31,32,33], and later refined in [14,17].Consequently, skipping the iteration index ℓ, let us consider the linear system to be solved: and its equivalent formulation, derived considering that matrix (52) is nonsingular, and with ξ > 0 a parameter to be later specified, Further, we consider the blending of the previous two equivalent formulations with weights θ(q) and I s − θ(q), where, by setting O ∈ R s×s the zero matrix, In so doing, one obtains the linear system with the coefficient matrix, such that [14]: This naturally induces the splitting matrix defining the blended iteration This latter iteration converges iff the spectral radius of the iteration matrix, The iteration is said to be A-convergent if (58) holds true for all q ∈ C − , the left-half complex plane, and L-convergent if, in addition, ρ(q) → 0, as q → ∞.
Since [14] θ(0)M (0) = I s , θ(q)M (q) → I s , q → ∞, the blended iteration is L-convergent iff it is A-convergent.For this purpose, we shall look for a suitable choice of the positive parameter ξ > 0. Considering that θ(q) is well defined for all q ∈ C − , the following statement easily follows from the maximum modulus theorem.

Theorem 3
The blended iteration is L-convergent iff the maximum amplification factor, The following result also holds true.Theorem 4 The eigenvalues of the iteration matrix I s − θ(q)M (q) are given by: Consequently, the maximum amplification factor is given by: Proof See [14, Theorem 2 and Equation ( 25)], by considering that is obtained at x = ξ −1 , so that (59) follows.
⊓ ⊔ Slightly generalizing the arguments in [14], we then consider the following choice of the parameter ξ, which is computed once forall, and always provides, in our experiments, an L-convergent iteration.In particular, the code fhbvm uses, at the moment, k = 22 and s = 20: the corresponding maximum amplification factor (59) is depicted in Figure 1, w.r.t. the order α of the fractional derivative, thus confirming this.
Coming back to the original problem (53), starting from the initial guess ∆γ = 0, and updating the r.h.s. as soon as a new approximation to γ is available, one has that the iteration (55)-(57) simplifies to: with ξ chosen according to (60), and Consequently, only the factorization of one matrix, having the same dimension m of the problem, is needed.Moreover, the initial guess γ 0 = 0 can be conveniently considered in (61).
Remark 6 It is worth mentioning that, due to the properties of the Kronecker product, the iteration (61) can be compactly cast in matrix form, thus avoiding an explicit use of the Kronecker product.This implementation has been considered in the code fhbvm used in the numerical tests.
Actually, according to Theorem 2, in the code fhbvm we automatically switch between the fixed-point iteration (50) or the blended iteration (61), depending on the fact that with tol < 1 a suitable tolerance.

Numerical Tests
In this section we report a few numerical tests using the Matlab © code fhbvm: the code implements a FHBVM (22,20) method using all the strategies discussed in the previous section.The calling sequence of the code is: [t,y,stats,err] = fhbvm( fun, y0, T, M ) where: In input: fun is the identifier (or the function handling) of the function evaluating the r.h.s. of the equation (also in vector mode), its Jacobian, and the order α of the fractional derivative (see help fhbvm for more details); -y0 is the initial condition; -T is the final integration time; -M is the parameter in (36) (it should be as small as possible); In output: t,y contain the computed mesh and solution; stats (optional) is a vector containing the following time statistics: 1. the pre-processing time for computing the parameters h 1 , r, and N (see Section 3.1) and the fractional itegrals (47), and (48) or (49); 2. the time for solving the problem; 3. the pre-processing time for computing the fractional itegrals (47), and (48) or (49) for the error estimation; 4. the time for solving the problem on the doubled mesh, for the error estimation; err (optional), if specified, contains the estimate of the absolute error.
This estimate, obtained on a doubled mesh, is relatively costly: for this reason, when the parameter is not specified, the solution on the doubled mesh is not computed.
For the first two problems, we shall also make a comparison with the Matlab © code flmm2 [27], 7 in order to emphasize the potentialities of the new code.All numerical tests have been done on a M2-Silicon based computer with 16GB of shared memory, using Matlab © R2023b.
The comparisons will be done by using a so called Work Precision Diagram (WPD), where the execution time (in sec) is plotted against accuracy.The accuracy, in turn, is measured through the mixed error significant computed digits (mescd) [39], defined, by using the same notation seen in ( 37), as being t i , i = 0, . . ., N , the computational mesh of the considered solver, and y(t i ) and ȳi the corresponding values of the solution and of its approximation.
Figure 2 contains the obtained results: as one may see, flmm2 reaches less than 12 mescd, since by continuing reducing the stepsize, at the 15-th mesh doubling the error starts increasing.On the other hand, the execution time essentially doubles at each new mesh doubling.Conversely, fhbvm can achieve full machine accuracy by employing a uniform mesh with stepsize h 1 = 1/M , with M very small, thus using very few mesh points.As a result, fhbvm requires very short execution times.
Figure 3 contains the obtained results, from which one deduces that flmm2 achieves about 5 mescd (with a run time of about 85 sec), whereas fhbvm has approximately 13 mescd, with an execution time of about 1 sec.Further, in Figure 4, we plot the true and estimated (absolute) errors for fhbvm in the case M = 10 (corresponding to a computational mesh made up of 251 meshpoints, with an initial stepsize h 1 ≈ 7.3 • 10 −12 , and a final stepsize h 250 ≈ 2): as one may see from the figure, there is a substantial agreement between the two errors.

Example 3
We now consider the following nonlinear problem [8]: having solution This problem is relatively simple and, in fact, both flmm2 and fhbvm solve it accurately.We use it to show the estimated error by using fhbvm with parameter M = 2, which produces a graded mesh with 41 mesh-points, with h 1 ≈ 1.8 • 10 −12 and h 40 ≈ 0.49.The absolute errors (true and estimated) for each component are depicted in Figure 5, showing a perfect agreement for both of them.In this case, the evaluation of the solution requires ≈ 0.04 sec, and the error estimation requires ≈ 0.11 sec.

Example 4
At last, we consider the following fractional Brusselator model:  By solving this problem using fhbvm with parameter M = 5, a graded mesh of 46 points is produced, with h 1 ≈ 6.1 • 10 −5 and h 45 ≈ 0.98.The maximum estimated error in the computed solution is less than 3.5 • 10 −13 , whereas the phase-plot of the solution is depicted in Figure 6.
In this case, the evaluation of the solution requires ≈ 0.04 sec, and the error estimation requires ≈ 0.14 sec.

Conclusions
In this paper we have described in full details the implementation of the Matlab © code fhbvm, able to solving systems of FDE-IVPs.The code is based on a FHBVM (22,20) method, as described in [8].We have also provided comparisons with another existing Matlab © code, thus confirming its potentialities.In fact, due to the spectral accuracy in time of the FHBVM (22,20) method, the generated computational mesh, which can be either a uniform or a graded one, depending on the problem at hand, requires relatively few mesh points.This, in turn, allows to reduce the execution time due to the evaluation of the memory term required at each step.We plan to further develop the code fhbvm, in order to provide approximations at prescribed mesh points, as well as to allow selecting different FHBVM(k, s), k ≥ s, methods [8].At last, we plan to extend the code to cope with values of the fractional derivative, α, greater than 1.
Declarations.The authors declare no conflict of interests, nor competing interests.No grants were received for conducting this study.

2
Here, P (a,b) j (x) denotes the j-th Jacobi polynomial with parameters a and b, in [−1, 1].