A Fractional B-spline Collocation Method for the Numerical Solution of Fractional Predator-Prey Models

: We present a collocation method based on fractional B-splines for the solution of fractional differential problems. The key-idea is to use the space generated by the fractional B-splines, i.e., piecewise polynomials of noninteger degree, as approximating space. Then, in the collocation step the fractional derivative of the approximating function is approximated accurately and efﬁciently by an exact differentiation rule that involves the generalized ﬁnite difference operator. To show the effectiveness of the method for the solution of nonlinear dynamical systems of fractional order, we solved the fractional Lotka-Volterra model and a fractional predator-pray model with variable coefﬁcients. The numerical tests show that the method we proposed is accurate while keeping a low computational cost.


Introduction
Fractional Calculus generalizes to real order the well known concept of derivative and integral of integer order. Even if the birth of fractional calculus can be set at the beginning of the 18th century, its use for the modeling of real-life problems is very recent. Nevertheless, in the last decades the interest in fractional differential and integral equations is rapidly growing and models involving fractional derivatives and/or integrals are widely used in several fields, from physics to continuum mechanics, from biophysics to electro-chemistry, from signal processing to control theory. For an introduction to fractional calculus see, for instance, [1][2][3][4] while a survey on its applications can be found in [5][6][7][8][9].
The reason why the fractional derivative is more suitable to model real-world problems is related to its non-local behavior that is able to introduce in a elegant way into the model, either memory effects in time or non-locality in space. In this paper we refer to the traditional Caputo fractional derivative introduced in [10] and extensively used to model a great variety of mechanical and biological phenomena. More recently, new fractional derivatives having non-singular kernel were introduced to better describe heterogeneous materials or structures with different scales. We refer the interested reader to the literature (see, for instance, [11][12][13][14][15] and references therein).
It is exclusively the property of being non-local that makes the solution of integro-differential equations of fractional order an exciting challenge. Analytical methods to solve fractional problems are mainly based on Laplace or Fourier transforms, series solutions, Laguerre's integral formula, direct solution using the Grunwald-Letnikov definition (see, for instance, [3,[16][17][18]). These methods could be difficult to use when solving complicated problems; moreover, the solution is given in a closed form that requires the evaluation of special functions with an involved expression, as the Mittag-Leffler function [19]. More recently, the Adomian decomposition method [20] and the homotopy perturbation method [21] have been used to approximate the solution of fractional differential equations of different types (see [18,[22][23][24][25][26][27][28] and references therein). In these methods the solution is expressed as a convergent series that is approximated by truncation. Since truncation affects the accuracy of the solution in large intervals, a large number of terms in the series is required to reach a good accuracy far away from the initial point. Moreover, the series does not have a general form so the expression of its terms has to be evaluated for each specific differential problem making the evaluation of accurate solution rather cumbersome.
For these reasons, the development of numerical methods to solve fractional differential problems is of paramount importance and the literature in this field is growing rapidly. When constructing numerical methods to approximate a fractional derivative, particular attention should be given to keep good accuracy and high efficiency while preserving the non-local behavior in the approximate derivative.
A first approach is to extend to the fractional case finite difference methods commonly used to numerically solve classical differential problems of integer order. Starting from the seminal works by Lubich [29,30] and Podlubny [3], several finite difference methods have been constructed and used to solve many fractional differential problems. Among the finite difference methods, here we want to cite the method introduced in [31], based on the Grunwald-Letnikov definition of the fractional derivative, and the Adams methods in [32][33][34]. For a review on finite difference methods and further references see [35]. The advantage of finite difference methods is in their easy implementation. On the other hand, the finite difference is a local operator that requires a high order to reproduce the non-locality of the fractional derivative so that these methods are computationally demanding when high accuracy is required [36].
In this regard, global methods that approximate the solution in all the discretization intervals by expansion in global bases are becoming increasingly popular. A first attempt in this direction was the Fourier spectral method introduced in [37]. Subsequently, spectral collocation and/or tau methods were addressed, for instance, in [38][39][40][41] while spectral Galerkin and Petrov-Galerkin methods were addressed in [42][43][44][45][46]. The drawback of spectral methods lies in the fact that the coefficient matrix of the resulting linear system is dense so that its numerical solution could be computationally demanding.
To overcome this problem the solution can instead be approximated by expansion in a local basis, usually a piecewise polynomial basis. Collocation methods based on polynomial splines were introduced for the first time in [47] and then considered in [48,49], while the Galerkin and Petrov-Galerkin methods based on finite elements were considered, for instance, in [50][51][52][53][54]. In wavelet methods as well, the solution is expanded in a compactly supported basis that has compression properties. Wavelet methods based on Chebyshev or spline wavelets were constructed, for instance, in [55][56][57][58][59].
All these methods were applied to solve a variety of fractional differential equations, i.e., linear and nonlinear differential problems, ordinary and partial differential equations, multi-term and multi-order problems. While Galerkin and Petrov-Galerkin methods are difficult to apply to nonlinear problems since the corresponding variational formulation cannot be easily obtained, collocation methods are particularly attractive for their easy implementation also in the case of nonlinear differential problems. Remark 1. The references listed above could not be exhaustive of the huge amount of literature in the field but are only intended to give a brief review of commonly used methods to approximate the solution of differential problems of fractional order. Further references can be found in the cited papers and books (see also [60,61]).
The aim of this paper is to solve nonlinear fractional differential systems by using the collocation method introduced in [62]. The method, described in Section 2.3, was successfully used to solve multi-term fractional differential equations [62], nonlinear fractional differential equations [63] and a time-fractional diffusion equation [64]. The key-idea is to use the space generated by the fractional B-splines introduced in [65] as approximating space. Fractional B-splines are piecewise polynomials of noninteger degree. Differently from classical polynomial B-splines, fractional B-splines do not have compact support but they decay rather fast toward infinity. More importantly, fractional derivatives of both classical and fractional B-splines have an exact expression that involves the generalized finite differences of B-splines of lower degree (for details and properties of fractional B-splines see Section 2.2). Thus, in the collocation step, the fractional derivative of the approximating function is approximated accurately and efficiently by an exact differentiation rule. Another advantage of this method is in that the same differentiation rule can be used to easily approximate derivatives of any order, both fractional and integer, without the need to construct ad hoc approximating formulas for differential problems of different types. Since fractional B-splines have a fast decay toward infinity, differential operators have a sparse representation in the functional space they generate. Thus, the proposed method ultimately results in the solution of a sparse linear system.
To show the effectiveness of the proposed collocation method for the solution of nonlinear fractional problems, we solved the fractional Lotka-Volterra model and a more general predator-prey model with variable coefficients (see Section 2.1 for details on the models). The numerical results displayed in Section 3 show that the method produces approximations with good accuracy also in a large interval while keeping the computational cost low. A discussion on the numerical results can be found in Section 4 while possible issues and improvements of the fractional B-spline collocation method are pointed out in Section 5.

The Fractional Order Predator-Prey Model
Among the various nonlinear fractional differential models used in the literature to describe natural phenomena, in the present paper we focus our interest on the fractional predator-prey model.
The classical predator-prey model was introduced separately by the mathematicians A.J. Lotka [66] and V. Volterra [67] at the beginning of the 20th century to describe the dynamics of ecological systems with predator-prey interactions. Starting from the archetypal Lotka-Volterra model, several generalizations were considered to describe predator-prey interactions in different fields, such as population biology, chemical kinetics, control systems, neural networks and epidemiology [68]. Recently, fractional order Lotka-Volterra models were introduced to model real-life growth phenomena since the non-local behavior of the fractional derivative can describe the memory properties of biological systems better than the usual derivative of integer order [9,22,23,[69][70][71].
In this paper, we consider a more general fractional predator-prey problem, where x(t) and y(t) are the densities of the prey and predator populations, respectively, and a(t), b(t), c(t), d(t) are positive functions representing the growth/decay rate of the populations and the competition coefficients [9,22]. In particular, a(t) is the growth rate of the prey, b(t) is the death rate of the prey, which is related to the efficiency of the predator to capture the prey, c(t) is the growth rate of the predator, which is proportional to the prey density, d(t) is the death rate of the predator. The symbol D β t f (t) denotes the Caputo derivative of fractional order of a function f having absolute integrable m-th derivative, defined as where Γ is the Euler's gamma function For details on fractional calculus and fractional derivatives see, for instance, [1,3,4,7,8].
The existence and uniqueness of the solution of the predator-prey model (1) were proved in [1] while its behavior was analyzed in [70]. When the coefficients a(t), b(t), c(t), d(t) are constant functions we retrieve the fractional Lotka-Volterra model, i.e., The stability analysis of (4) was studied in [72,73] (see also [74]).
In the following, we will construct a collocation method based on fractional B-splines to numerically solve the predator-prey models (1) and (4).

The Fractional B-splines
Let us denote by t n + , n ∈ N, the truncated power function, i.e., Then, the classical polynomial B-splines B n of degree n can be defined by the finite differences of t n + as where is the finite difference operator. Thus, the classical B-splines are piecewise polynomials of integer degree that are compactly supported in [0, n + 1], positive in (0, n + 1) and belong to C n−1 (R). Moreover, they reproduce polynomials up to degree n. (For details and further properties on polynomial B-splines see, for instance, [75]). Fractional derivatives of B-splines can be evaluated from definition (6) and the derivation rule We get In [65] the classical polynomial B-splines were generalized to noninteger power by introducing the truncated fractional power function Thus, the fractional B-splines are defined by where is the generalized finite difference operator and are the generalized binomial coefficients that extend to noninteger values, the usual definition of binomial coefficients of integer order. The fractional B-spline B α does not have compact support but it decays toward infinity as Moreover, B α is α-Hölder continuous, belongs to L 2 (R) and reproduces polynomials up to degree α [65]. The derivation rule (9) can be easily generalized to the fractional B-spline as Now, using the definition (11) we can write a unique differentiation rule for B n and B α , i.e., which holds for any B-spline both of integer and noninteger degrees, and for derivatives of any real or integer orders. We observe that the above formula expresses the derivative of a B-spline as a spline of lower degree so that it is easy to evaluate. In fact, even if the generalized difference operator ∆ β does not have compact support, the generalized binomial coefficients (13) decrease rather fast as k increases so that the series in (12) can be approximated by a truncated sum for computational purposes.
We observe that from (10)- (11) it follows that B α is a causal function with B (n) α (0) = 0 for n ∈ N\{0}, so that the Caputo fractional derivative coincides with the Riemann-Liouville fractional derivative [3] defined as As a consequence, the Riemann-Liouville derivative of a fractional B-spline can also be evaluated using the differentiation rule (16). In Figure 1a the polynomial and fractional B-splines are displayed for different values of the degree α. The picture shows that even if the fractional B-splines do not have compact support, they decay rather fast. Moreover, we observe that while polynomial B-splines are positive functions, fractional B-splines are not, even if the nonnegative part becomes smaller and smaller as α increases. In Figure 1b

The Fractional B-spline Collocation Method
In the fractional B-spline collocation method, the approximating function is a fractional spline, i.e., it is a linear combination of integer translates and dilates of a given fractional B-spline.
Let ∆ δ = {k δ, 0 ≤ k ≤ N}, N = T/δ, be a partition of equidistant knots in the interval [0, T]. We choose as a function basis for the fractional spline space a suitable set of integer translates and dilates of a fractional B-splines B α , i.e., The set of indexes N δ is chosen in order to take into account all the translates and dilates of B α that have support in [0, T]. When α = n ∈ N, due to the compact support of B n we get N δ = {k, −n ≤ k ≤ N − 1}. When α is not an integer, B α does not have compact support so that the number of basis functions belonging to B α,δ is infinite. Nevertheless, since B α decays rather fast toward infinity (cf. Figure 1), for computational purposes we can assume supp B α ≈ [0, n 0 + 1], where n 0 is a suitable truncation parameter. Thus, N δ = {k, −n 0 ≤ k ≤ N − 1}. The basis B α,δ has a number of boundary functions, that is the functions whose support is not all contained in [0, T]. More precisely, the number of boundary functions is 2n when α has an integer value while it is 2n 0 when α is not an integer. Here, the boundary functions are obtained just restricting B α (t/δ − k) to the interval [0, T]. In the following, we will denote by N δ the number of functions in the basis B α,δ . For a given partition ∆ δ , we look for two approximating functions of the form To determine the unknown coefficients {γ α,k } and {λ α,k } we ask x α,δ and y α,δ to satisfy the differential problem (1) on a set of collocation points. For the sake of simplicity, we choose the collocation points to be equidistant nodes. Denoting by η the step-size of the nodes, the set of collocation points is Thus, the discretized version of (1) is Using (19) in (20), we get where {γ α,k , λ α,k , k ∈ N δ } are the unknown coefficients. The nonlinear system (21) has 2N δ unknowns and 2(P + 1) > 2N δ equations and can be solved by the Gauss-Newton method. Let Γ = [γ α,k , k ∈ N δ ] T and Λ = [λ α,k , k ∈ N δ ] T be the unknown vectors and L = b α,k (t p ), 0 < p ≤ P, k ∈ N δ and M β = D β t b α,k (t p ), 0 < p ≤ P, k ∈ N δ be the collocation matrices of the basis functions and of their derivatives, respectively, evaluated in the collocation points t p , 0 < p ≤ P. Then, the nonlinear system (21) can be written in matrix form as The symbol • denotes the entrywise product of vectors. The Gauss-Newton method starts by two initial approximations Γ (0) , Λ (0) and approximates the solution of the nonlinear system by the iterative procedure where the vector with J (Γ, Λ) the Jacobian matrix of F, Here, with abuse of notation, the entrywise product V • L between the vector V and the matrix L means that V has to be intended as a matrix having as many columns as L, each column being a replica of the vector V itself. In practice, the iterative procedure is stopped when V ( ) is lower than a given tolerance. It is known that if the iterative procedure converges, its limit value is a solution of the nonlinear system. Nevertheless, the local convergence of the Gauss-Newton method is not guaranteed and particular attention should be paid in implementing the method.

Results
In this section, we will show the performances of the fractional collocation method we proposed by solving the predator-prey model described in Section 2.1. We used fractional derivatives of the same order for both predator and prey equations and solved the dynamical system for three different choices of the coefficients. In the numerical tests we used a polynomial spline of degree 4 as approximating function. Denoting by [0, T] the discretization interval, the partition for the spline space is formed by N + 1 = T/δ + 1 knots. The collocation points are P + 1 = T/η + 1. In the numerical tests we set δ = 1/2 8 and η = 1/2 9 .

The Lotka-Volterra Model
First of all, we considered the Lotka-Volterra model of fractional order, i.e., the predator-prey model with constant coefficients given in (4). We set a = 2 for the prey growth rate, b = 1 for the prey death rate, c = 1 for the predator growth rate, d = 3 for the predator death rate, and set x 0 = 1 and y 0 = 2 for the initial values of the prey and predator populations, respectively. We used the proposed collocation method to solved the problem for β = 1, 0.9, 0.8, 0.5, 0. 25

The Predator-Prey Model with Variable Coefficients
In the second set of tests we solved the predator-prey problem (1) for two sets of variable coefficients: and In the first case, the growth rates of both prey and predator increase with time while the death rates are constant. In the second case, the growth rates are constant while the death rates increase with time. In the numerical tests we set x 0 = 1.3, y 0 = 0.6 and β = 1, 2/3, 1/2, 1/3 (cf. [22]). The numerical solutions in the interval [0,2] are displayed in Figures 7 and 8.

Discussion
The numerical simulations in the previous section show that, as expected, the fractional order influences the population growth. In particular, the solution of the fractional Lotka-Volterra model is not periodic and the predator and prey populations reach the equilibrium point after an intial oscillatory behavior: the lower the derivative order is, the faster the equilibrium point is reached. As for the predator-prey model with variable coefficients, in Case (i) since the growth rate of both predator and prey increases with time, the two populations increase for a prolonged amount of time. Before increasing, the prey population has an initial decrease while the predator population reaches a local maximum. The order of the fractional derivative influences the time after which the two populations start to increase: the lower the derivative order, the shorter the time. In Case (ii) since the predator and prey death rates increase, after an initial increase both populations decrease for long times. In addition, in this case the initial increasing is influenced by the order of the fractional derivative being shorter when the fractional derivative is lower.
Since the predator-prey model can not be solved exactly, to evaluate the performance of the proposed method we compare the numerical results in the previous section with the approximate solutions obtained by the methods proposed in [22,72]. In particular, in [72] the classical Lotka-Volterra model was solved using the the predictor-corrector method proposed in [32,33] while the predator-prey model with variable coefficients was solved in [22] by the homotopy perturbation method [21,25]. Figures 2-4 show that the approximate solutions of the fractional Lotka-Volterra model produced by our collocation method are in good agreement with those ones displayed in [72] where just the values β = 1, 0.9, 0.8 were considered. In the case of derivative of integer order, the numerical solution in Figure 2 is periodic and does not suffer from instabilities as in the case of the predictor-corrector method. In fact, the phase state diagram in Figure 2b is a closed curve while the diagram in ( Figure 6 [72]) is a slightly increasing orbit. As stated by the theory in [72], when the derivative has fractional order, the numerical solution tends to the equilibrium point having coordinates (3, 2) (cf. Figures 3  and 4). The little inaccuracy-less than 0.7 × 10 −1 -in the approximation of the y coordinate of the equilibrium point is of the same order as in [72]. Figures 5 and 6 show that the error in approximating the equilibrium point increases when β decreases.
For the predator-prey model with variable coefficients, the numerical solutions displayed in Figure 7, which correspond to Case (i), are very similar to the approximate solutions obtained in [22]. Instead, the numerical solutions for Case (ii) shown in Figure 8 are rather different from the approximate solution displayed in [22]. The reason is that just four terms of the expansion series are used to evaluate the approximations while the homotopy perturbation method requires more terms to well approximate the solution in large interval (cf. [70] where a modification of the homotopy perturbation method is used and the solution is approximated using eight terms in the series). In fact, our numerical solution and the approximate solution in [22] are comparable near the origin.
To give an idea of the accuracy of the numerical solutions displayed in Section 3, we solved the predator-prey problem for decreasing values of the distance δ between the knots of the partition ∆ δ and evaluate the infinity norm of their difference. The values of the norm, listed in Tables 1-3, show that the difference decreases when δ decreases reaching a value in the order of 10 −12 when β = 1 and a value in the order of 10 −3 or less when β is of fractional order. As a consequence, we can argue that the numerical approximations obtained by the collocation method is sufficiently accurate also in the case of derivatives of fractional order and, in particular, of small order, which are more difficult to approximate.   Then, we solved the three models for different values of η. Several numerical tests not displayed here show that the numerical solution depends very little on η so that the accuracy does not increase sufficiently fast when reducing the time-step.

Conclusions
We used a new collocation method to numerically solve some predator-prey models of fractional order. The numerical results show that the method is accurate and efficient. A theoretical proof of the convergence of the method when δ and η go to zero follows from classical arguments in approximation theory (cf. [62]).
Some issues should be further analyzed to improve the method. First of all, the truncation of the basis functions {B α (t/δ − k), k ∈ N δ }, when restricted to the interval [0, T], produces few instabilities near the boundaries of the discretization interval; this could increase the ill-conditioning of the nonlinear system (21). A different approach should be used to generate stable bases on the interval (see, for instance, [76]). Another issue could arise from the truncation of the fractional B-spline that is necessary for numerical purposes. Finally, a different distribution of the collocation nodes could improve the accuracy of the method. These details will be analyzed in a forthcoming paper.

Conflicts of Interest:
The author declares no conflict of interest.