Two Collocation Type Methods for Fractional Differential Equations with Non-Local Boundary Conditions∗

A class of non-local boundary value problems for linear fractional differential equations with Caputo-type differential operators is considered. By using integral equation reformulation of the boundary value problem, we study the existence and smoothness of the exact solution. Using the obtained regularity properties and spline collocation techniques, we construct two numerical methods (Method 1 and Method 2) for finding approximate solutions. By choosing suitable graded grids, we derive optimal global convergence estimates and obtain some super-convergence results for Method 2 by requiring additional assumptions on equation and collocation parameters. Some numerical illustrations for verification of theoretical results is also presented.


Introduction
Fractional differential equations arise in various areas of science and engineering and have been proven to be a valuable tool in modelling various phenomena in physics, astrophysics, chemistry, geology, bioengineering, medicine, atmospheric science, material science, optics, mechanics and many other fields. For a more comprehensive list of applications of fractional calculus to science and engineering, we refer to [7] and the references cited therein. The mathematical aspects of fractional differential equations and various numerical methods for such equations are studied in the monographs [7,11,16,27]. A great deal of papers have been devoted to the numerical solution of initial value problems for fractional differential equations -some more recent results can be found in [6,18,21,23]. A smaller but growing number of papers concern the numerical solution of boundary value problems (see, e.g., [1,4,13,14,15,22,24,25,26,30]). However, these papers are only concerned with local boundary value conditions. While non-local boundary conditions are widely researched for differential equations with integer order (see, e.g., [9,10]), less attention has been paid for fractional differential equations with non-local boundary conditions. We refer to papers [2,3,5,29], which are concerned with various existence results for fractional boundary value problems with non-local conditions. It is known that we usually cannot expect the solution of a fractional differential equation to be smooth on the whole interval of integration (see, e.g., [22,23]). In collocation methods the singular behaviour of the exact solution can be taken into account by using polynomial splines and suitable graded grids (see, e.g., [21,24]). In the case of integral equations and integer order differential equations, a similar approach has been exploited by many authors (see, e.g., [12,17,19,20]).
In the present paper we construct two high-order collocation type methods for the numerical solution of the fractional differential equation with non-local boundary values in the following form: where b 1 , b 2 , b 3 ∈ (0, b], γ 0 , γ 1 , γ 2 , γ 3 , γ ∈ R := (−∞, ∞) and D α * is the Caputo differential operator of order α. We assume that 0 < β ≤ α < 1 and h, f ∈ C[0, b]. By C[0, b] we denote the Banach space of continuous functions u : [0, b] → R with the norm u ∞ = sup{|u(t)| : 0 ≤ t ≤ b}. The Caputo differential operator D δ * of order δ ∈ (0, 1) can be defined by formula (see, e.g, [11]) (D δ * y)(t) := (D δ [y − y(0)])(t), 0 ≤ t ≤ b. Here D δ y is the Riemann-Liouville fractional derivative of y : with J δ , the Riemann-Liouville integral operator, defined by where I is the identity mapping and Γ is the Euler gamma function. In the definition of D δ * y we assume that y ∈ C[0, b] and J 1−δ [y − y(0)] ∈ C 1 [0, b]. It is well known (see, e.g. [8]) that J δ , δ > 0, is linear, bounded and compact as an operator from L ∞ (0, b) into C[0, b], and we have for any y ∈ L ∞ (0, b) that (see, e.g. [16]) The rest of our paper is arranged as follows. In Section 2 we reformulate the problem (1.1)-(1.2) into two different integral equations. In Section 3, with the help of results obtained from the previous section, we study the existence and regularity of the exact solution to problem (1.1)-(1.2). Next, by using the integral equations obtained in Section 2, we construct two modified piecewise polynomial collocation schemes on graded grids for finding approximations to the exact solution. In Section 5, we study the attainable order of proposed algorithms. In particular, we show that, by using specific collocation parameters and graded grids, it is possible for one method to attain a global super-convergence rate. Finally, in Section 6 we test the theoretical results of both methods by a numerical example. The main results of the paper are given by Theorem 1, Theorem 2 and Theorem 3. In Theorem 1, the existence and regularity results of the exact solution of problem (1.1)-(1.2) are presented. In Theorems 2 and 3 the convergence rates of the proposed algorithms are given.
In an analogous way we find that Using these expressions we can see that y(t) is determined by the formula We now show two alternative integral equation reformulations of problem (1.1)-(1.2).

Integral equation for y
Let y ∈ C[0, b] be a solution to problem (1.1)-(1.2) so that z = D α * y ∈ C[0, b]. From (1.1) we see that z(t) = f (t) − h(t)y(t) and by substituting this equation into (2.2), we obtain that y is also a solution of an integral equation of the form

Integral equation for z
Let y ∈ C[0, b] be a solution to problem (1.1)-(1.2) so that D α * y ∈ C[0, b]. By substituting (2.2) into (1.1) and using (1.5), we obtain that z = D α * y is a solution of an integral equation in the form For given q ∈ N and ν ∈ R, ν < 1, by C q,ν (0, b] we denote the set of continuous functions y : [0, b] → R which are q times continuously differentiable in (0, b] and such that for all t ∈ (0, b] and i = 1, . . . , q the following estimates hold (cf. e.g. [8]): Here c = c(y) is a positive constant. In other words, where, for t > 0, Equipped with the norm y C q,ν (0,b] := y ∞ +|y| q,ν , the set C q,ν (0, b] becomes a Banach space. Note that In particular, a function of the form The two following lemmas are based on the corresponding results of [8].
Lemma 1. If y 1 , y 2 ∈ C q,ν (0, b], q ∈ N, ν < 1, then y 1 y 2 ∈ C q,ν (0, b], and with a constant c which is independent of y 1 and y 2 . The existence and regularity of a solution to (1.1)-(1.2) can be characterized by the following theorem.
Proof. (i) We first note that due to f ∈ C[0, b] and (1.3) the forcing function g y of equation y = T y y + g y (see (2.3) and (2.5)) belongs to C[0, b]. Further, due to (2.4) operator T y can be rewritten in the form with H, T 1 , T 2 and T 3 defined by the following formulas: Clearly, H is bounded as an operator from and, due to (3.2), operator T y is compact as an operator from We first note that g y , the forcing function of equation y = T y y + g y , belongs to C q,ν (0, b]. This clearly follows from f ∈ C q,µ (0, b], (2.5), (3.1) and Lemma 2.
With the help of Lemma 1 we obtain that H is bounded as an operator from C q,ν (0, b] into C q,ν (0, b]. Further, due to (3.1) we have 1 − α ≤ ν. Therefore, it follows from Lemma 2 that J α is compact as an operator from C q,ν (0, b] into C q,ν (0, b]. Thus, J α H is linear and compact as an operator from C q, Linear operators (functionals) T 1 , T 2 , T 3 : C q,ν (0, b] → R are bounded and consequently compact in C q,ν (0, b]. Since H is bounded as an operator C q,ν (0, b] into C q,ν (0, b], we obtain that T 1 H, T 2 H and T 3 H are linear and compact as operators from C q,ν (0, b] into C q,ν (0, b]. Thus, T y defined by (3.2) is linear and compact as an operator from C q,ν (0, b] into C q,ν (0, b]. Since the homogeneous equation y = T y y has in C q,ν (0, b] ⊂ C[0, b] only the trivial solution y = 0, it follows from the Fredholm alternative theorem that equation y = T y y + g y has a unique solution y ∈ C q,ν (0, b]. This yields that prob- , then we may in Theorem 1 set ν = 1 − α.

Numerical method
Let N ∈ N and let Π N := {t 0 , . . . , t N } be 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 (4.1) are distributed uniformly; for r > 1 the points (4.1) are more densely clustered near the left endpoint of the interval [0, b]. For given integer k ≥ 0 by S is denoted the standard space of piecewise polynomial functions: Here and π k denotes the set of polynomials of degree not exceeding k. Note that the elements of S In every interval [t j−1 , t j ], j = 1, . . . , N , we define m ∈ N collocation points t j1 , . . . , t jm by formula where η 1 . . . , η m are some fixed (collocation) parameters which do not depend on j and N and satisfy In the following we describe two different numerical methods for solving the problem (1.1)-(1.2), based on the integral equations found in Section 2.

Method 1 (method based on the integral equation for y)
We look for an approximate solution y N ∈ S where y N is determined by the following collocation conditions: y N (t jk ) = (T y y N )(t jk ) + g y (t jk ), k = 1, . . . , m, j = 1, . . . , N. (4.4) Here T y , g y and t jk are defined by (2.4),(2.5) and (4.2), respectively. If η 1 = 0, then by y N (t j1 ) we denote the right limit lim t→tj−1,t>tj−1 y N (t). If η m = 1, then y N (t jm ) denotes the left limit lim t→tj ,t<tj y N (t). The collocation conditions (4.4) form a system of equations whose exact form is determined by the choice of a basis in S (−1) m−1 (Π N ). If η 1 > 0 or η m < 1 then we can use the Lagrange fundamental polynomial representation: where ϕ λµ (t) := 0 for t ∈ [t λ−1 , t λ ] and We look for an approximate solutionȳ N to (1.1)-(1.2) in the form (see (2.2)) with T z , g z and t jk defined by (2.7),(2.8) and (4.2), respectively. If η 1 = 0, then by z N (t j1 ) we denote the right limit lim t→tj−1,t>tj−1 z N (t). If η m = 1, then z N (t jm ) denotes the left limit lim t→tj ,t<tj z N (t).
Similarly to previous subsection, the collocation conditions (4.8) form a system of equations whose form we determine by using the Lagrange fundamental polynomial representation, this time for z N : Note that both methods can also be used in the case equation (4.3) has η 1 = 0 and η m = 1. We then have t jm = t j+1,1 = t j , c jm = c j+1,1 (j = 1, . . . , N − 1), and hence in the systems (4.6) and (4.10) there are (m − 1)N + 1 equations and unknowns.

Convergence estimates
In this section we study the convergence and convergence order of the proposed Method 1 and Method 2. For this we need Lemmas 3-6 presented below. The proofs of Lemmas 3-5 follow from the results of [8,28]. The proof of Lemma 6 can be found in [21].
Theorem 2. (i) Let the assumptions introduced in the part (i) of Theorem 1 be fulfilled. Moreover, let m ∈ N and assume that the collocation points (4.2) with arbitrary parameters η 1 , . . . , η m satisfying (4.3) and grid points (4.1) are used.
Then problem (1.1)-(1.2) has a unique solution y ∈ C[0, b] such that D α * y ∈ C[0, b]. Moreover, there exists an integer N 0 such that, for N ≥ N 0 , equations (4.4) and (4.7) possess unique solutions y N andȳ N , respectively. We have for y, the solution of (1.1)-(1.2), that (ii) Assume that (i) holds and let h, f ∈ C m,µ (0, b], µ ∈ R, µ < 1. Then for all N ≥ N 0 the following error estimates hold: Here ν is defined by (3.1), r ∈ [1, ∞) is the grid parameter in (4.1) and c is a positive constant which does not depend on N .
Proof. We present the proof only for Method 1. In a similar way we can construct the proof for Method 2.
(i) First we prove the convergence (5.3). Consider equation y = T y y+g y (see (2.3)), with T y and g y given by (2.4) and (2.5), respectively. We already know that equation y = T y y + g y is uniquely solvable in C[0, b]. We now show that it is also uniquely solvable in L ∞ (0, b). Proceeding as in the proof of Theorem 1, operator T y can be rewritten in the form (3.2). It follows from Lemma 2 that J α , T 1 , T 2 and T 3 are compact operators from L ∞ (0, b) into C[0, b]. Therefore, due to (3.2), operator T y is compact as an operator from L ∞ (0, b) into C[0, b], thus also from L ∞ (0, b) into L ∞ (0, b). Further, g y ∈ C[0, b] ⊂ L ∞ (0, b) and the homogeneous equation y = T y y has in C[0, b] only the trivial solution y = 0. This together with T y ∈ L(L ∞ (0, b), C[0, b]) yields that y = T y y possesses also in L ∞ (0, b) only the trivial solution y = 0. Consequently, by Fredholm alternative theorem, equation y = T y y + g y with g y ∈ L ∞ (0, b) possesses a unique solution y ∈ L ∞ (0, b). In other words, operator I − T y is invertible in L ∞ (0, b) and its inverse (I − T y ) −1 is bounded:

Conditions (4.4) have an operator equation representation
with the interpolation operator P N defined in (5.1). Since T y as operator from L ∞ (0, b) into L ∞ (0, b) is compact, it follows from Lemma 4 and from the boundedness of ( for all sufficiently large N , say N ≥ N 0 , and where c is a constant not depending on N . Thus, for N ≥ N 0 , equation (5.5) provides a unique solution y N ∈ S −1 m−1 (Π N ). We have for it and y, the solution of equation y = T y y + g y , that (I − P N T y )(y − y N ) = y − P N y, N ≥ N 0 . Therefore, by (5.6), where c is a positive constant not depending on N . This, together with y ∈ C[0, b] and Lemma 3, yields the convergence (5.3).
It follows from Theorem 2 that in the case of sufficiently smooth h and f , using sufficiently large values of the grid parameter r in (4.1), for both first and second method by every choice of collocation parameters 0 ≤ η 1 < · · · < η m ≤ 1 a convergence of order O(N −m ) can be expected.
The following result shows that by a careful choice of parameters η 1 , . . . , η m and with a slightly higher requirement for the smoothness of functions h and f it is possible to establish a faster convergence rate for Method 2.
Proof. We follow the approach given in [25]. From Theorem 2 it follows that problem (1.1)-(1.2) possesses a unique solution y ∈ C q,ν (0, b] and thus equation (2.6) has a unique solution z = D α * y ∈ C q,ν (0, b], where q = m + 1 and ν = max{1 − α, µ}. We note that conditions (4.8) have an operator equation representation z N = P N T z z N + P N g z (5.9) with the interpolation operator P N defined in (5.1). It follows from Theorem 2 that there exists an integer N 1 > 0 such that for N ≥ N 1 equation (5.9) has a unique solution z N ∈ S (−1) m−1 (Π N ). We denotê where T z and g z are defined by (2.7) and (2.8), respectively. By substituting z N = P NẑN into (5.10) we see thatẑ N is a solution of the equation It can be shown (see [25]) that there is an integer N 0 ≥ N 1 such that, for N ≥ N 0 , equation (5.11) is uniquely solvable and that Therefore, we have with some positive constants c i (i = 0, . . . , 4), which are independent of N . We observe that for some positive constants c 1 and c 2 , which are independent of N . Therefore, , we obtain with the help of Lemma 6 for N ≥ N 0 the following estimate: with some constants c 0 , c 1 and c 2 which are independent of N . Further, for N ≥ N 0 we have z − z N = (z − P N z) + P N (z −ẑ N ), and thus, by (2.2) and (4.7) where c 0 and c 1 are positive constants not depending on N . This together with (5.12) and Lemma 6 yields the estimate (5.8).

Numerical illustration
Let us consider the following boundary value problem: , 0 ≤ t ≤ 1, (6.1) y(0) + y(1) + The exact solution y(t) and its Caputo derivative z(t) := (D 1 3 * y)(t) are given by the following formulas: Clearly, h, f ∈ C q,µ (0, 1] with µ = 2 3 and arbitrary q ∈ N. Therefore, by (3.1), ν = max{1 − α, µ} = 2/3. For a numerical comparison between Method 1 and Method 2 we find the approximate solutions y N andȳ N using collocation points (4.2) with the collocation parameters and Gauss points We point out that collocation parameters (6.4) satisfy the assumptions set in Theorem 3, while collocation parameters (6.3) do not. In Tables 1 and 2 we display some results from numerical experiments for Method 1 with different values of N and r. The results in Table 1 and Table 2 correspond to collocation parameters (6.3) and (6.4), respectively. Similarly, in Tables 3 and 4 we display the numerical results for Method 2 with different values of N and r, using collocation parameters (6.3) and (6.4). The errors ε N in Tables 1 and 2 and errorsε N in Tables 3 and 4 are calculated as follows: where τ jk := t j−1 + k(t j − t j−1 )/10, k = 0, . . . , 10, j = 1, . . . , N , y N is the approximation to y given by (4.4),ȳ N is the approximation to y given by (4.8), and the gridpoints t j are defined by (4.1). The ratios Table 2. Numerical results for problem (6.1)-(6.2) using Method 1, with Gauss collocation points (6.4) and m = 2. characterizing the observed convergence rate, are also presented.
In the case of collocation points (6.3), it follows from (5.4) with α = 1 3 and ν = 2 3 that, for sufficiently large N , if r ≥ 6, (6.5) where c 0 is a positive constant not depending on N . Due to (6.5), the ratios N and¯ N for r = 1, r = 3 and r ≥ 6 ought to be approximatively 2 1 3 ≈ 1.26, 2 1 = 2 and 2 2 = 4, respectively. These values are given in the last row of Table 1, Table 2 and Table 3. As the results show, for both methods the collocation points (6.3) are in agreement with the theoretical estimates given by Theorem 2. Actually, we can see that for small values of r the actual convergence rate is faster.
In the case of collocation points (6.4), it follows from the estimates in (5.8) of Theorem 3 with α = 1 3 and ν = 2 3 that, for sufficiently large N , where c 0 is a positive constant not depending on N . Due to (6.6), the ratios N for r = 1, r = 3 and r ≥ 6 ought to be approximatively 2 2 3 ≈ 1.58, 2 2 = 4 and 2 7 3 ≈ 5.04, respectively. These values are given in the last row of Table 4. As we can see, the numerical results are in good agreement with theoretical estimates. We note that only for Method 2 a global super-convergence for special collocation parameters has been observed. Moreover, it follows from Table 4 that, in general, the attainable order of global convergence of Method 2 on the conditions of Theorem 3 cannot be improved. Finally, from Tables 1-3 we see that for problem (6.1)-(6.2) the actual convergence rate with small values of r for both methods is better than the theoretical convergence rate given by Theorem 2.

Conclusions
In the present work we have studied a class of non-local boundary value problems for linear fractional differential equations involving Caputo-type fractional derivatives. Using two different integral equation reformulations of the boundary value problem, the existence and regularity of the exact solution has been investigated. With the help of the obtained integral equations, we have constructed two numerical schemes which are based on collocation techniques and graded grids. We have studied the attainable order of prepared algorithms and shown the convergence of these schemes. In particular, we have shown that by using specific collocation parameters and graded grids, one method attains a global super-convergence rate.