A computational approach with residual error analysis for the fractional-order biological population model

ABSTRACT In this study, a fractional Bernstein series solution method has been submitted to solve the fractional-order biological population model with one carrying capacity. The numerical method has been implemented by an effective algorithm written on the computer algebraic system Maple 15. An error-bound analysis is performed by using a process similar to the RK45 method. An error estimation technique relating to residual function is presented. In the numerical application, the variations in the population of prey and predator with respect to time and situations of these two species relative to each other are plotted. The outputs obtained from our method are compared with the homotopy perturbation Sumudu transform method and reproducing kernel Hilbert space method. The approximate solutions gained from the Bernstein series method are consistent with those of other methods. The advantage of our method is that it requires less computational cost compared with methods involving more complex operations.


Introduction
The interaction between species is one of the most interesting dynamics in ecology and biology. Mathematical models which include these relations are categorized as mutualism, competition and predator-prey models [1]. In predator-prey relations, one species feeds on another species. This interaction can take the form of resource-consumer, plant-herbivore and parasite-host systems [2]. There have been many theoretical and numerical studies on these dynamics over the years [3][4][5][6].
In recent years, fractional differential equations (FDEs) are one of the more favoured tools; it models complicated behaviours of dynamical systems from different fields. These include mechanics [7], electricity [8], electrochemistry [9], economy [10], biology [11,12] and epidemiology [13]. Since non-integer order derivatives contain integer orders, FDEs are a generalization of ordinary differential equations (ODEs). Systems modelled with FDEs reveal real-life phenomena more accurately. This is because they give more detailed information about the after-effects and memory of the system [14][15][16]. FDEs have gained great importance and popularity due to this feature, particularly in the field of biology and ecology modelling [17][18][19][20]. Recently, many studies have been published on systems modelled by using distinct fractional derivatives definitions. A new fractional model was proposed for the human liver involving the Caputo-Fabrizio derivative [21]. Boukhouima et al. [22] examined a fractional HIV population model and fractional cellular model theoretically by using the Caputo derivative, Ghanbari et al. [23], Bahaa and Hamiaz [24] analyzed the behaviour of a dynamic system that includes the Atangana-Baleanu fractional derivative.
In this work, a two-dimensional fractional-order biological population model with one carrying capacity is examined. The model is described as: with the initial conditions Here, 0 ≤ t ≤ b and D α t represent the Caputo fractional derivative of order α where 0 < α ≤ 1. x(t) and y(t) are respectively the population densities of the prey and the predator, a 1 is the prey's intrinsic growth and K 1 is the carrying capacity. a 2 denotes the predator's growth rate, β 1 and β 2 indicate the positive competition coefficients.
For the research, it is assumed that all variables and parameters of the model are non-negative at a random point in time because the system (1)-(2) is defined as the population of the species. Population densities x(t)and y(t) which are the solutions of the system will be restricted to the non-negative R 2 + . For case α = 1, the system has analytical solutions. From the two equations which give the solutions the positivity of x(t) and y(t) is observed (see Theorem 1 in [31]).
The main aim of this article is to introduce the Bernstein series solution method which comprises Bernstein polynomials and collocation method [32,33] for model (1)-(2) which was previously solved by homotopy perturbation Sumudu transform method (HPSTM) [31], and which gives a semi-analytic solution and reproducing kernel Hilbert space method (RKHSM) [34].
The other parts of the paper are arranged as follows. In section 2, some basic and necessary definitions are given about fractional derivatives. In Section 3, the form of Bernstein polynomials, matrix relations of expressions in the model (1)-(2) are given and the procedure is constituted. Error analysis of the method is presented in Section 4. In this section, first, the residual correction procedure is clarified. This process estimates the absolute error and also obtains new approximate solutions. Additionally, for the second analysis, another procedure that includes similar steps to RK45 methods is performed. The advantage of these methods is that both can be used even if the exact solution is not known. In Section 5, numerical outcomes are submitted by graphics and tables to show the effectiveness and simple applicability of the method. In the last section, the results are interpreted.

Some basic information regarding fractional calculus
In this section, some necessary definitions and some properties of the Caputo fractional derivatives are reminded.

Description of the method
In this section, we introduce the Bernstein series solution method to get the approximate solutions of the model (1)- (2). This technique depends on Bernstein polynomials and collocation methods. For this purpose, each term of the system (1)-(2) and the initial conditions (3) are formulated as matrix representations through the Bernstein polynomials.
The Bernstein polynomials of degree n defined on the interval [0, L] are given by If t β is substituted into t in the relation (6), the fractional Bernstein polynomials of degree n are attained as Thus, the fractional Bernstein series solutions are represented in the following forms where 0 < β < 1, c k and h k , k = 0, 1, . . . , n are the Bernstein coefficients.
If the solutions (8), (9) wanted to be specified in the matrix form, the following relations are given. where The expression of φ(x) can be written as where Hence, the fundamental matrix relations of solutions are in the form as Similarly, α order Caputo derivative of solutions stated in Equations (10) and (11) are written as where Here components of a matrix X * (t) are calculated by using property 3 in Definition 2.2. Hence, the fractionalorder derivative of x n,β (t) and y n,β (t) are converted into a matrix form: If the relations (13), (14), (16), (17) are substituted into the system (1)-(2), the following system of matrix equations are acquired as: +b 1 (C S X(t))(H S X(t)) = 0 HSX * (t) + a 2 HSX(t) − b 2 (C S X(t))(H S X(t)) = 0. (18) To find the matrix form of the initial conditions (13) and (14) are used in Equation (3) and therefore, matrix forms of conditions are gained as: Finally, by substituting the collocation points defined by into Equation (18), the system of 2(n + 1) equations with the unknown coefficients c k and h k, k = 0, 1, . . . , n is attained. By solving the system, the fractional Bernstein series solutions x n,β (t) and y n,β (t) are obtained.

Error analysis
In this part of the study, two error analyses for the method are given. One of these processes is the residual correction procedure. The second is similar to the error analysis of the RK45 method.
If a similar process is applied for R 2n and equation (2), the relation is obtained as: where e 2n,β (t) = y(t) − y n,β (t).
Since the exact solutions and the fractional Bernstein series solutions satisfy the initial conditions, the conditions for the errors are If e n,β (t) − e m n,β (t) < x(t) − x n,β (t) , then x m n,β is a better approach than x n,β . On the other hand, we can estimate the error e n,β (t) by e m n,β whenever e n,β (t) − e m n,β (t) < . In practice, the absolute errors can be estimated by e m n,β for m > n. As the second analysis, we get x n,β and x s,β which are any two Bernstein series solutions of (1)- (2). Suppose By using triangle inequality, we obtain (29) Note that e n,β (t) can be bounded by the difference between any two approximate solutions. Thus, if e n,β (t) is a monotone sequence, then x n+1,β (t) − x n,β (t) may nearly bound each of consecutive absolute errors.

Numerical application
In this section, the fractional biological model described in (1)-(2) is examined for different values of fractional order α. In the computations, we use Maple 15 computer algebraic system. The following norms are employed to measure the size of the error: We apply the method to the model for L = 1 and the dimensionless parameters a 1 = a 2 = 0.05, b 1 = 0.04, b 2 = 0.01, K 1 = 20, ξ 1 = 20, ξ 2 = 15. As α parameters vary, it can be observed how the population densities of prey (x) and predator (y) change in time and how their situations are observed in relation to each other in the same period. In Figures 1 and 2, for α = 0.65, 0.75, 0.90, 1 the variations of population densities of x(t) and y(t) are plotted respectively. It can be inferred from Figure 1 that the prey population decreases with respect to time t while α decreases. Figure 2 displays that the predator population increases as time progresses and α decreases. The interactions of two species for distinct α parameters are represented in Figures 3-8. Also in these figures, reproducing kernel Hilbert space method (RKHSM), homotopy perturbation Sumudu transform method (HPSTM) and the Bernstein series method are compared. From these illustrations, it can be concluded that the results of all numerical meth-    ods are consistent with each other. It is also inferred that the numerical results depend continuously on α.
The results for different values of n are tabulated in Table 1 to show the dependency of the method on the value of n. In Table 1, we give the norms of the difference between consecutive approximate solutions and the estimations of the absolute errors by residual corrections. From this table, it can be said that both error estimation methods support each other.
The accuracies of the approximate solutions are checked by measuring R in , i = 1, 2. Since the exact solutions of the model are unknown, we give the accuracies of the approximate solutions, estimated absolute errors in Tables 2 and 3      increase, residues decrease except n = m. Also, the running time of the programme gets bigger as n increases.

Conclusions
In this paper, we propose a numerical method to solve the fractional biological model with carrying capacity. This is a procedure that depends on Bernstein polynomials, collocation method and Caputo fractional derivative, also named as the Bernstein series solution method. Then, we give two techniques to estimate and bound the absolute error. One of these techniques is the residual error correction and the other is the difference between any approximate solutions. It should be stated that these two techniques can be applied to the problem even if the exact solution is unknown. The mentioned method and error analysis procedures are applied to a test example. Error estimations obtained from two different procedures overlap as can be seen in Table 1. From Tables 2-4, we can observe that increasing node numbers yield decreasing residues but it results in more computational costs.   When Figures 1 and 2 are analyzed, we can see the population of prey x(t) decrease and the population of predator y(t) increase with respect to time t. In these two figures, the numerical solutions converge to analytical solutions as α converges to 1. Figures 3-8 show the interactions of the prey and predator species for different α values. The approximate solutions attained from the mentioned method are also compared with the results obtained from HPSTM and RKHSM for the same α values. We can deduce from the computations and analysis that our method gives compatible approximate solutions with the semi-analytical methods. Also, the method includes simpler equations in comparison with the other methods.