The Dynamical Properties of a Class of Discrete Smith Diffusion Model

Yu-Ming Yan, Min-Hong Yang, Ying Lin, Lu-Di Chen, Xiao-Liang Zhou^*

School of Mathematics and Statistics, Lingnan Normal University, Zhanjiang, China.

DOI: 10.4236/ijmnta.2024.131001                

Abstract

In this paper, the dynamical properties of Smith type diffusion model with Dirichlet boundary conditions are studied. The properties of hyperbolic fixed points and non-hyperbolic fixed points of the model are analyzed. By using the central manifold theorem, the bifurcation phenomenon of the model is studied. The results show that flip, transcritical, pitchfork and Fold-flip bifurcations exist at non-hyperbolic fixed points.

Keywords

Hyperbolicfixed Point, Non-Hyperbolic Fixed Point, Central Manifold Theorem, Bifurcation

1. Introduction

For many years, the relationship between growth rate and density for population has been the object of discussion and experiment, and researchers have proposed many mathematical models to describe this relationship. In 1963, Smith [1] used continuous culture techniques to study a relatively complex adaptation of the metazoan daphnia population. According to the experimental data, the relationship between the specific growth rate and population density of daphnia was observed to be inconsistent with the prediction of the Logistic differential equation, that is, the relationship is not linear. Smith proposed the following model to describe what he observed

$\frac{\text{d}x}{\text{d}t}=xH\left(x\right)$ , (1)

where $H\left(x\right)$ is average growth rate of population and

$H\left(x\right)=r\left(\frac{k-x}{k+\frac{r}{c}x}\right)$

k is the mass of x at saturation, r is the growth rate of x without food restriction, c is the replacement rate of x at saturation. In 2020, Liu et al. based on Smith model (1) and considering linear yield, proposed the following model (see [2] )

$\frac{\text{d}x}{\text{d}t}=x\left(\frac{r\left(k-x\right)}{k+bx}-h\right)$ , (2)

where h (>0) is a linear harvest rate and b = r/c. They studied phenomena such as harvest behavior and equilibrium bifurcation caused by Allee effects.

Since diffusion can significantly change the spatial distribution of species, in recent years, many researchers have paid attention to the model with diffusion effects (see [3] [4] [5] ). For example, Meng et al. discussed the discrete population diffusion model

${\mu }_i^{t+1}=\frac{p{\mu }_i^t}{1+q{\mu }_i^t}+d{\Delta }^2{\mu }_{i-1}^t$

under the Dirichlet boundary condition (see [4] ) of

${\mu }_0^t=0={\mu }_{m+1}^t$ (3)

and studied the existence of steady-state solutions and bifurcation of the model.

In this paper, inspired by papers [2] and [4] , we will consider the diffusion model of model (2)

${\mu }_i^{t+1}={\mu }_i^t\left(\frac{r\left(k-{\mu }_i^t\right)}{k+b{\mu }_i^t}-h\right)+d{\Delta }^2{\mu }_{i-1}^t$ (4)

where $d>0$ is the diffusion coefficient, $t\in Z^{+}=\left\{0,1,2,\cdots \right\}$ , ${\Delta }^2$ is the second order difference operator, and ${\Delta }^2{\mu }_{i-1}^t={\mu }_{i+1}^t-2{\mu }_i^t+{\mu }_{i-1}^t$ , $i\in \left\{1,2,\cdots ,m\right\}$ . Supposing that the Dirichlet boundary condition (3) is as ${\mu }_0^t=0={\mu }_3^t$ and binary variables as

$x_t={\mu }_2^t$ , $y_t={\mu }_1^t$ , $t\in Z^{+}=\left\{0,1,2,\cdots \right\}$ ,

then, we rewrite model (4) (with replacing n with t) as the following two-dimensional discrete time element population model

$\{\begin{array}{l}{x}_{n+1}=x_n+\delta x_n\left(\frac{r\left(k-x_n\right)}{k+b{x}_n}-h\right)+d\left(-2x_n+y_n\right),\\ y_{n+1}=y_n+\delta y_n\left(\frac{r\left(k-y_n\right)}{k+b{y}_n}-h\right)+d\left(x_n-2y_n\right),\end{array}$ (5)

where $n\in Z^{+}=\left\{0,1,2,\cdots \right\}$ and $\delta >0$ . Obviously, the system (5) has three fixed points: zero fixed point $E_0=\left(0,0\right)$ and two positive fixed points

$E_1=\left(1,2-\frac{\delta }{d}\left(\frac{r\left(k-1\right)}{k+b}-h\right)\right),E_2=\left(2-\frac{\delta }{d}\left(\frac{r\left(k-1\right)}{k+b}-h\right),1\right)$ .

In the following, the hyperbolic and non-hyperbolic properties of three fixed points are analyzed, and the flip, transcritical and pitchfork bifurcations generated at E₀ are also studied.

2. Analysis of Hyperbolic and Non-hyperbolic Cases

We write the model (5) as a plane map on ${ℝ}^2$ as follows

$\left(\begin{array}{c}x\\ y\end{array}\right)\mapsto \left(\begin{array}{c}\left(1-\delta h-2d\right)x+\delta x\frac{r\left(k-x\right)}{k+bx}+dy\\ \left(1-\delta h-2d\right)y+\delta y\frac{r\left(k-y\right)}{k+by}+dx\end{array}\right)$ . (6)

For any fixed point $\left(x_i,y_i\right)$ , the Jacobian matrix of (6) is

$J\left(E_0\right)=\left(\begin{array}{cc}{A}_0& d\\ d& A_0\end{array}\right)$

where $A_0=1-\delta \left(h-r\right)-2d$ , and

$J\left(E_i\right)=\left(\begin{array}{cc}e+A_i& d\\ d& e+B_i\end{array}\right)$

where

$e=1-\delta h-2d$ , $A_i=\frac{\delta r\left(k-x_i\right)}{k+b{x}_i}-\frac{\delta rk\left(b+1\right)x_i}{{\left(k+b{x}_i\right)}^2}$ , $B_i=\frac{\delta r\left(k-y_i\right)}{k+b{y}_i}-\frac{\delta rk\left(b+1\right)y_i}{{\left(k+b{y}_i\right)}^2}$

for $i=1,2$ .

2.1. The Property of the Fixed Point E₀

In order to discuss the properties of E₀, we define

$\beta =h-r,{\mathcal{L}}_1=\left\{d|d=2-\delta \beta \right\},{\mathcal{L}}_2=\left\{d|d=\frac{1}{3}\left(2-\delta \beta \right)\right\},$

${\mathcal{L}}_3=\left\{d|d=-\delta \beta ,\beta <0\right\},{\mathcal{L}}_4=\left\{d|d=-\frac{1}{3}\delta \beta ,\beta <0\right\},$

${\mathcal{D}}_1=\left\{\left(d,\beta \right)|0<d<-\frac{1}{3}\delta \beta ,\beta <0\right\},$

${\mathcal{D}}_2\left\{\left(d,\beta \right)|0<\frac{1}{3}\left(2-\delta \beta \right)<d<-\delta \beta ,\beta <0\right\},$

${\mathcal{D}}_3=\left\{\left(d,\beta \right)|d>2-\delta \beta >0\right\},$

${\mathcal{K}}_1=\left\{\left(d,\beta \right)|0<-\delta \beta <d<\frac{1}{3}\left(2-\delta \beta \right),\beta \le 0\right\},$

${\mathcal{K}}_2=\left\{\left(d,\beta \right)|0<d<\min \left\{-\delta \beta ,\frac{1}{3}\left(2-\delta \beta \right)\right\}\right\},$

${\mathcal{K}}_3=\left\{\left(d,\beta \right)|\max \left\{0,-\delta \beta ,\frac{1}{3}\left(2-\delta \beta \right)\right\}<d<2-\delta \beta \right\}.$

Theorem 2.1 The fixed point E₀ has the following properties:

1) when $\left(d,\beta \right)\in {\mathcal{D}}_{1,2,3}$ , E₀ isan unstable node;

2) when $\left(d,\beta \right)\in {\mathcal{K}}_1$ , E₀ isa stable node;

3) when $\left(d,\beta \right)\in {\mathcal{K}}_{2,3}$ , E₀ is a saddle;

4) when $\left(d,\beta \right)\in {\mathcal{L}}_i\left(i=1,2,3,4\right)$ , E₀ is non-hyperbolic.

Proof. The characteristic equation corresponding to the Jacobian matrix of fixed point E₀ is

$F\left(\lambda \right)={\left[\lambda -\left(1-\delta \beta -2d\right)\right]}^2-d^2$ .

Thus, we get

${\lambda }_1=1-\delta \beta -d,{\lambda }_2=1-\delta \beta -3d,$

and

${\lambda }_1>{\lambda }_2\text{.}$

Therefore, we easily prove that

1) when $\left(d,\beta \right)\in {\mathcal{D}}_{1,2,3}$ , we have $\left|{\lambda }_1\right|>1$ , $\left|{\lambda }_2\right|>1$ . Thus E₀ is an unstable node;

2) when $\left(d,\beta \right)\in {\mathcal{K}}_1$ , we have $\left|{\lambda }_1\right|<1$ , $\left|{\lambda }_2\right|<1$ . Thus E₀ is a stable node;

3) when $\left(d,\beta \right)\in {\mathcal{K}}_{2,3}$ , we have $\left|{\lambda }_1\right|>1$ , $\left|{\lambda }_2\right|<1$ or $\left|{\lambda }_1\right|<1$ , $\left|{\lambda }_2\right|>1$ . Thus E₀ is astable node saddle point;

4) when $\left(d,\beta \right)\in {\mathcal{L}}_i\left(i=1,2,3,4\right)$ , we have $\left|{\lambda }_1\right|=1$ or $\left|{\lambda }_2\right|=1$ . Thus E₀ is a non-hyperbolic fixed point.

2.2. The Properties of Fixed Points $E_i\left(i=1,2\right)$

In order to discuss the stabilitiesof fixed point $E_i\left(i=1,2\right)$ , we give the following lemma of which the proof is obvious and will be omitted.

Lemma 2.1 Suppose $F\left(\lambda \right)={\lambda }^2+Q\lambda +S$ and ${\lambda }_1,{\lambda }_2$ are two roots of $F\left(\lambda \right)=0$ . Then

1) when $-2<Q<2$ , $Q^2\ge 4S$ and $F\left(-1\right)>0$ , then $0-1<{\lambda }_{1,2}<1$ ;

2) when $Q<-2$ and $Q^2\ge 4S$ , then ${\lambda }_{1,2}>1$ ;

3) when $Q>2$ and $Q^2\ge 4S$ , then ${\lambda }_{1,2}<-1$ ;

4) when $F\left(1\right)>0$ and $F\left(-1\right)<0$ , then $\left|{\lambda }_1\right|>1$ and $\left|{\lambda }_2\right|<1$ ;

5) when $F\left(1\right)=0$ and $Q\ne 0,-2$ , then ${\lambda }_1=1$ and ${\lambda }_2\ne 1$ ;

6) when $F\left(1\right)=0$ and $Q=0$ or −2, then ${\lambda }_{1,2}=1$ ;

7) when $F\left(-1\right)=0$ and $Q\ne 2$ , then ${\lambda }_1=-1$ and ${\lambda }_2\ne -1$ ;

8) when $F\left(-1\right)=0$ and $Q=2$ , then ${\lambda }_{1,2}=-1$ ;

9) when $0<S<1$ and $Q^2<4S$ , then ${\lambda }_{1,2}$ is a pair of complex roots and $\left|{\lambda }_{1,2}\right|<1$ ;

10) when $S>1$ and $Q^2<4S$ , then ${\lambda }_{1,2}$ is a pair of complex roots and $\left|{\lambda }_{1,2}\right|>1$ ;

11) when $S=1$ and $-2<Q<2$ , then ${\lambda }_{1,2}$ is a pair of complex roots and $\left|{\lambda }_{1,2}\right|=1$ .

Let

${\mu }_i=-\frac{1}{\delta r}\left(A_i+B_i\right)$ , ${\omega }_i=\left(e+A_i\right)\left(e+B_i\right)$ ,

for $i=1,2$ .

Theorem 2.2 All of topology types for the fixed points $E_i\left(i=1,2\right)$ are listed in Table 1.

Proof. Thecharacteristic equation corresponding to the Jacobian matrix of fixed point $E_{1,2}$ is

$F\left(\lambda \right)={\lambda }^2+Q\lambda +S=0$ ,

where

$Q=-\text{tr}\left(J\left(E_i\right)\right)=\delta r{\mu }_i-2e$ , $S=\det \left(J\left(E_i\right)\right)={\omega }_i-d^2$ ,

for $i=1,2$ .

Therefore, from lemma 2.1, we easily obtain the conclusion of Table 1 of the theorem.

3. Analysis of Bifurcation at Fixed Point E₀

Theorem 3.1 When $\left(d,\beta \right)\in {\mathcal{L}}_1=\left\{d|d=2-\delta \beta \right\}$ , system(5) undergoes a supercritical flip bifurcation at point $E_0\left(0,0\right)$ , i.e., as $\delta \beta$ goes from less than 2 to more than 2, system (5) bifurcates out a stable period-2 orbit at the fixed point $E_0\left(0,0\right)$ .

Proof. When $\left(d,\beta \right)\in {\mathcal{L}}_1=\left\{d|d=2-\delta \beta \right\}$ , characteristic values ${\lambda }_1=-1$ , ${\lambda }_2=-5+2\delta \beta$ and $\left|{\lambda }_2\right|\ne 1$ . Denote $\epsilon =d+\delta \beta -2$ and select $\epsilon$ as the bifurcation parameter. Therefore the mapping (6) becomes

Table 1. The topology types for the fixed points $E_i\left(i=1,2\right)$

$\left(\begin{array}{c}x\\ y\end{array}\right)\mapsto \left(\begin{array}{c}-\left(\epsilon +d+\delta r+1\right)x+\delta x\frac{r\left(k-x\right)}{k+bx}+dy\\ -\left(\epsilon +d+\delta r+1\right)y+\delta y\frac{r\left(k-y\right)}{k+by}+dx\end{array}\right)$ . (7)

The Taylor expansion of mapping (7) is

$\left(\begin{array}{c}x\\ y\end{array}\right)\mapsto \left(\begin{array}{cc}-\epsilon -d-1& d\\ d& -\epsilon -d-1\end{array}\right)\left(\begin{array}{c}x\\ y\end{array}\right)+\left(\begin{array}{c}\frac{-\delta r\left(b+1\right)}{k}{x}^2\\ \frac{-\delta r\left(b+1\right)}{k}{y}^2\end{array}\right)+O\left({\left|\left(x,y\right)\right|}^3\right)$ . (8)

We get the Jacobian matrix

$J\left(E_0\right)=\left(\begin{array}{cc}-\epsilon -d-1& d\\ d& -\epsilon -d-1\end{array}\right)$ ,

characteristic values

${\lambda }_1=-\epsilon -1,{\lambda }_2=-\epsilon -2d-1$

and the corresponding eigenvectors

${\left(1,1\right)}^{\tau },{\left(1,-1\right)}^{\tau }$ . (9)

From the eigenvectors the following transformation is obtained

$\left(\begin{array}{c}x\\ y\end{array}\right)=\left(\begin{array}{cc}1& 1\\ 1& -1\end{array}\right)(\; u\; v\; )$

and it can transform (8) into the following mapping (ε is treatedas an independent variable)

$\left(\begin{array}{c}u\\ v\\ \epsilon \end{array}\right)\mapsto \left(\begin{array}{ccc}-\epsilon -1& 0& 0\\ 0& -\epsilon -2d-1& 0\\ 0& 0& 1\end{array}\right)\left(\begin{array}{c}u\\ v\\ \epsilon \end{array}\right)+\left(\begin{array}{c}\frac{-\delta r\left(b+1\right)}{k}\left(u^2+v^2\right)\\ \frac{-2\delta r\left(b+1\right)}{k}uv\\ 0\end{array}\right)+O\left({\left|\left(u,v\right)\right|}^3\right)$ . (10)

From center manifold theorem, the stability of mapping (10) in small neighborhood of $\left(u,v\right)=\left(0,0\right)$ can be determined by single parameter mapping, which satisfies

$W_{loc}\left(0,0\right)=\left\{\left(u,v,\epsilon \right)\in R^3|v=h\left(u,\epsilon \right),h\left(0,0\right)=0,Dh\left(0,0\right)=0\right\}$ . (11)

Suppose the central manifold is as follows

$v=h\left(u,\epsilon \right)=A{u}^2+Bu\epsilon +C{\epsilon }^2+O\left({\left|\left(u,\epsilon \right)\right|}^3\right)$ . (12)

Then

$\begin{array}{l}h\left(\left(-\epsilon -1\right)u+\frac{-\delta r\left(b+1\right)}{k}\left(u^2+h^2\left(u,\epsilon \right)\right),\epsilon \right)\\ =\left(-\epsilon -2d-1\right)h\left(u,\epsilon \right)+\frac{-2\delta r\left(b+1\right)}{k}uh\left(u,\epsilon \right)+O\left({\left|\left(u,\epsilon \right)\right|}^3\right)\end{array}$ (13)

From (10), (12) and (14), we get

$A=B=C=0$ .

Then

$v=h\left(u,\epsilon \right)=O\left({\left|\left(u,\epsilon \right)\right|}^3\right)$ . (14)

The mapping (10) restricted to the central manifold (14) is

$T:u\mapsto f\left(u,\epsilon \right)=\left(-\epsilon -1\right)u+\frac{-\delta r\left(b+1\right)}{k}{u}^2+O\left({\left|\left(u,\epsilon \right)\right|}^3\right)$ . (15)

Thus

$f\left(0,0\right)=0,\frac{\partial f}{\partial u}\left(0,0\right)=-1,{\left(\frac{\partial f}{\partial \epsilon }\frac{{\partial }^{2}f}{\partial u^2}+2\frac{{\partial }^{2}f}{\partial u\partial \epsilon }\right)}_{\left(0,0\right)}=-2\ne 0,$

${\left[\frac{1}{2}{\left(\frac{{\partial }^{2}f}{\partial u^2}\right)}^2+\frac{1}{3}\left(\frac{{\partial }^{3}f}{\partial u^3}\right)\right]}_{\left(0,0\right)}=\frac{2{\delta }^2{r}^2{\left(b+1\right)}^2}{k^2}\ne 0.$

Therefore, from [6] we know that system (5) undergoes a supercritical flip bifurcation at point $E_0\left(0,0\right)$ .

Theorem 3.2 When $\left(d,\beta \right)\in {\mathcal{L}}_2=\left\{d|d=\frac{1}{3}\left(2-\delta \beta \right)\right\}$ and $d\ne 1$ , $\delta \beta \ne -1$ , system(5) undergoes a subcritical flip bifurcation at point $E_0\left(0,0\right)$ .

Proof. When $\left(d,\beta \right)\in {\mathcal{L}}_2=\left\{d|d=\frac{1}{3}\left(2-\delta \beta \right)\right\}$ and $d\ne 1$ , $\delta \beta \ne -1$ , we have

${\lambda }_1=\frac{1}{3}\left(1-2\delta \beta \right),{\lambda }_2=-1$

and $\left|{\lambda }_1\right|\ne 1$ . Similar to Theorem 3.1, the proof of this theorem can be obtained and will be omitted.

Theorem 3.3 When $\left(d,\beta \right)\in {\mathcal{L}}_3=\left\{d|d=-\delta \beta ,\beta <0\right\}$ and $d\ne 1$ , system(5) undergoes atranscritical bifurcation at point $E_0\left(0,0\right)$ .

Proof. When $\left(d,\beta \right)\in {\mathcal{L}}_3=\left\{d|d=-\delta \beta ,\beta <0\right\}$ and $d\ne 1$ , we have

${\lambda }_1=1,\left|{\lambda }_2\right|=\left|1+2\delta \beta \right|\ne 1$ .

Set $d=-\delta \beta +\epsilon$ and chose ε as the bifurcation parameter, thus mapping (6) is written in the following form

$\left(\begin{array}{c}u\\ v\\ \epsilon \end{array}\right)\mapsto \left(\begin{array}{ccc}1-\epsilon & 0& 0\\ 0& 1-\epsilon -2d& 0\\ 0& 0& \epsilon \end{array}\right)\left(\begin{array}{c}u\\ v\\ \epsilon \end{array}\right)+\left(\begin{array}{c}\frac{-\delta r\left(b+1\right)}{k}\left(u^2+v^2\right)\\ \frac{-2\delta r\left(b+1\right)}{k}uv\\ 0\end{array}\right)$ . (16)

Similar to the proof of theorem 3.1, one-dimensional equations under the restriction of a central manifold is obtained

$T:u\mapsto f\left(u,\epsilon \right)=\left(1-\epsilon \right)u+\frac{-\delta r\left(b+1\right)}{k}{u}^2+O\left({\left|\left(u,\epsilon \right)\right|}^3\right)$ . (17)

We have

$f\left(0,0\right)=0,\frac{\partial f}{\partial u}\left(0,0\right)=1,\frac{\partial f}{\partial \epsilon }\left(0,0\right)=0,$

$\frac{{\partial }^{2}f}{\partial u^2}\left(0,0\right)=\frac{-2\delta r\left(b+1\right)}{k}\ne 0,\frac{{\partial }^{2}f}{\partial u\partial \epsilon }\left(0,0\right)=-1\ne 0.$

Therefore, from [7] we know that system (5) undergoes atranscritical bifurcation at point $E_0\left(0,0\right)$ .

Theorem 3.4 When $\left(d,\beta \right)\in {\mathcal{L}}_4=\left\{d|d=-\frac{1}{3}\delta \beta ,\beta <0\right\}$ , system(5) undergoes a pitchfork bifurcation at point E₀.

Proof. When $\left(d,\beta \right)\in {\mathcal{L}}_4=\left\{d|d=-\frac{1}{3}\delta \beta ,\beta <0\right\}$ , we have

${\lambda }_1=1-\frac{2}{3}\delta \beta ,{\lambda }_2=1$ .

Set $d=\frac{1}{3}\left(-\delta \beta +\epsilon \right)$ and chose ε as the bifurcation parameter, thus mapping (6) is written in the following form

$\left(\begin{array}{c}u\\ v\\ \epsilon \end{array}\right)\mapsto \left(\begin{array}{ccc}2d-\epsilon +1& 0& 0\\ 0& -\epsilon +1& 0\\ 0& 0& \epsilon \end{array}\right)\left(\begin{array}{c}u\\ v\\ \epsilon \end{array}\right)+\left(\begin{array}{c}\frac{-\delta r\left(b+1\right)}{k}\left(u^2+v^2\right)\\ \frac{-2\delta r\left(b+1\right)}{k}uv\\ 0\end{array}\right)$ . (18)

Similar to the proof of theorem 3.1, one-dimensional equations under the restriction of a central manifold is obtained

$T:u\mapsto f\left(u,\epsilon \right)=\left(-\epsilon +1\right)v+-\frac{{\delta }^2{r}^2{\left(b+1\right)}^2}{k^{2}d}{v}^3+O\left({\left|\left(u,\epsilon \right)\right|}^3\right)$ . (19)

Then we have

$f\left(0,0\right)=0,\frac{\partial f}{\partial u}\left(0,0\right)=1,\frac{\partial f}{\partial \epsilon }\left(0,0\right)=0,$

$\frac{{\partial }^{2}f}{\partial v^2}\left(0,0\right)=0,\frac{{\partial }^{2}f}{\partial u\partial \epsilon }\left(0,0\right)=-1\ne 0,\frac{{\partial }^{3}f}{\partial v^3}\left(0,0\right)=-\frac{6{\delta }^2{r}^2{\left(b+1\right)}^2}{k^{2}d}\ne 0.$

Therefore, from [8] we know that system (5) undergoes a pitchfork bifurcation at point $E_0\left(0,0\right)$ .

Theorem 3.5 When $d=1$ and $\delta \beta =-1$ , system (5) undergoes fold-flip bifurcation at point E₀.

Proof. When $d=1$ and $\delta \beta =-1$ , we have ${\lambda }_1=1$ , ${\lambda }_2=-1$ . Let

$q=d-1,r=\delta \beta +1$

and chose q and r as the bifurcation parameter. Then mapping (6) may be

$\left(\begin{array}{c}x\\ y\end{array}\right)\mapsto \left(\begin{array}{c}\left(-2q-r-\delta r\right)x+\delta x\frac{r\left(k-x\right)}{k+bx}+\left(q+1\right)y\\ \left(-2q-r-\delta r\right)y+\delta y\frac{r\left(k-y\right)}{k+by}+\left(q+1\right)x\end{array}\right)$ . (20)

The Taylor expansion of mapping (20) is

$\left(\begin{array}{c}x\\ y\end{array}\right)\mapsto \left(\begin{array}{cc}-2q-r& q+1\\ q+1& -2q-r\end{array}\right)\left(\begin{array}{c}x\\ y\end{array}\right)+\left(\begin{array}{c}\frac{-\delta r\left(b+1\right)}{k}{x}^2\\ \frac{-\delta r\left(b+1\right)}{k}{y}^2\end{array}\right)+O\left({\left|\left(x,y\right)\right|}^3\right)$ . (21)

We get the Jacobian matrix

$J\left(E_0\right)=\left(\begin{array}{cc}-2q-r& q+1\\ q+1& -2q-r\end{array}\right)$ .

The characteristic values of $J\left(E_0\right)$ are

${\lambda }_1=-q-r+1,{\lambda }_2=-3q-r-1$

and the corresponding eigenvectors are

$q_1={\left(1,1\right)}^{\text{T}},q_2={\left(-1,1\right)}^{\text{T}},p_1={\left(\frac{1}{2},\frac{1}{2}\right)}^{\text{T}},p_2={\left(-\frac{1}{2},\frac{1}{2}\right)}^{\text{T}}$

which satisfy

$\{\begin{array}{l}J\left(E_0\right)q_1={\lambda }_1{q}_1,J\left(E_0\right)q_2={\lambda }_2{q}_2,\\ J{\left(E_0\right)}^{\text{T}}{p}_1={\lambda }_1{p}_1,J{\left(E_0\right)}^{\text{T}}{p}_2={\lambda }_2{p}_2,\\ \langle p_1,q_1\rangle =\langle p_2,q_2\rangle =1,\langle p_1,q_2\rangle =\langle p_2,q_1\rangle =0,\end{array}$

where $\langle \cdot ,\cdot \rangle$ is scalar product. Therefore, any vector $X={\left(x,y\right)}^{\text{T}}$ can be uniquely expressed as

$X={\xi }_1{q}_1+{\xi }_2{q}_2$

where ${\xi }_1$ and ${\xi }_2$ can be calculated by the following equation

$\{\begin{array}{l}{\xi }_1=\langle p_1,x\rangle ,\\ {\xi }_2=\langle p_2,x\rangle .\end{array}$

Then the mapping (21) can be rewritten in the following form with the new coordinates ${\xi }_1$ and ${\xi }_2$

$\left(\begin{array}{c}{\xi }_1\\ {\xi }_2\end{array}\right)\mapsto \left(\begin{array}{cc}{\beta }_1+1& 0\\ 0& {\beta }_2-1\end{array}\right)\left(\begin{array}{c}{\xi }_1\\ {\xi }_2\end{array}\right)+\left(\begin{array}{c}\frac{1}{2!}{g}_{20}{\xi }_1{}^2+\frac{1}{2!}{g}_{02}{\xi }_2^2\\ h_{11}{\xi }_1{\xi }_2\end{array}\right)+O\left({\left|\left({\xi }_1,{\xi }_2\right)\right|}^3\right)$ , (22)

where

${\beta }_1=-q-r,{\beta }_2=-3q-r,$

$g_{20}=\frac{-2\delta r\left(b+1\right)}{k},g_{02}=\frac{-2\delta r\left(b+1\right)}{k},h_{11}=\frac{-\delta r\left(b+1\right)}{k}\ne 0,$

$g_{11}=h_{20}=h_{02}=g_{ij}=h_{ij}=0,\text{for}i+j=3.$

Then from the result of [9] we know that system (5) undergoes Fold-flip bifurcation at point E₀.

Acknowledgements

This work has been supported by Guangdong Basic and Applied Basic Research Foundation (Grant No. 2022A1515010964, 2022A1515010193), the Key Project of Science and Technology Innovation of Guangdong College Students (Grant No. pdjh2023b0325).

Conflicts of Interest

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

References

[1]	Smith, F.E. (1963) Population Dynamics in Daphnia Magna and a New Model for Population Growth. Ecology, 44, 651-663. https://doi.org/10.2307/1933011
[2]	Liu, Y., Zhang, T. and Liu, X. (2020) Investigating the Interactions Between Allee Effect and Harvesting Behavior of a Single Species Model: An Evolutionary Dynamics Approach. Physica A: Statistical Mechanics and Its Applications, 549, Article 124323. https://doi.org/10.1016/j.physa.2020.124323
[3]	Yu, X., Yuan, S. and Zhang, T. (2018) Persistence and Ergodicity of a Stochastic Single Species Model with Allee Effect under Regime Switching. Communications in Nonlinear Science & Numerical Simulations, 59, 359-374. https://doi.org/10.1016/j.cnsns.2017.11.028
[4]	Meng, L., Li, X. and Zhang, G. (2017) Simple Diffusion Can Support the Pitchfork, the Flip Bifurcations, and the Chaos. Communications in Nonlinear Science and Numerical Simulation, 53, 202-212. https://doi.org/10.1016/j.cnsns.2017.04.025
[5]	Li, M.S., Zhou, X.L. and Xu, J.M. (2020) Dynamic Properties of a Discrete Population Model with Diffusion. Advances in Difference Equations, 2020, Article No. 580. https://doi.org/10.1186/s13662-020-03033-w
[6]	Kuznetsov, Y.A. (1998) Elements of Applied Bifurcation Theory. Springer-Verlag, New York.
[7]	Wiggins, S. (2003) Introduction to Applied Nonlinear Dynamical Systems and Chaos. SpringerVerlag, New York.
[8]	Guckenheimer, H.J. (1983) Nonlinear Oscillations, Dynamical Systems and Bifurcations of Vector Fields. Springer, New York. https://doi.org/10.1007/978-1-4612-1140-2
[9]	Kuznetsov, Y.A., Meijer, H.G.E. and Veen, L. (2004) The Fold-Flip Bifurcation. International Journal of Bifurcation and Chaos, 14, 2253-2282. https://doi.org/10.1142/S0218127404010576