Two-step collocation methods for fractional differential equations

We propose two-step collocation methods for the numerical solution of fractional differential equations. These methods increase the order of convergence of one-step collocation methods, with the same number of collocation points. Moreover, they are continuous methods, i.e. they furnish an approximation of the solution at each point of the time interval. We describe the derivation of two-step collocation methods and analyse convergence. Some numerical experiments confirm theoretical expectations.


1.
Introduction. In the last decades scientists paid much attention on Fractional Differential Equations (FDEs), since they are suitable to model a wide class of natural phenomena with memory, like the behavior of viscoelastic materials in mechanics [14,31,38], anomalous diffusion in transport dynamics of complex systems [2,27,33], Brownian motion [32,34], rheology [19].
Several numerical methods have been proposed so far for the solution of (1.1), as for example the product-integration methods (mainly in the predictor corrector framework) [17,18,20,21], Adomian decomposition methods [12], spectral methods [25], fractional linear multistep methods [22,30], collocation methods [1,35,36,37]. Here we focus our attention on collocation methods, whose main properties in the numerical solution of functional equation are high order of convergence, good stability, uniform approximation of the solution, the possibility of choosing a suitable collocation function basis. For these reasons, they have been largely applied for ordinary differential equations [23,24] and Volterra integral equations [6,7]. The first naive formulation of collocation methods for FDEs was given by Blank [1]. Later on Rawashdeh in [37] applied collocation methods to fractional integro-differential equations, without analyzing convergence. Pedas and Tamme gave the first rigorous foundation of collocation methods for linear and nonlinear FDEs [35,36]. In [36] it was proved that the error of collocation methods based on piecewise polynomials of degree m − 1 decay as O(N −m ), when a specially suited mesh of N points is adopted, called graded mesh. The use of graded mesh (already introduced for the solution of Abel-Volterra integral equations [5,6]) aims to recover the order of convergence of collocation methods when the solution is not enough smooth at the origin.
Our purpose is to obtain higher order collocation methods. A successful strategy to increase the order of convergence consists of applying multivalue methods, as done in the case of ordinary differential equations [13,28,29], where the multistep collocation methods were modified in [13] by also using the inherent quadratic stability technique [3,4]. Multivalue methods were also adopted in Volterra integral equations [11] and Volterra integro-differential equations [8]. Multivalue methods introduce additional information from the past, to compute the solution in the current interval and therefore the computational cost is not increased. Following this idea, we construct two step collocation methods based on piecewise polynomials of degree 2m − 1, which achieve an error decay of type O(N −2m ), on a suitable graded mesh.
The paper is structured as follows. In Sec. 2 we recall the main results on the existence, uniqueness and smoothness of the solution of (1.1). In Sec. 3 we illustrate two-step collocation methods on a graded mesh and give some details about the computational aspects. The analysis of convergence of the proposed method is carried out in Sec. 4. Some numerical experiments are reported in Sec. 5. Finally some concluding remarks complete the paper.
2. Existence, uniqueness and smoothness of solution. For sake of completeness, we report the following results on the existence, uniqueness and smoothness of solution of problem (1.1). Here α is the least integer greater or equal to α.
   and f : D → R be continuous and fulfil a Lipschitz condition with respect to y. Then there exists a unique solution of the IVP fractional problem (1.1) y : To analyze convergence of collocation methods, it is necessary to investigate thoroughly the behavior of y and of its derivatives, which may have singularities at the initial point t = 0. Therefore, it is useful to introduce the space C q,ν (0, b] (compare [7,35,36]).
Let C q,ν (0, b] with q ∈ N and ν ∈ (−∞, 1), be the space of continuous functions y : [0, b] → R, which are q times continuously differentiable in (0, b] such that: Under suitable hypotheses on f (t, y), the solution y of the problem (1.1) and its fractional derivative D α y belong to the space C q,ν (0, b].
The following theorem deals with the smoothness of the solution in the weighted space C q,ν (0, b].
such that ∀i, j ∈ N with i + j ≤ q and the following estimation holds: For α ∈ (0, 1) assume also that The function φ : [0, ∞) → R is assumed to be monotonically increasing. Suppose that the fractional IVP (1.1) possesses a solution y ∈ C[0, b] with D α y ∈ C[0, b]. Then y ∈ C q,ν (0, b] and D α y ∈ C q,ν (0, b]. 3. Two-step collocation methods. We observe that from (1.1) we have Then, by setting z = D α y, we have and J 0 y := y. Therefore the problem (1.1) is equivalent to with the grading exponent r ∈ R, r ≥ 1. Let We will look for an approximation v of the solution of (3.3) in the space of piecewise polynomial functions S where π k is the space of algebraic polynomials of degree k. From the definition, we note that v ∈ S We introduce the collocation parameters with (η 1 , η m ) = (0, 1), which lead to the set of collocation points We define the approximate solution of (3.3) as the piecewise polynomial v ∈ S where the functions v j (t) are constructed as follows. For j = 2, . . . , N , v j is a polynomial of degree 2m − 1 defined by the following conditions ) is a collocation condition at the collocation points t jk ∈ σ j , while equation (3.10b) is an interpolation condition at the collocation points of the previous interval where {L k (·)} k=1,...,2m are the fundamental Lagrange polynomials corresponding to the nodes {t j−1,k , t j,k |k = 1, . . . , m}, i.e.
We introduce a starting procedure to construct a function v 1 , which preserves the accuracy of the numerical scheme and vanishes outside the interval σ 1 . In particular, witht k =η k h 1 and 0 ≤η 1 < · · · <η 2m ≤ 1. Other choices of v 1 can be made, see Remark 2.
It is easy to verify that We observe that the interpolation condition (3.10b) is automatically satisfied by (3.11). Then, by imposing the collocation condition (3.10a), we obtain the nonlinear system Since (J α L λµ )(t jk ) = 0 when λ > j, the system can be written as follows Now we describe how to compute the fractional integrals appearing in (3.15). We define By simple calculations we get Thus {ϕ jµ (τ )} µ=1,...,2m coincide with the fundamental Lagrange polynomials corresponding to the nodes It will be useful to express these polynomial in this way In (3.15) the fractional integrals can be computed by using these equalities As a matter of fact, to compute (J α L jµ )(t jk ), we change the integration variable s = t j−1 + τ h j , use (3.14), (3.18) and (3.19), so we get Then, in the integral (J α L λµ )(t jk ), we change the integration variable s = t λ−1 + τ h λ , use (3.14) , (3.18) and (3.20), and obtain In (3.12) and (3.15), (J α v 1 )(t jk ) can be computed analogously, since we assumed that v 1 is a piecewise polynomial. For other choices of v 1 , a numerical quadrature might be necessary.
Here we illustrated two-step collocation methods for a scalar problem (1.1). In a multidimensional case, some algorithmic strategies used in [10] may be adapted, for an efficient implementation on a distributed-memory architecture.
4. Convergence analysis. In this section we investigate convergence and order of convergence of two-step collocation method as h → 0. To carry out this analysis we will need the definition of Fréchet differentiability and the following variant of a lemma on nonlinear equations, compare [36,Lemma 4.1].
Let (X, · ) be a Banach space and L(X, X) the Banach space of bounded linear operators A : X → X with the norm A L(X,X) := {sup Az : z ∈ X, z < 1}. T : B → X is said to be Fréchet differentiable at z 0 ∈ B if there exists a linear operator T (z 0 ) ∈ L(X, X) such that In this case T (z 0 ) is the (unique) Fréchet derivative of T at z 0 .
Lemma 4.1. Let the following conditions be fulfilled: 1. the equation (4.1) has a solution z * ∈ B, and the operator T is Fréchet differentiable at z * ; 2. there exists δ > 0 such that the operators T N , N ∈ N, are Fréchet differentiable in the ball z − z * ≤ δ, which is assumed to be contained in B; for any > 0 there exists a δ ∈ (0, δ] such that for every N ∈ N, compact;the homogeneous equation z = T (z * )z has in X only the trivial solution.
Then there exist N 0 ∈ N and δ 0 ∈ (0, δ] such that the equation (4.2) has for N ≥ N 0 a unique solution z N in the ball z − z * ≤ δ 0 . Thereby where c > 0 does not depend on N .
In the following lemma we prove that the error of the interpolation operator can be bounded by E N (2m, ν, r), where Proof. We have The application of [36,Lemma 4.2] for α = 0 yields ||u − P N u|| σ1 ≤ c 1 E N (2m, ν, r). If v is a polynomial of degree 2m − 1, Therefore: If Λ m is the Lebesgue constant corresponding to the Lagrange fundamental polynomials {ϕ jµ } 2m µ=0 , we have: This inequality and (4.13) yield: (4.16) Let v be the Taylor polynomial of degree 2m − 1 of function u. Then, by Lemma 4.2 applied to (4.16), we get: Now the thesis follows from (4.12), (4.17) and from the following simple inequalities: It can be verified that Λ m = Λ m2 , which is independent from N . As a matter of fact the maximum of {Λ mj } j is obtained when the ratio h j−1 /h j is minimum, and moreover We recall that L(X, Y ), is the Banach space of the linear bounded operators Proof. From (4.14) we have that Moreover whereφ k (τ ), k = 1, . . . , 2m, are the fundamental Lagrange polynomials corresponding to the nodesη 1 , . . . ,η 2m . Then whereΛ m is the Lebesgue constant corresponding to {φ k (τ )} 2m k=1 . Thus the thesis is proved.
If v is the solution of the two-step method defined by (3.9), (3.10) and (3.12), then v = T N v. As a matter of fact, for j = 2, . . . , N , it results that T N v |σj = p j |σj where (cfr. (4.10)) . . , m. Now, by applying the interpolation and collocation conditions (3.10), we have By applying the starting collocation conditions (3.12), it follows p 1 | σ1 = v 1 | σ1 . Therefore, the two-step collocation method may be written in an operator equation representation v = T N v.
If, in addition, the assumptions of Th. 2.4 with q := 2m and ν ∈ [1 − α, 1) are fulfilled, then for all N ≥ N 0 the following error estimate holds: N (2m, ν, r), with y N given by the formula (3.6).Here c is a constant not depending on N , r ∈ [1, ∞) is the grading exponent in (3.4) and E N is defined by (4.11).
Proof. We recall that the Fréchet derivative of T has been computed in [36,Th. 4.1] and T (z 0 ) ∈ L(L ∞ (0, b), P N is linear, therefore we have From Lemma 4.4 we obtain The right hand side expression goes to zero as z−z 0 → 0, thus we have From (3.1), (3.6) and (4.19), finally we get Remark 2. Theorem 4.5 holds also for starting procedures different from the one proposed in Sec. 3 for the computation of v 1 . As a matter of fact, it is sufficient that z − v 1 σ1 ≤ cE N (2m, ν, r) for some constant c > 0, and that the operator P N is defined in such a way that (P N v)| σ1 = v 1 and Lemma 4.3 holds.
For the one-step collocation methods proposed in [36], the numerical solution y N is obtained by collocation conditions similar to (3.10a), and require the solution of a nonlinear system of dimension m at each time-step. The error of this method is y N −y ∞ ≤ cE N (m, ν, r), for an arbitrary set of the collocation parameters. Thus, for a suitable choice of the grading exponent r, the error of the one-step collocation methods goes to zero as N −m , while for the two-step methods the error decays as N −2m , while for both the dimension of the nonlinear system to solve at each time-step is m.

Numerical experiments.
In this section we illustrate the performances of the two-step methods on some test problems from the literature, both to confirm the predicted order of convergence and for comparison with one-step collocation methods. The numerical experiments are carried out in Matlab and the nonlinear systems arising from the collocation conditions are solved by the Matlab routine fsolve.
In the first set of experiments we verify the order of convergence of two-step methods for m = 2 and m = 3. In Tables 1 and 2   where cd is the number of correct digits of the solution (the error is written as 10 −cd ). From the Tables 1 and 2 it is seen that the theoretical order of convergence is confirmed in any case.
In the second set of experiments we compare the two-step collocation method having m collocation parameters with the one-step collocation methods with m and 2m collocation parameters, the last one having the same order of convergence. In Figures 1 and 2 we show the work-precision diagrams, where on the x-axis there is the number of correct digits cd of the solution, while on the y-axis there is the number of evaluations of function f . The two-step methods have the best performances as expected, except when low accuracy is required. 6. Conclusions. We introduced two-step collocation methods for FDEs, which have an higher convergence order with respect to one-step collocation methods by using at each time-interval the known information from the previous time-interval, without increasing the dimension of the nonlinear system to solve. We carried out convergence analysis of the methods, considering a suitable function space which takes into account the behavior of the solution derivatives at the origin. Then we showed the performances of the methods by means of numerical experiments on some test equations from the literature, which also prove that the two-step methods compare favourably with respect to the one-step collocation methods.
Future developments may regard the investigation of a different collocation basis, as for example a trigonometric polynomial basis, possibly adapting some techniques from the exponential fitting [26], as done in [9].