Applications of Optimal Spline Approximations for the Solution of Nonlinear Time-Fractional Initial Value Problems

: Nonlinear fractional differential equations are widely used to model real-life phenomena. For this reason, there is a need for efﬁcient numerical methods to solve such problems. In this respect, collocation methods are particularly attractive for their ability to deal with the nonlocal behavior of the fractional derivative. Among the variety of collocation methods, methods based on spline approximations are preferable since the approximations can be represented by local bases, thereby reducing the computational load. In this paper, we use a collocation method based on spline quasi-interpolant operators to solve nonlinear time-fractional initial value problems. The numerical tests we performed show that the method has good approximation properties.


Introduction
Fractional differential equations are a well-established tool to model the nonlocal behavior of real systems. In the last few decades, fractional differential models have been widely used in many fields, such as continuum mechanics, biology, control theory [1][2][3][4]. At the same time, the quest for efficient numerical methods to solve fractional differential equations has become more and more intriguing (see [4,5] and references therein). Whereas in the beginning, the main effort was usually in adapting numerical methods used to solve ordinary differential equations to the fractional case, the challenge now is to construct numerical methods tailored for dealing with fractional derivatives. In this respect, collocation methods have received considerable attention for their ability to deal with the nonlocal behavior of the fractional derivative and are used to solve a variety of problems. For instance, a spectral collocation method based on a new family of fractional bases called Jacobi polyfractomials was used in [6] to solve steady-state and time-dependent fractional differential problems. The method was extended in [7] where initial value problems were considered. A collocation method based on Chebyshev polynomials was introduced in [8] and used to solve nonlinear multiterm initial value problems. These latter equations were also addressed in [9] in the context of collocation methods based on Laguerre polynomials. All these methods use polynomial bases, that is, bases having support in all the discretization interval. In order to maintain a low computational load, we look for an approximating solution that can be represented by a local basis. This goal can be achieved by collocation methods based on spline functions. These methods have a long history dating back to the 1970s (see [10][11][12]). In [13], spline collocation methods were used to solve fractional differential equations for the first time. Since then, they have been successfully used to solve a variety of problems (see, for instance, [14][15][16][17][18][19][20] and references therein). In these papers, the solution of the differential problem was approximated by an interpolating spline function, which requires us to impose restrictive conditions both on the spline space and on the interpolation points in order to guarantee convergence. These drawbacks can be overcome using quasi-interpolating splines since quasi-interpolation only requires us to interpolate polynomials up to a given degree, leaving more freedom in the choice of the spline space and of the collocation nodes [21,22]. In the literature, collocation methods based on spline quasi-interpolant operators were used to solve both integral and differential equations (see, for instance, [23][24][25]). More recently, a collocation method based on spline quasi-interpolant operators was proposed in [26] and then used to solve linear fractional differential equations in [27]. The method was proved to be accurate and efficient. Our aim is to generalize this method to the solution of nonlinear initial value problems having fractional derivative in time. In [26,27], we used truncated spline bases to represent the spline approximation. It is well known that truncated bases suffer from numerical instability at the initial point, which could produce instability in the approximate solution.
To overcome this problem, in this paper we use optimal spline bases, i.e., bases that are sufficiently smooth and satisfy the interpolation condition at the initial point. We show that the proposed method has good approximation properties when used to solve nonlinear differential problems.
The paper is organized as follows. In Section 2, we describe the differential problem we are interested in and recall some basic results on spline functions and spline quasiinterpolant operators. The numerical method is described in Section 3, while in Section 4 we show some numerical results. A discussion on the performance of the method is offered in Section 5.

Materials and Methods
We list here some symbols that will be used in the following sections. N denotes the set of natural numbers, Z denotes the set of integers and R denotes the set of real numbers while R + denotes the set of positive reals. The nonnegative real line will be denoted by I + = [0, ∞).

Initial Value Fractional Differential Problems
We are interested in solving the nonlinear initial value problem where D γ t y denotes the fractional derivative operator of order γ and y 0,ρ ∈ R, ρ ∈ N, are the initial conditions. Here, we consider the Caputo fractional derivative: where m = γ and Γ(t) is the Euler gamma function. The Caputo derivative is more suitable to model real-world phenomena because it retains many properties of the ordinary derivative. In particular, the Caputo derivative of a constant function is zero, which is not true for other types of fractional derivatives. Moreover, the initial conditions of differential problem (1) involve only ordinary derivatives.
Wellposedness of the solution to the fractional differential problem (1) was analyzed in [28] §6. In particular, if f (t, y) : I + → R is continuous and fulfills a Lipschitz condition with respect to y, then there exists a unique solution y ∈ C(I + ) that depends in a continuous way from the given data.

Optimal Spline Bases
We approximate the solution to Equation (1) by a spline linear operator, i.e., by a spline function constructed in order to satisfy special properties. We recall that a spline is a piecewise polynomial with prescribed smoothness. Spline functions can be represented as a linear combination of suitable spline bases. Here, we will use the optimal B-spline basis for its good approximation properties. Let S n be the space of splines of degree n. The starting point to construct an optimal basis for S n is the cardinal B-spline B n (t), which is a piecewise polynomial of degree n having knots on the integers. For our purposes, it is convenient to express B n by the truncated power function t n + := max(0, t) n : where [0, 1, . . . , n + 1] f denotes the divided difference operator on the integers. B n is compactly supported on I n := [0, n + 1], positive in (0, n + 1) and belongs to C n−1 (R). For further details on spline functions and cardinal B-splines see [29,30]. The integer translates {B n (t − k), k ∈ Z} form a basis for S n on the real line. In the semi-infinite interval I + the optimal spline basis is formed by the integer translates {B n (t − k), k ∈ N}, i.e., the cardinal B-splines having support all contained in I + , plus n left edge functions {B nk (t), 0 ≤ k ≤ n − 1}, which are cardinal B-splines with a knot of multiplicity n + 1 at the initial point: For details on the construction of optimal spline bases see [30,31].
The function system B n can be adapted to any sequence of equidistant knots by dilation. Let h denote the distance between the knots. The optimal basis on the knots {hi, i ∈ N} is: For easy reading, in the following we will use the notation We recall that optimality refers to the shape preserving properties of the basis B nh [32]. In particular, the basis B nh enjoys the variation diminishing property [33], i.e., for any where S − (·) denotes the number of sign changes of its argument. Moreover, the functions B nhk , k ∈ N, satisfy the initial conditions B nhk (0 + ) = δ k0 , k ∈ N. As a consequence, spline approximations preserve the shape of the data, as sign, monotonicity and convexity, and satisfy the left-edge interpolation condition, i.e., ∑ k∈N c k B nhk (0 + ) = c 0 . Moreover, we recall that optimal spline bases are numerically stable, i.e., where is the condition number of the basis B nh . This means that numerical errors are not amplified when evaluating spline approximations.
Finally, the optimal basis B nh reproduces polynomials up to degree n: where ξ (r) hk are real numbers that can be evaluated explicitly (see [30,34]). In particular, for r = 1, ξ (1) hk = 1 for all k ∈ N, so that the function system B nh forms a partition of unity: For r = 2, we obtain hk , k ∈ N} are the Greville points defined as

Discrete Spline Quasi-Interpolant Operators
To approximate the solution to the differential Problem (1) we need a suitable approximating operator. To this end, we choose quasi-interpolant operators that are more flexible than interpolating operators and can be designed to satisfy special properties of the function to be approximated. In particular, in this paper we consider the discrete spline quasi-interpolant operator introduced in [34] (see also [30,35]). Using the optimal basis B nh we can construct a sequence of refinable quasi-interpolant operators with good approximation properties.

The Numerical Method
We approximate the solution to the fractional differential Problem (1) by the discrete quasi-interpolant operator (9): where are the unknown coefficients whose expression depends on r (see (11) and (12)). We notice that, since B nhk (t) has compact support, for any value t ∈ I + the sum on k has a finite number of terms. Thus, considering a finite discretization interval I = [0, T], there is a finite number of unknown coefficients that can be evaluated by collocation as described below. Let {t p , 0 ≤ p ≤ M} be a set of collocation points belonging to I. We require the approximating function y nh , n ≥ γ , to solve the differential problem on the collocation points: Substituting (14) in the previous equations, we obtain the nonlinear system: where N h + 1 is the number of functions belonging to B nh such that supp(B nhk ) ∩ I = {0}.
The nonlinear system (16) has M + γ + 1 equations and (N h + 1)r unknowns with M + γ + 1 ≥ (N h + 1)r. The existence of a unique solution to (16) depends on the existence of a unique solution to (1). In particular, if f is sufficiently smooth in some neighborhood I of the exact solution to (1), the approximate solution y nh exists and is unique in I for sufficiently small h (see [12,15,36]). When M + γ + 1 > (N h + 1)r, Equations (16) form an overdetermined nonlinear system that can be solved by the nonlinear least squares method.
We recall that from (7), it follows that the approximation y nh is exact on polynomials up to degree n.

Numerical Results
In this section, we show the results we obtained when solving some fractional nonlinear differential equations by the collocation method (16). In the numerical tests we used cubic and quartic optimal bases. The explicit expression of the basis functions and of their fractional derivative can be found in [31]. In Figure 1, we show the cubic optimal basis B 3h and the quartic optimal basis B 4h in the interval [0, 1] when h = 1 8 . Their fractional derivative for different values of γ is shown in Figures 2 and 3.   For the sake of simplicity, in the tests, we set h = 1 2 j and choose equidistant nodes as collocation nodes, {t p = 2hp, 0 ≤ p ≤ M} with M = 2 j+1 T − 1. We obtain an overdetermined nonlinear system that we solved with the nonlinear least squares method. The accuracy of the numerical solution is evaluated by computing the infinity norm of the error evaluated on a refined grid:

Example 1
The aim of this example is to check the accuracy of the proposed method. To this end, we solve the nonlinear differential problem whose exact solution is y(t) = t µ − 1.
In the first test, we set µ = 2 so that the exact solution is a polynomial of degree 2. Due to the polynomial reproduction property (7), the approximation error y − y nh vanishes when n ≥ 2. In this test, we choose h = 1 8 , 1 16 , 1 32 so that the final overdetermined nonlinear system has N h = 1 h + n unknowns and M = 2 h + 1 equations. In Tables 1 and 2, we list the infinity norm of the error for decreasing values of h when using the quasi-interpolant operator (13) with n = 3 and n = 4. Since both approximations y 3h e y 4h reproduce quadratic polynomials, the error is very small; in particular, when h = 1 8 , the error is in the order of the machine precision. We notice that the error increases slightly when h increases. This is due to the numerical instabilities arising when the dimensions of the nonlinear system increase. To give an idea of the stability of the nonlinear system, in Tables 1 and 2, we also list the condition number of the Jacobian matrix. In the second test, we check the theoretical convergence order. To this end, we solve the differential problem (17) when µ = 1.9.
Since y nh is an approximating operator in the spline space, the convergence of the collocation method applied to the differential problem (17) is guarantee with convergence order at least µ if y is sufficiently smooth [29]. As a consequence, lim h→0 y(·) − y nh (·) ∞ = O(h −µ ).
In Figure 4, we show the semi-log plot of the infinity norm of the error for decreasing values of h = 2 −j when using as approximating function the quasi-interpolant operator (13) with n = 3 (red line) and n = 4 (blue line). The plots show that, as expected, the error decreases when h decreases. Moreover, the slope of the red and blue lines in the semilog plots are the same as the theoretical slope −µ. We observe that the accuracy of the approximation obtained with the optimal basis of degree n = 4 is higher than the accuracy obtained for n = 3. This is due to the higher smoothness of the approximating function.

Example 2
In this section, we use the proposed method to solve two nonlinear differential problems arising in real-world phenomena. First of all, we consider the fractional order logistic equation which is a nonlinear equation used to model population growth in biology and social studies. The existence and uniqueness of the solution to the logistic Equation (18) were proved in [37]. Its analytical solution is: where E γ (z) is the Mittag-Leffler function: (cf. [38]). In Figure 5, we show the numerical solutions y 3h (t) and y 4h (t) and the errors y(t) − y 3h (t) and y(t) − y 4h (t) for different values of γ when µ = 0.7, y 0 = 0.8 and h = 1 32 .  In the second example, we solve a nonlinear multiterm fractional differential equation, i.e., a differential equation characterized by the presence of both ordinary and fractional derivatives. This kind of equations is used to model viscoelastic materials. In particular, we consider the differential equation whose exact solution is y(t) = t 1.9 − 1. In Figure 6, we show the numerical solutions y 3h (t) and y 4h (t) and the errors y(t) − y 3h (t) and y(t) − y 4h (t) for different values of γ when h = 1 32 . The plots in Figures 5 and 6 show that the approximate solution has a good accuracy even in case of low values of γ. The accuracy increases as the order of the fractional derivative increases because the smoothness of the known term of the differential equation increases. Finally, we notice that, due to the interpolation property of the optimal basis, the error in the initial point is very small.

Conclusions
We presented a collocation method based on spline quasi-interpolant operators suitable to solve nonlinear differential problems having fractional derivative in time. The numerical results in Section 4 show that the method is convergent and accurate. We notice that the polynomial reproduction property (7) is a key ingredient to obtain high accuracy. Moreover, the approximations obtained by the quasi-interpolant operator (13) have shapepreserving properties. This means that the approximating function y nh preserves the shape of the function to be approximated. This can be useful in applications where we need to preserve the sign or the monotonicity of the approximation.
Some issues are worth analyzing in detail. First, in this paper, we used equidistant collocation nodes. It is well known that different distributions of nodes could improve the accuracy of spline collocation methods [11,14]. The behavior of the present method when using other kinds of nodes, such as Gaussian points or graded meshes, should be analyzed. Secondly, the stability and convergence of the method and the behavior of the conditioning of the nonlinear system should be studied in detail. Some preliminary results in the case when 0 < γ < 1 can be found in [36]. A thorough analysis of the general case is under study. Finally, we observe that the method can be easily extended to other kinds of fractional differential equations, such as boundary value problems (see [39]). This will be the subject of a forthcoming paper.
Author Contributions: E.P. and F.P. contributed equally to the manuscript. All authors have read and agreed to the published version of the manuscript.