On the Numerical Solution of Fractional Boundary Value Problems by a Spline Quasi-Interpolant Operator

: Boundary value problems having fractional derivative in space are used in several ﬁelds, like biology, mechanical engineering, control theory, just to cite a few. In this paper we present a new numerical method for the solution of boundary value problems having Caputo derivative in space. We approximate the solution by the Schoenberg-Bernstein operator, which is a spline positive operator having shape-preserving properties. The unknown coefﬁcients of the approximating operator are determined by a collocation method whose collocation matrices can be constructed efﬁciently by explicit formulas. The numerical experiments we conducted show that the proposed method is efﬁcient and accurate.


Introduction
Boundary value problems having fractional derivative in space are used to describe physical phenomena in which nonlocality effect are peculiar. For instance, they are used to model anomalous diffusion in biological tissues, viscoplastic materials in mechanical engineering or control of dynamical systems (see References [1][2][3][4][5] and references therein).
In particular, in this paper we are concerned with the Caputo fractional derivative [6]. The theoretical analysis of fractional boundary value problems (FBVPs) having Caputo derivative in space was addressed, for instance, in References [6][7][8][9][10][11]. We refer also to References [4,6,12,13] for the foundations of fractional calculus and details on fractional derivatives. We want to mention that in recent years other types of fractional derivatives were introduced, like the He's derivative [14] or the Fabrizio-Caputo derivative [15]. These derivatives are used to model physical phenomena characterized by the presence of structures with different scales.
In the literature several analytical and numerical methods were proposed for the solution of FBVPs. Analytical methods based on the homotopy perturbation technique [16] and the variational iteration method [17] were used, for instance, in References [18][19][20]. As for the numerical methods, several methods were proposed in the literature. In Reference [21] the authors proposed a Galerkin finite element approach to solve the one-dimensional steady state fractional advection dispersion equation. In Reference [22] the authors used finite difference methods to solve a nonlinear FBVP while in Reference [10] a collocation method based on spline functions was proposed to solve linear FBVPs. Spectral methods based on generalized Jacobi polynomials were used, for instance, in References [23,24]. In Reference [11] a Gegenbauer-based Nyström method was proposed to solve one-dimensional fractional-Laplacian boundary-value problems. In Reference [25] the authors solved linear and nonlinear FBVPs by a wavelet method. For an overview on numerical methods to solve fractional differential problems see, for instance, References [18,20,[26][27][28][29][30] and references therein.
In this paper we present a collocation method, based on a spline quasi-interpolant operator, for the solution of boundary value problems having Caputo derivative in space. Quasi-interpolant operators are approximating operators that reproduce polynomials up to a given degree. They have a greater flexibility with respect to interpolation operator. This freedom can be used to preserve special properties of the function to be approximated, like the sign, the shape or the area of its graph [31][32][33][34][35]. In particular, in this paper we propose a numerical method based on the Schoenberg-Bernstein operator [36], which is a positive operator having shape preserving properties. We determine the unknown coefficients of the approximating operator by a collocation method derived from References [37,38] and show through some numerical experiments that the method is efficient and accurate.
The paper is organized as follows. In Section 2 we describe the FBVPs we are interested in and the spline basis we use to construct approximating functions in this space. The main properties of the Schoenberg-Bernstein operator are also recalled. The details on the numerical method we propose are described in Section 3. The results of the numerical experiments are shown in Section 4 while some conclusions are given in Section 5.

Materials and Methods
In this section we describe the differential problem we are interested in (Section 2.1), the B-spline basis (Sections 2.2-2.4) used to construct the approximating function, and the main properties of the Schoenberg-Bernstein operator (Section 2.5).
Here, D γ x y denotes the Caputo fractional derivative in space defined as [6] D γ x y(x) := where Γ is the Euler's gamma function We assume y is a sufficient smooth function and the boundary conditions are linearly independent so that the differential problem has a unique solution [6,39,40].

The Cardinal B-Splines through the Truncated Power Function
The cardinal B-splines are compactly supported piecewise polynomials having breakpoints at the integers. They can be used to construct a function basis for the spline spaces. For details on spline functions see [41,42].
The system of the integer translates {B n (x − k), k ∈ Z} forms a basis for the n-degree spline space on the real line. Moreover, it reproduces polynomials up to degree n, has approximation order n + 1, and is a partition of unity, that is, for all x ∈ R .

B-Spline Bases on the Finite Interval
On the finite interval [0, L] a suitable basis for the spline space is the optimal basis, which is constructed using knots of multiplicity n + 1 at the endpoints of the interval [41,42].
For the sake of simplicity, we assume L is an integer greater than n. Thus, on the finite interval [0, L] we consider the sequence of integer knots I 0 = {x 0 , x 1 , . . . , x N }, where N = L + 2n + 1, and The optimal basis on the interval [0, L] is formed by L + n basis functions, that is, where the functions B k,n and B L+n−1−k,n , 0 ≤ k ≤ n − 1, are left and right boundary functions, respectively, while the functions B k,n , n ≤ k ≤ L − 1, are internal functions.
The analytical expression of the right boundary functions can be obtained using the central symmetry property, that is, We notice that by construction the following endpoint conditions hold where δ k0 denotes the Kronecker symbol. The optimal basis B n can be refined by using any sequence of equidistant knots on the interval [0, L] [41,42]. Denoting by h the refinement step, that is, the distance between the refined knot sequence, the refined basis is where N h = L/h + n. Once again, the basis N h,n has boundary and internal functions.
We notice that the refinement of the knots increases just the number of the internal functions while the number of the boundary functions remains the same. The optimal basis N h,3 for the cubic spline space on the interval [0, 1] with refinement step h = 1/8 is displayed in Figure 1.

Fractional Derivatives of Cardinal B-Splines
Since the internal functions are the basis functions having support all contained in the interval [0, L], their fractional derivatives can be easily evaluated by the differentiation rule where is the backward finite difference operator [37,44]. We notice that the fractional derivatives of the classical polynomial B-splines are fractional splines, that is, splines having noninteger degree (cf. Reference [44]). The fractional derivative D γ x B 3 of the cubic B-spline, evaluated for different values of γ, is shown in Figure 2. The plots show that the shape of the fractional derivative varies continuously with γ so that the order of the derivative acts as a tension parameter. The explicit expression of the fractional derivative of the left boundary functions can be obtained by applying the fractional differentiation operator to their analytical expression (5) (cf. Reference [43]), that is, where The fractional derivative of the right boundary functions can be obtained by the symmetry property, that is, Thus, the fractional derivative of the boundary functions is a linear combination of the fractional derivative of the translates of the truncated power function whose derivative has expression [43]

The Schoenberg-Bernstein Operator
A spline quasi-interpolant operator is a spline approximation of a given function that reproduces polynomials up to a given degree. There are several types of quasi-interpolant operators depending on which properties of the function to be approximated we require to preserve (see, for instance, References [31][32][33][34][35]). In this paper we consider the Schoenberg-Bernstein operator [36] S n y(x) = where θ k , 0 ≤ k ≤ L + n − 1, are the Greville nodes, that is, the coefficients that guarantee the reproduction of linear functions. They satisfy the property Even if the Schoenberg-Bernstein operator reproduces only polynomials of degree not greater than 1, it has many properties useful in applications. In particular, the operator is a positive operator and has shape preserving properties, that is, for any linear function Λ(x) it holds where S − (y) denotes the number of strict sign changes of its argument. Moreover, the operator satisfies the endpoint conditions, that is, S n y(0) = y(0) , S n y(L) = y(L) .
The Schoenberg-Bernstein operator is refinable meaning that we can construct a refined version of the operator using the refined basis N h,n , that is, where θ k,h , 0 ≤ k ≤ N h , are the refined Greville nodes satisfying Finally, the Schoenber-Bernstein operator is convergent with approximation order 1 [45], that is, We notice that usually the limit is taken either for h → 0 and n held fix or for n → ∞ and h held fix.

The Collocation Method
We solve the fractional differential problem (1) by a collocation method based on the refinable Schoenberg-Bernstein operator (13), that is, To determine the unknown values y(θ k,h ), 0 ≤ k ≤ N h , we choose a set of collocation points and collocate the differential problem in those points. For the sake of simplicity, here we assume the collocation points are a set of equidistant nodes on the interval [0, L] having space step δ, that is, Thus, collocating (1) on the the nodes X δ and using (13) we get This is a linear system that can be written in matrix form as where is the unknown vector, are the collocation matrices of the refinable basis N h,n and of its fractional derivative, respectively. We notice that the entries of the matrices A h,δ and D h,δ can be efficiently evaluated by the formulas given in Sections 2.2-2.4. The vectors are the know terms, the vectors contain the parameters, and the vectors contain the boundary values of the basis functions and of their first derivative, respectively. Here, the symbol F • A denotes the entrywise product between matrices. In the case when F is a vector, F has to be intended as a matrix having as many columns as A, each column being a replica of the vector F itself.
The linear system (17) has N δ − 1 + γ equations and N h + 1 unknowns. To guarantee the existence of a unique solution the refinement step h, the distance of the collocation points δ and the degree of the B-spline n have to be chosen such that N δ − 1 + γ ≥ N h + 1. When N δ − 1 + γ > N h + 1 we obtain an overdetermined linear system that can be solved by the least squares method [38].

Numerical Results
In this section we show the performance of the proposed method by solving some FBVPs. In the following tests we approximate the solution of the differential problem by the cubic Schoenber-Bernstein operator S h,3 y.

Example 1
In the first example we consider the fractional differential equation where f (x) = 1 and g(x) = The analytical solution is y(x) = x. We approximate the solution by the Schoenberg-Bernstein operator S h,3 y with h = 1/4. We choose δ = h/2 = 1/8 so that the final linear system has 9 equations and 7 unknowns. Since the operator S h,3 y reproduces any linear functions, in this case the approximation is exact. In Table 1 the infinity norm of the approximation error, evaluated as is listed for different values of γ. As expected, the error is in the order of the machine precision. We notice that when γ = 0 we are solving the approximation problem The unknown values y(θ k,h ) ≡ θ k,h , 0 ≤ k ≤ 6, are listed in Table 2. For all the values of γ they coincide at the machine precision with the Greville nodes of the interval [0, 1] [45].

Example 2
In the second example we solve the fractional differential equation where f (x) = 1 and g(x) = Γ(ν+1)x ν−γ Γ(ν−γ+1) + x ν . The analytical solution is y(x) = x ν . We approximate the solution by the Schoenberg-Bernstein operator S h,3 using different values of h. In all the tests we set δ = h/2. The infinity norm of the approximation error when ν = 4 for different values of γ is shown in Table 3. The analytical solution, the numerical solution and the approximation error e h,3 = y(x) − y h,3 (x) evaluated at the collocation nodes are shown in Figure 4 in the case when h = 1/128.

Conclusions
We have presented a collocation method based on the Schoenberg-Bernstein quasi-interpolant operator and used the method to efficiently solve boundary value problems having Caputo fractional derivative. The numerical results shown in Section 4 show that the method is accurate and exact on linear functions. We notice that we can increase the approximation order using different kind of quasi-interpolant operators, like the discrete operators introduced in References [32,35]. Some first results in this direction can be found in Reference [47]. Finally, we notice that even if the B-spline basis is centrally symmetric in the interval [0, L], its Caputo fractional derivative is not. The symmetry could be recovered replacing the Caputo derivative with the Riesz derivative defined as [4] d γ d|x| γ y(x) = − 1 2 cos( πγ 2 ) (D γ 0,x + D γ x,L )y(x) , which is centrally symmetric in the interval [0, L] . The solution of boundary value problems having Riesz fractional derivative in space will be the subject of a forthcoming paper.