Numerical analysis of Bazykin–Berezovskaya model

Abstract

In this manuscript, a Bazykin–Berezovskaya model with diffusion by strong Allee effects is studied. Neumann boundary conditions are used to see the positive solution of a diffusion system. Local stability analyses are discussed for all the equilibrium points. The analysis of stability for the proposed scheme is also given. Implicit finite difference schemes like: Euler, Crank–Nicolson (CN) and non-standard finite difference (NSFD) are used to verify the simulation by numerically. A comparison reveals that NSFD method is unconditionally stable for any temporal step-size.

Keywords:

• Bazykin–Berezovskaya (BB) model
• backward Euler method
• CN method
• equilibrium points
• NSFD method
• stability analysis

1. Introduction

In biomathematics, the Predator-Pray (PP) model is making an impact to discuss and understand the dynamics of populations. In 1925 and 1926, Alferd J. Lotka [1] and Vito Volterra [2] initially developed the model for a polynomial second-degree differential equation known as the Lotka-Volterra (LV) system. Lotka developed and studied the application of insectivore animal and plant family. On the other hand, the same model discussed by Volterra to explain the percentage of caught predictory fish. Volterra considered the system of linear ordinary differential equation by taking the PP model (1) $\begin{array}{r}{x}^{\prime }=x(\lambda -\mathrm{\mu y}),y^{\prime }=y(\xi -\mathrm{\zeta x}),\end{array}$(1) where $\lambda ,\mu ,\xi$ and ζ are positive numbers. Later on, the generalized form of LV system for the multiple dimensions was (2) $\begin{array}{r}{x}^{\prime }=x_l({\lambda }_{l0}+\sum _{k=1}^m{\lambda }_{\mathrm{lk}}{x}_k),l=0,1,2,\dots ,n.\end{array}$(2) Many researchers extended this model for a variety of applications especially in two dimensions, i.e. one predator and other one is pray species [3–7]. The theoretical aspect of PP model is very interesting for the development of population dynamics. In recent years, the authors focused on the phenomenon of the Allee effect. The correlation between the density (population size) and rate of growth is discussed on the Allee effect. Depending upon the nature of population growth, there are two types of the Allee effect, i.e. Strong Allee effect and Weak Allee effect. A weak Allee effect means, the population growth rate has reduced as per capita growth rate. Moreover, a population size rate becomes negative in case of a strong Allee effect under zero density. This complete procedure was developed by Allee in Ref. [8]. The Allee effect can be written as, (3) $\begin{array}{r}{x}^{\prime }=\delta (1-\frac{x}{p})(x-d)x,\end{array}$(3) where δ represents growth rate. The Allee effect depends on the value of d, if d>0, it behaves as a strong Allee effect and its becomes weak if $d\le 0$. In literature, there are different works that exist to analyse the outcomes of Allee effect on different population models and conclude that it can have more impact on the dynamical system. The stabilizing or destabilizing of a dynamical system by replacing the stability of some singularities in Allee effects, and the time taken to reach the stable solution of that dynamical system to take much longer.

In 1988, Bazykin and Berezovskaya introduced the new PP model which is known as Bazykin–Berezovskaya (BB) model [9]. This model related to the strong Allee effect. Voorn et al. [10] considered the following Bazykin-Berezovskaya model with strong Allee effects; (4) $\begin{array}{r}\frac{\mathrm{d}X}{\mathrm{d}t}=-X(L-X)(K-X)-\mathrm{XY},\frac{\mathrm{d}Y}{\mathrm{d}t}=-C(M-X)Y.\end{array}$(4) The model is in dimensionless and X, Y are the sizes of prey and predator population, respectively. The parameter K is the carrying capacity, L is the strong Allee effect threshold, C is the feeding efficiency, M is the predator mortality rate of Lotka-Voltrra model [10,11]. Bifurcation review in (4(4) $\begin{array}{r}\frac{\mathrm{d}X}{\mathrm{d}t}=-X(L-X)(K-X)-\mathrm{XY},\frac{\mathrm{d}Y}{\mathrm{d}t}=-C(M-X)Y.\end{array}$(4) ) with weak Allee effects discussed by Saad and Boubaker [12]. Further, the same authors studied the fractional order qualitative analysis. The noise-induced eradication in (4(4) $\begin{array}{r}\frac{\mathrm{d}X}{\mathrm{d}t}=-X(L-X)(K-X)-\mathrm{XY},\frac{\mathrm{d}Y}{\mathrm{d}t}=-C(M-X)Y.\end{array}$(4) ) investigated by Bashkirtseva and Ryashk [13]. In engineering, ecology, biology, and physics, the discrete-time population (DTP) models are applicable [14–21]. In addition, Din [22], Pal et al. [23], and Kartal et al. [24] investigated that (DTP) models may be more complex than continuous time models. Recently different authors published in planar systems [25–29].

The authors of Ref. [30] have used the fractional derivative of polynomials of the first kind in the Caputo sense to obtain a spectral solution for delay differential equations with the help of an explicit Chebyshev tau method. The authors of [31] has used the shifted Gegenbauer-Gauss collocation method for solving fractional differential equations with delay. Youssri et al. [32] has used the shifted Chebeyshev polynomials of the third kind in the Caputo sense to obtain fractional differential equation and fractional delay differential equations. In another article, Youssri et al. [33] has applied generalized Lucas polynomial sequence treatment to obtain the solution of fractional pantograph differential equation. The authors of the papers [34–42] considered various kinds of ordinary and functional differential equations of higher order and Volterra integro-differential equations. They investigated the stability and some other qualitative analysis of solutions of that equations using suitable tools and techniques.

2. Development of diffusion model

Over here, we consider 1-D coupled BB population model: (5) $\begin{array}{r}\{\begin{array}{rl}{\mathrm{\partial }}_{\tau }X(\tau ,x)& =-\mathrm{XY}+(L-X)X(X-K)+\mathrm{\Delta X}{\kappa }_1,\\ & x\in \mathfrak{U},\tau >0,\\ {\mathrm{\partial }}_{\tau }Y(\tau ,x)& =C(X-M)Y+{\kappa }_2\mathrm{\Delta Y}x\in \mathfrak{U},\tau >0,\\ X(0,x)& =X_0;Y(0,x)=Y_{0}x\in \mathfrak{U},\\ {\mathrm{\partial }}_{x}X(\tau ,x)& ={\mathrm{\partial }}_{x}Y(\tau ,x)=0,x\in \mathrm{\partial }\mathfrak{U},\tau >0,\end{array}\end{array}$(5) where $\mathfrak{U}$ is a bounded domain in ${\mathbb{R}}^N$ with sufficiently smooth boundary $\mathrm{\partial }\mathfrak{U}$. Here $X=X(x,\tau )$ and $Y=Y(x,\tau )$ symbolize the sizes of prey and predator respectively at time $\tau >0$, and a point $x\in \mathfrak{U}$. The parameters ${\kappa }_1>0$ and ${\kappa }_2>0$ are diffusion rates for each X and Y. Certainly due to biological reasons X and Y have to be non-negative functions. $X_0$ and $Y_0$ are also non-negative. We denote ${\mathrm{\partial }}_{\tau }$ as ordinary derivative with respect to τ.

2.1. Steady states of model without diffusion

The system (4(4) $\begin{array}{r}\frac{\mathrm{d}X}{\mathrm{d}t}=-X(L-X)(K-X)-\mathrm{XY},\frac{\mathrm{d}Y}{\mathrm{d}t}=-C(M-X)Y.\end{array}$(4) ) which is given in the absence of diffusion has four equilibrium points and these are $E_1=(0,0),E_2=(L,0),E_3=(K,0)$ and $E_4=(M,(M-L)(K-M))$.

2.2. Diffusionless case

In this case, we have ${\kappa }_1={\kappa }_2=0$.

1. For equilibrium point $E_1$, the Jacobian matrix is: (6) $\begin{array}{r}J=\left[\begin{array}{cc}-\mathrm{LK}& 0\\ 0& -\mathrm{CM}\end{array}\right].\end{array}$(6) This means the two eigenvalues are ${\lambda }_1=-\mathrm{LK}$ and ${\lambda }_2=-\mathrm{CM}$. Since all the parameters of the system are positive, so $0>{\lambda }_1$ and $0>{\lambda }_2$. Thus according to the RH criterion, if all the eigenvalues are negative or have negative real parts, then the model is locally asymptotically stable at that equilibrium point. So the model under consideration is locally asymptotically stable at $E_1$.

2. For equilibrium point $E_2$, the Jacobian matrix is: (7) $\begin{array}{r}J=\left[\begin{array}{cc}L(K-L)& -L\\ 0& C(L-M)\end{array}\right].\end{array}$(7) This means the two eigenvalues are ${\lambda }_1=-L(L-K)$ and ${\lambda }_2=-C(M-L)$. Since all the parameters of the system are positive, so ${\lambda }_1<0$ if K<L and ${\lambda }_2<0$ if L<M. Thus according to the RH criterion, if all the eigenvalues are negative or have negative real parts, then the model is locally asymptotically stable at that equilibrium point. So the model is locally asymptotically stable at $E_2$ if K<L<M.

3. For equilibrium point $E_3$, the Jacobian matrix is: (8) $\begin{array}{r}J=\left[\begin{array}{cc}-K(K-L)& -K\\ 0& -C(M-K)\end{array}\right].\end{array}$(8) This means the two eigenvalues are ${\lambda }_1=-K(K-L)$ and ${\lambda }_2=-C(M-K)$. Since the all the parameters of the system are positive, so ${\lambda }_1<0$ if L<K and ${\lambda }_2<0$ if K<M. Thus according to the RH criterion, if all the eigenvalues are negative or have negative real parts, then the model is locally asymptotically stable at that equilibrium point. So the model is locally asymptotically stable at $E_3$ if L<K<M.

4. For equilibrium point $E_4$, the Jacobian matrix is: (9) $\begin{array}{r}J=\left[\begin{array}{cc}M(K+L-2M)& -M\\ (\mathrm{CK}-\mathrm{CM})(M-L)& 0\end{array}\right].\end{array}$(9) The characteristic equation of it is (10) $\begin{array}{r}{\lambda }^2+M(2M-K-L)\lambda +(C{M}^2-\mathrm{CLM})(K-M)=0.\end{array}$(10) It can be written as (11) $\begin{array}{r}{\lambda }^2+a_1\lambda +a_2=0,\end{array}$(11) where $a_1=M(2M-K-L)$ and $a_2=\mathrm{MC}(M-L)(K-M)$. Then according to the RH criterion if $a_1>0$ and $a_2>0$, then the system is locally asymptotically stable at that equilibrium point. It could be possible only if 2M−K−L>0 or $M>(K+L)/2$, M>L and K>M. In other words, it can be written as $K>M>(K+L)/2$ and M>L otherwise this model is unstable.

2.3. Diffusion model with steady-states

We have to discuss the basic properties of the Bazykin–Berezovskaya population system of non-homogeneous steady states (12) $\begin{array}{r}\begin{array}{rl}0& ={\kappa }_1\mathrm{\Delta }X+X(X-L)(K-X)-\mathrm{XY},\\ 0& ={\kappa }_2\mathrm{\Delta Y}+C(X-M)Y,\end{array}\end{array}$(12) with the suitable conditions $X_x(x,\tau )=Y_x(x,\tau )=0,x\in \mathrm{\partial }\mathfrak{U}$.

Definition 2.1

If a constant solution is stable in the absence of diffusion but becomes unstable when diffusion is present, the constant solution is said to be turing unstable.

Proposition 2.1

[43]

Suppose that $h\in C(\overline{\mathfrak{U}}\times \mathbb{R})$ and $U\in C^2(\mathfrak{U})\cap C^1(\overline{\mathfrak{U}})$. It follows that

(a)	If (13) $\begin{array}{r}\mathrm{\Delta }U(Z)+h(Z,U(Z))\ge 0,\end{array}$(13) with $U_Z\le 0$ on $\mathrm{\partial }\mathfrak{U}$ and $U(Z_0)=\mathrm{ma}{x}_{\mathfrak{U}}(U(Z))$, then (14) $\begin{array}{r}g(Z_0,U(Z_0))\ge 0.\end{array}$(14)
(b)	If (15) $\begin{array}{r}\mathrm{\Delta }U(Z)+h(Z,U(Z))\le 0,\end{array}$(15) with $U_Z\ge 0$ on $\mathrm{\partial }\mathfrak{U}$ and $U(Z_0)=\mathrm{mi}{n}_{\mathfrak{U}}(U(Z))$, then (16) $\begin{array}{r}h(Z_0,U(Z_0))\ge 0.\end{array}$(16)

3. Numerical investigation

We used the following numerical analysis based on backward Euler (BE) method, CN method and NSFD method to investigate the solutions and their stability.

3.1. The BB model

In this section, the system (5(5) $\begin{array}{r}\{\begin{array}{rl}{\mathrm{\partial }}_{\tau }X(\tau ,x)& =-\mathrm{XY}+(L-X)X(X-K)+\mathrm{\Delta X}{\kappa }_1,\\ & x\in \mathfrak{U},\tau >0,\\ {\mathrm{\partial }}_{\tau }Y(\tau ,x)& =C(X-M)Y+{\kappa }_2\mathrm{\Delta Y}x\in \mathfrak{U},\tau >0,\\ X(0,x)& =X_0;Y(0,x)=Y_{0}x\in \mathfrak{U},\\ {\mathrm{\partial }}_{x}X(\tau ,x)& ={\mathrm{\partial }}_{x}Y(\tau ,x)=0,x\in \mathrm{\partial }\mathfrak{U},\tau >0,\end{array}\end{array}$(5) ) will discretized: (17) $\begin{array}{r}\{\begin{array}{rl}{X}_{\tau }& =X(X-L)(K-X)-\mathrm{XY}+{\kappa }_1\mathrm{\Delta X},\\ Y_{\tau }& =C(X-M)Y+{\kappa }_2\mathrm{\Delta Y}.\end{array}\end{array}$(17) Using FD for Equation (17(17) $\begin{array}{r}\{\begin{array}{rl}{X}_{\tau }& =X(X-L)(K-X)-\mathrm{XY}+{\kappa }_1\mathrm{\Delta X},\\ Y_{\tau }& =C(X-M)Y+{\kappa }_2\mathrm{\Delta Y}.\end{array}\end{array}$(17) ), firstly $[0,P{]}^2\times [0,T]$ are shared into $A^2\times B$ with space and time step sizes like, $\mathrm{d}x=P/A$ and $\mathrm{d}\tau =T/B$. Then the mesh nodes are $x_p=p\mathrm{d}x$ and ${\tau }_q=q\mathrm{d}\tau$, where $p=0,1,2,\dots ,A$ and $q=0,1,2,\dots ,B$.

We can write $X_p^q$ and $Y_p^q$ as the FD approximations of $X(p\mathrm{d}x,q\mathrm{d}\tau )$ and $Y(p\mathrm{d}x,q\mathrm{d}\tau )$. First- and second-order temporal and spatial derivative difference formulas are: (18) $\begin{array}{rl}(C{)}_{\tau }{|}_u^v& ={\delta }_{\tau }(C{)}_u^v=\frac{(C{)}_u^{v+1}-(C{)}_u^v}{\mathrm{d}\tau }+O(\mathrm{d}\tau ),\end{array}$(18) (19) $\begin{array}{rl}(C{)}_{\mathrm{xx}}{|}_u^v& ={\delta }_x^2(C{)}_u^v=\frac{(C{)}_{u-1}^v-2(C{)}_u^v+(C{)}_{u+1}^v}{(\mathrm{d}x{)}^2}\\ & +O((\mathrm{d}x{)}^2),\end{array}$(19) (20) $\begin{array}{rl}(C{)}_{\mathrm{xx}}{|}_u^{v+1}& ={\delta }_x^2(C{)}_u^{v+1}=\frac{(C{)}_{u-1}^{v+1}-2(C{)}_u^{v+1}+(C{)}_{u+1}^{v+1}}{(\mathrm{d}x{)}^2}\\ & +O((\mathrm{d}x{)}^2).\end{array}$(20)

3.2. BE method

Applying the values of $X_{\tau }{|}_p^q$ and $X_{\mathrm{xx}}{|}_p^q$ in the first equation of (17(17) $\begin{array}{r}\{\begin{array}{rl}{X}_{\tau }& =X(X-L)(K-X)-\mathrm{XY}+{\kappa }_1\mathrm{\Delta X},\\ Y_{\tau }& =C(X-M)Y+{\kappa }_2\mathrm{\Delta Y}.\end{array}\end{array}$(17) ) and the values of $Y_{\tau }{|}_p^q$ and $Y_{\mathrm{xx}}{|}_p^q$ in the second equation of (17(17) $\begin{array}{r}\{\begin{array}{rl}{X}_{\tau }& =X(X-L)(K-X)-\mathrm{XY}+{\kappa }_1\mathrm{\Delta X},\\ Y_{\tau }& =C(X-M)Y+{\kappa }_2\mathrm{\Delta Y}.\end{array}\end{array}$(17) ), we have (21) $\begin{array}{rl}& {\delta }_{\tau }{X}_p^q={\kappa }_1{\delta }_x^2{X}_p^{q+1}+X_p^q(X_p^q-L)(K-X_p^q)-X_p^q{Y}_p^q,\end{array}$(21) (22) $\begin{array}{rl}& {\delta }_{\tau }{Y}_p^q={\kappa }_2{\delta }_x^2{Y}_p^{q+1}+C(X_p^q-M)Y_p^q.\end{array}$(22) After some simple calculation, we have (23) $\begin{array}{rl}& -{\mathcal{R}}_1(X_{p-1}^{q+1}+X_{p+1}^{q+1})+(1+2{\mathcal{R}}_1)X_p^{q+1}\\ & =X_p^q+\mathrm{d}\tau (X_p^q(X_p^q-L)(K-X_p^q)-X_p^q{Y}_p^q)\end{array}$(23) (24) $\begin{array}{rl}& -{\mathcal{R}}_2(Y_{p-1}^{q+1}+Y_{p+1}^{q+1})+(1+2{\mathcal{R}}_2)Y_p^{q+1}\\ & =Y_p^q+\mathrm{d}\tau (C(X_p^q-M)Y_p^q),\end{array}$(24) where ${\mathcal{R}}_1=\frac{{\kappa }_1\mathrm{d}\tau }{(\mathrm{d}x{)}^2}$ and ${\mathcal{R}}_2=\frac{{\kappa }_2\mathrm{d}\tau }{(\mathrm{d}x{)}^2}$.

This technique is stable for Equations (23(23) $\begin{array}{rl}& -{\mathcal{R}}_1(X_{p-1}^{q+1}+X_{p+1}^{q+1})+(1+2{\mathcal{R}}_1)X_p^{q+1}\\ & =X_p^q+\mathrm{d}\tau (X_p^q(X_p^q-L)(K-X_p^q)-X_p^q{Y}_p^q)\end{array}$(23) ) and (24(24) $\begin{array}{rl}& -{\mathcal{R}}_2(Y_{p-1}^{q+1}+Y_{p+1}^{q+1})+(1+2{\mathcal{R}}_2)Y_p^{q+1}\\ & =Y_p^q+\mathrm{d}\tau (C(X_p^q-M)Y_p^q),\end{array}$(24) ).

Figure 1. The initial distribution of the numbers of Prey and Predator for all cases of the system.

3.3. CN method

Applying the values of $X_{\tau }{|}_p^q$, $X_{\mathrm{xx}}{|}_p^q$ and $X_{\mathrm{xx}}{|}_p^{q+1}$ in the first equation of (4(4) $\begin{array}{r}\frac{\mathrm{d}X}{\mathrm{d}t}=-X(L-X)(K-X)-\mathrm{XY},\frac{\mathrm{d}Y}{\mathrm{d}t}=-C(M-X)Y.\end{array}$(4) ) and the values of $Y_{\tau }{|}_p^q$, $Y_{\mathrm{xx}}{|}_p^q$ and $Y_{\mathrm{xx}}{|}_p^{q+1}$ in Equation (5(5) $\begin{array}{r}\{\begin{array}{rl}{\mathrm{\partial }}_{\tau }X(\tau ,x)& =-\mathrm{XY}+(L-X)X(X-K)+\mathrm{\Delta X}{\kappa }_1,\\ & x\in \mathfrak{U},\tau >0,\\ {\mathrm{\partial }}_{\tau }Y(\tau ,x)& =C(X-M)Y+{\kappa }_2\mathrm{\Delta Y}x\in \mathfrak{U},\tau >0,\\ X(0,x)& =X_0;Y(0,x)=Y_{0}x\in \mathfrak{U},\\ {\mathrm{\partial }}_{x}X(\tau ,x)& ={\mathrm{\partial }}_{x}Y(\tau ,x)=0,x\in \mathrm{\partial }\mathfrak{U},\tau >0,\end{array}\end{array}$(5) ), we have (25) $\begin{array}{rl}{\delta }_{\tau }{X}_p^q& =\frac{{\kappa }_1}{2}({\delta }_x^2{X}_p^{q+1}+{\delta }_x^2{X}_p^q)+X_p^q(X_p^q-L)(K-X_p^q)\\ & -X_p^q{Y}_p^q,\end{array}$(25) (26) $\begin{array}{rl}{\delta }_{\tau }{Y}_p^q& =\frac{{\kappa }_2}{2}({\delta }_x^2{Y}_p^{q+1}+{\delta }_x^2{Y}_p^q)+C(X_p^q-M)Y_p^q.\end{array}$(26) After some simple calculation, we have (27) $\begin{array}{rl}& -\frac{{\mathcal{R}}_1}{2}(X_{p-1}^{q+1}+X_{p+1}^{q+1})+(1+{\mathcal{R}}_1)X_p^{q+1}\\ & =\frac{{\mathcal{R}}_1}{2}(X_{p-1}^q+X_{p+1}^q)+(1-{\mathcal{R}}_1)X_p^{q+1}\\ & +\mathrm{d}\tau (X_p^q(X_p^q-L)(K-X_p^q)-X_p^q{Y}_p^q),\end{array}$(27) (28) $\begin{array}{rl}& -{\mathcal{R}}_2(Y_{p-1}^{q+1}+Y_{p+1}^{q+1})+2(1+{\mathcal{R}}_2)Y_p^{q+1}\\ & ={\mathcal{R}}_2(Y_{p-1}^q+Y_{p+1}^q)+2(1-{\mathcal{R}}_2)Y_p^{q+1}\\ & +2\mathrm{d}\mathrm{\tau C}(X_p^q-M)Y_p^q,\end{array}$(28) where ${\mathcal{R}}_1=\frac{{\kappa }_1\mathrm{d}\tau }{(\mathrm{d}x{)}^2}$ and ${\mathcal{R}}_2=\frac{{\kappa }_2\mathrm{d}\tau }{(\mathrm{d}x{)}^2}$.

This technique is unconditionally stable for Equations (17(17) $\begin{array}{r}\{\begin{array}{rl}{X}_{\tau }& =X(X-L)(K-X)-\mathrm{XY}+{\kappa }_1\mathrm{\Delta X},\\ Y_{\tau }& =C(X-M)Y+{\kappa }_2\mathrm{\Delta Y}.\end{array}\end{array}$(17) ).

3.4. Proposed NSFD method

We will design the NSFD method for system (17(17) $\begin{array}{r}\{\begin{array}{rl}{X}_{\tau }& =X(X-L)(K-X)-\mathrm{XY}+{\kappa }_1\mathrm{\Delta X},\\ Y_{\tau }& =C(X-M)Y+{\kappa }_2\mathrm{\Delta Y}.\end{array}\end{array}$(17) ). Applying the values of $X_{\tau }{|}_p^q$, $X_{\mathrm{xx}}{|}_p^q$, $Y_{\tau }{|}_p^q$, $Y_{\mathrm{xx}}{|}_p^q$ in Equation (4(4) $\begin{array}{r}\frac{\mathrm{d}X}{\mathrm{d}t}=-X(L-X)(K-X)-\mathrm{XY},\frac{\mathrm{d}Y}{\mathrm{d}t}=-C(M-X)Y.\end{array}$(4) ), we have (29) $\begin{array}{rl}\frac{X_p^{q+1}-X_p^q}{\mathrm{d}\tau }& ={\kappa }_1\left(\frac{X_{p-1}^q-2X_p^{q+1}+X_{p+1}^q}{(\mathrm{d}x{)}^2}\right)\\ & +(K+L)(X_p^q{)}^2-X_p^{q+1}(\mathrm{KL}+(X_p^q{)}^2)\\ & -X_p^{q+1}{Y}_p^q\end{array}$(29) (30) $\begin{array}{rl}\frac{Y_p^{q+1}-Y_p^q}{\mathrm{d}\tau }& ={\kappa }_2\left(\frac{Y_{p-1}^q-2Y_p^{q+1}+Y_{p+1}^q}{(\mathrm{d}x{)}^2}\right)\\ & +C(X_p^q{Y}_p^q-M{Y}_p^{q+1})\end{array}$(30) or (31) $\begin{array}{rl}{X}_p^{q+1}& =X_p^q+{\mathcal{R}}_1(X_{p-1}^q-2X_p^{q+1}+X_{p+1}^q)\\ & +\mathrm{d}\tau (K+L)(X_p^q{)}^2-\mathrm{d}\tau X_p^{q+1}(\mathrm{KL}+(X_p^q{)}^2)\\ & -\mathrm{d}\tau X_p^{q+1}{Y}_p^q,\end{array}$(31) (32) $\begin{array}{rl}{Y}_p^{q+1}& =Y_p^q+{\mathcal{R}}_2(Y_{p-1}^q-2Y_p^{q+1}+Y_{p+1}^q)\\ & +\mathrm{d}\mathrm{\tau C}(X_p^q{Y}_p^q-M{Y}_p^{q+1}).\end{array}$(32) After some calculation, we have the following explicit NSFD scheme (33) $\begin{array}{r}\{\begin{array}{rl}{X}_p^{q+1}& =\frac{X_p^q+{\mathcal{R}}_1(X_{p-1}^q+X_{p+1}^q)+\mathrm{d}\tau (K+L)(X_p^q{)}^2}{(1+2{\mathcal{R}}_1+\mathrm{d}\tau (\mathrm{KL}+(X_p^q{)}^2)+\mathrm{d}\tau Y_p^q)},\\ Y_p^{q+1}& =\frac{Y_p^q+{\mathcal{R}}_2(Y_{p-1}^q+Y_{p+1}^q)+\mathrm{d}\mathrm{\tau C}(X_p^q{Y}_p^q)}{(1+2{\mathcal{R}}_2+\mathrm{d}\mathrm{\tau CM})},\end{array}\end{array}$(33) where ${\mathcal{R}}_1=\frac{{\kappa }_1\mathrm{d}\tau }{(\mathrm{d}x{)}^2}$ and ${\mathcal{R}}_2=\frac{{\kappa }_2\mathrm{d}\tau }{(\mathrm{d}x{)}^2}$.

This approach is completely stable for Equations (17(17) $\begin{array}{r}\{\begin{array}{rl}{X}_{\tau }& =X(X-L)(K-X)-\mathrm{XY}+{\kappa }_1\mathrm{\Delta X},\\ Y_{\tau }& =C(X-M)Y+{\kappa }_2\mathrm{\Delta Y}.\end{array}\end{array}$(17) ).

3.4.1. Analysis of stability NSFD scheme

To obtain the steadiness norms, Von-Neumann steadiness scheme is used here. After linearizing (31(31) $\begin{array}{rl}{X}_p^{q+1}& =X_p^q+{\mathcal{R}}_1(X_{p-1}^q-2X_p^{q+1}+X_{p+1}^q)\\ & +\mathrm{d}\tau (K+L)(X_p^q{)}^2-\mathrm{d}\tau X_p^{q+1}(\mathrm{KL}+(X_p^q{)}^2)\\ & -\mathrm{d}\tau X_p^{q+1}{Y}_p^q,\end{array}$(31) ) and then applying the Von-Neumann steadiness norms, we get (34) $\begin{array}{rl}& Z(\mathrm{\Delta t}+t){\mathrm{e}}^{\mathrm{\iota \beta x}}\\ & =Z(t){\mathrm{e}}^{\mathrm{\iota \beta x}}+{\mathcal{R}}_1(Z(t){\mathrm{e}}^{\mathrm{\iota \beta }(x-\mathrm{\Delta x})}-2Z(t+\mathrm{\Delta t}){\mathrm{e}}^{\mathrm{\iota \beta x}}\\ & +Z(t){\mathrm{e}}^{\mathrm{\iota \beta }(x+\mathrm{\Delta x})})-\mathrm{d}\mathrm{\tau KLZ}(t+\mathrm{\Delta t}){\mathrm{e}}^{\mathrm{\iota \beta x}},\end{array}$(34) (35) $\begin{array}{rl}& (1+2{\mathcal{R}}_1+\mathrm{KL}\mathrm{d}\tau )Z(t+\mathrm{\Delta t}){\mathrm{e}}^{\mathrm{\iota \beta x}}\\ & =(1+2{\mathcal{R}}_1\cos (\beta \mathrm{\Delta x}))Z(t){\mathrm{e}}^{\mathrm{\iota \beta x}},\end{array}$(35) (36) $\begin{array}{rl}& \frac{Z(t+\mathrm{\Delta t})}{Z(t)}=\frac{(1+2{\mathcal{R}}_1\cos (\beta \mathrm{\Delta x}))}{1+2{\mathcal{R}}_1+\mathrm{LK}\mathrm{d}\tau }\\ & =\frac{(1+2{\mathcal{R}}_1-4{\mathcal{R}}_1{\sin }^2(\frac{\beta \mathrm{\Delta x}}{2}))}{1+2{\mathcal{R}}_1\mathrm{KL}\mathrm{d}\tau },\end{array}$(36) (37) $\begin{array}{rl}& \left|\frac{Z(t+\mathrm{\Delta t})}{Z(t)}\right|=\left|\frac{1+2{\mathcal{R}}_1-4{\mathcal{R}}_1{\sin }^2(\frac{\beta \mathrm{\Delta x}}{2})}{1+2{\mathcal{R}}_1\mathrm{KL}\mathrm{d}\tau }\right|\\ & \le \left|\frac{1-2{\mathcal{R}}_1}{1+2{\mathcal{R}}_1+\mathrm{KL}\mathrm{d}\tau }\right|<1.\end{array}$(37) Similarly, (38) $\begin{array}{rl}\left|\frac{T(t+\mathrm{\Delta t})}{T(t)}\right|& =\left|\frac{1+2{\mathcal{R}}_1-4{\mathcal{R}}_1{\sin }^2(\frac{\beta \mathrm{\Delta x}}{2})}{1+2{\mathcal{R}}_2}\right|\\ & \le \left|\frac{-2{\mathcal{R}}_2+1}{2{\mathcal{R}}_1+1+\mathrm{CM}\mathrm{d}\tau }\right|<1.\end{array}$(38) From (37(37) $\begin{array}{rl}& \left|\frac{Z(t+\mathrm{\Delta t})}{Z(t)}\right|=\left|\frac{1+2{\mathcal{R}}_1-4{\mathcal{R}}_1{\sin }^2(\frac{\beta \mathrm{\Delta x}}{2})}{1+2{\mathcal{R}}_1\mathrm{KL}\mathrm{d}\tau }\right|\\ & \le \left|\frac{1-2{\mathcal{R}}_1}{1+2{\mathcal{R}}_1+\mathrm{KL}\mathrm{d}\tau }\right|<1.\end{array}$(37) ) and (38(38) $\begin{array}{rl}\left|\frac{T(t+\mathrm{\Delta t})}{T(t)}\right|& =\left|\frac{1+2{\mathcal{R}}_1-4{\mathcal{R}}_1{\sin }^2(\frac{\beta \mathrm{\Delta x}}{2})}{1+2{\mathcal{R}}_2}\right|\\ & \le \left|\frac{-2{\mathcal{R}}_2+1}{2{\mathcal{R}}_1+1+\mathrm{CM}\mathrm{d}\tau }\right|<1.\end{array}$(38) ), it is clear that the proposed nonstandard FD scheme for system (17(17) $\begin{array}{r}\{\begin{array}{rl}{X}_{\tau }& =X(X-L)(K-X)-\mathrm{XY}+{\kappa }_1\mathrm{\Delta X},\\ Y_{\tau }& =C(X-M)Y+{\kappa }_2\mathrm{\Delta Y}.\end{array}\end{array}$(17) ) is completely stable.

3.4.2. Consistency of NSFD method

To observe the uniformity of NSFD technique, Taylor series is used for $X{|}_p^{q+1}$, $X{|}_{p+1}^q$ and $X{|}_{p-1}^q$ are given below: (39) $\begin{array}{rl}X{|}_p^{q+1}-X{|}_p^q& =(\mathrm{d}\tau )X_{\tau }{|}_p^q+\frac{(\mathrm{d}\tau {)}^2}{2}{X}_{\mathrm{\tau \tau }}{|}_p^q\\ & +\frac{(\mathrm{d}\tau {)}^3}{6}{X}_{\mathrm{\tau \tau \tau }}{|}_p^q+\dots .\end{array}$(39) (40) $\begin{array}{rl}X{|}_{p+1}^q-X{|}_p^q& =(\mathrm{d}x)X_x{|}_p^q+\frac{(\mathrm{d}x{)}^2}{2}{X}_{\mathrm{xx}}{|}_p^q\\ & +\frac{(\mathrm{d}x{)}^3}{6}{X}_{\mathrm{xxx}}{|}_p^q+\dots .\end{array}$(40) (41) $\begin{array}{rl}X{|}_{p-1}^q-X{|}_p^q& =-(\mathrm{d}x)X_x{|}_p^q+\frac{(\mathrm{d}x{)}^2}{2}{X}_{\mathrm{xx}}{|}_p^q\\ & -\frac{(\mathrm{d}x{)}^3}{6}{X}_{\mathrm{xxx}}{|}_p^q+\dots .\end{array}$(41) (42) $\begin{array}{rl}X{|}_p^{q+1}& =X{|}_p^q+{\mathcal{R}}_1(X{|}_{p-1}^q-2X{|}_p^{q+1}+X{|}_{p+1}^q)\\ & +\mathrm{d}\tau (K+L)(X{|}_p^q{)}^2-\mathrm{d}\mathrm{\tau X}{|}_p^{q+1}\\ & \times (\mathrm{KL}+(X{|}_p^q{)}^2)-\mathrm{d}\tau (X{|}_p^{q+1})(Y{|}_p^q).\end{array}$(42) Putting the values of $X{|}_p^{q+1}$, $X{|}_{p+1}^q$ and $X{|}_{p-1}^q$ in above equation, we have (43) $\begin{array}{rl}& (X_{\tau }{|}_p^q+\frac{(\mathrm{d}\tau )}{2}{X}_{\mathrm{\tau \tau }}{|}_p^q+\frac{(\mathrm{d}\tau {)}^2}{6}{X}_{\mathrm{\tau \tau \tau }}{|}_p^q+\dots ..)\\ & \times (1+2\frac{{\kappa }_1\mathrm{d}\tau }{(\mathrm{d}x{)}^2}+\mathrm{d}\tau (\mathrm{KL}+(X_p^q{)}^2)+Y_p^q)\\ & ={\kappa }_1(X_{\mathrm{xx}}{|}_p^q+\frac{(\mathrm{d}x{)}^2}{12}{X}_{\mathrm{xxxx}}{|}_p^q+\dots )+(K+L)\\ & \times (X{|}_p^q{)}^2-\mathrm{d}\tau (\mathrm{KL}+(X{|}_p^q{)}^2)-\mathrm{d}\tau (X{|}_p^q)(Y{|}_p^q).\end{array}$(43) By substituting the values $\mathrm{d}\tau =(\mathrm{d}x{)}^3$ and then taking $\mathrm{d}x\to 0$, the above equation gives (44) $\begin{array}{r}{X}_{\tau }={\kappa }_1\mathrm{\Delta }X+X(X-L)(K-X)-\mathrm{XY}.\end{array}$(44) Similarly, Taylor series is used for $Y{|}_p^{q+1}$, $Y{|}_{p+1}^q$ and $Y{|}_{p-1}^q$ are given below: (45) $\begin{array}{rl}Y{|}_p^{q+1}-Y{|}_p^q& =(\mathrm{d}\tau )Y_{\tau }{|}_p^q+\frac{(\mathrm{d}\tau {)}^2}{2}{Y}_{\mathrm{\tau \tau }}{|}_p^q\\ & +\frac{(\mathrm{d}\tau {)}^3}{6}{Y}_{\mathrm{\tau \tau \tau }}{|}_p^q+\dots ,\end{array}$(45) (46) $\begin{array}{rl}Y{|}_{p+1}^q-Y{|}_p^q& =(\mathrm{d}x)Y_x{|}_p^q+\frac{(\mathrm{d}x{)}^2}{2}{Y}_{\mathrm{xx}}{|}_p^q\\ & +\frac{(\mathrm{d}x{)}^3}{6}{Y}_{\mathrm{xxx}}{|}_p^q+\dots .,\end{array}$(46) (47) $\begin{array}{rl}Y{|}_{p-1}^q-Y{|}_p^q& =-(\mathrm{d}x)Y_x{|}_p^q+\frac{(\mathrm{d}x{)}^2}{2}{Y}_{\mathrm{xx}}{|}_p^q\\ & -\frac{(\mathrm{d}x{)}^3}{6}{Y}_{\mathrm{xxx}}{|}_p^q+\dots .\end{array}$(47) Putting the values of $Y{|}_p^{q+1}$, $Y{|}_{p+1}^q$ and $Y{|}_{p-1}^q$ in above equation, we have (48) $\begin{array}{rl}& (Y_{\tau }{|}_p^q+\frac{(\mathrm{d}\tau )}{2}{Y}_{\mathrm{\tau \tau }}{|}_p^q+\frac{(\mathrm{d}\tau {)}^2}{6}{Y}_{\mathrm{\tau \tau \tau }}{|}_p^q+\cdots )\\ & \times (1+2\frac{{\kappa }_2\mathrm{d}\tau }{(\mathrm{d}x{)}^2}+\mathrm{CM}\mathrm{d}\tau )\\ & ={\kappa }_2(Y_{\mathrm{xx}}{|}_p^q+\frac{(\mathrm{d}x{)}^2}{12}{Y}_{\mathrm{xxxx}}{|}_p^q+\cdots )\\ & -\mathrm{CM}(Y{|}_p^q)+C(X_p^p)(Y{|}_p^q).\end{array}$(48) Putting $\mathrm{d}\tau =(\mathrm{d}x{)}^3$ and taking $\mathrm{d}x\to 0$, the above equation simplifies to the equation given below (49) $\begin{array}{r}{Y}_{\tau }={\kappa }_2\mathrm{\Delta }Y+C(X-M)Y.\end{array}$(49)

4. Result and discussion

The initial distribution of the number of Prey and Predator are given in Figure 1. 3D mesh simulations and 2D solutions for prey and predator using backward Euler scheme are given in Figures 2– 9. Simulations 2, 4, 6 and 8 converge to all the four equilibrium points verifying the theoretical part, whereas in the simulations 3, 5, 7 and 9 one or both simulations fails to converge the true equilibrium points when increasing the stepsize $\mathrm{d}\tau$. When using Crank–Nicolson scheme, 3D mesh plots and 2D solutions for $X(x,\tau )$ and $Y(x,\tau )$ are given in Figures 10– 17. Simulations 10,12,14, and 16 converge to true equilibrium points, whereas in Figures 11, 13, 15 and 17, one or both populations fails to converge to true equilibrium points when taking temporal stepsize higher than the stepsize taken for 10, 12 14 and 16. When using NSFD scheme for the system (17(17) $\begin{array}{r}\{\begin{array}{rl}{X}_{\tau }& =X(X-L)(K-X)-\mathrm{XY}+{\kappa }_1\mathrm{\Delta X},\\ Y_{\tau }& =C(X-M)Y+{\kappa }_2\mathrm{\Delta Y}.\end{array}\end{array}$(17) ), the simulations for both populations converge to the stable solutions even for $\mathrm{d}\tau =1,5,10,100$ (Figures 18–25). All the simulations are done on Intel Corei5, 7th generation Dell Laptop.

Figure 2. Numerical solutions of the system (17(17) $\begin{array}{r}\{\begin{array}{rl}{X}_{\tau }& =X(X-L)(K-X)-\mathrm{XY}+{\kappa }_1\mathrm{\Delta X},\\ Y_{\tau }& =C(X-M)Y+{\kappa }_2\mathrm{\Delta Y}.\end{array}\end{array}$(17) ) using backward Euler scheme with temporal step-size $\mathrm{d}\tau =1$. Here, $X_0=0.8x$ and $Y_0=0.02x$ for 0<x< = 0.5 and $X_0=0.8(1-x)$ and $Y_0=0.02(1-x)$ for 0.5<x< = 1, with L = 0.2, M = 0.1, K = 1, C = 1, ${\kappa }_1=0.1$, and ${\kappa }_2=0.1$. Top: The left side of the plot shows the profile of prey $X(\tau ,x)$, while the right hand side shows the profile of predator $Y(\tau ,x)$. Bottom: The solution of $X(\tau ,x)$ and $Y(\tau ,x)$ tends to the first steady state i.e. $E_1=(0,0)$. The elapsed time is 0.074853 s.

Figure 3. Numerical solutions of the system (17(17) $\begin{array}{r}\{\begin{array}{rl}{X}_{\tau }& =X(X-L)(K-X)-\mathrm{XY}+{\kappa }_1\mathrm{\Delta X},\\ Y_{\tau }& =C(X-M)Y+{\kappa }_2\mathrm{\Delta Y}.\end{array}\end{array}$(17) ) using backward Euler scheme with temporal step-size $\mathrm{d}\tau =3$. Here, $X_0=0.8x$ and $Y_0=0.02x$ for 0<x< = 0.5 and $X_0=0.8(1-x)$ and $Y_0=0.02(1-x)$ for 0.5<x< = 1, with $L=0.2,M=0.1,K=1,C=1,{\kappa }_1=0.1$, and ${\kappa }_2=0.1$. Top: The left side of the plot shows the profile of prey $X(\tau ,x)$, while the right hand side shows the profile of predator $Y(\tau ,x)$. Bottom: The solution of $X(\tau ,x)$ fails to converge the steady state, whereas $Y(\tau ,x)$ tends to the true steady state. The elapsed time is 0.062071 s.

Figure 4. Numerical solutions of the system (17(17) $\begin{array}{r}\{\begin{array}{rl}{X}_{\tau }& =X(X-L)(K-X)-\mathrm{XY}+{\kappa }_1\mathrm{\Delta X},\\ Y_{\tau }& =C(X-M)Y+{\kappa }_2\mathrm{\Delta Y}.\end{array}\end{array}$(17) ) using backward Euler scheme with temporal step-size $\mathrm{d}\tau =1$. Here, $X_0=0.8x$ and $Y_0=0.02x$ for 0<x< = 0.5 and $X_0=0.8(1-x)$ and $Y_0=0.02(1-x)$ for 0.5<x< = 1, with L = 0.7, M = 1.2, K = 0.1, C = 1, ${\kappa }_1=0.1$ and ${\kappa }_2=0.1$. Top: The left side of the plot shows the profile of prey $X(\tau ,x)$, while the right hand side shows the profile of predator $Y(\tau ,x)$. Bottom: The solution of $X(\tau ,x)$ and $Y(\tau ,x)$ tends to the second steady state i.e. $E_2=(L,0)$. The elapsed time is 0.145165 s.

Figure 5. Numerical solutions of the system (17(17) $\begin{array}{r}\{\begin{array}{rl}{X}_{\tau }& =X(X-L)(K-X)-\mathrm{XY}+{\kappa }_1\mathrm{\Delta X},\\ Y_{\tau }& =C(X-M)Y+{\kappa }_2\mathrm{\Delta Y}.\end{array}\end{array}$(17) ) using backward Euler scheme with temporal step-size $\mathrm{d}\tau =60/59$. Here, $X_0=0.8x$ and $Y_0=0.02x$ for 0<x< = 0.5 and $X_0=0.8(1-x)$ and $Y_0=0.02(1-x)$ for 0.5<x< = 1, with L = 0.7, M = 1.2, K = 0.1, C = 1, ${\kappa }_1=0.1$ and ${\kappa }_2=0.1$. Top: The left side of the plot shows the profile of prey $X(\tau ,x)$, while the right hand side shows the profile of predator $Y(\tau ,x)$. Bottom: The solution of $X(\tau ,x)$ and $Y(\tau ,x)$ fails to converge the true steady state $E_2$. The elapsed time is 0.153060 s.

Figure 6. Numerical solutions of the system (17(17) $\begin{array}{r}\{\begin{array}{rl}{X}_{\tau }& =X(X-L)(K-X)-\mathrm{XY}+{\kappa }_1\mathrm{\Delta X},\\ Y_{\tau }& =C(X-M)Y+{\kappa }_2\mathrm{\Delta Y}.\end{array}\end{array}$(17) ) using backward Euler scheme with temporal step-size $\mathrm{d}\tau =1$. Here, $u_0=0.8x$ and $v_0=0.02x$ for 0<x< = 0.5 and $u_0=0.8(1-x)$ and $v_0=0.02(1-x)$ for 0.5<x< = 1, with L = 0.1, M = 1, K = 0.2, C = 1, ${\kappa }_1=0.1$ and ${\kappa }_2=0.1$. Top: The left side of the plot shows the profile of prey $X(\tau ,x)$, while the right hand side shows the profile of predator $Y(\tau ,x)$. Bottom: The solution of $U(\tau ,x)$ and $V(\tau ,x)$ tends to the second steady state i.e. $E_3=(K,0)$. The elapsed time is 0.143194 s.

Figure 7. Numerical solutions of the system (17(17) $\begin{array}{r}\{\begin{array}{rl}{X}_{\tau }& =X(X-L)(K-X)-\mathrm{XY}+{\kappa }_1\mathrm{\Delta X},\\ Y_{\tau }& =C(X-M)Y+{\kappa }_2\mathrm{\Delta Y}.\end{array}\end{array}$(17) ) using backward Euler scheme with temporal step-size $\mathrm{d}\tau =75/60$. Here, $X_0=0.8x$ and $Y_0=0.02x$ for 0<x< = 0.5 and $U_0=0.8(1-x)$ and $V_0=0.02(1-x)$ for 0.5<x< = 1, with L = 0.1, M = 1, K = 0.2, C = 1, ${\kappa }_1=0.1$ and ${\kappa }_2=0.1$. Top: The left side of the plot shows the profile of prey $X(\tau ,x)$, while the right hand side shows the profile of predator $Y(\tau ,x)$. Bottom: The solution of $X(\tau ,x)$ and $Y(\tau ,x)$ fails to converge the true steady state $E_3$. The elapsed time is 0.130104 s.

Figure 8. Numerical solutions of the system (17(17) $\begin{array}{r}\{\begin{array}{rl}{X}_{\tau }& =X(X-L)(K-X)-\mathrm{XY}+{\kappa }_1\mathrm{\Delta X},\\ Y_{\tau }& =C(X-M)Y+{\kappa }_2\mathrm{\Delta Y}.\end{array}\end{array}$(17) ) using backward Euler scheme with temporal step-size $\mathrm{d}\tau =1$. Here, $X_0=0.8x$ and $Y_0=0.02x$ for 0<x< = 0.5 and $X_0=0.8(1-x)$ and $Y_0=0.02(1-x)$ for 0.5<x< = 1, with L = 0.2, M = 1, K = 1.2, C = 1, ${\kappa }_1=0.1$ and ${\kappa }_2=0.1$. Top: The left side of the plot shows the profile of prey $X(\tau ,x)$, while the right hand side shows the profile of predator $Y(\tau ,x)$. Bottom: The solution of $X(\tau ,x)$ and $Y(\tau ,x)$ tends to the second steady state i.e. $E_4=(M,(M-L)(K-M))$. The elapsed time is 0.124455 s.

Figure 9. Numerical solutions of the system (17(17) $\begin{array}{r}\{\begin{array}{rl}{X}_{\tau }& =X(X-L)(K-X)-\mathrm{XY}+{\kappa }_1\mathrm{\Delta X},\\ Y_{\tau }& =C(X-M)Y+{\kappa }_2\mathrm{\Delta Y}.\end{array}\end{array}$(17) ) using backward Euler scheme with temporal step-size $\mathrm{d}\tau =79/60$. Here, $X_0=0.8x$ and $Y_0=0.02x$ for 0<x< = 0.5 and $X_0=0.8(1-x)$ and $Y_0=0.02(1-x)$ for 0.5<x< = 1, with L = 0.2, M = 1, K = 1.2, C = 1, ${\kappa }_1=0.1$ and ${\kappa }_2=0.1$. Top: The left side of the plot shows the profile of prey $X(\tau ,x)$, while the right hand side shows the profile of predator $Y(\tau ,x)$. Bottom: The solution of $X(\tau ,x)$ and $Y(\tau ,x)$ fails to converge the true steady states. The elapsed time is 0.031224 s.

Figure 10. Numerical solutions of the system (17(17) $\begin{array}{r}\{\begin{array}{rl}{X}_{\tau }& =X(X-L)(K-X)-\mathrm{XY}+{\kappa }_1\mathrm{\Delta X},\\ Y_{\tau }& =C(X-M)Y+{\kappa }_2\mathrm{\Delta Y}.\end{array}\end{array}$(17) ) using Crank–Nicolson Method with temporal step-size $\mathrm{d}\tau =22/30$. Here, $X_0=0.8x$ and $Y_0=0.02x$ for 0<x< = 0.5 and $X_0=0.8(1-x)$ and $Y_0=0.02(1-x)$ for 0.5<x< = 1, with L = 0.2, M = 0.1, K = 1, C = 1, ${\kappa }_1=0.1$, and ${\kappa }_2=0.1$. Top: The left side of the plot shows the profile of prey $X(\tau ,x)$, while the right hand side shows the profile of predator $Y(\tau ,x)$. Bottom: The solution of $X(\tau ,x)$ and $Y(\tau ,x)$ tends to the first steady state i.e. $E_1=(0,0)$. The elapsed time is 0.100409 s.

Figure 11. Numerical solutions of the system (17(17) $\begin{array}{r}\{\begin{array}{rl}{X}_{\tau }& =X(X-L)(K-X)-\mathrm{XY}+{\kappa }_1\mathrm{\Delta X},\\ Y_{\tau }& =C(X-M)Y+{\kappa }_2\mathrm{\Delta Y}.\end{array}\end{array}$(17) ) using Crank–Nicolson method with temporal step-size $\mathrm{d}\tau =25/30$. Here, $X_0=0.8x$ and $Y_0=0.02x$ for 0<x< = 0.5 and $X_0=0.8(1-x)$ and $Y_0=0.02(1-x)$ for 0.5<x< = 1, with L = 0.2, M = 0.1, K = 1, C = 1, ${\kappa }_1=0.1$, and ${\kappa }_2=0.1$. Top: The left side of the plot shows the profile of prey $X(\tau ,x)$, while the right hand side shows the profile of predator $Y(\tau ,x)$. Bottom: The solution of $X(\tau ,x)$ fails to preserve the positivity for the system, whereas $Y(\tau ,x)$ tends to the true steady state. The elapsed time is 0.098213 s.

Figure 12. Numerical solutions of the system (17(17) $\begin{array}{r}\{\begin{array}{rl}{X}_{\tau }& =X(X-L)(K-X)-\mathrm{XY}+{\kappa }_1\mathrm{\Delta X},\\ Y_{\tau }& =C(X-M)Y+{\kappa }_2\mathrm{\Delta Y}.\end{array}\end{array}$(17) ) using Crank–Nicolson method with temporal step-size $\mathrm{d}\tau =195/600$. Here, $X_0=0.8x$ and $Y_0=0.02x$ for 0<x< = 0.5 and $X_0=0.8(1-x)$ and $Y_0=0.02(1-x)$ for 0.50<x< = 1, with L=0.7, M=1.2, K=0.1, C=1, ${\kappa }_1=0.1$, and ${\kappa }_2=0.1$. Top: The left side of the plot shows the profile of prey $X(\tau ,x)$, while the right hand side shows the profile of predator $Y(\tau ,x)$. Bottom: The solution of $X(\tau ,x)$ and $Y(\tau ,x)$ tends to the second steady state i.e. $E_2=(L,0)$. The elapsed time is 0.154392 s.

Figure 13. Numerical solutions of the system (17(17) $\begin{array}{r}\{\begin{array}{rl}{X}_{\tau }& =X(X-L)(K-X)-\mathrm{XY}+{\kappa }_1\mathrm{\Delta X},\\ Y_{\tau }& =C(X-M)Y+{\kappa }_2\mathrm{\Delta Y}.\end{array}\end{array}$(17) ) using Crank–Nicolson method with temporal step-size $\mathrm{d}\tau =196/600$. Here, $X_0=0.8x$ and $Y_0=0.02x$ for 0<x< = 0.5 and $X_0=0.8(1-x)$ and $Y_0=0.02(1-x)$ for 0.5<x< = 1, with L = 0.7, M = 1.2, K = 0.1, C = 1, ${\kappa }_1=0.1$, and ${\kappa }_2=0.1$. Top: The left side of the plot shows the profile of prey $X(\tau ,x)$, while the right hand side shows the profile of predator $Y(\tau ,x)$. Bottom: The solution of $X(\tau ,x)$ tends to the true steady state, whereas $Y(\tau ,x)$ fails to preserve the positivity for the system. The elapsed time is 0.159403 s.

Figure 14. Numerical solutions of the system (17(17) $\begin{array}{r}\{\begin{array}{rl}{X}_{\tau }& =X(X-L)(K-X)-\mathrm{XY}+{\kappa }_1\mathrm{\Delta X},\\ Y_{\tau }& =C(X-M)Y+{\kappa }_2\mathrm{\Delta Y}.\end{array}\end{array}$(17) ) using Crank–Nicolson method with temporal step-size $\mathrm{d}\tau =192/600$. Here, $X_0=0.8x$ and $Y_0=0.02x$ for 0<x< = 0.5 and $X_0=0.8(1-x)$ and $Y_0=0.02(1-x)$ for 0.5<x< = 1, with L = 0.1, M = 1, K = 0.2, C = 1, ${\kappa }_1=0.1$, and ${\kappa }_2=0.1$. Top: The left side of the plot shows the profile of prey $X(\tau ,x)$, while the right hand side shows the profile of predator $Y(\tau ,x)$. Bottom: The solution of $X(\tau ,x)$ and $Y(\tau ,x)$ tends to the second steady state i.e. $E_3=(K,0)$. The elapsed time is 0.162833 s.

Figure 15. Numerical solutions of the system (17(17) $\begin{array}{r}\{\begin{array}{rl}{X}_{\tau }& =X(X-L)(K-X)-\mathrm{XY}+{\kappa }_1\mathrm{\Delta X},\\ Y_{\tau }& =C(X-M)Y+{\kappa }_2\mathrm{\Delta Y}.\end{array}\end{array}$(17) ) using Crank–Nicolson method with temporal step-size $\mathrm{d}\tau =193/600$. Here, $X_0=0.8x$ and $Y_0=0.02x$ for 0<x< = 0.5 and $u_0=0.8(1-x)$ and $v_0=0.02(1-x)$ for 0.5<x< = 1, with L = 0.1, M = 1, K = 0.2, C = 1, ${\kappa }_1=0.1$, and ${\kappa }_2=0.1$. Top: The left side of the plot shows the profile of prey $X(\tau ,x)$, while the right hand side shows the profile of predator $Y(\tau ,x)$. Bottom: The solution of $X(\tau ,x)$ tends to the true steady state, whereas $Y(\tau ,x)$ fails to preserve the positivity for the system. The elapsed time is 0.144537 s.

Figure 16. Numerical solutions of the system (17(17) $\begin{array}{r}\{\begin{array}{rl}{X}_{\tau }& =X(X-L)(K-X)-\mathrm{XY}+{\kappa }_1\mathrm{\Delta X},\\ Y_{\tau }& =C(X-M)Y+{\kappa }_2\mathrm{\Delta Y}.\end{array}\end{array}$(17) ) using Crank–Nicolson method with temporal step-size $\mathrm{d}\tau =67/200$. Here, $X_0=0.8x$ and $Y_0=0.02x$ for 0<x< = 0.5 and $X_0=0.8(1-x)$ and $Y_0=0.02(1-x)$ for 0.5<x< = 1, with L = 0.2, M = 1, K = 1.2, C = 1, ${\kappa }_1=0.1$, and ${\kappa }_2=0.1$. Top: The left side of the plot shows the profile of prey $X(\tau ,x)$, while the right hand side shows the profile of predator $Y(\tau ,x)$. Bottom: The solution of $X(\tau ,x)$ and $Y(\tau ,x)$ tends to the second steady state i.e. $E_4=(M,(M-L)(K-M))$. The elapsed time is 0.064987 s.

Figure 17. Numerical solutions of the system (17(17) $\begin{array}{r}\{\begin{array}{rl}{X}_{\tau }& =X(X-L)(K-X)-\mathrm{XY}+{\kappa }_1\mathrm{\Delta X},\\ Y_{\tau }& =C(X-M)Y+{\kappa }_2\mathrm{\Delta Y}.\end{array}\end{array}$(17) ) using Crank–Nicolson method with temporal step-size $\mathrm{d}\tau =70/200$. Here, $X_0=0.8x$ and $Y_0=0.02x$ for 0<x< = 0.5 and $X_0=0.8(1-x)$ and $Y_0=0.02(1-x)$ for 0.5<x< = 1, with L = 0.2, M = 1, K = 1.2, C = 1, ${\kappa }_1=0.1$, and ${\kappa }_2=0.1$. Top: The left side of the plot shows the profile of prey $X(\tau ,x)$, while the right hand side shows the profile of predator $Y(\tau ,x)$. Bottom: The solution of $X(\tau ,x)$ tends to the true steady state, whereas $Y(\tau ,x)$ fails to preserve the positivity for the system. The elapsed time is 0.067225 s.

Figure 18. Numerical solutions of the system (17(17) $\begin{array}{r}\{\begin{array}{rl}{X}_{\tau }& =X(X-L)(K-X)-\mathrm{XY}+{\kappa }_1\mathrm{\Delta X},\\ Y_{\tau }& =C(X-M)Y+{\kappa }_2\mathrm{\Delta Y}.\end{array}\end{array}$(17) ) using nonstandard FD method with temporal step-size $\mathrm{d}\tau =1$. Here, $X_0=0.8x$ and $Y_0=0.02x$ for 0<x< = 0.5 and $X_0=0.8(1-x)$ and $Y_0=0.02(1-x)$ for 0.5<x< = 1, with L = 0.2, M = 0.1, K = 1, C = 1, ${\kappa }_1=0.1$, and ${\kappa }_2=0.1$. Top: The left side of the plot shows the profile of prey $X(\tau ,x)$, while the right hand side shows the profile of predator $Y(\tau ,x)$. Bottom: The solution of $X(\tau ,x)$ and $Y(\tau ,x)$ tends to the first steady state $E_1=(0,0)$. The elapsed time is 0.044335 s.

Figure 19. Numerical solutions of the system (17(17) $\begin{array}{r}\{\begin{array}{rl}{X}_{\tau }& =X(X-L)(K-X)-\mathrm{XY}+{\kappa }_1\mathrm{\Delta X},\\ Y_{\tau }& =C(X-M)Y+{\kappa }_2\mathrm{\Delta Y}.\end{array}\end{array}$(17) ) using nonstandard FD method with temporal step-size $\mathrm{d}\tau =5$. Here, $X_0=0.8x$ and $Y_0=0.02x$ for 0<x< = 0.5 and $X_0=0.8(1-x)$ and $Y_0=0.02(1-x)$ for 0.5<x< = 1, with L = 0.2, M = 0.1, K = 1, C = 1, ${\kappa }_1=0.1$, and ${\kappa }_2=0.1$. Top: The left side of the plot shows the profile of prey $X(\tau ,x)$, while the right hand side shows the profile of predator $Y(\tau ,x)$. Bottom: The solution of $X(\tau ,x)$ and $Y(\tau ,x)$ tends to the first steady state. The elapsed time is 0.050633 s.

Figure 20. Numerical solutions of the system (17(17) $\begin{array}{r}\{\begin{array}{rl}{X}_{\tau }& =X(X-L)(K-X)-\mathrm{XY}+{\kappa }_1\mathrm{\Delta X},\\ Y_{\tau }& =C(X-M)Y+{\kappa }_2\mathrm{\Delta Y}.\end{array}\end{array}$(17) ) using nonstandard FD method with temporal step-size $\mathrm{d}\tau =1$. Here, $X_0=0.8x$ and $Y_0=0.02x$ for 0<x< = 0.5 and $X_0=0.8(1-x)$ and $Y_0=0.02(1-x)$ for 0.5<x< = 1, with L = 0.7, M = 1.2, K = 0.1, C = 1, ${\kappa }_1=0.1$, and ${\kappa }_2=0.1$. Top: The left side of the plot shows the profile of prey $X(\tau ,x)$, while the right hand side shows the profile of predator $Y(\tau ,x)$. Bottom: The solution of $X(\tau ,x)$ and $Y(\tau ,x)$ tends to the second steady state $E_2=(L,0)$. The elapsed time is 0.039919 s.

Figure 21. Numerical solutions of the system (17(17) $\begin{array}{r}\{\begin{array}{rl}{X}_{\tau }& =X(X-L)(K-X)-\mathrm{XY}+{\kappa }_1\mathrm{\Delta X},\\ Y_{\tau }& =C(X-M)Y+{\kappa }_2\mathrm{\Delta Y}.\end{array}\end{array}$(17) ) using nonstandard FD method with temporal step-size $\mathrm{d}\tau =5$. Here, $X_0=0.8x$ and $Y_0=0.02x$ for 0<x< = 0.5 and $X_0=0.8(1-x)$ and $Y_0=0.02(1-x)$ for 0.5<x< = 1, with L = 0.7, M = 1.2, K = 0.1, C = 1, ${\kappa }_1=0.1$, and ${\kappa }_2=0.1$. Top: The left side of the plot shows the profile of prey $X(\tau ,x)$, while the right hand side shows the profile of predator $Y(\tau ,x)$. Bottom: The solution of $X(\tau ,x)$ and $Y(\tau ,x)$ converges the true steady state $E_2$. The elapsed time is 0.037491  s.

Figure 22. Numerical solutions of the system (17(17) $\begin{array}{r}\{\begin{array}{rl}{X}_{\tau }& =X(X-L)(K-X)-\mathrm{XY}+{\kappa }_1\mathrm{\Delta X},\\ Y_{\tau }& =C(X-M)Y+{\kappa }_2\mathrm{\Delta Y}.\end{array}\end{array}$(17) ) using nonstandard FD method with temporal step-size $\mathrm{d}\tau =1$. Here, $X_0=0.8x$ and $Y_0=0.02x$ for 0<x< = 0.5 and $X_0=0.8(1-x)$ and $Y_0=0.02(1-x)$ for 0.5<x< = 1, with L = 0.1, M = 1, K = 0.2, C = 1, ${\kappa }_1=0.1$, and ${\kappa }_2=0.1$. Top: The left side of the plot shows the profile of prey $X(\tau ,x)$, while the right hand side shows the profile of predator $Y(\tau ,x)$. Bottom: The solution of $X(\tau ,x)$ and $Y(\tau ,x)$ tends to the second steady state $E_3=(K,0)$. The elapsed time is 0.034039 s.

Figure 23. Numerical solutions of the system (17(17) $\begin{array}{r}\{\begin{array}{rl}{X}_{\tau }& =X(X-L)(K-X)-\mathrm{XY}+{\kappa }_1\mathrm{\Delta X},\\ Y_{\tau }& =C(X-M)Y+{\kappa }_2\mathrm{\Delta Y}.\end{array}\end{array}$(17) ) using nonstandard FD method with temporal step-size $\mathrm{d}\tau =5$. Here, $X_0=0.8x$ and $Y_0=0.02x$ for 0<x< = 0.5 and $X_0=0.8(1-x)$ and $Y_0=0.02(1-x)$ for 0.5<x< = 1, with L = 0.1, M =1, K = 0.2, C = 1, ${\kappa }_1=0.1$, and ${\kappa }_2=0.1$. Top: The left side of the plot shows the profile of prey $X(\tau ,x)$, while the right hand side shows the profile of predator $Y(\tau ,x)$. Bottom: The solution of $X(\tau ,x)$ and $Y(\tau ,x)$ converges the true steady state $E_3$. The elapsed time is 0.044540 s.

Figure 24. Numerical solutions of the system (17(17) $\begin{array}{r}\{\begin{array}{rl}{X}_{\tau }& =X(X-L)(K-X)-\mathrm{XY}+{\kappa }_1\mathrm{\Delta X},\\ Y_{\tau }& =C(X-M)Y+{\kappa }_2\mathrm{\Delta Y}.\end{array}\end{array}$(17) ) using nonstandard FD method with temporal step-size $\mathrm{d}\tau =1$. Here, $X_0=0.8x$ and $Y_0=0.02x$ for 0<x< = 0.5 and $X_0=0.8(1-x)$ and $Y_0=0.02(1-x)$ for 0.5<x< = 1, with L = 0.2, M = 1, K = 1.2, C = 1, ${\kappa }_1=0.1$, and ${\kappa }_2=0.1$. Top: The left side of the plot shows the profile of prey $X(\tau ,x)$, while the right hand side shows the profile of predator $Y(\tau ,x)$. Bottom: The solution of $X(\tau ,x)$ and $Y(\tau ,x)$ tends to the second steady state i.e. $E_4=(M,(M-L)(K-M))$. The elapsed time is 0.041376 s.

Figure 25. Numerical solutions of the system (17(17) $\begin{array}{r}\{\begin{array}{rl}{X}_{\tau }& =X(X-L)(K-X)-\mathrm{XY}+{\kappa }_1\mathrm{\Delta X},\\ Y_{\tau }& =C(X-M)Y+{\kappa }_2\mathrm{\Delta Y}.\end{array}\end{array}$(17) ) using nonstandard FD method with temporal step-size $\mathrm{d}\tau =5$. Here, $X_0=0.8x$ and $Y_0=0.02x$ for 0<x< = 0.5 and $X_0=0.8(1-x)$ and $Y_0=0.02(1-x)$ for 0.5<x< = 1, with L = 0.2, M = 1, K = 1.2, C = 1, ${\kappa }_1=0.1$, and ${\kappa }_2=0.1$, Top: The left side of the plot shows the profile of prey $X(\tau ,x)$, while the right hand side shows the profile of predator $Y(\tau ,x)$. Bottom: The solution of $X(\tau ,x)$ and $Y(\tau ,x)$ converges the true steady states. The elapsed time is 0.034553 s.

5. Conclusions

In this manuscript, we have proposed a BB model by strong Allee effects. The main objective is to investigate the effects of diffusion term in the dynamics of the BB model. We have developed three finite difference schemes to solve system (4(4) $\begin{array}{r}\frac{\mathrm{d}X}{\mathrm{d}t}=-X(L-X)(K-X)-\mathrm{XY},\frac{\mathrm{d}Y}{\mathrm{d}t}=-C(M-X)Y.\end{array}$(4) ). We used Taylor's theorem to check our proposed nonstandard FD scheme is accurate in first order and consistent in spatiotemporal. The analysis of nonstandard FD is done by using the Von-Neumann stability. The graphical portraits exhibits the fact that the proposed nonstandard FD is stable for both small or large temporal step-size. On the other hand, both well-known methods fail to preserve the positivity property and show divergence on different step size $\mathrm{d}\tau$. Hence the proposed scheme is more effective, accurate, reliable and consistent for under observed model. Finally, the authors have concluded that NSFD scheme has an edge on the other two schemes.

Acknowledgments

Acknowledgements

All authors contributed equally and significantly in writing this article. All authors read and approved the final manuscript.

Disclosure statement

No potential conflict of interest was reported by the author(s).