Diffusive predator-prey model with Beddington-DeAngelis functional response: theoretical analysis and numerical results

1 Introduction

In natural ecosystems, the classical Lotka–Volterra model has been introduced to describe interactions between species [1, 2]. By integrating logistic growth for the prey population and the functional response of the predator population (which connects predation to predator–prey interactions), the Lotka–Volterra model has been extended into a prey–predator population model, where the functional response defines the feeding rate per predator on the prey population [1, 2]. A Holling-type functional response effectively captures the prey-dependent rate within a typical range but fails to account for interactions among predators and has been challenged by some biology and physiology communities [3, 4]. However, Beddington [5], Crowley and Martin [6], and Hassell and Varley [7] functional responses are of a predator-dependent type. A number of natural mechanisms can generate the Beddington–DeAngelis functional response. Beddington [5] and DeAngelis et al. [8] proposed that individuals in a population of more than two predators not only allocate time to search and handle their prey but also to encounter other predators. The ratio αuva+bu+cv represents the Beddington–DeAngelis functional response, while mua+bu+cv represents the predator's growth rate, where u and v are the prey and predator densities, respectively [5, 8].

The Beddington–DeAngelis functional response is similar to the Holling type II functional response but contains an extra term, cv, describing mutual interference by predators. This model has the potential to protect the predator from extinction, but under some conditions, it may also lead to the extinction of the prey [9]. The following model is considered from Song et al. [10].

∂u(x,t)∂t=d1∂2u∂x2+ru(1-uρ)-αuva+bu+cv,(1)

∂v(x,t)∂t=d2∂2v∂x2+v(-d+mua+bu+cv),(2)

with initial conditions

u(x,0)=u0(x), v(x,0)=v0(x),

where x ∈ Ω.

The boundary conditions are

∂u∂x(x,t)=0, ∂v∂x(x,t)=0,

where x ∈ ∂Ω, t > 0.

We note that u(x, t) and v(x, t) represent the prey and predator densities at the spatial location x and time t, respectively. The domain Ω is bounded with a smooth boundary ∂Ω. Additionally, a, b, c, d, m, and α are positive constants. The diffusion coefficients are denoted by d₁ and d₂, respectively. The parameter r represents the intrinsic growth rate, while ρ is the carrying capacity of the prey in the absence of predation.

There are three types of bifurcation [11]:

• (1) Turing
• (2) Hopf
• (3) Saddle-node

Turing and Hopf bifurcation are two important kinds of transitions giving birth to inhomogeneous solution. On a disk, these two bifurcations can lead to equivalent Turing-Hopf bifurcations [11, 12]. Zhang et al. [13] investigated a discrete prey-predator system with harvesting of both species and the Beddington–DeAngelis functional response.

Poincaré [14] laid the rigorous geometric groundwork for the qualitative study of differential equations and bifurcation theory. Later, Hopf [15] provided the mathematical framework for the birth of periodic orbits from an equilibrium, a phenomenon now known as the Hopf bifurcation. A bifurcation in dynamical systems occurs when a small, continuous change in a system's parameter causes a sudden qualitative change in its long-term behavior. These transitions define how systems lose stability, create new equilibrium points, or transition into cyclic and chaotic motion [12].

The main types of bifurcations are categorized according to how they alter system behavior [12]:

• Local bifurcation.
• Global bifurcation.
• Bifurcation in maps.

Local bifurcation occurs at a specific equilibrium point and is analyzed using Taylor series expansion. Some examples of local bifurcations are saddle–node, transcritical, pitchfork, and Hopf bifurcations.

A pitchfork bifurcation refers to a symmetry breaking event where a stable equilibrium loses its stability and gives rise to two stable equilibria that are symmetric [16]. In reaction-diffusion systems, a Turing instability is associated with an instability driven by diffusion that affects a spatially uniform steady state, resulting in the formation of stationary spatial patterns [11]. A Hopf bifurcation takes place when a set of complex conjugate eigenvalues crosses the imaginary axis, leading to the creation or annihilation of a periodic orbit [16].

Li et al. [17] combined theoretical analysis with numerical simulations to explore Turing instability and bifurcating periodic solutions. Yan and Zhang [18] focused on the local (global) asymptotic stability and Turing instability of the positive constant steady state, analyzing the impact of diffusion on stability. Tripathi et al. [19] studied the dynamical analysis of a prey predator model, where the predator's feeding rate followed the Beddington–DeAngelis functional response. This response, similar to Holling type II, included an extra term for predator interference. The study examined the effects of prey refuge and predator interference on system dynamics, identifying conditions for persistence. Stability analysis using Lyapunov functions showed that a Hopf bifurcation occurred when the prey refuge parameter crossed a critical value. The findings indicated that predator survival depended on this threshold, beyond which extinction occurred. Numerical simulations supported the results [19].

Xu and Fu [20] mainly studied the stability and Turing instability of the positive equilibrium. They also investigated the existence and nonexistence of nonconstant positive steady states.

Wen and Fu [21] demonstrated that strong cross-diffusion effects can destabilize uniform steady states and gradually induce Turing patterns driven by interspecies interaction pressures. Tan et al. [22] studied a Leslie-Gower predator-prey model with Beddington–DeAngelis functional response. They mainly studied stability, Turing instability and the Hopf bifurcation and showed the existence of Hopf bifurcation steady state solution and Turing-Hopf bifurcation via numerical results.

Song et al. [10] primarily examined the Turing instability and Hopf bifurcation of the coexisting equilibrium (u^*, v^*). Based on the model given by Equations 1, 2, they studied the Turing-Hopf bifurcation through bifurcation parameters m and d₂. Additionally, they obtained the normal form for the Turing-Hopf bifurcation to unfold the complex dynamics near the Turing-Hopf bifurcation point.

Farshid and Jalilian [23] investigated a cross-diffusion Leslie–Gower predator–prey model with the Beddington–DeAngelis functional response and showed that variations in cross-diffusion parameters can trigger both Turing instability and Hopf bifurcation, leading to rich spatial patterns. Along similar lines, Yan and Zhang [24] studied a reaction–diffusion predator–prey system under Neumann boundary conditions and derived stability criteria, conditions for Turing instability, and non-constant steady-state solutions, further showing that the Beddington–DeAngelis response increases spatiotemporal complexity. Together, these studies highlight the significant role of the Beddington–DeAngelis response and diffusion mechanisms in shaping complex dynamical behaviors in predator–prey systems.

Myint et al. [25] investigated the dynamical properties and positive steady states of a diffusive predator–prey system with fear effect and Beddington–DeAngelis functional response.

In this study, we perform a coexistence equilibrium analysis and derive the mathematical conditions under which a positive equilibrium exists for a diffusive predatory-prey system with Beddington–DeAngelis functional response. We then investigate the local stability of the equilibrium point and provide detailed conditions for the onset of Turing instability. To complement the theoretical findings, we construct the Forward Time Central Space (FTCS) numerical scheme and present some numerical profiles to validate the theoretical results.

There are many ways in which this work is different to that of Song et al. [10]. Firstly, the theoretical analysis is done in slightly different way. In Song et al. [10], numerical results are shown without explaining the numerical methods used and derivation of the methods is not given as well as the analysis of the methods. Thirdly, in our work, we construct Forward Time Central Space method. Derivations are provided. Stability and convergence are studied. The numerical rate of convergence is computed to verify the suitability of the method to solve the problem considered.

2 Equilibrium analysis

Letting α = βr and d = mδ, the system of equations given by Equations 1, 2 can be rewritten as

∂u(x,t)∂t=d1∂2u∂x2+ru(1-uρ-βva+bu+cv),(3)

∂v(x,t)∂t=d2∂2v∂x2+mv(-δ+ua+bu+cv),(4)

where x ∈ Ω and t ∈ (0, T]. A coexisting equilibrium point is an equilibrium point where all of its coordinates are positive [26]. We analyze the existence of a coexisting equilibrium. The system of equations without diffusion terms is given as follows:

ru(1-uρ-βva+bu+cv)=0,(5)

mv(-δ+ua+bu+cv)=0.(6)

The system of equations given by Equations 5, 6 has boundary equilibria (0, 0) and (ρ, 0).

Assuming (u^*, v^*) is the coexisting equilibrium, we discuss the existence of (u^*, v^*).

Using Equation 6, we get

ua+bu+cv=δ,

and on rearrangement, we obtain

u=δ(cv+a)1-bδ.(7)

For the solution u to be valid, we must have

δ(cv+a)1-bδ>0.

This implies that

bδ<1.(8)

Substituting Equation 7 into Equation 5 gives

rδ(cv+a)1-bδ[1-1ρδ(cv+a)(1-bδ)-βva+bδ(cv+a)1-bδ+cv]=0.

On expansion, we obtain

rδ(cv+a)ρ(1-bδ)-rδ2(cv+a)2ρ(1-bδ)2=rδβv.(9)

Expanding Equation 9, we get the following equation

-rδ2c2v2+(rδρc-rδ2ρbc+2rδ2βρb-2rδ2ca-rδβρ-rδ3βρb2)v+rδρa-rδ2ρab-rδ2a2=0.(10)

Assumming r ≠ 0 and δ ≠ 0 and dividing Equation 10 by −rδ, we get

c2δv2+v(δρbc-ρc-2δβρb+2δca+βρ+δ2βρb2)-ρa+δρab+δa2=0.

Let h(v) represent the quadratic expression given by Equation 11:

h(v)=c2δv2+v(2δca+cρ(-1+bδ)+ρβ(bδ-1)2)+aρ(-1+bδ)+a2δ.(11)

We solve h(v) = 0 to obtain Equations 12, 13:

v=-(2δca+cρ(bδ-1)+ρβ(bδ-1)2)+G2c2δ,(12)

and

v=-(2δca+cρ(bδ-1)+ρβ(bδ-1)2)-G2c2δ,(13)

where G is given by Equation 14:

G=(2δca+cρ(bδ-1)+ρβ(bδ-1)2)2 -4c2δ(a2δ+aρ(bδ-1)).(14)

Since v must be greater than 0, we reject the negative solution of v.

For v to exist the following conditions given by Equation 15 and Equation 16 must be satisfied:

(2δca+cρ(bδ-1)+ρβ(bδ-1)2)2-4c2δ(a2δ+aρ(bδ-1))>0,(15)

and

-(2δca+cρ(bδ-1)+ρβ(bδ-1)2)+G>0.(16)

Using Equation 16, we have inequalities Equations 17, 18:

(2δca+cρ(bδ-1)+ρβ(bδ-1)2)2-4c2δ(a2δ+aρ(bδ-1))>[2δa+cρ(bδ-1)+ρβ(bδ-1)2]2,(17)

which can be simplified as

-4c2δ(a2δ+aρ(bδ-1))>0.(18)

Now from Equation 18, we have that

δ<ρa+ρb.

Therefore, for v to exist, we need the inequality given by Equation 19:

δ<min(1b,ρa+ρb).(19)

3 Stability of equilibrium point

Using dynamical system theory, the stability of (u^*, v^*) is determined by the sign of the trace of the Jacobian matrix of (fg) [11, 27].

The equilibrium point (u^*, v^*) is locally asymptotically stable when the trace of Jacobi matrix J for (fg)is negative and determinant of Jacobian matrix is positive [27].

A positive determinant means that a linear transformation preserves the orientation of space. It indicates that the matrix is invertible and maps space in a right-handed manner.

We will examine the local asymptotic stability of the positive constant steady state of the system of Equations 3, 4 without diffusion, i.e.

0=ru(1-uρ-βva+bu+cv),0=mv(-δ+ua+bu+cv).

Let

f=ru-ru2ρ-βvrua+bu+cv,

and

g=mv(-δ+ua+bu+cv),

then the Jacobian matrix for (fg) is

J=(∂f∂u∂f∂v∂g∂u∂g∂v).

We therefore have

∂f∂u=r-2ruρ-(βvr(a+bu+cv)-βvrub(a+bu+cv)2)=r-2ruρ -(βvra+βv2rc)(a+bu+cv)2,(20)

∂f∂v=-βru(a+bu+cv)+βvruc(a+bu+cv)2=-βrua-βru2b(a+bu+cv)2,(21)

∂g∂u=mv(a+bu+cv)-mvub(a+bu+cv)2=mva+mv2c(a+bu+cv)2,(22)

∂g∂v=-mδ+mu(a+bu+cv)-mvuc(a+bu+cv)2=-mδ +mua+mu2b(a+bu+cv)2.(23)

The Jacobian matrix for (fg) at (u^*, v^*) is

J=(a11a12a21a22),(24)

where

a11=r-2ruρ-βvra+β(v)2rc(a+bu+cv*)2,(25)

a12=-βrua-βr(u)2b(a+bu+cv)2,(26)

a21=mva+m(v)2c(a+bu+cv)2,(27)

a22=-mδ+mua+m(u)2b(a+bu+cv)2.(28)

The determinant of the matrix J is det(J)=a11a22-a12a21, which gives

det(J)=-rmδ+2rumδρ +rmua+rm(u)2b+βvramδ+β(v)2rcmδ(a+bu+cv)2 -(2r(u)3mb+2r(u)2maρ(a+bu+cv*)2),(29)

and the trace of the matrix J is

tr(J)=a11+a22=r-2ruρ+mua+m(u)2b(a+bu+cv)2 -βrva+βr(v)2c(a+bu+cv*)2-mδ.

We conclude section 3 by highlighting the importance of Equations 20–29.

4 Turing instability of positive equilibrium point

We derive the conditions for Turing instability of the model given by the system of equations:

∂u∂t=d1∂2u∂x2+f(u,v),(30)

∂v∂t=d2∂2v∂x2+g(u,v).(31)

We use the approach used by Somathilake and Burrage [28] to obtain conditions for Turing instability for a space-fractional reaction-diffusion model for pattern formation in coral reefs.

Letting u = u₀+u₁ and v = v₀+v₁ and linearizing the system of equations about a steady state (u₀, v₀), using first order expansion, we get

∂u∂t=d1∂2u∂x2+f(u0,v0)+u1fu(u0,v0)+v1fv(u0,v0).(32)

Since (u₀, v₀) is a steady state, we get

∂∂t(u1)=d1∂2∂x2(u1)+u1fu(u0,v0)+v1fv(u0,v0).(33)

Similarly, we get

∂v∂t=d2∂2v∂x2+g(u0,v0)+u1gu(u0,v0)+v1gv(u0,v0).(34)

Also,

∂(v0+v1)∂t=d2∂2∂x2(v0+v1)+g(u0,v0)+u1gu(u0,v0) +v1gv(u0,v0),(35)

and since (u₀, v₀) is a steady state, this gives

∂v1∂t=d2∂2v1∂x2+u1gu(u0,v0)+v1gv(u0,v0).(36)

We can therefore express the system as

∂∂t(u1v1)=(d1Δ+a11a12a21d2Δ+a22)(u1v1),(37)

where

a11=fu(u0,v0), a12=fv(u0,v0), a21=gu(u0,v0),a22=gv(u0,v0).(38)

Using perturbations u1(x,t)=eσteIθx and v1(x,t)=eσteIθx [29], where I=-1, we obtain

∂u1∂t=σu1, ∂2u1∂x2=-θ2u1,∂v1∂t=σv1, ∂2v1∂x2=-θ2v1.

The perturbations yield the following set of equation in matrix form as

[σu1σv1]=[a11-d1θ2a12a21a22-d2θ2][u1v1].(39)

We therefore have

σu1=(a11-d1θ2)u1+a12v1,(40)

and

σv1=a21u1+(a22-d2θ2)v1.(41)

From Equation 40, we can express u₁ as

u1=a12v1σ-a11+d1θ2.(42)

Substituting Equation 42 into Equation 41, we get

(σ-a22+d2θ2)(σ-a11+d1θ2)=a21a12.

Letting μ = θ², we obtain

σ2-G(μ)σ+h(μ)=0,

where

G(μ)=a11-d1μ-d2μ+a22,

and

h(μ)=d1d2μ2-(d1a22+a11d2)μ+a22a11-a21a12.

Solving for σ gives

σ=G(μ)±[G(μ)]2-4h(μ)2.(43)

The trace of the matrix from Equation 39 is a₁₁−d₁μ+a₂₂−d₂μ and we need the trace to be negative for u₀, v₀ to be locally asymptotically stable. Hence, we need

a11+a22-μ(d1+d2)<0i.e. G(μ)<0(44)

For σ to take real values, we need [G(μ)]² − 4h(μ) > 0, i.e., [G(μ)]² > 4h(μ). Since G(μ) < 0, it follows that we require h(μ) < 0.

The system undergoes a Turing bifurcation if and only if h(μ) < 0.

We note that

h(μ)=d1d2μ2-(d1a22+a11d2)μ+a22a11-a21a12.

Differentiating h(μ) with respect to μ, we get

dhdμ=2d1d2μ-(d1a22+a11d2).

Solving dhdμ=0, we obtain

μc=d1a22+a11d22d1d2.(45)

h(μ) attains a minimum value which is

hmin=d1d2(d1a22+a11d22d1d2)2 -(d1a22+a11d2)(d1a22+a11d22d1d2)+a22a11-a21a12 =4d1d2(a22a11-a21a12)-(d1a22+a11d2)24d1d2.(46)

For spatial patterns to exist, we require μ_*c* > 0 and h_min < 0.

Thus, the two additional conditions for Turing instability are

d1a22+a11d2>0,(47)

4d1d2(a11a22-a21a12)-(d1a22+a11d2)24d1d2<0.(48)

The steps in section 4 are described by Equations 30–48.

5 Problem from Song et al.

Song et al. [10] solved the system of Equations 1, 2 subject to the initial conditions:

u(x,0)=2.0-0.1cos4x,v(x,0)=2.0+0.1cos4x,

The boundary conditions are

∂u∂x(x=0,t)=∂u∂x(x=4π,t)=0.

where x ∈ [0, 4π]. They use the following parameters: β = 0.8, ρ = 10, r = 5, m = 9, a = 0.2, b = 0.7, c = 0.2, d₁ = 1, d₂ = 9.Since 1b≈1.429,ρa+ρb≈1.389, they chose δ = 1, as δ<min(b,ρa+ρb). Using Equation 12, they obtain v^* = 2 and using Equation 7, u^* = 2. The values of a₁₁, a₁₂, a₂₁, and a₂₂ are 1.8, −3.2, 2.7, and −1.8, respectively.

The trace of the Jacobian matrix is a₁₁+a₂₂ = 0, which is one of the conditions needed for a Hopf bifurcation. The determinant is calculated as det(J)=5.4 and also we have

d1a22+d2a11=14.47>0,

4d1d2(a11a22-a21a12)-(d1a22+d2a11)24d1d2=-0.36<0.

6 Problem considered and choice of parameters

Section 6 is presented using Equations 49–54.

6.1 Conditions for Turing instability

The four conditions required for the occurrence of Turing instability are:

C1=det(J)=a11a22-a12a21>0,(49)

C2=tr(J)=a11+a22<0,(50)

C3=d1a22+d2a11>0,(51)

C4=4d1d2(a11a22-a12a21)-(d1a22+d2a11)24d1d2<0.(52)

Apart from the four conditions, a fifth condition is obtained by considering the initial conditions and the spatial domain of the problem.

On the interval [0, 4π], the eigenvalues of -∂2u∂x2 with homogeneous Neumann boundary conditions are

μN=(Nπ4π)2 =N216,

where N = 0, 1, 2, ...

Since

u(x,0)=2.0-0.1cos(4x),and v(x,0)=2.0+0.1cos(4x),

the wave number θ = 4.

Now, μ = θ² from Section 4. Therefore, we obtain μ = 4² = 16. Solving N216=16 gives N = 16.

Using Equation 45,

μ=d1a22+a11d22d1d2.

Therefore, d1a22+a11d22d1d2=16. The fifth condition for Turing instability is therefore

C5=d1a22+a11d2-32d1d2=0.(53)

We attempt in Section 6 to use some of the parameters in Song et al. [10].

6.2 Problem considered

We solve the same problem considered by Song et al. [10], but consider five different cases.

We solve Equations 3, 4 subject to the following initial conditions

u(x,0)=2.0-0.1cos(4x),v(x,0)=2.0+0.1cos(4x),

with x ∈ [0, 4π]

The boundary conditions are homogeneous Neumann boundary conditions, i.e., ∂u∂x(x,t)=0, and ∂v∂x(x,t)=0, for x ∈ ∂Ω, t ≠ 0.

For all five cases, the parameter values are fixed as

β=0.8,ρ=10, r=5, a=0.2, b=0.7, c=0.2.

The values of d₁, d₂, and m differ.

The parameter δ is calculated using Equation 18

δ<min(10.7,100.2+10×0.7)

that is,

δ<min(1.428571,1.388889).

We can choose δ = 1. Then the equilibrium values u^* and v^* are obtained as

u=2.000000, v=2.000000.

This means that for all five cases, we have

δ=1, u=2.000000, v=2.000000.

The coefficients are therefore

C1=0.6m,C2=1.8-0.2m,C3=1.8d2-0.2d1m,C4=14[2.4d1d2m-(1.8d2-0.2d1m)2d1d2].

The coefficient C₅ is obtained as

C5=-0.2d1m+1.8d2-32d1d2.(54)

6.3 Case 1 (Turing instability)

We need C₂ < 0, we solve

1.8-0.2m<0,

which gives

m>9.

We can choose m = 9.2 for example. Next solving C₅ = 0, gives

d1=22.5d223+400d2.

We next obtain 2D plots of C₃ and −C₄ vs. d₂ ∈ (0, 20) in Figure 1 to find a possible value of d₂. We need C₃, −C₄ > 0.

We can choose d₂ = 9 (similar value chosen in Song et al. [10]). After obtaining d₂ = 9, we calculate the corresponding value of d₁ as 22.5(9)23+400(9)=0.055893.

This case corresponds to Turing instability with the possibility of a Hopf bifurcation pattern, as C₂ = −0.04 (very close to zero).

6.4 Case 2 (Turing instability)

Solving C₂ = 1.8 − 0.2m < 0 gives m > 9. Hence, we can take m = 15.

Then, solving C₅ = 0 gives

d2=15d1-9+160d1.

To determine a suitable value of d₁, we plot the coefficients C₃ and −C₄ vs. d₁ ∈ (0, 0.055), as shown in Figure 2. We find a suitable value of d₁ to be 0.05.

We note that this 2D plot shows a discontinuity near d₁ = 0.06 for the curve −C₄. Therefore, we restrict the range of d₁ to (0, 0.055). Using d₁ = 0.05, the corresponding value of d₂ is 0.75.

6.5 Case 3 (Hopf bifurcation)

Here, we need C₂ = 0. Solving 1.8 − 0.2m = 0 gives m = 9. Then, solving C₅ = 0 gives

d1=9d29+160d2.

To determine a suitable value of d₂, we plot C₃ and −C₄ vs. d₂ ∈ (0, 20), as shown in Figure 3. From the figure, we choose d₂ = 9.

Substituting d₂ = 9 into the above expression gives the corresponding value of d₁ as d₁ = 0.055901.

6.6 Case 4

In case 4, we consider a situation where the first four conditions for Turing instability holds but the last condition is not satisfied, i.e.,

C1>0, C2<0, C3>0, C4<0, C5 ≠ 0.

Solving C₂ < 0 gives m > 9. Hence, we choose m = 15. We also choose d₂ = 9. The 2D plots of C₃ and −C₄ vs. d₁ ∈ (0, 20) are shown in Figure 4, and from the plots, we can choose d₁ = 0.2.

6.7 Case 5

In Case 5, we consider a situation where only three conditions for Turing instability are satisfied.

In this case,

C1>0, C3>0, C4<0,

while the remaining two conditions,

C2<0 and C5=0,

are not satisfied.

Solving C₂ > 0 gives m < 9. Hence, we choose m = 8. We also choose d₂ = 9, as in Case 4.

The plots of −C₃ and C₄ vs. d₁ are shown in Figure 5, from which we choose d₁ = 0.2.

We note from Table 1 that we expect Turing instability for Case 2. A Hopf bifurcation applies to Case 3. For Case 1, we expect Turing instability with a possible Hopf bifurcation. No Turing instability or Hopf bifurcation is expected for Cases 4 and 5.

Table 1. Values of parameters for the five cases.

7 Forward Time Central Space Scheme

7.1 Derivation

We construct the Forward Time Central Space Scheme (FTCS) scheme in order to discretise the system of equations given by Equations 3, 4.

We present the discretisation as follows:

Uin+1-Uink=d1(Ui+1n-2Uin+Ui-1nh2) +rUin(1-Uinρ-βVina+bUin+cVin)(55)

and

Vin+1-Vink=d2(Vi+1n-2Vin+Vi-1nh2) +mVin(-δ+Uina+bUin+cVin).(56)

Rewriting Equations 55, 56, we obtain

Uin+1=Uin+kd1h2(Ui+1n-2Uin+Ui-1n) +rkUin(1-Uinρ-βVina+bUin+cVin),(57)

Vin+1=Vin+kd2h2(Vi+1n-2Vin+Vi-1n) +mkVin(-δ+Uina+bUin+cVin).(58)

7.2 Consistency of the FTCS scheme

To check for consistency and to determine the order of accuracy of the FTCS scheme, we obtain the Taylor series of Equations 57, 58 about the point (t_*n*, x_*i*) as follows:

U+kUt+k2Utt2!+k3Uttt3!+k4Utttt4!+Θ(k5)=U +kd1h2(U+hUx+h2Uxx2!+h3Uxxx3! +h4Uxxxx4!+Θ(h5))-2Ukd1h2+kd1h2(U-hUx +h2Uxx2!-h3Uxxx3!+h4Uxxxx4!+Θ(h5)) +rUk(1-Uρ-βVa+bU+cV),(59)

and

V+kVt+k2Vtt2!+k3Vttt3!+k4Vtttt4!+Θ(k5)=V +kd2h2(V+hVx+h2Vxx2!+h3Vxxx3!+h4Vxxxx4!) -2Vkd2h2+kd2h2(V-hVx+h2Vxx2!-h3Vxxx3! +h4Vxxxx4!+Θ(h5)) +kmV(-δ+Ua+bU+cV).(60)

After some rearrangement of Equation 59, we obtain

k(Ut+kUtt2!+k2Uttt3!+k3Utttt4!+Θ(k4)) =kd1(2Uxx2!+2h2Uxxxx4!+Θ(h4)) +rUk(1-Uρ-βVa+bU+cV),(61)

which can be rewritten as

Ut+kUtt2!+k2Uttt3!+k3Utttt4!+Θ(k4) =d1(2Uxx2!+2h2Uxxxx4!+Θ(h4)) +rU(1-Uρ-βVa+bU+cV).(62)

Hence,

Ut-d1Uxx-rU(1-Uρ-βVa+bU+cV)=-k2Utt +d1h212Uxxxx+Θ(k2)+Θ(h4).(63)

Using Equation 60, we have

k(Vt+kVtt2!+k2Vttt3!+k3Vtttt4!+Θ(k4)) =kd2(2Vxx2!+2h2Vxxxx4!+Θ(h4)) +kmV(-δ+Ua+bU+CV).

Hence

Vt-d2Vxx-mV(-δ+Ua+bU+CV)=-kVtt2 +d2h2Vxxxx12+Θ(k2)+Θ(h4).(64)

As k, h → 0, Equations 63, 64 give Equations 65 and 66:

Ut-d1Uxx-rU(1-Uρ-βVa+bU+cV)=0,(65)

and

Vt-d2Vxx-mV(-δ+Ua+bU+cV)=0.(66)

We conclude that the FTCS scheme is consistent with the PDEs given by the system of Equations 3, 4. FTCS is first-order accurate in time and second-order accurate in space.

7.3 Stability

This is explained using Equations 67–70.

We consider Equation 57, i.e.,

Uin+1=Uin+kd1h2(Ui+1n-2Uin+Ui-1n) +rkUin(1-Uinρ-βVina+bUin+cVin).(67)

Using the freezing coefficient technique, we have

Uin+1=(1-2kd1h2)Uin+kd1h2(Ui+1n+Ui-1n) +rkUin-rkρUmaxUin-rβkVmaxUina+bUmax+cVmax,(68)

where

Umax=max|U(x,t)|

and

Vmax=max|V(x,t)|.

If we leave out the very complicated term, we get

Uin+1=(1-2kd1h2)Uin+kd1h2(Ui+1n+Ui-1n)+rkUin-rkUmaxρUin.(69)

Using Equation 69, we obtain the amplification factor as

ξ1=1+kd1h2(2cosω-2)+rk-rkUmaxρ.(70)

We observe that the amplification factor depends on d₁, r, ρ, k, h, U_max, and ω.

We need different values of d₁ and k for the five cases.

The range of values of k such that

with h=4π100, for the five cases is displayed in Table 2.

Table 2. Range of values of k satisfying |ξ₁| ≤ 1.000000 for the five cases.

We consider Equation 58, i.e.,

Vin+1=Vin+kd2h2(Vi+1n-2Vin+Vi-1n) +mkVin(-δ+Vina+bVi-1n+cVin).

Using freezing coefficient technique and Von Neumann analysis, we get the amplification factor

ξ2=1+2kd2h2(cosω-1)-δmk.

and we observe that it is dependent on d₂, k, h, ω, δ, and m.

We need different values of m, d₂ for the 5 cases.

The range of values of k such that

with h=4π100 for the five cases is displayed inTable 3.

Table 3. Range of values of k satisfying |ξ₂| ≤ 1.000000 for the five cases.

From Tables 2, 3, we now obtain the range of values of k when the FTCS scheme is stable for the five different cases when h=4π100. We next present some profiles for the 5 cases when FTCS scheme is used with k = 10⁻⁴ and h=4π100 in Figures 6–10.

Table 4. Range of values of k for which the FTCS scheme is stable for the five cases when h=4π100.

8 Numerical results for problem 1

Results using FTCS for case 1 (Figure 6).

Results using FTCS for case 2 (Figure 7).

Results using FTCS for case 3 (Figure 8).

Results using FTCS for case 4 (Figure 9).

Results using FTCS for case 5 (Figure 10).

9 Numerical rate of convergence

Since we do not possess a known exact solution, we determine the numerical rate of convergence over time by utilizing [30, 31].

RT=ln (EkEk/2)ln 2,

where E_*k* = ||U_*k*−U_2*k*||_∞ and E_*k*/2 = ||U_*k*/2−U_*k*||_∞ represent the discrete maximum norm errors.

To obtain the rate of convergence for FTCS, we start with k = 0.0001 and h=4π100 and the time step is reduced by a factor of two. The numerical rate of convergence in time is close to 1 as shown in Table 5.

Table 5. Numerical rate of convergence in time for the FTCS with h=4π100 at t = 1 for Case 1.

10 Conclusion

In this work, we perform a theoretical analysis of the predator-prey model with the Beddington-DeAngelis functional response. We determine the conditions for a unique coexistence equilibrium and obtain the conditions for Turing instability and Hopf bifurcation. We construct the Forward Time Central Space Scheme and study some of its properties. The Scheme is used to obtain plots of prey and predator density vs. time vs. space for five different cases describing Turing instabilities, Hopf bifurcation, and neither.

Generative AI statement

The author(s) declared that Generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

Edited by: Haitao Song, Shanxi University, China

Reviewed by: Xun Cao, Harbin Institute of Technology, China

Adisak Seesanea, Thammasat University, Thailand

Footnotes

k, temporal step size; h, spatial step size in the direction of x; ξ, amplification factor; t, time; u, the prey density; v, the predator density; Uin, numerical solution for the prey density at the grid point (t_*n*, x_*i*); Vin, numerical solution for the predator density at the grid point (t_*n*, x_*i*); ρ, carrying capacity for prey density.