Optimal B-Spline Bases for the Numerical Solution of Fractional Differential Problems

Efficient numerical methods to solve fractional differential problems are particularly required in order to approximate accurately the nonlocal behavior of the fractional derivative. The aim of the paper is to show how optimal B-spline bases allow us to construct accurate numerical methods that have a low computational cost. First of all, we describe in detail how to construct optimal B-spline bases on bounded intervals and recall their main properties. Then, we give the analytical expression of their derivatives of fractional order and use these bases in the numerical solution of fractional differential problems. Some numerical tests showing the good performances of the bases in solving a time-fractional diffusion problem by a collocation–Galerkin method are also displayed.


Introduction
In the last few decades, fractional calculus has proved to be a powerful tool for describing real-world phenomena.Differential problems of fractional order, initially introduced to model anomalous diffusion in viscoelastic materials, are now used in several fields, from physics to population dynamics, from signal processing to control theory [1][2][3][4][5].For an introduction of fractional calculus, we refer to [6][7][8][9].
As fractional differential models have become widespread, the development of efficient numerical methods to approximate their solution has become of primary interest.In fact, the nonlocality of the fractional derivative poses a severe challenge in its approximation, and, to face this problem, several numerical methods were proposed in the literature.For a review on numerical methods, see, for instance, [10][11][12][13][14][15].
In particular, we are interested in the solution of differential problems on bounded intervals.In this case, it is important to have a function basis available that naturally fulfills boundary and/or initial conditions and does not show numerical instabilities at the boundaries.From this point of view, B-spline bases are especially suitable and their use in the numerical solution of classical differential problems is widely diffused (see, for instance, [16][17][18]).Despite the indubitable success of these methods, the use of B-spline bases for the solution of differential problems of fractional order is not yet very common and limited to few examples.For instance, in [19], the authors used a linear B-spline basis to solve a multi-order fractional differential problem by the operational matrix method, while, in [20], a collocation method based on cubic B-spline wavelets was constructed to solve ordinary differential equations of fractional order.
The aim of this paper is to show how optimal B-spline bases can be profitably used for the solution of fractional differential problems.To this end, we not only describe in detail the construction of optimal B-spline bases of any degree but also give the analytical expression of the basis functions and of their fractional derivatives.These expressions are an efficient tool to construct numerical methods based on spline functions.
The paper is organized as follows.In Section 2, we show how to construct polynomial B-spline bases and recall their main properties.In particular, the construction of cardinal B-splines on the real line through the divided difference operator is described in Sections 2.1-2.2,while, in Section 2.3, we show how to construct optimal B-spline bases on bounded intervals.Their analytical expression in case of equidistant nodes is given in Section 2.4.Section 2.5 is devoted to the evaluation of the fractional derivatives of the B-spline bases.Finally, in Section 2.6, we show how to use optimal bases for the numerical solution of a time-fractional diffusion problem by a collocation-Galerkin method.In Section 3, we use the optimal B-spline basis of degree 3 to solve some test problems; some numerical results are also displayed.Section 4 contains a discussion on the results and highlights possible extensions.

Materials and Methods
In this section, we define the B-spline bases through the divided difference operator and recall their main properties.Then, we give the analytical expression of the basis functions and of their fractional derivatives.Finally, B-spline bases are used to solve a fractional differential problem by a collocation-Galerkin method.

The Divided Difference Operator
be a sequence of simple, i.e., non coincident, knots and let be the determinant of the collocation matrix associated with the function system U = {u 1 (x), u 2 (x), . . ., u n (x)} and the knot sequence (1).
The divided difference operator [x 1 , x 2 , . . ., x n ] f of a function f over the knots X is defined as The formula above can be generalized to the case of coincident knots as follows.Let the knots occur more than once, i.e., and let be the multiplicities of the coincident knots, i.e., where η 1 < η 2 < . . .< η r are the non coincident knots of the sequence X .Assuming the functions in U are sufficiently smooth, the collocation matrix associated with the function system U over the knot sequence (7) is given by Thus, if the function f has sufficient derivatives, the divided difference operator of f over the sequence of multiple knots X can be defined by using the collocation matrix (8) in Definition (4), i.e., where D denotes the determinant of the collocation matrix as defined in Equation (8).

The Cardinal B-Splines
The polynomial splines are piecewise polynomials having a certain degree of smoothness which is related to the degree of the polynomial pieces.Spline functions can be represented as a linear combination of B-splines that form a local basis for the spline spaces [21].
The cardinal B-splines, i.e., the polynomial B-splines having break points on integer knots, can be defined by applying the divided difference operator (4) to the truncated power function defined as Then, the cardinal B-spline B n of degree n on the integer knots where ∆ n is the finite difference operator From the above definition, it follows that the cardinal B-spline B n is compactly supported on [0, n + 1], positive in (0, n + 1) and belongs to C n−1 (R).Moreover, the system of the integer translates forms a function basis that is a partition of unity, i.e., and reproduces polynomials up to degree n [21].
The basis B n can be generalized to any set of equidistant knots X h = {x k = k h, k ∈ Z}, where h > 0 is the space step, by scaling, i.e., For details and further properties on polynomial B-splines, see, for instance, [21].
It is always possible to construct bases on bounded intervals by restricting the B-spline basis B nh on an interval [a, b].However, in this way, we obtain bases that could be numerically unstable since the functions at the boundaries are obtained by truncation.Moreover, boundary or initial conditions are not easy to satisfy.In the following section, we will show how to construct stable bases on the interval that easily satisfy boundary or initial conditions.

Optimal B-Spline Bases
Optimal B-spline bases on bounded intervals can be defined by introducing multiple knots at the boundaries of the given interval [21].
Let L ≥ n + 1 be an integer and let I = {0, 1, . . ., L − 1, L} be the sequence of the integer knots of the interval I = [0, L].We extend I to a sequence I of L + 1 + 2n knots by introducing knots of multiplicity n + 1 at the boundaries of the interval, i.e., The functions of the optimal B-spline basis of degree n with knots I are given by the divided difference so that each function N kn has compact support with From (19), it follows that the basis N n has L − n interior functions, i.e., the functions N kn , n ≤ k ≤ L − 1 having all simple knots, and 2n edge functions, i.e., the functions N kn for 0 ≤ k ≤ n − 1 and L ≤ k ≤ L − 1 + n that have a multiple knot at the left or right boundary, respectively.
The basis N n is centrally symmetric, i.e., Moreover, the functions N kn satisfy the boundary conditions where D r x denotes the usual derivative of integer order r.Finally, the basis N n forms a partition of unity, i.e., and is stable, i.e., for any sequence c = {c 0 , c 1 , . . ., where is the condition number of the basis N n .From (25), it follows that where ∆c denotes a perturbation of the sequence c.As a consequence, optimal B-spline bases are numerically stable so that numerical errors are not amplified when evaluating spline approximations.The basis N n can be generalized to any sequence of equidistant knots on a bounded interval [a, b] by mapping x → (x − a)/h, i.e., where L = (b − a)/h.The interior functions are the L − n functions while the 2n edge functions are All the properties above make the basis N n optimal when used in approximation problems [22,23].

Analytical Expression of the Optimal B-Spline Bases
Let N n be the optimal basis (18) associated with the knot sequence (17).The interior functions The left edge functions , can be evaluated by using the finite difference operator (9) for coincident knots.For 0 ≤ k ≤ n − 1, we get By Definition (8), we get where D kn is the (n − k + 1) order diagonal matrix and T kn is the (k + 1) order collocation matrix Thus, for the numerator in the last equality of (31), we have For the denominator, a similar calculation gives where P kn is the k + 1 order collocation matrix Theorem 1.Let I be the sequence of multiple knots (17) on the interval I = [0, L].The analytical expression of the left edge functions is given by where |T nk | and |P kn | are the determinants of the collocation matrices ( 35) and (38), respectively.
The right edge functions , can be easily obtained recalling the symmetry property (21).
In the following corollary, we give the analytical expression of the first three left boundary functions.

Corollary 1.
Let I be the sequence of multiple knots (17) on the interval I = [0, L].For 0 ≤ k ≤ 2, the left edge functions N kn are given by It is interesting to observe that In particular, when n = 1 (linear case), there is just one left edge function N 01 (x) = B 1 (x + 1) and the basis N 1 coincides with the basis B 1 restricted to the interval I.We recall that the basis N 1 is the unique interpolatory B-spline basis.

Fractional Derivatives of the Optimal B-Spline Bases
In this section, we give the explicit expression of the derivatives of fractional order of the functions in the basis N n .
Let Γ be the Euler's gamma function Definition 1. Assuming f is a sufficiently smooth function, the Caputo fractional derivative of f is defined as Definition 2. The Riemann-Liouville fractional derivative of a function f is defined as We recall that the Caputo derivative and the Riemann-Liouville derivative can be obtained from each other by If f satisfies homogeneous initial conditions, i.e., f (k) (0 + ) = 0, 0 ≤ k ≤ m − 1, the Riemann-Liouville derivative coincides with the Caputo derivative.
To use the basis N n for the solution of differential problems of fractional order, we need the expression of the fractional derivatives of the functions N kn .
The derivatives of the interior functions can be easily evaluated by the differentiation rule where ∆ n is the finite difference operator (13) [15,24].Since B n satisfies homogeneous boundary conditions, the Riemann-Liouville derivative is equal to the Caputo derivative.
Theorem 2. The Caputo derivative of the left edge functions are given by where Proof.First of all, we write the explicit expression of the matrix T kn , i.e., The claim follows by expanding the determinant of T kn along the last column.
Thus, to evaluate the fractional derivatives of the edge functions, we need the fractional derivatives of the translates of the truncated power function.By Definition (44), we get Integration rules for rational functions give [25] x where Θ(z) is the Heaviside function and 2 F 1 (a, b, c, z) is the hypergeometric function with (q) k denoting the rising Pochhammer symbol Thus, we get the analytical expression of the Caputo derivative (cf.also [20]).
Theorem 3.For 0 < β < n with m − 1 < β < m and m a positive integer, the Caputo derivative of the translates of the truncated power function is given by The Riemann-Liouville derivatives of the functions N kn can be obtained by the Caputo derivatives using the initial conditions (22) in Equation ( 46).

The Collocation-Galerkin Method
In this section, we show how to use the optimal basis introduced in the previous sections to solve a fractional differential problem, i.e., the time-fractional diffusion problem (cf.[26][27][28]) where D β t u, 0 < β < 1, denotes the partial Caputo derivative with respect to the time t.
where L = η/h is the approximating space generated by the optimal B-spline basis of the interval [0, η] on equidistant knots with space step h > 0 (cf.( 27)).We observe that the basis N 0 nh naturally satisfies the homogeneous boundary conditions For any h held fix, the spline Galerkin method looks for an approximating function that solves the variational problem where ( f , g) = η 0 f g.Writing (62) in a weak form and using (61), we get the system of fractional ordinary differential equations with initial condition Here, C(t) = (c k (t)) 1≤k≤L−2+n is the unknown function vector, Q = (q kj ) 1≤k,j≤L−2+n is the mass matrix, L = ( kj ) 1≤k,j≤L−2+n is the stiffness matrix and F = ( f k (t)) 1≤k≤L−2+n is the load vector whose entries are The entries of Q and L can be evaluated explicitly by using the integration and differentiation rules for B-splines [21].The entries of F can be evaluated by quadrature formulas especially designed for Galerkin methods [29].Now, we introduce the sequence of temporal knots.Let where T = τ/s with s > 0 the time step, be a set of equidistant nodes in the interval [0, τ] having a knots of multiplicity m + 1 at the left boundary.We assume the unknown functions c k (t) belong to the spline space We observe that just the functions N rms , 0 ≤ r ≤ n, are boundary functions while all the other functions in the basis are cardinal B-splines.
To solve the fractional differential system (63), we use the collocation method on equidistant nodes introduced in [28,30].Let P = {t p , 1 ≤ p ≤ P} with P ≥ T − 1 + m be a set of non coincident collocation points in the interval (0, τ].Then, collocating (63) on the nodes t p , we get the linear system where Λ = (λ rk ) 1≤r≤T−1+m,1≤k≤L−2+n is the unknown vector, are collocation matrices, and F = ( f k (t p )) 1≤k≤L−2+n,1≤p≤P .When P = T − 1 + m, (68) is a square linear system and the collocation matrix A is non singular if and only if [21], while B is always non singular [31].Otherwise, ( 68) is an overdetermined linear system that can be solved in the least squares sense.
We notice that we must pay a special attention to the evaluation of the entries of A since they involve the values of the fractional derivative D β t N rms in the collocation points.This can be done by the differentiation rules given in Section 2.5.

Results
To test the performance of the optimal B-spline basis when used to solve fractional differential problems, we solved the time-fractional diffusion problem (57) by the collocation-Galerkin method described in Section 2.6.In the tests, we used the optimal basis N 3 having degree 3 since the cubic B-spline has a small support but is sufficiently smooth.In fact, the cubic B-spline belongs to C 2 (I) and its support has length 4. In particular, the basis N 3 has 3 right (left) edge functions with support In the following section, we give the analytical expression of the functions N k3 and of their fractional derivatives; then, we show some numerical results.

The Optimal B-Spline Basis of Degree n = 3
When n = 3, we get the cubic B-spline basis.The optimal B-spline basis N 3 of the interval [0, L], with L an integer greater than 3, is associated with the sequence of integer knots L − 3 interior functions and six edge functions, three for each boundary.
The interior functions N 3k , 3 ≤ k ≤ L − 1, are the translates of the cubic cardinal B-spline, i.e., where From Corollary 1, it follows that the three left edge functions N kn , 0 ≤ k ≤ 2, are given by Now, using Theorems 3, we obtain the analytical expression of the Caputo derivative of fractional order 0 < β < 1 of the translates of the truncated power function Thus, substituting (78) in (48), we get The optimal basis and its fractional derivatives in case of equidistant knots on an interval [a, b] can be evaluated by scaling as shown in ( 27)- (29).The optimal B-spline basis in the interval [0, 2] for h = 0.25 is shown in Figure 1a.The Caputo derivatives of fractional order β = 0.25, 0.5, 0.75 are shown in Figure 1b-d.To give an idea of the condition number of the discretization matrix (Q ⊗ A + L ⊗ B) (cf.(68)), in Tables 1-3, we list the condition number of the matrix when using the cubic B-spline basis for different values of the order of the fractional derivative and different values of the space and time steps.

Numerical Tests
In this section, we show the numerical results we obtained when solving the differential problem (57) for three different expressions of the known term f (x, t) and different values of the fractional derivative.In the tests, we set τ = 1 and η = 2 and chose as collocation nodes the equidistant points of the interval [0, 1] with time step δ.All the tests were performed in the Mathematica environment [32].

Test 1
To test the accuracy of the method, we first solved the time-fractional diffusion Equation (57) when In this case, the exact solution is the bivariate polynomial Since the cubic B-spline basis reproduces polynomials up to degree 3, the approximation error of the collocation-Galerkin method is zero so that we expect the numerical solution coincides with the exact solution.
We evaluated the numerical solution for δ = 0.25, s = 0.5 and h = 0.25 and different values of the order of the time fractional derivative β.
Let us define the error by e h,δ (x, t) = u(x, t) − u h,δ (x, t) , where u h,δ (x, t) denotes the numerical solution of the fractional differential system (63) obtained by the collocation method.
The tests show that, as expected, the error is on the order of the machine precision even using a large step size.The numerical solution and the error for β = 0.25 are shown in Figure 2.

Test 2
To test the accuracy of the method when approximating time fractional derivatives, we solved the time-fractional diffusion Equation (57) for where 1 F 1 (a, b, z) is the Kummer's confluent hypergeometric function where N 0 = N\{0} (cf.[9]).In this case, the exact solution is so that the error for the spatial approximation is negligible.The L 2 -norm of the error and the numerical convergence order are listed in Table 4.The numerical solution in the case when β = 0.25 obtained for δ = 0.03125, s = 2δ and h = 0.25 is shown in Figure 3.The numerical results show that the method converges with convergence order approximatively equal to 4 and gives a good approximation even using few collocation points.Moreover, the accuracy of the numerical solution does not depend on the order of the fractional derivative.

Test 3
Finally, we solved the time-fractional diffusion equation (57) for In this case, the exact solution is The L 2 -norm of the error and the numerical convergence order are listed in Table 5.The numerical solution in the case when β = 0.25 obtained for δ = h = 0.03125 and s = 2δ is shown in Figure 4. Also in this case, the numerical solution has a good accuracy that is not affected by the order of the fractional derivative.The numerical convergence order is approximatively 4. Table 5. Test 3: The L 2 -norm of the error and the numerical convergence order ρ h,δ as a function of the time step δ and the space step h for different values of the order β of the fractional derivative.

Discussion
In this paper, we showed how optimal B-spline bases have good approximation properties that make them particularly suitable to be used in the solution of fractional differential problems.First of all, boundary and initial conditions can be satisfied easily in view of properties ( 22)- (23) so that the numerical solution we obtained does not suffer from the numerical instabilities near the boundaries that appear when truncated functions at the boundaries are used (cf.[28]).Moreover, optimal B-spline bases are stable (cf.(25)) meaning that numerical errors are not amplified.Nevertheless, Tables 1-3 show that the condition number of the discretization matrix increases as the dimension of the matrix increases.The conditioning gets worse when the order of the fractional derivative becomes smaller.In this case, the linear system (68) can be accurately solved by Krylov methods [33].
As for the convergence order, we observe that several methods proposed in the literature to solve fractional differential problems have convergence order that depends on the order of the fractional derivative (see, for instance, [26,27]).In the case of the cubic B-spline basis, the numerical tests show that the numerical convergence order does not depend on the order of the fractional derivative being close to 4 for any value of β.This result is in line with the Strang-Fix theory for classical differential problems [34].We expect that the convergence order for the optimal B-spline basis N n is n + 1 so that it can be increased easily by increasing the degree of the basis.In fact, optimal bases of high degree and their fractional derivatives can be evaluated by using the analytical expression given in Sections 2.4-2.5.As far as we know, the analytical expressions (39) and ( 48) are new and, together with ( 27)- (29), are an easy tool to evaluate optimal B-spline bases on any set of equidistant knots on bounded intervals.Finally, we observe that the accuracy of the numerical solution we obtained is higher than the accuracy obtained in [26] where a quadrature formula was used to approximate the time fractional derivative and a finite element method was used for the spatial approximation.
The optimal bases we described can be used to solve other fractional differential problems, for instance problems involving the Riesz derivative [2], or can be used in other numerical methods, for instance in the operational matrix method or in wavelet methods.We notice that adaptive wavelet methods can be used to solve differential problems whose solutions are non smooth functions.Optimal bases generating multiresolution analyses and wavelet spaces can be obtained starting from a special family of refinable functions [35,36] so that the procedure described above can be generalized to other refinable bases.

Table 1 .
The condition number K (rounded to the nearest integer) and the dimensions n rows × n cols of the discretization matrix (Q ⊗ A + L ⊗ B) for α = 0.25 as a function of the space step h and the time step δ.

Table 2 .
The condition number K (rounded to the nearest integer) and the dimensions n rows × n cols of the discretization matrix (Q ⊗ A + L ⊗ B) for α = 0.5 as a function of the space step h and the time step δ.

Table 3 .
The condition number K (rounded to the nearest integer) and the dimensions n rows × n cols of the discretization matrix (Q ⊗ A + L ⊗ B) for α = 0.75 as a function of the space step h and the time step δ.