Spline Collocation for Multi-Term Fractional Integro-Differential Equations with Weakly Singular Kernels

We consider general linear multi-term Caputo fractional integro-differential equations with weakly singular kernels subject to local or non-local boundary conditions. Using an integral equation reformulation of the proposed problem, we first study the existence, uniqueness and regularity of the exact solution. Based on the obtained regularity properties and spline collocation techniques, the numerical solution of the problem is discussed. Optimal global convergence estimates are derived and a superconvergence result for a special choice of grid and collocation parameters is given. A numerical illustration is also presented.

However, when working with problems stemming from real-world applications, it is only rarely possible to find the solution of a given fractional differential or integral equation in closed form, and even if such an analytic solution is available, it is typically too complicated to be used in practice. Therefore, in general, numerical methods are required for solving fractional differential and integral equations. As a consequence, the last few decades have witnessed a steadily increasing development and analysis of numerical methods for fractional differential equations, of which a good deal are concerned with the numerical solution of initial or boundary value problems with one fractional derivative in the equation, see, for example, the works [21][22][23][24][25] for initial value problems and [26][27][28][29][30][31] for boundary value problems. A comprehensive survey of the most important methods for fractional initial value problems, along with a detailed introduction to the subject and a brief summary about the convergence behaviour of the methods is given in the monograph [32]. Less attention has been paid to the numerical solution of equations with multiple fractional derivatives (the so-called multi-term equations) [33][34][35][36][37][38][39][40][41], fractional differential equations with non-local boundary conditions [42][43][44][45][46] and fractional integro-differential equations with weakly singular kernels [47][48][49]. This motivated us in the present paper to focus on constructing effective numerical methods for fractional weakly singular integro-differential equations with local or non-local boundary conditions. In order to construct high-order numerical methods for solving fractional differential and integro-differential equations, one needs some information about the regularity of the exact solution of the underlying problem. This becomes even more significant if we wish to study the optimal order of convergence of the proposed algorithms. However, fractional differential and integral equations pose an extra challenge compared to integer order differential equations. For example, it is well known that, in the case of integer order differential equations, the smoothness properties of a solution are determined by certain assumptions on the given data. Very simple examples show (see also Theorem 1 below) that, in general, we can not expect this to be true for fractional differential equations-a solution of a fractional differential equation will generally exhibit non-smooth behaviour even in the case of smooth data. Thus, when constructing high order numerical methods for fractional differential and integral equations, one should take into account, in some way, the possible non-smooth behaviour of an exact solution. Numerical methods which assume smooth solutions for fractional differential equations are valid only for a tiny subclass of problems, as is made clear in [50,51]. In the numerical solution of weakly singular integral equations of the second kind by collocation methods the possible non-smooth behaviour of the exact solution of the underlying problem can be taken into account by using special non-uniform grids (see, for example, [52][53][54][55]). Similar ideas for solving certain types of fractional differential equations have been successfully used in [28,36,39,48,56].
The main purpose of the present paper is to extend the results of [39,48] and the corresponding results of [56] to a much larger class of linear multi-term fractional equations, where the terms can be Caputo fractional derivatives of arbitrary order and weakly singular integrals involving Caputo fractional derivatives (for the exact problem setting see Section 3 below). To the best of our knowledge, up to now there has been no discussion about the numerical solution of such equations. In the present paper, using special non-uniform grids reflecting the possible singular behaviour of the exact solution, we construct a high order collocation-based numerical method for these equations. Our analysis hinges on a smoothness result for the exact solution of the underlying problem, given by Theorem 1. The theoretical convergence results for the proposed algorithm are proved by Theorems 2 and 3. These theorems show how to choose grid and collocation parameters so that the best possible order of convergence is attained. Notably, we show that with a careful choice of collocation parameters and with a suitably graded grid, the constructed method will obtain a superconvergence rate on the whole interval of integration [0, b]. We also provide numerical results to verify the theoretical analysis.
Without loss of generality, in this paper we shall assume that the starting point for fractional operators is located at the origin 0.

Basic Notations and Definitions
By N we denote the set of all positive integers {1, 2, . . . }, by N 0 the set of all nonnegative integers {0, 1, 2, . . . } and by R the set of all real numbers (−∞, ∞). By α we denote the smallest integer greater or equal to a real number α.
Let b > 0 be a fixed real number. By L 1 (0, b) we denote the Banach space of measurable functions u : where µ(Ω) is the Lebesgue measure of the set Ω. By C[0, b] we denote the Banach space of continuous functions u : By C m (∆) we denote the set of m times (m ∈ N 0 , for m = 0 we set C 0 (∆) = C(∆)) continuously differentiable functions on ∆, with Next, we recall the definitions and some properties of Riemann-Liouville integral and Caputo fractional differential operators, see [1,2].
For given δ ∈ (0, ∞) by J δ we denote the Riemann-Liouville integral operator of order δ, defined by where Γ(δ) := ∞ 0 e −s s δ−1 ds is the Euler gamma function. For δ = 0 we set J 0 := I, where I is the identity operator. If δ > 0, then (J δ u)(t) exists for almost all t ∈ [0, b] and the function J δ u is also an element of L 1 (0, b). We have for any u ∈ L ∞ (0, b) that Moreover, it is well known (see, e.g., [57]) that the operator J δ (δ > 0) is linear, bounded and compact from L we denote the Taylor polynomial of degree m − 1 for the function u ∈ C m−1 [0, b] at the point 0: By D δ Cap we denote the Caputo fractional differential operator of order δ > 0, defined by the formula In the definition (4) we assume that u ∈ C n−1 [0, b] is such that the integer order derivative However, this is not a necessary condition. In [58], Vainikko gives a comprehensive description of In particular, he has derived necessary and sufficient conditions for the existence of D δ where z := D δ Cap u and c λ ∈ R (λ = 0, . . . , δ − 1) are some constants. Note also that for where t ≥ 0.

Problem Setting
In the present paper, we will consider a general class of linear multi-term fractional weakly singular integro-differential equations in the form subject to the conditions Here y = y(t) is the unknown function and by D α i Cap y (i = 0, . . . , p) and D θ j Cap y (j = 0, . . . , q) are denoted the Caputo fractional derivatives of y of orders α i and θ j , respectively. We assume that the following conditions (10) are fulfilled: n 0 , n 1 ∈ N 0 , n 0 < n, n 1 < n, β ij0 , β ijk , β i , γ i ∈ R, Note that for certain values of coefficients β ij0 , β ijk and β i the problem (8)-(10) takes the form of an initial value problem or a multi-point boundary value problem.
Following some ideas of [39,48,56], we will construct a high-order method for the numerical solution of problem (8) and (9). We first introduce an integral equation reformulation of the underlying problem (Section 4) and prove some results about the existence, uniqueness and regularity of the exact solution of (8) and (9) (Section 5). Using this information and spline collocation techniques, the numerical solution of the problem is discussed (Section 6). Optimal global convergence estimates are derived and a superconvergence result for a special choice of collocation parameters is established (Section 7). Finally, numerical illustrations confirming the convergence estimates are given (Section 8).
In the sequel, we assume that the matrix M is regular. Observe that the matrix M is regular if and only if from all polynomials y of degree n − 1 only y = 0 satisfies the homogeneous conditions corresponding to the conditions (9) by γ i = 0, i = 0, . . . , n − 1. Indeed, substituting (11) with z = 0 into (14) we obtain a homogeneous system of algebraic equations with respect to c 0 , c 1 , . . . , c n−1 . This system coincides with (13) by γ i = 0 (i = 0, . . . , n − 1) and z = 0. Therefore, the homogeneous system corresponding to (13) has only the trivial solution c 0 = c 1 = · · · = c n−1 = 0 (and thus M is regular) if and only if from all polynomials y of degree n − 1 only y = 0 satisfies (14). Let be the inverse of M. Using M −1 , the solution of the system (13) can be written in the form where Therefore, a function y in the form (11) satisfies the conditions (9) if and only if it can be expressed by the formula where Suppose now that y ∈ C n−1 [0, b] is a solution of problem (8) and (9) (18) and (19), respectively. Inserting (17) into (8), we see that Cap y satisfies the integral equation where t ∈ [0, b]. Conversely, it turns out that if z ∈ C[0, b] is a solution of Equation (20) then y defined by (17) belongs to C n−1 [0, b] and is a solution to (8) and (9). In this sense, Equation (20) is equivalent to the problem (8) and (9).

Existence, Uniqueness and Smoothness of the Solution
In order to characterize the possible singular behaviour of a solution of a fractional differential equation, we introduce a weighted space C m,ν (0, b] of smooth functions on (0, b] (cf. [57,59]).
For given b ∈ R, b > 0, m ∈ N and ν ∈ R, ν < 1, by C m,ν (0, b] we denote the set of continuous functions u : where c is a positive constant. In other words, where, for t > 0, Equipped with the norm Moreover, a function of the form Observe that as ν increases so do the possible singularities of the derivatives of the functions in C q,ν (0, b]. Next we introduce some auxiliary results. Their proofs can be found in [57].
, m ∈ N, ν < 1, then y 1 y 2 ∈ C m,ν (0, b], and with a constant c which is independent of y 1 and y 2 .

Lemma 2.
Let η ∈ R, η < 1 and let K ∈ C(∆). Then operator S defined by is compact as an operator from L The following theorem characterizes the existence, uniqueness and regularity properties of the solution of (8) and (9). Theorem 1. (i) Suppose that assumptions (10) hold. Moreover, assume that problem (8) and (9) with f = 0 and γ i = 0 (i = 0, . . . , n − 1) has in C[0, b] only the trivial solution y = 0, and from all polynomials y of degree n − 1 only y = 0 satisfies the conditions (14).

Numerical Method
In order to take into account the potential non-smoothness of the exact solution y = y(t) of (8) and (9) at the origin t = 0, we introduce on the interval [0, b] a graded grid Π N (N ∈ N). More precisely, let N ∈ N, then Π N := {t 0 , . . . , t N } is a partition (a graded grid) of the interval [0, b] with the grid points where the grading exponent r ∈ R, r ≥ 1. If r = 1, then the grid points (32) are distributed uniformly; for r > 1 the points (32) are more densely clustered near the left endpoint of the interval [0, b].
Next, for a given integer k ≥ 0, by S (−1) k (Π N ) we denote the standard space of piecewise polynomial functions: and π k denotes the set of polynomials of degree not exceeding k. Note that the elements of S (−1) k (Π N ) may have jump discontinuities at the interior points t 1 , . . . , t N−1 of the grid Π N . Let m ∈ N. In every subinterval [t j−1 , t j ] (j = 1, . . . , N) we define m collocation points t j1 , . . . , t jm by where η 1 . . . , η m are some fixed parameters which do not depend on j and N and satisfy Approximations to the solution z of (20) we find by collocation conditions z N (t jk ) = (Tz N )(t jk ) + g(t jk ), k = 1, . . . , m, j = 1, . . . , N, where T, g and {t jk } are defined by (25), (26) and (33), respectively. Note that conditions (35) for finding z N ∈ S (−1) m−1 (Π N ) have an operator equation representation where P N = P N,m : If η 1 = 0, then by (P N v)(t j1 ) we denote the right limit lim t→t j−1 ,t>t j−1 (P N v)(t). If η m = 1, then by (P N v)(t jm ) we denote the left limit lim t→t j ,t<t j (P N v)(t). The collocation conditions (35) form a system of equations whose exact form is determined by the choice of a basis in the space S (−1) m−1 (Π N ). If η 1 > 0 or η m < 1 then we can use the Lagrange fundamental polynomial representation: where, for µ = 1, . . . , m, λ = 1, . . . , N, Then, z N ∈ S (−1) m−1 (Π N ) and z N (t jk ) = c jk , k = 1, . . . , m, j = 1, . . . , N. Substituting z N in the form (38) to (35), we obtain a system of linear algebraic equations with respect to the coefficients {c jk }: Solving this system of equations, we obtain the coefficients {c jk } and thus have found the approximation z N in the form (38). Note that for the computation of (TΦ λµ )(t jk ) we need to find the weakly singular integrals J δ Φ λµ (t jk ) (δ > 0), which can be found exactly. Approximation y N to y we find by the formula (cf. (17)) where G and Q are defined by (18) and (19), respectively. By substituting z N in the form (38) into (41), we obtain the following expression for the approximate solution y N :

Convergence Analysis
In this section, we study the convergence and convergence order of the proposed algorithms. For this we need Lemmas 3-6 below. The proofs of Lemmas 3-5 follow from the results of [57,59]. The proof of Lemma 6 can be found in [60]. In what follows, for Banach spaces E and F, by L(E, F) is denoted the Banach space of linear bounded operators A : E → F with the norm A L(E,F) = sup{ Au F : u ∈ E, u E ≤ 1}.  m−1 (Π N ) (N ∈ N) be defined by (37). Then where r ∈ [1, ∞) is the grading exponent in (32) and c is a positive constant independent of N.
Then, problem (8) and (9) has a unique solution y ∈ C n−1 [0, b] such that D m−1 (Π N ), determining by (41) a unique approximation y N to y, the solution of (8) and (9), and (ii) If, in addition to (i), we assume that . . . , q), where µ ∈ R, µ < 1, then for all N ≥ N 0 the following error estimate holds: where ν is determined by the Formula (30), r ≥ 1 is the grading parameter in (32) and c is a constant independent of N.
Proof. (i) First, we prove the convergence (46). To this end, we need to show that equation z = Tz + g (see (20)), with T and g given by (21) and (22), is uniquely solvable in L ∞ (0, b). We observe that T is compact as an operator from where c is a constant independent of N. Thus, for N ≥ N 0 , Equation (36) provides a unique solution z N ∈ S −1 m−1 (Π N ). Note that for z, the solution of equation z = Tz + g, it holds P N z = P N Tz + P N g. We have for z and z N that Therefore, by (48), where c is a positive constant independent of N. Using (41) we obtain that From this and (49) it follows that where c and c 1 are some positive constants independent of N. This together with z ∈ C[0, b] and Lemma 3 yields the convergence (46). (ii) If K ∈ C m (∆), h, f ∈ C m,µ (0, b], m ∈ N, µ ∈ R, µ < 1, then it follows from the part (ii) of Theorem 1 that z ∈ C m,ν (0, b], with ν given by (30). This together with (51) and Lemma 5 yields the estimate (47).
It follows from Theorem 2 that in case of sufficiently smooth f , d i (i = 0, . . . , p − 1) and K i (i = 0, . . . , q), by using sufficiently large values of r for every choice of collocation parameters 0 ≤ η 1 < · · · < η m ≤ 1 a convergence of order O(N −m ) can be expected. From Theorem 3 below, we see that by a careful choice of parameters η 1 , . . . , η m and by a slightly more restrictive smoothness requirement on f , d i (i = 0, . . . , p − 1) and K i (i = 0, . . . , q) it is possible to establish a faster convergence.
Proof. From Theorem 2, we know that problem (8) and (9) has a unique solution y ∈ Cap y ∈ C[0, b] and there exists an integer N 0 such that for all N ≥ N 0 Equation (36) has a unique solution z N for which (49) is valid. Denotê where T and g are defined by (21) and (22), respectively. From (36) we see that z N = P NẑN . Substituting this expression of z N into (55) we obtain that From (20) and (56) follows the identity we obtain, with the help of (48), that This together with (25) yields where with L i defined by (23). Here and below c, c 1 , c 2 and c 3 are generic positive constants which are independent of N. With the help of (3) and from the boundedness of the Riemann-Liouville integral operator we obtain for N ≥ N 0 the following estimates: It follows from (57)-(61) that where α * = min{α p − α p−1 , α p − n 1 }. Since z N = P NẑN , we obtain with the help of (17) and (41) that Using (18), (59) and (60) we obtain This together with (62) and (63) yields Because of Theorem 1 we have z ∈ C m+1,ν (0, b] and due to (3) for α * ≥ 1. Therefore, it follows from (64) and Lemma 6 that the estimates (53) and (54) are true.
In the case m = 3, it follows from (53) with α 1 = 1 2 and ν = 3 4 that, for sufficiently large N, where c 1 is a positive constant independent of N. Due to (70), the ratios Θ N for r = 3, r = 4 and r ≥ 14 3 ought to be approximately 2 2.25 ≈ 4.76, 2 3 = 8.00 and 2 3.5 ≈ 11.31, respectively. These values are given in the last row of Table 2.
As we can see from Tables 1 and 2, the numerical results are in good accord with the theoretical estimates given by Theorem 3 (by estimates (69) and (70)).
for which the condition (ii) in Theorem 3 is not satisfied. The results of the numerical experiments are presented in Table 3, with the ratios (68) corresponding to the estimate (69) given in the last row of the table. As we can see, although the numerical convergence rate for r = 1 and r = 2 is still in accord with Theorem 3, for r = 3 and r = 10 3 the numerical results do not attain the predicted superconvergence rate. Instead, the highest attained numerical convergence rate is close to 2 2 = 4, which is the maximal convergence rate predicted by Theorem 2. Table 3. Numerical results for problem (65) and (66) with m = 2 and η 1 = 0.1, η 2 = 0.9.
In the case m = 2 it follows from (54) with α 1 = 11 10 and ν = 9 10 that, for sufficiently large N, where c 0 is a positive constant independent of N. Due to (73), the ratios Θ N for r = 1, r = 2 and r ≥ 30 11 ought to be approximately 2 1.1 ≈ 2.14, 2 2.2 ≈ 4.60 and 2 3 = 8, respectively. These values are given in the last row of Table 4.
As we can see from Tables 4 and 5, the numerical results are in accord with the theoretical estimates given by Theorem 3 (by estimates (73) and (74)).

Concluding Remarks
In this work, we have introduced and analyzed a high order numerical method for solving a wide class of linear multi-term fractional weakly singular integro-differential equations with Caputo fractional derivatives for local or non-local boundary conditions. For certain values of coefficients, the considered problem is an initial value problem or a multi-point boundary value problem. We have reformulated the proposed problem as an integral equation with respect to the highest order Caputo derivative in the fractional integro-differential equation. Using this reformulation, we have first studied the existence, uniqueness and regularity of the problem. We have shown that, in general, the exact solution of the problem is non-smooth, even if the initial data is smooth. On the basis of the results obtained regarding the smoothness of the exact solution, with the help of special graded grids and spline collocation techniques we have constructed an effective numerical method that recovers its optimal convergence order. Moreover, we have shown that, by a careful choice of grid and collocation parameters, the method will obtain a global superconvergence rate. Note that (see [45,56]) the superconvergence advantage might not hold if the reformulated integral equation is obtained with respect to the exact solution of the original problem. The main conclusions of the paper extend known ones and are formulated in Theorems 1-3. From these results, an optimal choice of grid and collocation parameters can be made for finding a numerical solution to any problem given in the form (8) and (9) that satisfies the conditions of Theorems 1-3. More precisely, in order to apply Theorems 2 or 3, we first of all have to calculate the smoothness parameter ν, characterizing the regularity of the exact solution of the underlying problem. After finding the value of ν and choosing an order m − 1 for the polynomials used in the piecewise polynomial method, one can select a grid parameter r large enough so that the corresponding conditions set for the convergence rates in Theorems 2 or 3 are fulfilled. Note that while the conditions given by Theorems 2 and 3 do not set an upper bound for the parameter r, for the optimal convergence rate the smallest possible value of r satisfying the corresponding inequality tends to give a more precise result.