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

Authors: Behrooz Basirat; Hamid Reza Elahi
Journal: International Journal of Differential Equations
 (2020)

---

## Abstract

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 simplifies the problem. For these models, the explicit formula determining the stability and the direction is given. Numerical examples illustrate the reliability and efficiency of the proposed scheme.


---

## Body

## 1. 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–6]. The first predator-prey model with delay was proposed by Volterra and Brelot [7]. Let p1t and p2t 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]:(1)p1′t=p1tr1−a11p1t−a12∫−∞tFt−sp2sds,p2′t=p2t−r2+a21∫∞tGt−sp1sds−a22p2t,where r1>0 is the growth rate of the prey in the absence of predators, a11>0 denotes the self-regulation constant of the prey, a12>0 describes the predation of the prey by predators, r2>0 is the death rate of predators in the absence of the prey, a21>0 is the conversion rate for the predators, and a22≥0 describes the transpacific competition among predators. Fs and Gs 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 Ft−s=δt−s−τ1andGt−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:(2)p1′t=p1tr1−a11p1t−a12p2t−τ1,p2′t=p2t−r2+a21p1t−τ2−a22p2t.We will find the approximate solutions of the system (2) with initial conditions:(3)p10=α1,p20=α2,where α1 and α2 are positive numbers.The 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–12] are proposed to handle an approximate solution of the predator–prey system, such as finite-difference, finite element, and spectral methods.The organization of this paper is as follows. In Section2, 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.


## 2. 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:(4)p1′t=p1tfp1t,p2t−τ1,p2′t=p2tgp1t−τ2,p2t,where τi≥0, i=1,2, is a constant. Denote C=C−τ,0×R, where τ=maxτ1,τ2. Assume that f:R×C⟶R and g:R×C⟶R satisfy the following assumptions:S1: there exists a point p1∗,p2∗ with p1∗,p2∗>0 for which fp1∗,p2∗=gp1∗,p2∗=0,S2: f and g are continuously differentiable such that ∂f/∂p1<0, ∂f/∂p2<0, ∂g/∂p1>0, and ∂g/∂p2>0.AssumptionS1 ensures that p1∗,p2∗ is a positive equilibrium point of system (4). Let ϕt=p1t−p1∗ and ψt=p2t−p2∗. Then, the linearized system at p1∗,p2∗ is(5)ϕ′t=∂f∂p1p1∗,p2∗p1∗ϕt+∂f∂p2p1∗,p2∗p1∗ψt−τ1,ψ′t=∂g∂p1p1∗,p2∗p2∗ϕt−τ2+∂g∂p2p1∗,p2∗p2∗ψt.Therefore, we have the following theorem.Theorem 1.
Suppose thatf and g satisfy the assumptions S1 and S2. Define(6)a=∂f∂p1p1∗,p2∗,b=∂f∂p2p1∗,p2∗,c=∂g∂p1p1∗,p2∗,d=∂g∂p2p1∗,p2∗.Then, the positive equilibriump1∗,p2∗ of the delayed predator-prey system (4) is absolutely stable if and only if ad+bc>0.For proof, see [13].By Theorem1, we have the following condition. If(7)r1a21−r2a11>0,then system (2) has a positive equilibrium E∗=p1∗,p2∗, where(8)p1∗=r1a22+r2a12a11a22+a12a21,p2∗=r1a21−r2a11a11a22+a12a21.The conditionad+bc>0 becomes a11a22−a12a21>0. Thus, by Theorem 1, we have the following result on the stability of E∗=p1∗,p2∗.Corollary 1.
If condition (7) is satisfied, i.e., if the positive equilibrium E∗=p1∗,p2∗ of system (2) exists, then it is absolutely stable if and only if a11a22−a12a21>0.He [14] and Lu and Wang [15] showed that the positive equilibrium E∗=p1∗,p2∗ of system (2) is globally stable. The above stability result depends on the assumption that a11a22−a12a21>0. If a11a22−a12a21<0, then Faria [16] proved the following result.Corollary 2.
Ifa11a22−a12a21<0 and a11a22≠0, then there is a critical value τ2,0>0, such that E∗ of system (2) is asymptotically stable when τ2<τ2,0 and unstable when τ2>τ2,0.


## 3. Matrix Relations for Euler Polynomials

Bernoulli polynomialsBnx are usually defined from the generating function for t<2π:(9)Gx,t≔textet−1=∑n=0∞Bnxtnn!,and Bernoulli numbers βn:=Bn0 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]:(10)tet−1=∑n=0∞βntnn!.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–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:(11)∑k=0n−1nkBkt=ntn−1,n=2,3,…,Ent+∑k=0nnkEkt=2tn,n=1,2,….Euler polynomialsEnt can also be defined as polynomials of degree n≥0 satisfying the conditions:(1)
Em′t=mEm−1t,m≥1(2)
Emt+1+Emt=2tm,m≥1The above conditions express a recurrence relation to Euler polynomials. An explicit formula forEnt is represented as(12)Ent=1n+1∑k=1n+12−2k+1n+1kβktn+1−k,where βk=Bk0 is the first Bernoulli numbers for each k=0,1,…,n [21].


### 3.1. Function Approximation

LetH=L20,1, and assume that E0t,E1t,…,ENt be the set of Euler polynomials, and PH=spanE0t,E1t,…,ENt, and ht be an arbitrary element in H. Since PH is a finite dimensional vector space, ht has the unique best approximation out of PH such as p˜t, that is,(13)ht−p˜t≤ht−pt,∀pt∈PH.Sincep˜t∈PH, there exist unique coefficients p0,p1,…,pN, such that(14)ht≃p˜t=∑n=0NpnEnt=PEt,h′t≃p˜′t=∑n=0NpnEn′t=PE′t,where(15)P=p0,p1,…,pN,Et=E0t,E1t,…,ENtT,E′t=DEt,such that(16)D=00⋯0010⋯0002⋯00⋮⋮⋱⋮⋮00⋯N0.Note that, from property (12) of Euler polynomials, we can represent Et in a matrix form as follows:(17)Et=MTt,where(18)M=100⋯02−23222β210⋯02−24333β32−23323β21⋯0⋮⋮⋮⋱⋮2−2N+2N+1N+1N+1βN+12−2N+1N+1NN+1βN2−2NN+1N−1N+1βN−1⋯1,is an N+1×N+1 matrix, and Tt=1,t,t2,…,tNT.


## 3.1. Function Approximation

LetH=L20,1, and assume that E0t,E1t,…,ENt be the set of Euler polynomials, and PH=spanE0t,E1t,…,ENt, and ht be an arbitrary element in H. Since PH is a finite dimensional vector space, ht has the unique best approximation out of PH such as p˜t, that is,(13)ht−p˜t≤ht−pt,∀pt∈PH.Sincep˜t∈PH, there exist unique coefficients p0,p1,…,pN, such that(14)ht≃p˜t=∑n=0NpnEnt=PEt,h′t≃p˜′t=∑n=0NpnEn′t=PE′t,where(15)P=p0,p1,…,pN,Et=E0t,E1t,…,ENtT,E′t=DEt,such that(16)D=00⋯0010⋯0002⋯00⋮⋮⋱⋮⋮00⋯N0.Note that, from property (12) of Euler polynomials, we can represent Et in a matrix form as follows:(17)Et=MTt,where(18)M=100⋯02−23222β210⋯02−24333β32−23323β21⋯0⋮⋮⋮⋱⋮2−2N+2N+1N+1N+1βN+12−2N+1N+1NN+1βN2−2NN+1N−1N+1βN−1⋯1,is an N+1×N+1 matrix, and Tt=1,t,t2,…,tNT.


## 4. 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 pit and pi′t, for i=1,2, by equations (14)–(17), as follows:(19)pit≃PiEt=PiMTt,pi′t≃PiE′t=PiDEt=PiDMTt,where Pi is the unknown N+1 vector of coefficients.Remark 1.
Recall that the Hadamard product of two matrices, sayA and B, of the same dimension m×n, is denoted by A∘B, defined to be an m×n-matrix, in which the i,jth entry is equal to the product of i,jth entries of A and B, i.e., A∘Bi,j=Ai,jBi,j.
Bearing in mind the previous remark, we can simplify the delay part of system (2) to the matrix form as(20)Et−τ=MTt−τ=MBτTt,where B0=I and for any τ≠0, Bτ is the Hadamard product C∘G (i.e., the product of entry by entry….) such that for 1≤j≤i≤N+1, Ci,j=i−1j−1 and Gi,j=−τi−j, and for i<j, Ci,j=Gi,j=0. In other words, we may evaluate entries of Bτ as follows:(21)Bτi,j=i−1j−1−τi−j,i≥j,0,i<j.Remark 2.
Note that IfA=aijm×n and Bp×q are two matrices, then the Kronecker product A⊗B is defined to be the following mp×nq block matrix:(22)A⊗B=a11B⋯a1nB⋮⋱⋮am1B⋯amnB.
By the mixed-product property of a Kronecker product, we know that for matricesA, B, C, and D, which are of the size that matrix products AC and BD are well-defined, the following relation holds true:(23)AC⊗BD=A⊗BC⊗D.
Now, we have to estimate all expressions consisting the products ofp1t−τ2 and p2t−τ1, existing in system (2) (note that, for τk=0, the products pit−τkpjt, pitpjt−τk, and pitpjt can be approximated in matrix form, where i,j∈1,2). In order to do this, let i,j⊆1,2, i′=3−i, and j′=3−j. Then, by using equations (19) and (20), we have(24)pit−τi′pjt−τj′≃PiEt−τi′PjEt−τj′=PiMBτi′Tt⊗PjMBτj′Tt=Pi⊗PjM⊗MBτi′⊗Bτj′Tt⊗Tt=Pi,jMBτi′,τj′Tt,where(25)Pi,j=Pi⊗Pj,M=M⊗M,Bτi′,τj′=Bτi′⊗Bτj′,Tt=Tt⊗Tt.
With the above notations, we can rewrite system (2) as(26)P1DMTt=r1P1MTt−a11P1,1MTt−a12P1,2MB0,τ1Tt,P2DMTt=−r2P2MTt+a21P2,1MB0,τ2Tt−a22P2,2MTt.
Suppose that(27)M1t=DMTt−r1MTt,M2t=DMTt+r2MTt,A1t=a12MB0,τ1Tt,A2t=−a21MB0,τ2Tt,C1t=a11MTt,C2t=a22MTt.
Then, the fundamental matrix equations of system (2) can be written in the following form:(28)PMt+P¯At+P^Ct=ft,where(29)P=P1,P2,Mt=M1t00M2t,P¯=P1,2,P2,1,At=A1t00A2t,P^=P1,1,P2,2,Ct=C1t00C2t,ft=0,01×2.


### 4.1. Implementing the Collocation Method

Equation (28) gives us 2 nonlinear equations. Since the number of unknowns for each vector P1 and P2 in (15) is N+1 and our system has 4 equations, the total number of unknowns is 2N+1; then, we collocate equation (28) at N Newton–Cotes points given as(30)ts=bNs,s=1,2,…,N,and hence we have equation (28) and initial conditions:(31)PMts+P¯Ats+P^Cts=fts;s=1,2,…,N,P1MT0=α1,P2MT0=α2.After solving linear system (31), we get P1 and P2. 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.


## 4.1. Implementing the Collocation Method

Equation (28) gives us 2 nonlinear equations. Since the number of unknowns for each vector P1 and P2 in (15) is N+1 and our system has 4 equations, the total number of unknowns is 2N+1; then, we collocate equation (28) at N Newton–Cotes points given as(30)ts=bNs,s=1,2,…,N,and hence we have equation (28) and initial conditions:(31)PMts+P¯Ats+P^Cts=fts;s=1,2,…,N,P1MT0=α1,P2MT0=α2.After solving linear system (31), we get P1 and P2. 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.


## 5. 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. Throughout the calculations, the absolute error is defined by(32)Errort=pit−y1t,and maximum absolute error is defined by(33)Error=Maxpiti−y1ti;ti=i1000,i=1,2,..,1000.The computations associated with these examples were performed using MATLAB.Since the truncated series (17) is approximate solutions of system (31), when the function p1,Nt, p2,Nt and its derivatives are substituted in system (2), the resulting equation must be satisfied approximately; that is, for ts,s=1,2,…,N,(34)R1,Nts=p1,N′ts−p1,Nts1−p2,Nts≅0,R2,Nts=p2,N′ts−p2,Nts−1+p1,Nts≅0,and MaxR1,Nts−R2,Nts≤10−rs, (rs positive integer). If max10−rs=10−r (r positive integer) is prescribed, then the truncation limit N is increased until the difference Ri,Nts 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:(35)p1′t=p1t1−p1t−p2t−1.7,p2′t=p2t−1+2p1t−1.8−p2t,which has also a positive equilibrium E∗=p1∗=2/3,p2∗=1/3[22]. The numerical simulations for Example 1 are shown in Figures 1 and 2.Figure 1
The stable behavior of the prey and predator populations for Example1 to the equilibrium population points. (a) The graph on p1 − t plane of Example 1. (b) The graph on p2 − t plane of Example 1.
(a)
(b)Figure 2
The phase graph of Example1 in t=500 with the equilibrium points p1∗=2/3 and p2∗=1/3.Example 2.
Consider the Lotka–Volterra predator-prey models with discrete delay given in [23] by(36)p1′t=p1t1−p2t,p2′t=p2t−1+p1t,with initial conditions(37)p10=1.3,p20=0.6.
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. This 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]. The 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.Table 1
Errors for Example2.
t
R
1
,
N
t
R
2
,
N
t
Present method
Method in [23]
Present method
Method in [23]
N
=
5
N
=
12
N
=
5
N
=
5
N
=
12
N
=
5
0.0
5.88378
e
−
8
3.56178
e
−
15
≃
10
−
5
,
10
−
10
3.415153
e
−
8
8.881278
e
−
15
≃
10
−
5
,
10
−
10
0.1
5.76570
e
−
7
1.76292
e
−
14
≃
10
−
5
,
10
−
10
1.45283
e
−
7
8.743266
e
−
14
≃
10
−
5
,
10
−
10
0.2
1.36470
e
−
6
2.87548
e
−
14
≃
10
−
5
,
10
−
10
2.87360
e
−
8
1.65978
e
−
14
≃
10
−
5
,
10
−
10
0.3
3.65207
e
−
8
4.27880
e
−
12
≃
10
−
5
,
10
−
10
3.61146
e
−
8
2.37310
e
−
13
≃
10
−
5
,
10
−
10
0.4
6.42077
e
−
6
3.62315
e
−
12
≃
10
−
5
,
10
−
10
1.54657
e
−
7
2.02147
e
−
13
≃
10
−
5
,
10
−
10
0.5
4.04666
e
−
6
5.97320
e
−
12
≃
10
−
5
,
10
−
10
5.76562
e
−
7
4.63123
e
−
13
≃
10
−
5
,
10
−
10
0.6
2.37656
e
−
6
7.28578
e
−
12
≃
10
−
5
,
10
−
10
4.63832
e
−
7
5.21534
e
−
13
≃
10
−
5
,
10
−
10
0.7
3.34357
e
−
6
7.51433
e
−
12
≃
10
−
5
,
10
−
10
5.66931
e
−
7
5.75760
e
−
13
≃
10
−
5
,
10
−
10
0.8
2.21441
e
−
5
1.99510
e
−
12
≃
10
−
5
,
10
−
10
3.65471
e
−
6
7.44423
e
−
12
≃
10
−
5
,
10
−
10
0.9
3.27256
e
−
4
5.44983
e
−
12
≃
10
−
5
,
10
−
10
2.44449
e
−
5
4.42886
e
−
12
≃
10
−
5
,
10
−
10
1.0
2.27723
e
−
5
4.08768
e
−
12
≃
10
−
5
,
10
−
10
4.33326
e
−
5
2.75555
e
−
10
≃
10
−
5
,
10
−
10Example 3.
Suppose in a closed ecosystem there are only 2 types of animals: the predator and the prey. They form a simple food-chain where the predator species hunts the prey species, while the prey grazes vegetation. The size of the 2 populations can be described by a simple system of2 nonlinear first-order differential equations:(38)p1′t=p1tA−Bp2t,p2′t=p2t−C+Dp1t,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 populationp1t increases at a rate dp1t/dt=Ap1t (proportional to the number of prey) but is simultaneously destroyed by predators at a rate dp1t/dt=−Bp1tp2t (proportional to the product of the number of prey and predators). A predator population p2t decreases at a rate dp2t/dt=−Cp2t (proportional to the number of predators), but increases at a rate dp2t/dt=Dp1tp2t (again proportional to the product of the numbers of prey and predators). By Theorem 1, we have the the positive equilibrium E∗=p1∗=C/D,p2∗=A/B.
For this example, the stability of each of the three steady states can be assessed more formally using the approach discussed. The absolute errors forN=7,9 are estimated by R1,N and R2,N and are presented in Tables 2 and 3. We see that if N increases, then the absolute errors decrease more rapidly.Table 2
Comparison of the absolute errors obtained by different rates forR1,N in Example 3.
t
A
=
0.45, B=0.15, C=0.45, D=0.9
A
=
0.8, B=0.4, C=0.4, D=0.6
A
=
0.16, B=0.48, C=0.3, D=0.15
N
=
7
N
=
9
N
=
7
N
=
9
N
=
7
N
=
9
0.2
3.75745
e
−
12
1.30547
e
−
14
2.34442
e
−
11
2.12354
e
−
13
1.65098
e
−
11
2.77903
e
−
15
0.4
1.00987
e
−
12
3.22244
e
−
14
1.02348
e
−
11
6.77926
e
−
14
1.30007
e
−
11
5.99001
e
−
13
0.6
6.43322
e
−
13
9.91158
e
−
15
4.82091
e
−
13
4.66258
e
−
15
8.20989
e
−
13
5.56193
e
−
15
0.8
4.49987
e
−
11
9.00302
e
−
14
3.35303
e
−
11
6.39937
e
−
15
8.96511
e
−
11
4.09824
e
−
15
1.0
2.66561
e
−
10
9.21367
e
−
14
6.12546
e
−
11
9.09012
e
−
15
5.36768
e
−
11
2.09891
e
−
15Table 3
Comparison of the absolute errors obtained by different rates forR2,N in Example 3.
t
A
=
0.45, B=0.15, C=0.45, D=0.9
A
=
0.8, B=0.4, C=0.4, D=0.6
A
=
0.16, B=0.48, C=0.3, D=0.15
N
=
7
N
=
9
N
=
7
N
=
9
N
=
7
N
=
9
0.2
3.86301
e
−
11
4.21008
e
−
15
8.23171
e
−
11
7.9878
e
−
15
5.76173
e
−
11
2.21378
e
−
15
0.4
2.45611
e
−
12
9.59823
e
−
15
8.06578
e
−
11
3.28935
e
−
15
7.20981
e
−
12
4.98678
e
−
15
0.6
1.87501
e
−
12
2.00256
e
−
15
1.63809
e
−
12
4.99876
e
−
15
2.65410
e
−
12
1.57792
e
−
15
0.8
2.12409
e
−
10
1.55729
e
−
14
8.28909
e
−
11
1.37811
e
−
15
2.00256
e
−
11
3.32209
e
−
14
1.0
9.12781
e
−
10
2.81009
e
−
14
3.12927
e
−
10
8.44606
e
−
15
1.34690
e
−
10
5.67813
e
−
14


## 6. 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. The 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.



---

*Source: 9176784-2020-04-01.xml*