An explicit fourth-order accurate compact method for the Allen-Cahn equation

1.   Introduction

The dynamic relationship among predators and prey has for some time been and will keep on being one of the predominant subjects in both mathematical ecology and ecology because of its widespread presence and significance ^[1]. There are a lot of research articles about the dynamic behavior of predator-prey system without and with various types of functional responses. It is important to note that outcomes of hiding behavior of prey on the dynamics of prey-predator interactions can be predictably significant. Although the effects of prey refuges on population dynamics are actually very complex, they can be divided into two categories for modelling purposes ^[2]: the first effect, which affects prey growth positively and predator growth negatively, is the reduction of prey mortality as a result of decreased predation success. The second may be the exchange and spin-off of hiding behavior of prey, which may or may not be beneficial for all interacting populations. It is anticipated that, in recent years, most of the work has been done on dynamical characteristics of continuous-time predator-prey system with a prey refuge designated by differential equation without time delay ^[3,4]. For instance, Leslie ^[5,6] has investigated the following predator-prey model where carrying capacity of the predator's environment is proportional to the number of prey:

where $H$ is prey density, $P$ represents predator density, $r_1, a_1, b_1, r_2, a_2$ are positive constants. Biologically $r_2$ and $r_1$ denote intrinsic growth rate of predator and prey, respectively; carrying capacity of prey is $\frac{r_1}{b_1}$ and carrying capacity of predator is $\frac{r_2H}{a_2}$. Further, Chen et al. ^[7] have extended the model of Leslie ^[5,6] by incorporating a refuge defending $mH$ of the prey, and a resulting continuous-time system designated by differential equations takes the form:

while $m\in[0, 1)$. This leaves $(1-m)H$ of the prey available to the predator. In contrast to continuous-time models, discrete-time models designated by maps or different equations are more applicable than differential models, if populations have non-overlapping generations, and these results also provide more efficient computational results for numerical simulations ^[8,9]. So, for non-overlapping generations, many mathematicians have investigated the dynamical behavior of discrete-time models as compared to the continuous-time model (see ^{[10,11,12,13,14,15,16,17,18,19,20,21,22]}). For instance, Zhuang and Wen ^[23] have studied the local dynamics of the following model which is a discrete form of (1.2) by forward Euler's formula:

where $h$ is the step size. More precisely, Zhuang and Wen ^[23] have proved that boundary and interior fixed points of the discrete model (1.3) is a sink, saddle, source and non-hyperbolic under certain parametric condition(s). Motivated by the mentioned studies, the purpose of this paper is to explore further complicate dynamical analysis of the discrete model (1.3), and so our key contributions in this regard include:

(ⅰ) Topological classifications at fixed points of the discrete model (1.3).

(ⅱ) Existence of bifurcations sets at fixed points.

(ⅲ) Detail bifurcation analysis at fixed points of the discrete model (1.3).

(ⅳ) Study of chaos by state feedback control strategy.

(ⅴ) Verification of theoretical results numerically.

(ⅵ) Influence of prey refuge.

The organization of the rest of the paper is as fellows: local stability with topological classifications at fixed points of a discrete model (1.3) is explored in Section 2. In Section 3, we have studied detailed bifurcation analysis at fixed points, whereas Section 4 is about the study of chaos control of the discrete model (1.3). The numerical simulation that verifies the correctness of theoretical results are presented in Section 5. In Section 6, we will discuss the influence of prey refuge, whereas the conclusion is given in Section 7.

2.   Local stability of fixed points

In this section, we will analyze the local stability of fixed points of discrete model (1.3). Here, due to biological meaning of discrete model (1.3), first we will pick up nonnegative fixed points. It is clear that if $E_{HP}\left(H, P\right)$ is a fixed point of discrete model (1.3), then

From (1.1), the straightforward calculation yields the following results.

Theorem 2.1. (ⅰ) $\forall\ r_1, r_2, a_1, a_2, b_1, m, h > 0$, $E_{H0}(\frac{r_1}{b_1}, 0)$ is a boundary fixed point of discrete model (1.3);

(ⅱ) $E^+_{HP}\left(\frac{r_1a_2}{a_2b_1+a_1r_2(1-m)^2}, \frac{r_1r_2(1-m)}{a_2b_1+a_1r_2(1-m)^2}\right)$ is the interior fixed point of discrete model (1.3) if $m < 1$.

Now using linearization, at $E_{HP}\left(H, P\right)$ following variational matrix $V|_{E_{HP}\left(H, P\right)}$ is constructed:

Now, we will study local stability at boundary and interior fixed points of model (1.3) by the stability theory ^[24,25,26]. It should be noted that at $E_{H0}(\frac{r_1}{b_1}, 0)$, (2.2) gives

The characteristic roots of $V|_{E_{H0}\left(\frac{r_1}{b_1}, 0\right)}$ are

Based on (2.4), one can easily obtain the following results.

Theorem 2.2. (A) $\forall\ r_1, r_2, a_1, a_2, b_1, m, h > 0$, $E_{H0}(\frac{r_1}{b_1}, 0)$ of discrete model (1.3) is never sink;

(B) $E_{H0}(\frac{r_1}{b_1}, 0)$ of discrete model (1.3) is a source if

(C) $E_{H0}\left(\frac{r_1}{b_1}, 0\right)$ of discrete model (1.3) is a saddle if

(D) $E_{H0}\left(\frac{r_1}{b_1}, 0\right)$ of discrete model (1.3) is non-hyperbolic if

Now at $E_{HP}^+\left(\frac{r_1a_2}{a_2b_1+a_1r_2(1-m)^2}, \frac{r_1r_2(1-m)}{a_2b_1+a_1r_2(1-m)^2}\right)$, (2.2) gives

The characteristic equation of $V|_{E_{HP}^+\left(\frac{r_1a_2}{a_2b_1+a_1r_2(1-m)^2}, \frac{r_1r_2(1-m)}{a_2b_1+a_1r_2(1-m)^2}\right)}$ is

with

Finally, roots of (2.9) are

where

Now, following two theorems are established based on sign of $\Delta$, i.e. $\Delta < 0$ and $\Delta\geq0$, respectively.

Theorem 2.3. If $\Delta = \left(2-h r_2-\frac{a_2b_1h r_1}{a_2b_1+a_1r_2(1-m)^2}\right)^2-4\left(1+h r_2(h r_1-1)-\frac{a_2b_1h r_1}{a_2b_1+a_1r_2(1-m)^2}\right) < 0$ then at $E_{HP}^+\left(\frac{r_1a_2}{a_2b_1+a_1r_2(1-m)^2}, \frac{r_1r_2(1-m)}{a_2b_1+a_1r_2(1-m)^2}\right)$ following results hold:

(A) $E_{HP}^+\left(\frac{r_1a_2}{a_2b_1+a_1r_2(1-m)^2}, \frac{r_1r_2(1-m)}{a_2b_1+a_1r_2(1-m)^2}\right)$ of discrete model (1.3) is a stable focus if

with

(B) $E_{HP}^+\left(\frac{r_1a_2}{a_2b_1+a_1r_2(1-m)^2}, \frac{r_1r_2(1-m)}{a_2b_1+a_1r_2(1-m)^2}\right)$ of discrete model (1.3) is an unstable focus if

(C) $E_{HP}^+\left(\frac{r_1a_2}{a_2b_1+a_1r_2(1-m)^2}, \frac{r_1r_2(1-m)}{a_2b_1+a_1r_2(1-m)^2}\right)$ of discrete model (1.3) is non-hyperbolic if

Theorem 2.4. If $\Delta = \left(2-h r_2-\frac{a_2b_1h r_1}{a_2b_1+a_1r_2(1-m)^2}\right)^2-4\left(1+h r_2(h r_1-1)-\frac{a_2b_1h r_1}{a_2b_1+a_1r_2(1-m)^2}\right)\geq0$ then at $E_{HP}^+\left(\frac{r_1a_2}{a_2b_1+a_1r_2(1-m)^2}, \frac{r_1r_2(1-m)}{a_2b_1+a_1r_2(1-m)^2}\right)$ following results hold:

(A) $E_{HP}^+\left(\frac{r_1a_2}{a_2b_1+a_1r_2(1-m)^2}, \frac{r_1r_2(1-m)}{a_2b_1+a_1r_2(1-m)^2}\right)$ of discrete model (1.3) is a stable node if

with

(B) $E_{HP}^+\left(\frac{r_1a_2}{a_2b_1+a_1r_2(1-m)^2}, \frac{r_1r_2(1-m)}{a_2b_1+a_1r_2(1-m)^2}\right)$ of discrete model (1.3) is an unstable node if

(C) $E_{HP}^+\left(\frac{r_1a_2}{a_2b_1+a_1r_2(1-m)^2}, \frac{r_1r_2(1-m)}{a_2b_1+a_1r_2(1-m)^2}\right)$ of discrete model (1.3) is non-hyperbolic if

3.   Analysis of bifurcation

In this section, we study the flip bifurcation at boundary fixed point $E_{H0}(\frac{r_1}{b_1}, 0)$, and flip and hopf bifurcations at $E_{HP}^+\left(\frac{r_1a_2}{a_2b_1+a_1r_2(1-m)^2}, \frac{r_1r_2(1-m)}{a_2b_1+a_1r_2(1-m)^2}\right)$ of discrete model (1.3) by center manifold theorem and bifurcation theory ^{[28,29,30,31,32]}. If condition (2.7) of Theorem 2.2 holds, then from (2.4) one gets $\lambda_1|_{(2.7)} = -1$ but $\lambda_2|_{(2.7)} = 1+h r_2\ne 1$ or $-1$, which implies that at $E_{H0}(\frac{r_1}{b_1}, 0)$ discrete model (1.3) may undergoes flip bifurcation if $(h, r_1, r_2, b_1, a_1, a_2, m)$ located in the set:

But following this theorem, it follows that if $(h, r_1, r_2, b_1, a_1, a_2, m)\in \mathcal{F}|_{E_{H0}\left(\frac{r_1}{b_1}, 0\right)}$ then discrete model (1.3) must undergo flip bifurcation.

Theorem 3.1. Discrete model (1.3) undergoes flip bifurcation if $(h, r_1, r_2, b_1, a_1, a_2, m)\in \mathcal{F}|_{E_{H0}\left(\frac{r_1}{b_1}, 0\right)}$.

Proof. Since, with respect to $P = 0$, discrete model (1.3) is invariant. So, if it is restricted on $P = 0$ then

Now from (3.2), we denote

Now if $r_1 = {r_1}^\ast = \frac{2}{h}$ and $H = H^\ast = \frac{r_1}{b_1}$ then from (3.3) one gets

and

From (3.4)–(3.6), and Table 4.4.1 of ^[30] one can concluded that if $(h, r_1, r_2, b_1, a_1, a_2, m)\in \mathcal{F}|_{E_{H0}\left(\frac{r_1}{b_1}, 0\right)}$ then flip bifurcation must exist at $E_{H0}(\frac{r_1}{b_1}, 0)$. Additionally $\Omega_1 = \left(\frac{\partial^2 f}{\partial H^2}\right)\left(\frac{\partial f}{\partial {r_1}}\right)+2\left(\frac{\partial^2 f}{\partial H\partial r_1}\right)\biggr|_{r_1 = {r_1}^\ast = \frac{2}{h}, \ H = H^\ast = \frac{r_1}{b_1}} = -2h$, $\Omega_2 = \left(\frac{\partial^3 f}{\partial H^3}\right)+3 \left(\left(\frac{\partial^2 f}{\partial H^2}\right)\right)^2 = 12h^2b_1^2$ and finally $\Omega_1\Omega_2 = -24b_1^2h^3 < 0$ which implies that at $E_{H0}(\frac{r_1}{b_1}, 0)$ discrete model (1.3) undergoes supercritical flip bifurcation.

Now if $\Delta = \left(2-h r_2-\frac{a_2b_1h r_1}{a_2b_1+a_1r_2(1-m)^2}\right)^2-4\left(1+h r_2(h r_1-1)-\frac{a_2b_1h r_1}{a_2b_1+a_1r_2(1-m)^2}\right) < 0$ then roots of (2.11) are pair of complex conjugate satisfying $\left|\lambda_{1, 2}\right|_{ (2.16)} = 1$ which implies that hopf bifurcation may exists if $(h, r_1, r_2, b_1, a_1, a_2, m)$ are locate in the set:

But following result follows that if $(h, r_1, r_2, b_1, a_1, a_2, m)\in \mathcal{N}|_{E_{HP}^+\left(\frac{r_1a_2}{a_2b_1+a_1r_2(1-m)^2}, \frac{r_1r_2(1-m)}{a_2b_1+a_1r_2(1-m)^2}\right)}$ then at $E_{HP}^+\left(\frac{r_1a_2}{a_2b_1+a_1r_2(1-m)^2}, \frac{r_1r_2(1-m)}{a_2b_1+a_1r_2(1-m)^2}\right)$ discrete model (1.3) undergoes hopf bifurcation.

Theorem 3.2. If $(h, r_1, r_2, b_1, a_1, a_2, m)\in \mathcal{N}|_{E_{HP}^+\left(\frac{r_1a_2}{a_2b_1+a_1r_2(1-m)^2}, \frac{r_1r_2(1-m)}{a_2b_1+a_1r_2(1-m)^2}\right)}$ then at interior equilibrium $E_{HP}^+\left(\frac{r_1a_2}{a_2b_1+a_1r_2(1-m)^2}, \frac{r_1r_2(1-m)}{a_2b_1+a_1r_2(1-m)^2}\right)$, discrete model (1.3) undergo hopf bifurcation.

Proof. Since $(h, r_1, r_2, b_1, a_1, a_2, m)\in \mathcal{N}|_{E_{HP}^+\left(\frac{r_1a_2}{a_2b_1+a_1r_2(1-m)^2}, \frac{r_1r_2(1-m)}{a_2b_1+a_1r_2(1-m)^2}\right)}$ and clearly $r_1$ is the bifurcation parameter. So, if $r_1$ varies in a small neighborhood of $r_1^*$, $i.e$, $r_1 = r_1^*+\epsilon$ where $\epsilon < < 1$ then model (1.3) becomes

with $E_{HP}^+\left(\frac{\left(r_1^*+\epsilon\right)a_2}{a_2b_1+a_1r_2(1-m)^2}, \frac{\left(r_1^*+\epsilon\right)r_2(1-m)}{a_2b_1+a_1r_2(1-m)^2}\right)$ is the interior fixed point. Now the roots of characteristic equation of $V|_{E_{HP}^+\left(\frac{\left(r_1^*+\epsilon\right)a_2}{a_2b_1+a_1r_2(1-m)^2}, \frac{\left(r_1^*+\epsilon\right)r_2(1-m)}{a_2b_1+a_1r_2(1-m)^2}\right)}$ at $E_{HP}^+\left(\frac{\left(r_1^*+\epsilon\right)a_2}{a_2b_1+a_1r_2(1-m)^2}, \frac{\left(r_1^*+\epsilon\right)r_2(1-m)}{a_2b_1+a_1r_2(1-m)^2}\right)$ of model (3.8) is

where

From (3.9) and (3.10) we have

with $\frac{d|\lambda_{1, 2}|}{d \epsilon}|_{\epsilon = 0} = \frac{1}{2}\left(h^2 r_2-\frac{a_2b_1h}{a_2b_1+a_1r_2(1-m)^2}\right)\not = 0$. Moreover, for occurrence of hopf bifurcation at $E_{HP}^+\left(\frac{r_1a_2}{a_2b_1+a_1r_2(1-m)^2}, \frac{r_1r_2(1-m)}{a_2b_1+a_1r_2(1-m)^2}\right)$ it is also require that $\lambda_{1, 2}^v\neq 1, \ v = 1, \cdots, 4$ if $\epsilon = 0$ that corresponds to $p(0)\neq -2, 0, 1, 2$. But $q(0) = 1$ if (2.16) holds, and so $p(0)\not = -2, 2$. Hence $p(0)\not = 0, 1$ which gives

Now using the following transformations in order to transform $E_{HP}^+\left(\frac{r_1a_2}{a_2b_1+a_1r_2(1-m)^2}, \frac{r_1r_2(1-m)}{a_2b_1+a_1r_2(1-m)^2}\right)$ of (3.8) to origin:

with $H^* = \frac{r_1a_2}{a_2b_1+a_1r_2(1-m)^2} \ and\ P^* = \frac{r_1r_2(1-m)}{a_2b_1+a_1r_2(1-m)^2}$. From (3.13) and (3.8) one writes

Now if $\epsilon = 0$ then we will explore normal form of (3.14). On expanding (3.14) up to order-$2^{nd}$ at $F_{00}\left(0, 0\right)$ we write

with

Now using the following transformation, we transform linear part of (3.15) to canonical form

with

From (3.17) and (3.15) we write

where

and

From (3.20) we have

Finally, following condition required to be non-zero for (3.19) undergo hopf bifurcation

where

The calculation yields

From (3.25) and (3.23) if we get $\ell\not = 0$ as $(h, r_1, r_2, b_1, a_1, a_2, m)\in \mathcal{N}|_{E_{HP}^+\left(\frac{r_1a_2}{a_2b_1+a_1r_2(1-m)^2}, \frac{r_1r_2(1-m)}{a_2b_1+a_1r_2(1-m)^2}\right)}$ then at $E_{HP}^+\left(\frac{r_1a_2}{a_2b_1+a_1r_2(1-m)^2}, \frac{r_1r_2(1-m)}{a_2b_1+a_1r_2(1-m)^2}\right)$ discrete model (1.3) undergoes hopf bifurcation. Further supercritical (respectively subcritical) hopf bifurcation take place if $\ell < 0$ (respectively $\ell > 0$).

Now if $\Delta = \left(2-h r_2-\frac{a_2b_1h r_1}{a_2b_1+a_1r_2(1-m)^2}\right)^2-4\left(1+h r_2(h r_1-1)-\frac{a_2b_1h r_1}{a_2b_1+a_1r_2(1-m)^2}\right) > 0$ and condition (2.20) of Theorem 2.4 holds then $\lambda_1|_{(2.20)} = -1$ but $\lambda_2|_{(2.20)} = 3-h r_2-\frac{2 a_2b_1(h r_2-2)}{a_1h(1-m)^2{r_2}^2+a_2b_1(h r_2-2)} \ne 1 \ {\mathrm{or}} \ -1$ that concludes the fact that at $E_{HP}^+\left(\frac{r_1a_2}{a_2b_1+a_1r_2(1-m)^2}, \frac{r_1r_2(1-m)}{a_2b_1+a_1r_2(1-m)^2}\right)$ discrete model (1.3) may undergoes flip bifurcation if $\left(h, r_1, r_2, b_1, a_1, a_2, m\right)$ are located in the set:

But following result follows that if $(h, r_1, r_2, b_1, a_1, a_2, m)\in \mathcal{F}|_{E_{HP}^+\left(\frac{r_1a_2}{a_2b_1+a_1r_2(1-m)^2}, \frac{r_1r_2(1-m)}{a_2b_1+a_1r_2(1-m)^2}\right)}$ then at $E_{HP}^+\left(\frac{r_1a_2}{a_2b_1+a_1r_2(1-m)^2}, \frac{r_1r_2(1-m)}{a_2b_1+a_1r_2(1-m)^2}\right)$ discrete model (1.3) undergoes flip bifurcation.

Theorem 3.3. If $(h, r_1, r_2, b_1, a_1, a_2, m)\in \mathcal{F}|_{E_{HP}^+\left(\frac{r_1a_2}{a_2b_1+a_1r_2(1-m)^2}, \frac{r_1r_2(1-m)}{a_2b_1+a_1r_2(1-m)^2}\right)}$ then discrete model (1.3) undergoes flip bifurcation.

Proof. If $r_1$ varies in a small neighborhood of $r_1^*$ then model (1.3) takes the form (3.8) which further transform into the following form

where

by (3.13). By

(3.27) becomes

where

Now center manifold $M^cF_{00}(0, 0)$ of (3.30) at $F_{00}(0, 0)$ is determined in a neighborhood of $\epsilon$, and therefore mathematical expression for $M^cF_{00}(0, 0)$ is

where the computation yields

Finally, (3.30) restricted to $M^cF_{00}(0, 0)$ is

where

So, following discriminatory quantities are non-zero for existence of the flip bifurcation:

The simplification yields

and

where the involved quantities $A = a_1h^3r_2^3(1-m)^2+a_2\left(-8+h r_2\left(6-2 h r_2+b_1h(h r_2-2)\right)\right)$ and $B = a_2(h r_2-2)^2\left(a_2b_1+a_1r_2(1-m)^2\right)\left(a_1h r_2^2(1-m)^2(h r_2-4)+a_2b_1(h r_2-2)^2\right)$. From (3.38) if one gets $\jmath_2\not = 0$ as $(h, r_1, r_2, b_1, a_1, a_2, m)\in \mathcal{F}|_{E_{HP}^+\left(\frac{r_1a_2}{a_2b_1+a_1r_2(1-m)^2}, \frac{r_1r_2(1-m)}{a_2b_1+a_1r_2(1-m)^2}\right)}$ then at equilibrium $E_{HP}^+\left(\frac{r_1a_2}{a_2b_1+a_1r_2(1-m)^2}, \frac{r_1r_2(1-m)}{a_2b_1+a_1r_2(1-m)^2}\right)$ discrete model (1.3) undergoes flip bifurcation. Additionally, period-$2$ points from $E_{HP}^+\left(\frac{r_1a_2}{a_2b_1+a_1r_2(1-m)^2}, \frac{r_1r_2(1-m)}{a_2b_1+a_1r_2(1-m)^2}\right)$ are stable (respectively unstable) if $\jmath_2 > 0$ (respectively $\jmath_2 < 0)$.

4.   Chaos control

In the literature, there are many techniques to control chaos for the discrete-time models like state feedback control, pole placement and hybrid control methods ^{[33,34,35,36,37]}. In our understudied discrete-time model (1.3), we will use the state feedback control strategy that stabilizes the chaotic orbits at an unstable equilibrium point $E_{HP}^+\left(\frac{r_1a_2}{a_2b_1+a_1r_2(1-m)^2}, \frac{r_1r_2(1-m)}{a_2b_1+a_1r_2(1-m)^2}\right)$. It is important here to mention that in this method, the chaotic discrete-time model (1.3) is transformed into a piecewise linear system to attain an optimal controller that minimizes the upper bound, and then solving the optimization problem under certain constraints. So, on adding control force $U_t$ to model (1.3), it becomes

where $U_t = -k_1\left(H_t-H^*\right)-k_2\left(P_t-P^*\right)$ and $k_{1, 2}$ are control parameters. Now for (4.1), the variational matrix $V^C|_{E_{HP}^+\left(\frac{r_1a_2}{a_2b_1+a_1r_2(1-m)^2}, \frac{r_1r_2(1-m)}{a_2b_1+a_1r_2(1-m)^2}\right)}$ is

where

Now if $\lambda_{1, 2}$ are roots of characteristic equation of $V^C|_{E_{HP}(H, P)}$ then

and

Since marginal stability determined from the conditions $\lambda_1 = \pm 1, \ \lambda_1\lambda_2 = 1$ which further implies the fact that $|\lambda_{1, 2}| < 1$. If $\lambda_1\lambda_2 = 1$ then from (4.5) we get

If $\lambda_1 = 1$ then from (4.4) and (4.5) we get

Finally, if $\lambda_1 = -1$ then from (4.4) and (4.5) we get

Therefore, from (4.6)–(4.8) lines $L_1, \ L_2$ and $L_3$ in $(k_1, k_2)$-plane gives the triangular region that further gives $|\lambda_{1, 2}| < 1$.

5.   Numerical simulations

In this section, we will give some numerical simulation to verify theoretical results.

5.1.   Flip bifurcation at $E_{H0}(\frac{r_1}{b_1}, 0)$

If $h = 1.4, \ b_1 = 0.18, \ a_1 = 0.3, \ m = 0.8, \ r_2 = 0.27, \ a_2 = 0.14$ then from (2.7) one gets: $r_1 = 1.4285714285714286$. From theoretical point of view, $E_{H0}\left(\frac{r_1}{b_1}, 0\right)$ of discrete model (1.3) undergoes a flip bifurcation if $r_1 = 1.4285714285714286$. So, if $r_1 = {r_1}^\ast = \frac{2}{h} = 1.4285714285714286$ and $H = H^\ast = \frac{r_1}{b_1} = 7.936507936507937$ then from (3.4)–(3.6) the computation yields:

and

Equations (5.1)–(5.3) indicate that non-degenerate conditions hold, and so at $E_{H0}\left(\frac{r_1}{b_1}, 0\right) = E_{H0}\left(7.936507936507937, 0\right)$ discrete model (1.3) undergoes flip bifurcation. In addition, the simple calculation also yields $\Omega_1 = \left(\frac{\partial^2 f}{\partial H^2}\frac{\partial f}{\partial {r_1}}+2\frac{\partial^2 f}{\partial H\partial r_1}\right)\biggr|_{{r_1}^\ast = 1.4285714285714286, H^\ast = 7.936507936507937} = -2.8$ and $\Omega_2 = \left(\frac{\partial^3 f}{\partial H^3}+3\left(\frac{\partial^2 f}{\partial H^2}\right)^2\right)\biggr|_{{r_1}^\ast = 1.4285714285714286, \ H^\ast = 7.936507936507937} = 0.7620480000000001$. Finally, $\Omega_1\Omega_2 = -2.1337344 < 0$ which shows that model (1.3) undergoes supercritical flip bifurcation. Hence Maximum Lyapunov exponents (M. L. E.) and flip bifurcation diagram at $E_{H0}(\frac{r_1}{b_1}, 0)$ are drawn in Figure 1.

5.2.   Hopf bifurcation at $E_{HP}^+\left(\frac{r_1a_2}{a_2b_1+a_1r_2(1-m)^2}, \frac{r_1r_2(1-m)}{a_2b_1+a_1r_2(1-m)^2}\right)$

If $h = 1.39$, $b_1 = 0.56$, $a_1 = 0.61$, $m = 0.48$, $r_2 = 1.05$, $a_2 = 1.65$ then from (2.16) we get $r_1 = 1.7008190945069108$, and so $E_{HP}^+\left(\frac{r_1a_2}{a_2b_1+a_1r_2(1-m)^2}, \frac{r_1r_2(1-m)}{a_2b_1+a_1r_2(1-m)^2}\right)$ of discrete model (1.3) is a stable (respectively, an unstable) focus if $0 < r_1 < 1.7008190945069108$ (respectively $r_1 > 1.7008190945069108$). For this if $r_1 = 1.686 < 1.7008190945069108$ then Figure 2a indicates that $E_{HP}^+\left(\frac{r_1a_2}{a_2b_1+a_1r_2(1-m)^2}, \frac{r_1r_2(1-m)}{a_2b_1+a_1r_2(1-m)^2}\right) = E_{HP}^+\left(2.5354742181672614, 0.8390114685571666\right)$ of discrete model (1.3) is a stable focus, that means that all orbits goes to the interior fixed point $E_{HP}^+\left(2.5354742181672614, 0.8390114685571666\right)$ and additionally Figure 2b–2h also indicates same nature of solution if $r_1$ respectively are $r_1 = 1.62$, $1.58$, $1.66$, $1.695$, $1.698$, $1.0$, $1.6 < 1.7008190945069108$. Furthermore, if $r_1 = 1.72 > 1.7008190945069108$ then Figure 3a indicates that $E_{HP}^+\left(2.586604777726981, 0.8559310355387465\right)$ of discrete model (1.3) changes the nature of solution and as a result stable curves take place. Hereafter it is shown numerically that under consideration model undergoes supercritical hopf bifurcation if $r_1 = 1.72 > 1.7008190945069108$, i.e., $\ell < 0$. Therefore, if $r_1 = 1.72$ then $\frac{d|\lambda_{1, 2}|}{d \epsilon}|_{\epsilon = 0} = 0.4255699494665087 > 0$, and additionally from (3.9) and (3.25) we get

and

On substituting (5.4) and (5.5) in (3.23) we get $\ell = -1.1066864173641395 < 0$ which confirms the correctness of theoretical results and so supercritical hopf bifurcation takes place. Similarly Figure 3b–3h also shown same nature of solution if $r_1$ respectively are $r_1 = 1.735$, $1.75$, $1.776$, $1.797$, $1.83 > 1.7008190945069108$ and so for $r_1 = 1.735$, $1.75$, $1.776$, $1.797$, $1.83$, $1.8$, $2.1 > 1.7008190945069108$ model (1.3) undergoes supercritical hopf bifurcation with $\ell < 0$ (see Table 1). The M. L. E. and bifurcation diagrams are plotted in Figure 4.

5.3.   Flip bifurcation at $E_{HP}^+\left(\frac{r_1a_2}{a_2b_1+a_1r_2(1-m)^2}, \frac{r_1r_2(1-m)}{a_2b_1+a_1r_2(1-m)^2}\right)$

If $h = 1.4, b_1 = 0.18, a_1 = 0.3, m = 0.8, r_2 = 0.27, a_2 = 0.14$ then from non-hyperbolic condition (2.20) we get $r_1 = 1.6620447594316734$. Theoretically, if $r_1 = 1.6620447594316734$ then at interior fixed point $E_{HP}^+\left(\frac{r_1a_2}{a_2b_1+a_1r_2(1-m)^2}, \frac{r_1r_2(1-m)}{a_2b_1+a_1r_2(1-m)^2}\right)$, discrete model (1.3) undergoes a flip bifurcation, i.e., if $r_1 = 1.6620447594316734$ then from (3.37) one gets $\jmath_1 = 1.455433013797801\not = 0$. Further from (3.38) we get $\jmath_2 = 0.025512710056662495 > 0$ which represent that stable period-2 points bifurcate from the interior fixed point $E^+_{HP}\left(\frac{r_1a_2}{a_2b_1+a_1r_2(1-m)^2}, \frac{r_1r_2(1-m)}{a_2b_1+a_1r_2(1-m)^2}\right)$, and hence M. L. E. and flip bifurcation diagram are drawn in Figure 5.

5.4.   Simulation for correctness of theoretical results in Section 4

If $h = 0.4, b_1 = 0.58, a_1 = 0.03, m = 0.15, r_2 = 1.1, a_2 = 0.14, r_1 = 0.04$ then from (4.6)–(4.8) we get

and

The lines (5.6)–(5.8) define a triangular region that gives $|\lambda_{1, 2}| < 1$ (See Figure 6). Finally, $t$ vs $H_t$ and $P_t$ have drawn for (4.1) with respectively $k_{1, 2} = -0.15083845871908988, 0.02073349564264731$, which predict that unstable trajectories are stabilized (See Figure 7).

6.   Influence of prey refuge

In this section, the following two cases are to be considered:

Case Ⅰ: Prey densities increase due to the influence of prey refuge

For this, from $P$-component of interior fixed point $E_{HP}^+\left(\frac{r_1a_2}{a_2b_1+a_1r_2(1-m)^2}, \frac{r_1r_2(1-m)}{a_2b_1+a_1r_2(1-m)^2}\right)$, one can observe that if $m\in(0, 1)$, then following inequality holds obviously:

From (6.1), the $H$-component of interior fixed point $E_{HP}^+\left(\frac{r_1a_2}{a_2b_1+a_1r_2(1-m)^2}, \frac{r_1r_2(1-m)}{a_2b_1+a_1r_2(1-m)^2}\right)$ give

which implies that for fixed refuge, the prey refuge can increase the prey densities. Furthermore, $\frac{d}{dm}\left(\frac{r_1a_2}{a_2b_1+a_1r_2(1-m)^2}\right) = \frac{2r_1r_2a_1a_2}{\left(a_2b_1+a_1r_2\left(1-m\right)^2\right)^2} > 0\ \forall \ m\in(0, 1)$. This shows the fact that $\frac{r_1a_2}{a_2b_1+a_1r_2(1-m)^2}$ is strictly increasing function of $m$, that is, increasing the amount of refuge results the increase of prey densities.

Case Ⅱ: Predator densities decreases due to the influence of prey refuge

Again from the $P$-component of interior fixed point $E_{HP}^+\left(\frac{r_1a_2}{a_2b_1+a_1r_2(1-m)^2}, \frac{r_1r_2(1-m)}{a_2b_1+a_1r_2(1-m)^2}\right)$, one has $\frac{d}{dm}\left(\frac{r_1r_2(1-m)}{a_2b_1+a_1r_2(1-m)^2}\right) = \frac{r_1r_2\left(a_1r_2(1-m)^2-a_2b_1\right)^2} {\left(a_2b_1+a_1r_2\left(1-m\right)^2\right)^2}$. Additionally if $a_1r_2-a_2b_1\le0$, that is, $a_1r_2\le a_2b_1$ then $\frac{d}{dm}\left(\frac{r_1r_2(1-m)}{a_2b_1+a_1r_2(1-m)^2}\right) = \frac{r_1r_2\left(a_1r_2(1-m)^2-a_2b_1\right)^2} {\left(a_2b_1+a_1r_2\left(1-m\right)^2\right)^2} < 0\ \forall\ m\in(0, 1)$. This implies that $\frac{r_1r_2(1-m)}{a_2b_1+a_1r_2(1-m)^2}$ is strictly non-increasing function of $m$, that is, predator densities decreases due to the influence of prey refuge. Moreover $\frac{r_1r_2(1-m)}{a_2b_1+a_1r_2(1-m)^2}$ has maximum value $\frac{r_1r_2}{a_2b_1+a_1r_2}$ at $m = 0$.

7.   Conclusions

In this paper, we have investigated local behavior at fixed points, chaos and bifurcation of a discrete time model (1.3). More precisely, it is shown that $\forall \ r_1, \ r_2, \ a_1, \ a_2, \ b_1, \ m, \ h > 0$, model (1.3) has boundary fixed point $E_{H0}\left(\frac{r_1}{b_1}, 0\right)$ and if $m < 1$ then it has interior fixed point $E_{HP}^+\left(\frac{r_1a_2}{a_2b_1+a_1r_2(1-m)^2}, \frac{r_1r_2(1-m)}{a_2b_1+a_1r_2(1-m)^2}\right)$. Further at $E_{H0}\left(\frac{r_1}{b_1}, 0\right)$ and $E_{HP}^+\left(\frac{r_1a_2}{a_2b_1+a_1r_2(1-m)^2}, \frac{r_1r_2(1-m)}{a_2b_1+a_1r_2(1-m)^2}\right)$ the local dynamical characteristics have been studied, and proved that $E_{H0}\left(\frac{r_1}{b_1}, 0\right)$ of discrete model (1.3) is never sink; source if $r_1 > \frac{2}{h}$; saddle if $0 < r_1 < \frac{2}{h}$ and non-hyperbolic if $r_1 = \frac{2}{h}$. Moreover interior fixed point $E_{HP}^+\left(\frac{r_1a_2}{a_2b_1+a_1r_2(1-m)^2}, \frac{r_1r_2(1-m)}{a_2b_1+a_1r_2(1-m)^2}\right)$ of discrete mathematical model (1.3) is a stable focus if $0 < r_1 < \frac{r_2a_2b_1+{r_2}^2a_1(1-m)^2}{h a_2b_1r_2-a_2b_1+h a_1{r_2}^2(1-m)^2}$ with (2.14) holds; an unstable focus if $r_1 > \frac{r_2a_2b_1+{r_2}^2a_1(1-m)^2}{h a_2b_1r_2-a_2b_1+h a_1{r_2}^2(1-m)^2}$; non-hyperbolic if $r_1 = \frac{r_2a_2b_1+{r_2}^2a_1(1-m)^2}{h a_2b_1r_2-a_2b_1+h a_1{r_2}^2(1-m)^2}$; stable node if $0 < \frac{2(h r_2-2)\left(a_2b_1+a_1r_2(1-m)^2\right)}{h^2r_2\left(a_2b_1+a_1r_2(1-m)^2\right)-2 a_2b_1h}$ with (2.18) holds; an unstable node if $r_1 > \frac{2(h r_2-2)\left(a_2b_1+a_1r_2(1-m)^2\right)}{h^2r_2\left(a_2b_1+a_1r_2(1-m)^2\right)-2 a_2b_1h}$ and non-hyperbolic if $r_1 = \frac{2(h r_2-2)\left(a_2b_1+a_1r_2(1-m)^2\right)}{h^2r_2\left(a_2b_1+a_1r_2(1-m)^2\right)-2 a_2b_1h}$. We have also explored existence of bifurcation scenarios at fixed points, and proved that flip bifurcation exists at $E_{H0}\left(\frac{r_1}{b_1}, 0\right)$ if $(h, r_1, r_2, b_1, a_1, a_2, m)\in \mathcal{F}|_{E_{H0}\left(\frac{r_1}{b_1}, 0\right)} = \left\{\left(h, r_1, r_2, b_1, a_1, a_2, m\right), \ r_1 = \frac{2}{h}\right\}$. It is proved that at $E_{HP}^+\left(\frac{r_1a_2}{a_2b_1+a_1r_2(1-m)^2}, \frac{r_1r_2(1-m)}{a_2b_1+a_1r_2(1-m)^2}\right)$ model (1.3) undergoes hopf bifurcation if $(h, r_1, r_2, b_1, a_1, a_2, m)\in \mathcal{N}|_{E_{HP}^+\left(\frac{r_1a_2}{a_2b_1+a_1r_2(1-m)^2}, \frac{r_1r_2(1-m)}{a_2b_1+a_1r_2(1-m)^2}\right)}$ $= \left\{\left(h, r_1, r_2, b_1, a_1, a_2, m\right): \Delta < 0 \ {\mathrm{and}} \ r_1 = \frac{r_2a_2b_1+{r_2}^2a_1(1-m)^2}{h a_2b_1r_2-a_2b_1+h a_1{r_2}^2(1-m)^2}\right\}$ and flip bifurcation if $(h, r_1, r_2, b_1, a_1, a_2, m)$ $\in\mathcal{F}|_{E_{HP}^+\left(\frac{r_1a_2}{a_2b_1+a_1r_2(1-m)^2}, \frac{r_1r_2(1-m)}{a_2b_1+a_1r_2(1-m)^2}\right)}$ $= \left\{\left(h, r_1, r_2, b_1, a_1, a_2, m\right): \Delta > 0 \ {\mathrm{and}} \ r_1 = \frac{2(h r_2-2)\left(a_2b_1+a_1r_2(1-m)^2\right)}{h^2r_2\left(a_2b_1+a_1r_2(1-m)^2\right)-2a_2b_1h}\right\}$. By state feedback control strategy, chaos at $E_{HP}^+\left(\frac{r_1a_2}{a_2b_1+a_1r_2(1-m)^2}, \frac{r_1r_2(1-m)}{a_2b_1+a_1r_2(1-m)^2}\right)$ in discrete model (1.3) is also investigated. Next numerically verified theoretical results. Our numerical simulation reveals that if parameter crosses $(h, r_1, r_2, b_1, a_1, a_2, m)\in \mathcal{N}|_{E_{HP}^+\left(\frac{r_1a_2}{a_2b_1+a_1r_2(1-m)^2}, \frac{r_1r_2(1-m)}{a_2b_1+a_1r_2(1-m)^2}\right)}$ then model (1.3) undergoes the supercritical Neimark-Sacker bifurcation at $E_{HP}^+\left(\frac{r_1a_2}{a_2b_1+a_1r_2(1-m)^2}, \frac{r_1r_2(1-m)}{a_2b_1+a_1r_2(1-m)^2}\right)$, and so biologically this implies that there exists a periodic or quasi-periodic oscillation between prey and predator populations. Furthermore if $(h, r_1, r_2, b_1, a_1, a_2, m)$ $\in\mathcal{F}|_{E_{HP}^+\left(\frac{r_1a_2}{a_2b_1+a_1r_2(1-m)^2}, \frac{r_1r_2(1-m)}{a_2b_1+a_1r_2(1-m)^2}\right)}$ then model undergoes the flip bifurcation which indicates that the prey population will not remain steady, resulting in a biological imbalance in the ecosystem. Finally, we have also discussed the influence of prey refuge in the understudied discrete model.

Conflict of interest

The authors declare that they have no conflicts of interest regarding the publication of this paper.