An Improved Collocation Approach of Euler Polynomials Connected with Bernoulli Ones for Solving Predator-Prey Models with Time Lag

This paper deals with an approach to obtaining the numerical solution of the Lotka–Volterra predator-prey models with discrete delay using Euler polynomials connected with Bernoulli ones. By using the Euler polynomials connected with Bernoulli ones and collocation points, this method transforms the predator-prey model into a matrix equation. The main characteristic of this approach is that it reduces the predator-prey model to a system of algebraic equations, which greatly simpliﬁes the problem. For these models, the explicit formula determining the stability and the direction is given. Numerical examples illustrate the reliability and eﬃciency of the proposed scheme.


Introduction
In recent years, population models in various fields of mathematical biology have been proposed and studied extensively. Predator-prey interaction is the fundamental structure in population dynamics. Understanding the dynamics of predator-prey models will be very helpful for investigating multiple species interactions [1][2][3][4][5][6]. e first predator-prey model with delay was proposed by Volterra and Brelot [7]. Let p 1 (t) and p 2 (t) denote the population density of prey and predator at time t. In this work, we first, consider the modified version of the predator-prey model with delay [8]: where r 1 > 0 is the growth rate of the prey in the absence of predators, a 11 > 0 denotes the self-regulation constant of the prey, a 12 > 0 describes the predation of the prey by predators, r 2 > 0 is the death rate of predators in the absence of the prey, a 21 > 0 is the conversion rate for the predators, and a 22 ≥ 0 describes the transpacific competition among predators. F(s) and G(s) are nonnegative continuous delay kernels defined and integrable on [0, ∞), which weigh the contribution of the predation occurred in the past to the change rate of the prey and predators, respectively. Detailed study on stability and bifurcation of system (1) can be found in Cushing [9]. When F(t − s) � δ(t − s − τ 1 ) and G(t − s) � δ(t − s − τ 2 )-where δ is the delta function and τ 1 ≥ 0 and τ 2 ≥ 0 are the hunting delay and the maturation time of the prey species, respectively-system (1) reduces to the following Lotka-Volterra predator-prey model with two discrete delays: We will find the approximate solutions of the system (2) with initial conditions: where α 1 and α 2 are positive numbers. e exact solution for the predator-prey models in the general form does not exist; therefore, numerical methods are needed to computing the population densities of the prey and predator. Some methods [10][11][12] are proposed to handle an approximate solution of the predator-prey system, such as finite-difference, finite element, and spectral methods. e organization of this paper is as follows. In Section 2, we study the stability of the predator-prey model. In Section 3, matrix relations for Euler polynomials connected with Bernoulli ones will be discussed. In Section 4, we shall present our method in detail. Finally, we give a few numerical examples in Section 5.

Stability
In this section, we review some results on the stability of the predator-prey model. To begin with, we first consider system (2) appearing in the interspecific interaction terms of both equations: Assume that f: R × C ⟶ R and g: R × C ⟶ R satisfy the following assumptions: : f and g are continuously differentiable such that zf/zp 1 < 0, zf/zp 2 < 0, zg/zp 1 > 0, and zg/zp 2 > 0.
erefore, we have the following theorem.
Theorem 1. Suppose that f and g satisfy the assumptions (S 1 ) and (S 2 ). Define en, the positive equilibrium (p * 1 , p * 2 ) of the delayed predator-prey system (4) is absolutely stable if and only if a d + bc > 0.
He [14] and Lu and Wang [15] showed that the positive equilibrium E * � (p * 1 , p * 2 ) of system (2) is globally stable. e above stability result depends on the assumption that a 11 a 22 − a 12 a 21 > 0. If a 11 a 22 − a 12 a 21 < 0, then Faria [16] proved the following result.

Matrix Relations for Euler Polynomials
Bernoulli polynomials B n (x) are usually defined from the generating function for |t| < 2π: and Bernoulli numbers β n : � B n (0) are contained in many mathematical formulas such as the Taylor expansion in a neighborhood of the origin of trigonometric and hyperbolic tangent and cotangent functions, and they can be obtained by the generating function [17]: Euler studied Bernoulli polynomials, and these polynomials are employed in the integral representation of differentiable periodic functions and play an important role in the approximation of such functions using polynomials. Many early Euler and Bernoulli polynomial implementations can be contained in [18][19][20]. Euler polynomials are strictly related to Bernoulli and are used in Taylor's expansion in the trigonometric and hyperbolic secant function districts. Recursive computation of Bernoulli and Euler polynomials can be obtained by using the following formulas: Euler polynomials E n (t) can also be defined as polynomials of degree n ≥ 0 satisfying the conditions: e above conditions express a recurrence relation to Euler polynomials. An explicit formula for E n (t) is represented as where β k � B k (0) is the first Bernoulli numbers for each k � 0, 1, . . . , n [21].
, and h(t) be an arbitrary element in H. Since P H is a finite dimensional vector space, h(t) has the unique best approximation out of P H such as p(t), that is, Since p(t) ∈ P H , there exist unique coefficients where P � p 0 , p 1 , . . . , p N , such that Note that, from property (12) of Euler polynomials, we can represent E(t) in a matrix form as follows: where International Journal of Differential Equations is

Numerical Solution of the Lotka-Voltra Predator-Prey Models with Discreet Delay Using Euler Operational Matrices
Consider the Lotka-Voltra predator-prey models with discreet delay (2). For implementing the Euler operational matrix on the predator-prey model, we first approximate p i (t) and p i ′ (t), for i � 1, 2, by equations (14)-(17), as follows: where P i is the unknown (N + 1) vector of coefficients.

Remark 1.
Recall that the Hadamard product of two matrices, say A and B, of the same dimension m × n, is denoted by A ∘ B, defined to be an m × n-matrix, in which the (i, j) th entry is equal to the product of (i, j) th entries of A and B, i.e., Bearing in mind the previous remark, we can simplify the delay part of system (2) to the matrix form as where B 0 � I and for any τ ≠ 0, B τ is the Hadamard product C ∘ G (i.e., the product of entry by entry. . ..) such that for In other words, we may evaluate entries of B τ as follows: By the mixed-product property of a Kronecker product, we know that for matrices A, B, C, and D, which are of the size that matrix products AC and B D are well-defined, the following relation holds true: Now, we have to estimate all expressions consisting the products of p 1 (t − τ 2 ) and p 2 (t − τ 1 ), existing in system (2) (note that, for τ k � 0, the products p i (t − τ k )p j (t), p i (t)p j (t − τ k ), and p i (t)p j (t) can be approximated in matrix form, where i, j ∈ 1, 2 { }). In order to do this, let en, by using equations (19) and (20), we have where P i,j � P i ⊗ P j ,

⎧ ⎨ ⎩ (26)
Suppose that en, the fundamental matrix equations of system (2) can be written in the following form: where P � P 1 , P 2 ,

Implementing the Collocation Method. Equation (28)
gives us 2 nonlinear equations. Since the number of unknowns for each vector P 1 and P 2 in (15) is (N + 1) and our system has 4 equations, the total number of unknowns is 2(N + 1); then, we collocate equation (28) at N Newton-Cotes points given as and hence we have equation (28) and initial conditions: PM t s + PA t s + PC t s � f t s ; s � 1, 2, . . . , N, After solving linear system (31), we get P 1 and P 2 . For collocating equation (31), we have used the Newton-Cotes points because of their simplicity and their good utility in our implementation regarding the speed and accuracy of answers. However, we could use other points such as the Gauss points, Clenshaw-Curtis points, and Lobatto points.

Numerical Examples
To test the method in terms of its precision and efficacy, we turn our attention to show three numerical performance results. Analytic approximate solutions of high-order accuracy are presented in most of the cases. roughout the calculations, the absolute error is defined by and maximum absolute error is defined by e computations associated with these examples were performed using MATLAB.
Since the truncated series (17) is approximate solutions of system (31), when the function p 1,N (t), p 2,N (t) and its derivatives are substituted in system (2), the resulting equation must be satisfied approximately; that is, for t s , s � 1, 2, . . . , N, and Max |R 1,N (t s ) − R 2,N (t s )| ≤ 10 − r s , (r s positive integer). If max10 − r s � 10 − r (r positive integer) is prescribed, then the truncation limit N is increased until the difference R i,N (t s ) at each of the points becomes smaller than the prescribed 10 − r . All computations were carried out using MATLAB. Example 1. We consider the following Lotka-Volterra predator-prey model with delays: which has also a positive equilibrium E * � (p * 1 � (2/3), p * 2 � (1/3)) [22]. e numerical simulations for Example 1 are shown in Figures 1 and 2. International Journal of Differential Equations

International Journal of Differential Equations
Example 2. Consider the Lotka-Volterra predator-prey models with discrete delay given in [23] by with initial conditions We selected our second example from [23], which solved this predator-prey interaction by the Bessel collocation approach. We compared our results with the Bessel collocation results, and the outcomes are tabulated in Table 1. For this example too, the improved Collocation approach of Euler polynomials connected with Bernoulli ones is incisive.
is is because for the same number of basis functions, it obtains better results. Table 1 shows the present method results for Example 2 in comparison with method of [23]. e superiority of Euler polynomials operational matrices method compared with Bessel collocation approach is clear here, which is because with the same number of basis functions, we get very better results.
Example 3. Suppose in a closed ecosystem there are only 2 types of animals: the predator and the prey. ey form a simple food-chain where the predator species hunts the prey species, while the prey grazes vegetation. e size of the 2 populations can be described by a simple system of 2 nonlinear first-order differential equations: which a set of fixed positive constants. A: the growth rate of the prey, B: the rate at which predators destroy the prey, C: the death rate of predators, and D: the rate at which predators increase by consuming the prey. A prey population p 1 (t) increases at a rate (dp 1 (t)/dt) � Ap 1 (t) (proportional to the number of prey) but is simultaneously destroyed by predators at a rate (dp 1 (t)/dt) � −Bp 1 (t)p 2 (t) (proportional to the product of the number of prey and predators). A predator population p 2 (t) decreases at a rate (dp 2 (t)/dt) � −Cp 2 (t) (proportional to the number of predators), but increases at a rate (dp 2 (t)/dt) � Dp 1 (t)p 2 (t) (again proportional to the product of the numbers of prey and predators). By eorem 1, we have the the positive equilibrium E * � (p * 1 � (C/D), p * 2 � (A/B)). For this example, the stability of each of the three steady states can be assessed more formally using the approach discussed. e absolute errors for N � 7, 9 are estimated by R 1,N and R 2,N and are presented in Tables 2 and 3. We see that if N increases, then the absolute errors decrease more rapidly.

Conclusion
In this paper, we have proposed a numerical approach for solving the Lotka-Volterra predator-prey models with discrete delay by utilizing Euler polynomials connected with Bernoulli ones. e cognate matrices and mixed-product property of Kronecker products, besides the collocation method, have been utilized for transforming a predator-prey system to a linear system of algebraic equations that can be solved facilely. To the best of our cognizance, this is the first work coalescing the Euler polynomial connected with the Bernoulli polynomial and collocation points for solving Lotka-Volterra predator-prey models. Finally, numerical examples reveal that the present method is very precise and convenient for solving predator-prey models.

Data Availability
No data were used to support this study.

Conflicts of Interest
e authors declare that they have no conflicts of interest.