Hopf Bifurcation Analysis of a Diffusive Nutrient–Phytoplankton Model with Time Delay

Abstract

In this paper, we studied a nutrient–phytoplankton model with time delay and diffusion term. We studied the Turing instability, local stability, and the existence of Hopf bifurcation. Some formulas are obtained to determine the direction of the bifurcation and the stability of periodic solutions by the central manifold theory and normal form method. Finally, we verify the above conclusion through numerical simulation.

1. Introduction

One of the most complex and difficult problems in water pollution treatment is the prevention and control of algal bloom. Due to the complexity of the pollution source and the difficulty factor of material removal, it takes a lot of energy, but it is not very effective. Therefore, scientists search for better methods to prevent and cure algal bloom, especially using mathematical models, in order to find reasonable prevention and cure measures [1,2,3,4,5,6,7]. In addition, many scholars further study the dynamics of the $N-P$ model by considering factors such as time delay and diffusion [8,9,10,11,12]. M. Rehim et al. studied a nutrient–plankton–zooplankton system with toxic phytoplankton and three delays, and showed the phenomenon of stability switches [8]. Y. Wang and W. Jiang considered a differential algebraic phytoplankton–zooplankton system with delay and harvesting, and indicated that the toxic liberation delay of phytoplankton may affect the stability of the coexisting equilibrium [10]. In particular, Huppert et al. [13] considered the following $N-P$ model

$$\begin{array}{c}\hfill \left\{\begin{array}{c}\frac{dN\left(t\right)}{dt}=a-bNP-eN,\hfill \\ \frac{dP\left(t\right)}{dt}=cNP-dP,\hfill \end{array}\right.\end{array}$$

(1)

where N is the nutrient level and P is the density of phytoplankton. a denotes the constant external nutrient inflow. b represents the maximal nutrient uptake rate. c represents the maximal conversion rate of nutrients into phytoplankton. d stands for the per capita mortality rate of phytoplankton. e denotes the per capita loss rate of nutrients. Relevant research work has analyzed the reasonable, deterministic, and empirical relationship between the abundance of toxin-producing phytoplankton and the diversity of plankton communities with large amounts of plankton but no toxins (called nontoxic plankton plants, NTP) [14]. In the case of toxic substances released by toxic phytoplankton (TPP), a simple model of vegetative phytoplankton was proposed and analyzed to understand the dynamic changes of the phenomenon of the seasonal mass reproductive cycle. The presence of chemical and toxic substances helps explain this phenomenon [15,16,17]. In [18], Chakraborty et al. considered the effect of toxins produced by toxic phytoplankton on the death of nontoxic phytoplankton, and produced the following equation

$$\begin{array}{c}\hfill \left\{\begin{array}{c}\frac{dN}{dt}=a-bNP-eN,\hfill \\ \frac{dP}{dt}=cNP-dP-\frac{\theta P^2}{{\mu }^2+P^2}.\hfill \end{array}\right.\end{array}$$

(2)

where $\theta$ is the release rate of toxic chemicals by the TPP population, and $\mu$ denotes the half-saturation constant.

Since the spatial distribution of nutrients and phytoplankton is inhomogeneous, there is diffusion. In addition, there is a time delay in the conversion from nutrients to phytoplankton. So, we incorporate reaction diffusion and time delay into the model (2), that is

$$\begin{array}{c}\hfill \left\{\begin{array}{c}\frac{\partial N}{\partial t}=d_1▵N+a-bNP-eN,\hfill \\ \frac{\partial P}{\partial t}=d_2▵P+cPN(t-\tau )-dP-\frac{\theta P^2}{{\mu }^2+P^2},\hfill \end{array}\right.\end{array}$$

(3)

where $d_1$ and $d_2$ are diffusion coefficients for N and P, respectively. ∆ is the Laplace operator. This is based on the assumption that the prey and predator are not stationary and will spread randomly. $\tau$ is the time delay that occurs for nutrients to be converted to phytoplankton. For analysis convenience, we have denoted

$$h=\frac{b}{a},s=\frac{e}{a},\alpha =\frac{d}{c},\beta =\frac{\theta }{c}.$$

The corresponding problem has the following form

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}\frac{\partial N}{\partial t}=d_1▵N+a(-hNP-Ns+1),\hfill & x\in (0,l\pi ),t>0,\hfill \\ \frac{\partial P}{\partial t}=d_2▵P+cP\left(-\alpha +N(t-\tau )-\frac{\beta P}{{\mu }^2+P^2}\right),\hfill & x\in (0,l\pi ),t>0,\hfill \\ N_x(0,t)=P_x(0,t)=0,N_x(l\pi ,t)=P_x(l\pi ,t)=0,\hfill & t>0,\hfill \\ N(x,t)=N_0(x,t)\ge 0,P(x,t)=P_0(x,t)\ge 0,\hfill & x\in [0,l\pi ],t\in [-\tau ,0].\hfill \end{array}\right.\end{array}$$

(4)

The content of the paper is arranged as follows. In Section 2, we study the stability and the existence of the Hopf bifurcation. In Section 3, we analyze the property of Hopf bifurcation. In Section 4, we provide a numerical simulation to verify the previous conclusions. Finally, we conclude this paper.

2. Stability Analysis

In [18], Chakraborty et al. studied the existence of equilibria. We cite the following result. The equilibrium points satisfy the following equation

$$\begin{array}{c}\hfill \left\{\begin{array}{c}1-hNP-sN=0,\hfill \\ -\alpha +N-\frac{\beta P}{{\mu }^2+P^2}=0,\hfill \end{array}\right.\end{array}$$

(5)

It can be calculated that trivial equilibrium $\left(\frac{1}{s},0\right)$ and interior equilibrium $(N_{*},P_{*})$, where $N_{*}=\frac{1}{h{P}_{*}+s},$ and $P_{*}$ is a root of the equation

$$h\alpha P^3+(h\beta +s\alpha -1)P^2+\left(h\alpha {\mu }^2+s\alpha \right)P-{\mu }^2(1-s\alpha )=0.$$

We provide the result from [18] as follows.

Lemma 1.

The existence of a positive equilibrium for the model (4) can be divided into the following cases.

(1)

If $1-s\alpha \le 0$, system (2) has no positive equilibrium.

(2)

If $0<1-s\alpha \le h\beta$, system (2) has one unique positive equilibrium.

(3)

If $1-s\alpha >h\beta$, then system (2) has either three or one positive equilibrium.

In what follows, we always assume that $0<1-s\alpha \le h\beta$, and we study the stability of problem (4) for $(N_{*},P_{*})$. Denote

$$N_1\left(t\right)=N(·,t)N_2\left(t\right)=P(·,t),N={(N_1,N_2)}^T,$$

$$X=C([0,l\pi ],{\mathbb{R}}^2),and{\mathcal{C}}_{\tau }:=C([-\tau ,0],X).$$

The linearized system of (4) at $(N^{*},P^{*})$ is

$$\dot{N}=(\mathbb{D}\mathsf{\Delta }+L)N,$$

(6)

where

$$\mathbb{D}=\left(\begin{array}{cc}{d}_1& 0\\ 0& d_2\end{array}\right),dom\left(\mathbb{D}\mathsf{\Delta }\right)=\{{(N,P)}^T:N,P\in C^2([0,l\pi ],{\mathbb{R}}^2),N_x,P_x=0,x=0,l\pi \},$$

and $L:{\mathcal{C}}_{\tau }\mapsto X$ is defined by

$$\begin{array}{c}\hfill L\varphi (·)=L_1\varphi (·)+L_2\varphi (·-\tau ),\end{array}$$

for $\varphi ={({\varphi }_1,{\varphi }_2)}^T\in {\mathcal{C}}_{\tau }$ with

$$L_1=\left(\begin{array}{cc}-aA& -aB\\ 0& D\end{array}\right),L_2=\left(\begin{array}{cc}0& 0\\ \tilde{C}& 0\end{array}\right),$$

$$A=h{P}_{*}+s,B=\frac{h}{h{P}_{*}+s},\tilde{C}=c{P}_{*},D=\frac{c\beta P_{*}\left(P_{*}^2-{\mu }^2\right)}{{\left({\mu }^2+P_{*}^2\right)}^2}.$$

(7)

The characteristic equations are

$${\lambda }^2+\lambda A_n+B_n+C_n{e}^{-\lambda \tau }=0,n\in {\mathcal{N}}_0,$$

(8)

where $A_n=(d_1+d_2)\frac{n^2}{l^2}+aA-D,$ $B_n=d_1{d}_2\frac{n^4}{l^4}+(aA{d}_2-D{d}_1)\frac{n^2}{l^2}-aAD,$ $C=aB\tilde{C}.$

2.1. Non-Delay Model

When $\tau =0$, the characteristic becomes

$${\lambda }^2-T_n\lambda +D_n=0,n\in {\mathcal{N}}_0,$$

(9)

where

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{T}_n=-(d_1+d_2)\frac{n^2}{l^2}+D-aA,\hfill \\ D_n=d_1{d}_2\frac{n^4}{l^4}+(aA{d}_2-D{d}_1)\frac{n^2}{l^2}+a(B\tilde{C}-AD),\hfill \end{array}\right.\end{array}$$

and the eigenvalues are given by

$${\lambda }_n=\frac{T_n\pm \sqrt{T_n^2-4D_n}}{2},n\in {\mathcal{N}}_0.$$

(10)

Then, make hypothesis

$$a>a_0:=\frac{D}{A},B\tilde{C}-AD>0.$$

(11)

Theorem 1.

Suppose $d_1=d_2=0$, $\tau =0$, and hypothesis (11) hold, then the equilibrium $(N_{*},P_{*})$ is locally asymptotically stable.

Proof. 

Suppose $d_1=d_2=0$, $\tau =0$, and hypothesis (11) hold, we can obtain $T_0<0$, $D_0>0$, so the real part of the roots of the characteristic equation is negative, then the equilibrium $(N_{*},P_{*})$ is locally asymptotically stable. □

It is calculated that the discriminant of $D_n$ is $\mathsf{\Gamma }=a^2{A}^2{d}_2^2+2a{d}_1{d}_2\left(AD-2B\tilde{C}\right)+D^2{d}_1^2$, and

$$a_{\pm }=\frac{d_1(2B\tilde{C}-AD)\pm d_1\sqrt{4B\tilde{C}(B\tilde{C}-AD)}}{A^2{d}_2}$$

(12)

$${\sigma }_{\pm }=\frac{-(aA{d}_2-D{d}_1)\pm \sqrt{{(aA{d}_2-D{d}_1)}^2-4d_1{d}_{2}a(B\tilde{C}-AD)}}{2d_1{d}_2}.$$

(13)

It is easy to verify that $a_{-}<\frac{d_1}{d_2}{a}_0<a_{+}$ under the hypothesis (11).

Theorem 2.

Suppose $d_1>0$, $d_2>0$, $\tau =0$, and hypothesis (11) hold. For the system (4), we have the following conclusion.

(1)

If $a\ge \frac{d_1}{d_2}{a}_0$, then the equilibrium $(N_{*},P_{*})$ is locally asymptotically stable.

(2)

If $a_{-}<a<\frac{d_1}{d_2}{a}_0$, then the equilibrium $(N_{*},P_{*})$ is locally asymptotically stable.

(3)

If $a_0<a<a_{-}$, and there is no $k\in \mathcal{N}$ such that $\frac{k^2}{l^2}\in ({\sigma }_{-},{\sigma }_{+})$, then the equilibrium $(N_{*},P_{*})$ is locally asymptotically stable.

(4)

If $a_0<a<a_{-}$, and there is a $k\in \mathcal{N}$ such that $\frac{k^2}{l^2}\in ({\sigma }_{-},{\sigma }_{+})$, then the equilibrium $(N_{*},P_{*})$ is Turing unstable.

Proof. 

We can obtain $T_n<0$ and $D_n>0$ for $a\ge \frac{d_1}{d_2}{a}_0$. It can be concluded that all the characteristic roots have a negative real part. Then, the equilibrium $(N_{*},P_{*})$ is locally asymptotically stable (so, statement (1) is true). In the same way, statements (1)–(3) are also correct. Suppose the conditions in statement (4) are true, then at least there is a positive real part of eigenvalue root. Then, the equilibrium $(N_{*},P_{*})$ is Turing unstable. □

2.2. Delay Model

Now, suppose $\tau >0$, one of the conditions (1)–(3) in Theorem 2 and hypothesis (11) hold. Assume $i\omega (\omega >0)$ is a solution of Equation (8), we can obtain

$$-w^2+i{A}_{n}w+B_n+Ccosw\tau -iCsinw\tau =0.$$

Then we have

$$\begin{array}{c}\hfill \left\{\begin{array}{c}-w^2+B_n+Ccosw\tau =0,\hfill \\ w{A}_n-Csinw\tau =0,\hfill \end{array}\right.\end{array}$$

(14)

which leads to

$$w^4+(A_n^2-2B_n)w^2+B_n^2-C^2=0.$$

(15)

Let $z={\omega }^2$, Equation (15) is

$$z^2+(A_n^2-2B_n)z+B_n^2-C^2=0.$$

(16)

By direct computation, we have

$$\begin{array}{cc}\hfill & A_n^2-2B_n={\left(aA+d_1\frac{n^2}{l^2}\right)}^2+{\left(D-d_2\frac{n^2}{l^2}\right)}^2>0,\hfill \\ \hfill & B_n+C=D_n>0,\hfill \\ \hfill & B_n-C=d_1{d}_2\frac{n^4}{l^4}+(aA{d}_2-D{d}_1)\frac{n^2}{l^2}-a(AD+B\tilde{C}).\hfill \end{array}$$

(17)

Define

$$\mathbb{M}=\{m\in {\mathbb{N}}_0|B_n-C<0withn=m\}.$$

(18)

Lemma 2.

Suppose one of the conditions (1)–(3) in Theorem 2 and hypothesis (11) hold. If $\mathbb{M}=\varnothing$, then Equation (16) has no positive root. If $\mathbb{M}\ne \varnothing$, then the equation has positive roots.

Proof. 

The roots of Equation (16) are

$$z_n^{\pm }=\frac{1}{2}[-(A_n^2-2B_n)\pm \sqrt{{(A_n^2-2B_n)}^2-4(B_n^2-C^2)}]$$

(19)

It is easy to verify that $z_n^{+}>0$ if and only if $n\in \mathbb{M}$, and $z_n^{-}$ is always negative or a non real number. □

Suppose one of the conditions (1)–(3) in Theorem 2 and hypothesis (11) hold, from Equation (14), we can obtain

$$sin\omega \tau =\frac{\omega A_n}{C}>0,cos\omega \tau =\frac{{\omega }^2-B_n}{C}.$$

For $n\in \mathbb{M}$, then Equation (8) has a pair of purely imaginary roots $\pm i{\omega }_n$ at ${\tau }_n^j,j\in {\mathbb{N}}_0$,

$${\omega }_n=\sqrt{z_n},{\tau }_n^j={\tau }_n^0+\frac{2j\pi }{{\omega }_n},{\tau }_n^0=\frac{1}{{\omega }_n}arccos\frac{{\omega }_n^2-B_n}{C}.$$

(20)

Lemma 3.

Suppose one of the conditions (1)–(3) in Theorem 2 and hypothesis (11) hold. Then

$$Re\left[\frac{d\lambda }{d\tau }\right]{|}_{\tau ={\tau }_n^j}>0forn\in \mathbb{M}andj\in {\mathbb{N}}_0.$$

Proof. 

From (8), we can obtain

$${\left(\frac{d\lambda }{d\tau }\right)}^{-1}=\frac{A_n+2\lambda }{\lambda C{e}^{-\lambda \tau }}-\frac{\tau }{\lambda }.$$

Then

$$Re{\left(\frac{d\lambda }{d\tau }\right)}_{\tau ={\tau }_n^j}^{-1}=\frac{A_n^2{w}^2+2w^2(w^2-B_n)}{A_n^2{w}^4+{(w^3-B_{n}w)}^2}=\frac{\sqrt{{(A_n^2-2B_n)}^2-4(B_n^2-C^2)}}{w^4+B_n^2+(A_n^2-2B_n)w^2}>0.$$

□

Denote $\mathcal{D}:=\{{\tau }_n^j:{\tau }_m^j\ne {\tau }_n^k,m\ne n,m,n\in \mathbb{M},j,k\in {\mathbb{N}}_0\},$ and ${\tau }_{*}=min\{\tau \in \mathcal{D}\}$.

Theorem 3.

For system (4), assume one of the conditions (1)–(3) in Theorem 2 and hypothesis (11) hold, then we have the following conclusion.

(1)

If $\mathbb{M}=\varnothing$, $(N_{*},P_{*})$ is locally asymptotically stable for $\tau \ge 0$.

(2)

If $\mathbb{M}\ne \varnothing$, $(N_{*},P_{*})$ is locally asymptotically stable for $\tau \in [0,{\tau }_{*})$ and unstable for $\tau >{\tau }_{*}$.

(3)

Hopf bifurcation occurs when $\tau ={\tau }_0^j$ ($j\in {\mathbb{N}}_0,n\in \mathbb{M}$).

Proof. 

If $\mathbb{M}=\varnothing$, then $B_n-C>0$ and $B_n^2-C^2>0,$ so Equation (16) has no positive root; then, the roots of Equation (8) all have negative real parts. Therefore, $(N_{*},P_{*})$ is locally asymptotically stable. Similarly, statement (2) is also correct. When $\tau ={\tau }_0^j$ ($j\in {\mathbb{N}}_0,n\in \mathbb{M}$) implying that $T_n=0$, then Hopf bifurcation occurs near $(N_{*},P_{*}).$ □

3. Property of Hopf Bifurcation

By the method [19,20,21], we study the property of Hopf bifurcation. For fixed $j\in {\mathbb{N}}_0$ and $n\in \mathbb{M}$, denote $\tilde{\tau }={\tau }_n^j$. Let $\tilde{N}(x,t)=N(x,\tau t)-N_{*}$, and $\tilde{P}(x,t)=P(x,\tau t)-P_{*}$. The system (4) (drop the tilde) is

$$\begin{array}{c}\hfill \left\{\begin{array}{c}\frac{\partial N}{\partial t}=\tau (d_1▵N+a(1-h(N+N_{*})(P+P_{*})-s(N+N_{*})),\hfill \\ \frac{\partial P}{\partial t}=\tau (d_2▵P+c(N(t-1)-N_{*})(P+P_{*})-\alpha (P+P_{*})-\frac{\beta {(P+P_{*})}^2}{{\mu }^2+{(P+P_{*})}^2}.\hfill \end{array}\right.\end{array}$$

(21)

Let

$$\tau =\tilde{\tau }+\mu ,N_1\left(t\right)=N(·,t),N_2\left(t\right)=P(·,t)andN={(N_1,N_2)}^T.$$

Then (21) is written as

$$\frac{dN\left(t\right)}{dt}=\tilde{\tau }D\mathsf{\Delta }N\left(t\right)+L_{\tilde{\tau }}\left(N_t\right)+F(N_t,\mu ),$$

(22)

where

$$\begin{array}{c}\hfill L_{\mu }\left(\varphi \right)=\mu \left(\begin{array}{c}-aA{\varphi }_1\left(0\right)-aB{\varphi }_2\left(0\right)\\ C{\varphi }_1(-1)+D{\varphi }_2\left(0\right)\end{array}\right),\end{array}$$

(23)

$$F(\varphi ,\mu )=\mu D\mathsf{\Delta }\varphi +L_{\mu }\left(\varphi \right)+f(\varphi ,\mu ),$$

(24)

with

$$\begin{array}{cc}\hfill f(\varphi ,\mu )& =(\tilde{\tau }+\mu ){(F_1(\varphi ,\mu ),F_2(\varphi ,\mu ))}^T,\hfill \\ \hfill F_1(\varphi ,\mu )& =a(1-h({\varphi }_1\left(0\right)+N_{*})({\varphi }_2\left(0\right)+P_{*})-s({\varphi }_1\left(0\right)+N_{*})+A{\varphi }_1\left(0\right)+B{\varphi }_2\left(0\right)),\hfill \\ \hfill F_2(\varphi ,\mu )& =c(({\varphi }_1(-1)+N_{*})({\varphi }_2\left(0\right)+P_{*})-\alpha ({\varphi }_2\left(0\right)+P_{*})-\frac{\beta \left(P_{*}+{\varphi }_2\left(0\right)\right){}^2}{{\mu }^2+\left(P_{*}+{\varphi }_2\left(0\right)\right){}^2}).\hfill \end{array}$$

Then

$$\frac{dN\left(t\right)}{dt}=\tilde{\tau }D\mathsf{\Delta }N\left(t\right)+L_{\tilde{\tau }}\left(N_t\right)$$

(25)

has characteristic roots ${\mathsf{\Lambda }}_n:=\{i{\omega }_n\tilde{\tau },-i{\omega }_n\tilde{\tau }\}$. Its linear functional differential equation is

$$\frac{dz\left(t\right)}{dt}=-\tilde{\tau }D\frac{n^2}{l^2}z\left(t\right)+L_{\tilde{\tau }}\left(z_t\right).$$

(26)

There exists a $2\times 2$ matrix function ${\eta }^n(\sigma ,\tilde{\tau })$ $-1\le \sigma \le 0$, such that

$$-\tilde{\tau }D\frac{n^2}{l^2}\varphi \left(0\right)+L_{\tilde{\tau }}\left(\varphi \right)={\int }_{-1}^{0}d{\eta }^n(\sigma ,\tau )\varphi \left(\sigma \right).$$

Choose

$${\eta }^n(\sigma ,\tau )=\left\{\begin{array}{cc}\tau E\hfill & \sigma =0,\hfill \\ 0\hfill & \sigma \in (-1,0),\hfill \\ -\tau F\hfill & \sigma =-1,\hfill \end{array}\right.E=\left(\begin{array}{cc}-aA-d_1\frac{n^2}{l^2}\hfill & -aB\hfill \\ 0\hfill & D-d_2\frac{n^2}{l^2}\hfill \end{array}\right),F=\left(\begin{array}{cc}0\hfill & 0\hfill \\ C\hfill & 0\hfill \end{array}\right).$$

(27)

Define

$$\begin{array}{cc}\hfill (\psi ,\varphi )& =\psi \left(0\right)\varphi \left(0\right)-{\int }_{-1}^0{\int }_{\xi =0}^{\sigma }\psi (\xi -\sigma )d{\eta }^n(\sigma ,\tilde{\tau })\varphi \left(\xi \right)d\xi \hfill \\ \hfill & =\psi \left(0\right)\varphi \left(0\right)+\tilde{\tau }{\int }_{-1}^0\psi (\xi +1)F\varphi \left(\xi \right)d\xi ,\hfill \end{array}$$

(28)

for $\varphi \in C([-1,0],{\mathbb{R}}^2)$, $\psi \in C([-1,0],{\mathbb{R}}^2)$. Choose $p_1\left(\theta \right)={(1,\xi )}^T{e}^{i{\omega }_n\tilde{\tau }\sigma }(\sigma \in [-1,0]),$$p_2\left(\sigma \right)=\overline{p_1\left(\sigma \right)}$ is a basis of $A\left(\tilde{\tau }\right)$ with ${\mathsf{\Lambda }}_n$ and $q_1\left(r\right)=(1,\eta )e^{-i{\omega }_n\tilde{\tau }r}(r\in [0,1]),q_2\left(r\right)=\overline{q_1\left(r\right)}$ is a basis of $A^{*}$ with ${\mathsf{\Lambda }}_n$, where

$$\xi =\frac{C{e}^{-i\tilde{\tau }{\omega }_n}}{i{\omega }_n+\frac{d_2{n}^2}{l^2}-D},\eta =\frac{B}{i{\omega }_n+D-\frac{d_2{n}^2}{l^2}}.$$

Let $\mathsf{\Phi }=({\mathsf{\Phi }}_1,{\mathsf{\Phi }}_2)$ and ${\mathsf{\Psi }}^{*}={({\mathsf{\Psi }}_1^{*},{\mathsf{\Psi }}_2^{*})}^T$ with

$${\mathsf{\Phi }}_1\left(\sigma \right)=\frac{p_1\left(\sigma \right)+p_2\left(\sigma \right)}{2},{\mathsf{\Phi }}_2\left(\sigma \right)=\frac{p_1\left(\sigma \right)-p_2\left(\sigma \right)}{2i},for\theta \in [-1,0].$$

In addition,

$${\mathsf{\Psi }}_1^{*}\left(r\right)=\frac{q_1\left(r\right)+q_2\left(r\right)}{2},{\mathsf{\Psi }}_2^{*}\left(r\right)=\frac{q_1\left(r\right)-q_2\left(r\right)}{2i},forr\in [0,1].$$

Then we can compute by (28)

$$D_1^{*}:=({\mathsf{\Psi }}_1^{*},{\mathsf{\Phi }}_1),D_2^{*}:=({\mathsf{\Psi }}_1^{*},{\mathsf{\Phi }}_2),D_3^{*}:=({\mathsf{\Psi }}_2^{*},{\mathsf{\Phi }}_1),D_4^{*}:=({\mathsf{\Psi }}_2^{*},{\mathsf{\Phi }}_2).$$

Define $({\mathsf{\Psi }}^{*},\mathsf{\Phi })=({\mathsf{\Psi }}_j^{*},{\mathsf{\Phi }}_k)=\left(\begin{array}{c}{D}_1^{*}{D}_2^{*}\\ D_3^{*}{D}_4^{*}\end{array}\right)$ and construct a new basis $\mathsf{\Psi }$ for $P_{*}$ by

$$\mathsf{\Psi }={({\mathsf{\Psi }}_1,{\mathsf{\Psi }}_2)}^T={({\mathsf{\Psi }}^{*},\mathsf{\Phi })}^{-1}{\mathsf{\Psi }}^{*}.$$

Then $(\mathsf{\Psi },\mathsf{\Phi })=I_2$. In addition, define $f_n:=({\beta }_n^1,{\beta }_n^2)$, where

$${\beta }_n^1=\left(\begin{array}{c}cos\frac{n}{l}x\\ 0\end{array}\right),{\beta }_n^2=\left(\begin{array}{c}0\\ cos\frac{n}{l}x\end{array}\right).$$

We also define

$$d·f_n=d_1{\beta }_n^1+d_2{\beta }_n^2,ford={(d_1,d_2)}^T\in {\mathcal{D}}_1,$$

and

$$<N,P>:=\frac{1}{l\pi }{\int }_0^{l\pi }{N}_1\overline{P_1}dx+\frac{1}{l\pi }{\int }_0^{l\pi }{N}_2\overline{P_2}dx$$

for $N=(N_1,N_2)$, $P=(P_1,P_2)$, $N,P\in X$, and $<\varphi ,f_0>={(<\varphi ,f_0^1>,<\varphi ,f_0^2>)}^T$. Equation (21) can be rewritten as

$$\begin{array}{c}\hfill \frac{dN\left(t\right)}{dt}=A_{\tilde{\tau }}{N}_t+R(N_t,\mu ),\end{array}$$

(29)

where

$$\begin{array}{c}\hfill R(N_t,\mu )=\left\{\begin{array}{cc}0,\hfill & \theta \in [-1,0),\hfill \\ F(N_t,\mu ),\hfill & \theta =0.\hfill \end{array}\right.\end{array}$$

(30)

By the decomposition of ${\mathcal{C}}_1$, the above solution is

$$\begin{array}{c}\hfill N_t=\mathsf{\Phi }\left(\begin{array}{c}{x}_1\hfill \\ x_2\hfill \end{array}\right)f_n+h(x_1,x_2,\mu ),\end{array}$$

(31)

with

$$\left(\begin{array}{c}{x}_1\hfill \\ x_2\hfill \end{array}\right)=(\mathsf{\Psi },<N_t,f_n>),$$

and

$$h(x_1,x_2,\mu )\in P_S{\mathcal{C}}_1,h(0,0,0)=0,Dh(0,0,0)=0.$$

The solution of (22) is

$$\begin{array}{c}\hfill N_t=\mathsf{\Phi }\left(\begin{array}{c}{x}_1\left(t\right)\hfill \\ x_2\left(t\right)\hfill \end{array}\right)f_n+h(x_1,x_2,0).\end{array}$$

(32)

Let $z=x_1-i{x}_2$, and notice that $p_1={\mathsf{\Phi }}_1+i{\mathsf{\Phi }}_2$. Then, we can obtain

$$\mathsf{\Phi }\left(\begin{array}{c}{x}_1\hfill \\ x_2\hfill \end{array}\right)f_n=({\mathsf{\Phi }}_1,{\mathsf{\Phi }}_2)\left(\begin{array}{c}\frac{z+\overline{z}}{2}\hfill \\ \frac{i(z-\overline{z})}{2}\hfill \end{array}\right)f_n=\frac{1}{2}(p_{1}z+\overline{p_1}\overline{z})f_n,$$

and

$$h(x_1,x_2,0)=h(\frac{z+\overline{z}}{2},\frac{i(z-\overline{z})}{2},0).$$

Hence, (32) is

$$\begin{array}{cc}\hfill N_t& =\frac{1}{2}(p_{1}z+\overline{p_1}\overline{z})f_n+h(\frac{z+\overline{z}}{2},\frac{i(z-\overline{z})}{2},0)\hfill \\ \hfill & =\frac{1}{2}(p_{1}z+\overline{p_1}\overline{z})f_n+W(z,\overline{z}),\hfill \end{array}$$

(33)

where

$$W(z,\overline{z})=h(\frac{z+\overline{z}}{2},\frac{i(z-\overline{z})}{2},0).$$

From [19], z meets

$$\begin{array}{c}\hfill \dot{z}=i{\omega }_n\tilde{\tau }z+g(z,\overline{z}),\end{array}$$

(34)

among them

$$\begin{array}{c}\hfill g(z,\overline{z})=({\mathsf{\Psi }}_1\left(0\right)-i{\mathsf{\Psi }}_2\left(0\right))<F(N_t,0),f_n>.\end{array}$$

(35)

Let

$$\begin{array}{c}\hfill W(z,\overline{z})=W_{20}\frac{z^2}{2}+W_{11}z\overline{z}+W_{02}\frac{{\overline{z}}^2}{2}+\cdots ,\end{array}$$

(36)

$$\begin{array}{c}\hfill g(z,\overline{z})=g_{20}\frac{z^2}{2}+g_{11}z\overline{z}+g_{02}\frac{{\overline{z}}^2}{2}+\cdots ;\end{array}$$

(37)

from Equations (33) and (36), we can obtain

$$N_t\left(0\right)=\frac{1}{2}(z+\overline{z})cos\left(\frac{nx}{l}\right)+W_{20}^{\left(1\right)}\left(0\right)\frac{z^2}{2}+W_{11}^{\left(1\right)}\left(0\right)z\overline{z}+W_{02}^{\left(1\right)}\left(0\right)\frac{{\overline{z}}^2}{2}+\cdots ,$$

$$P_t\left(0\right)=\frac{1}{2}(\xi +\overline{\xi }\overline{z})cos\left(\frac{nx}{l}\right)+W_{20}^{\left(2\right)}\left(0\right)\frac{z^2}{2}+W_{11}^{\left(2\right)}\left(0\right)z\overline{z}+W_{02}^{\left(2\right)}\left(0\right)\frac{{\overline{z}}^2}{2}+\cdots ,$$

$$N_t(-1)=\frac{1}{2}(z{e}^{-i{\omega }_n\tilde{\tau }}+\overline{z}{e}^{i{\omega }_n\tilde{\tau }})cos\left(\frac{nx}{l}\right)+W_{20}^{\left(1\right)}(-1)\frac{z^2}{2}+W_{11}^{\left(1\right)}(-1)z\overline{z}+W_{02}^{\left(1\right)}(-1)\frac{{\overline{z}}^2}{2}+\cdots ,$$

and

$$\begin{array}{cc}\hfill {\overline{F}}_1(N_t,0)& =\frac{1}{\tilde{\tau }}{F}_1={\alpha }_1{N}_t\left(0\right)P_t\left(0\right)+O\left(4\right),\hfill \end{array}$$

(38)

$$\begin{array}{cc}\hfill {\overline{F}}_2(N_t,0)& =\frac{1}{\tilde{\tau }}{F}_2=c{N}_t(-1)P_t\left(0\right)+{\beta }_1{P}_t^2\left(0\right)+{\beta }_2{P}_t^3\left(0\right)+O\left(4\right),\hfill \end{array}$$

(39)

with

$${\alpha }_1=-ah,{\beta }_1=\frac{\beta c\left(3{\mu }^2{P}_{*}^2-{\mu }^4\right)}{{\left({\mu }^2+P_{*}^2\right)}^2},{\beta }_2=-\frac{4\beta c{\mu }^2{P}_{*}\left(P_{*}^2-{\mu }^2\right)}{{\left({\mu }^2+P_{*}^2\right)}^4}.$$

Hence,

$$\begin{array}{cc}\hfill {\overline{F}}_1(N_t,0)=& {cos}^2\left(\frac{nx}{l}\right)(\frac{z^2}{2}{\chi }_{20}+z\overline{z}{\chi }_{11}+\frac{{\overline{z}}^2}{2}{\overline{\chi }}_{20})+\frac{z^2\overline{z}}{2}cos\frac{nx}{l}{\varsigma }_{11}+\frac{z^2\overline{z}}{2}{cos}^3\frac{nx}{l}{\varsigma }_{12}+\cdots ,\hfill \end{array}$$

(40)

$$\begin{array}{cc}\hfill {\overline{F}}_2(N_t,0)=& {cos}^2\left(\frac{nx}{l}\right)(\frac{z^2}{2}{\varrho }_{20}+z\overline{z}{\varrho }_{11}+\frac{{\overline{z}}^2}{2}{\overline{\varrho }}_{20})+\frac{z^2\overline{z}}{2}cos\frac{nx}{l}{\varsigma }_{21}+\frac{z^2\overline{z}}{2}{cos}^3\frac{nx}{l}{\varsigma }_{22}+\cdots ,\hfill \end{array}$$

(41)

$$\begin{array}{cc}\hfill <F(N_t,0),f_n>=& \tilde{\tau }({\overline{F}}_1(N_t,0)f_n^1+{\overline{F}}_2(N_t,0)f_n^2)\hfill \\ \hfill =& \frac{z^2}{2}\tilde{\tau }\left(\begin{array}{c}{\chi }_{20}\\ {\varsigma }_{20}\end{array}\right)\mathsf{\Gamma }+z\overline{z}\tilde{\tau }\left(\begin{array}{c}{\chi }_{11}\\ {\varsigma }_{11}\end{array}\right)\mathsf{\Gamma }+\frac{{\overline{z}}^2}{2}\tilde{\tau }\left(\begin{array}{c}{\overline{\chi }}_{20}\\ {\overline{\varsigma }}_{20}\end{array}\right)\mathsf{\Gamma }+\frac{z^2\overline{z}}{2}\tilde{\tau }\left(\begin{array}{c}{\kappa }_1\\ {\kappa }_2\end{array}\right)+\cdots \hfill \end{array}$$

(42)

with

$$\mathsf{\Gamma }=\frac{1}{l\pi }{\int }_0^{l\pi }{cos}^3\left(\frac{nx}{l}\right)dx,$$

$$\begin{array}{cc}\hfill & {\chi }_{20}=\frac{{\alpha }_1\xi }{2},{\chi }_{11}={\alpha }_1\left(\frac{\eta }{4}+\frac{\xi }{4}\right),{\varsigma }_{12}=0,\hfill \\ \hfill & {\varsigma }_{11}={\alpha }_1(\xi W_{11}^1\left(0\right)+W_{11}^2\left(0\right)+\frac{W_{20}^1\left(0\right)}{2}\eta +\frac{W_{20}^2\left(0\right)}{2}),\hfill \\ \hfill & {\varrho }_{20}=\frac{1}{2}\xi e^{-i\tau {\omega }_n}\left(c+{\beta }_1\xi e^{i\tau {\omega }_n}\right),\hfill \\ \hfill & {\varrho }_{11}=\frac{1}{4}{e}^{-i\tau {\omega }_n}\left(c\left(\eta +\xi e^{2i\tau {\omega }_n}\right)+2{\beta }_1\eta \xi e^{i\tau {\omega }_n}\right),\hfill \\ \hfill & {\varsigma }_{21}=W_{11}^2\left(0\right)\left(2{\beta }_1\xi +c{e}^{-i\tau {\omega }_n}\right)+W_{20}^2\left(0\right)\left({\beta }_1\eta +\frac{1}{2}c{e}^{i\tau {\omega }_n}\right)+c\xi W_{11}^1(-1)+\frac{c\eta W_{02}^1(-1)}{2},\hfill \\ \hfill & {\varsigma }_{22}=\frac{3}{4}{\beta }_2\eta {\xi }^2.\hfill \end{array}$$

(43)

$$\begin{array}{cc}\hfill & {\kappa }_1={\varsigma }_{11}\frac{1}{l\pi }{\int }_0^{l\pi }{cos}^2\left(\frac{nx}{l}\right)dx+{\varsigma }_{12}\frac{1}{l\pi }{\int }_0^{l\pi }{cos}^4\left(\frac{nx}{l}\right)dx,\hfill \end{array}$$

$$\begin{array}{cc}\hfill & {\kappa }_2={\varsigma }_{21}\frac{1}{l\pi }{\int }_0^{l\pi }{cos}^2\left(\frac{nx}{l}\right)dx\hfill \\ \hfill & +{\varsigma }_{22}\frac{1}{l\pi }{\int }_0^{l\pi }{cos}^4\left(\frac{nx}{l}\right)dx.\hfill \end{array}$$

Denote

$${\mathsf{\Psi }}_1\left(0\right)-i{\mathsf{\Psi }}_2\left(0\right):=\left({\gamma }_1{\gamma }_2\right).$$

Notice that

$$\frac{1}{l\pi }{\int }_0^{l\pi }{cos}^3\left(\frac{nx}{l}\right)dx=0,n\in \mathbb{N},$$

and we have

$$\begin{array}{cc}\hfill ({\mathsf{\Psi }}_1\left(0\right)& -i{\mathsf{\Psi }}_2\left(0\right))<F(N_t,0),f_n>=\hfill \\ \hfill & \frac{z^2}{2}({\gamma }_1{\chi }_{20}+{\gamma }_2{\varsigma }_{20})\mathsf{\Gamma }\tilde{\tau }+z\overline{z}({\gamma }_1{\chi }_{11}+{\gamma }_2{\varsigma }_{11})\mathsf{\Gamma }\tilde{\tau }+\frac{{\overline{z}}^2}{2}({\gamma }_1{\overline{\chi }}_{20}+{\gamma }_2{\overline{\varsigma }}_{20})\mathsf{\Gamma }\tilde{\tau }\hfill \\ \hfill & +\frac{z^2\overline{z}}{2}\tilde{\tau }[{\gamma }_1{\kappa }_1+{\gamma }_2{\kappa }_2]+\cdots ,\hfill \end{array}$$

(44)

From (35), (37) and (44), we have $g_{20}=g_{11}=g_{02}=0$, for $n\in \mathbb{N}$. If $n=0$, we obtain

$$g_{20}={\gamma }_1\tilde{\tau }{\chi }_{20}+{\gamma }_2\tilde{\tau }{\varrho }_{20},g_{11}={\gamma }_1\tilde{\tau }{\chi }_{11}+{\gamma }_2\tilde{\tau }{\varrho }_{11},g_{02}={\gamma }_1\tilde{\tau }{\overline{\chi }}_{20}+{\gamma }_2\tilde{\tau }{\overline{\varrho }}_{20}.$$

Furthermore, for $n\in {\mathbb{N}}_0$, $g_{21}=\tilde{\tau }({\gamma }_1{\kappa }_1+{\gamma }_2{\kappa }_2)$. Now, we compute $W_{20}\left(\theta \right)$ and $W_{11}\left(\theta \right)$ for $\theta \in [-1,0]$. From [19], we obtain

$$\dot{W}(z,\overline{z})=W_{20}z\dot{z}+W_{11}\dot{z}\overline{z}+W_{11}z\dot{\overline{z}}+W_{02}\overline{z}\dot{\overline{z}}+\cdots ,$$

$$A_{\tilde{\tau }}W(z,\overline{z})=A_{\tilde{\tau }}{W}_{20}\frac{z^2}{2}+A_{\tilde{\tau }}{W}_{11}z\overline{z}+A_{\tilde{\tau }}{W}_{02}\frac{{\overline{z}}^2}{2}+\cdots ,$$

$$\dot{W}(z,\overline{z})=A_{\tilde{\tau }}W+H(z,\overline{z}),$$

where

$$\begin{array}{cc}\hfill H(z,\overline{z})& =H_{20}\frac{z^2}{2}+W_{11}z\overline{z}+H_{02}\frac{{\overline{z}}^2}{2}+\cdots \hfill \\ \hfill & =X_{0}F(N_t,0)-\mathsf{\Phi }(\mathsf{\Psi },<X_{0}F(N_t,0),f_n>·f_n).\hfill \end{array}$$

(45)

Hence, we have

$$\begin{array}{c}\hfill (2i{\omega }_n\tilde{\tau }-A_{\tilde{\tau }})W_{20}=H_{20},-A_{\tilde{\tau }}{W}_{11}=H_{11},(-2i{\omega }_n\tilde{\tau }-A_{\tilde{\tau }})W_{02}=H_{02},\end{array}$$

(46)

that is,

$$\begin{array}{c}\hfill W_{20}={(2i{\omega }_n\tilde{\tau }-A_{\tilde{\tau }})}^{-1}{H}_{20},W_{11}=-A_{\tilde{\tau }}^{-1}{H}_{11},W_{02}={(-2i{\omega }_n\tilde{\tau }-A_{\tilde{\tau }})}^{-1}{H}_{02}.\end{array}$$

(47)

By (44), we have

$$\begin{array}{cc}\hfill H(z,& \overline{z})=-\mathsf{\Phi }\left(0\right)\mathsf{\Psi }\left(0\right)<F(N_t,0),f_n>·f_n\hfill \\ \hfill =& -(\frac{p_1\left(\theta \right)+p_2\left(\theta \right)}{2},\frac{p_1\left(\theta \right)-p_2\left(\theta \right)}{2i})\left(\begin{array}{c}{\mathsf{\Phi }}_1\left(0\right)\\ {\mathsf{\Phi }}_2\left(0\right)\end{array}\right)<F(N_t,0),f_n>·f_n\hfill \\ \hfill =& -\frac{1}{2}[p_1\left(\theta \right)({\mathsf{\Phi }}_1\left(0\right)-i{\mathsf{\Phi }}_2\left(0\right))+p_2\left(\theta \right)({\mathsf{\Phi }}_1\left(0\right)+i{\mathsf{\Phi }}_2\left(0\right))]<F(N_t,0),f_n>·f_n\hfill \\ \hfill =& -\frac{1}{2}[(p_1\left(\theta \right)g_{20}+p_2\left(\theta \right){\overline{g}}_{02})\frac{z^2}{2}+(p_1\left(\theta \right)g_{11}+p_2\left(\theta \right){\overline{g}}_{11})z\overline{z}+(p_1\left(\theta \right)g_{02}+p_2\left(\theta \right){\overline{g}}_{20})\frac{{\overline{z}}^2}{2}]+\cdots .\hfill \end{array}$$

Therefore by (45), for $\theta \in [-1,0)$,

$$\begin{array}{c}\hfill H_{20}\left(\theta \right)=\left\{\begin{array}{cc}0\hfill & n\in \mathbb{N},\hfill \\ -\frac{1}{2}(p_1\left(\theta \right)g_{20}+p_2\left(\theta \right){\overline{g}}_{02})·f_0\hfill & n=0,\hfill \end{array}\right.\end{array}$$

$$\begin{array}{c}\hfill H_{11}\left(\theta \right)=\left\{\begin{array}{cc}0\hfill & n\in \mathbb{N},\hfill \\ -\frac{1}{2}(p_1\left(\theta \right)g_{11}+p_2\left(\theta \right){\overline{g}}_{11})·f_0\hfill & n=0,\hfill \end{array}\right.\end{array}$$

$$\begin{array}{c}\hfill H_{02}\left(\theta \right)=\left\{\begin{array}{cc}0\hfill & n\in \mathbb{N},\hfill \\ -\frac{1}{2}(p_1\left(\theta \right)g_{02}+p_2\left(\theta \right){\overline{g}}_{20})·f_0\hfill & n=0,\hfill \end{array}\right.\end{array}$$

and

$$H(z,\overline{z})\left(0\right)=F(N_t,0)-\mathsf{\Phi }(\mathsf{\Psi },<F(N_t,0),f_n>)·f_n,$$

where

$$\begin{array}{c}\hfill H_{20}\left(0\right)=\left\{\begin{array}{cc}\tilde{\tau }\left(\begin{array}{c}{\chi }_{20}\\ {\varrho }_{20}\end{array}\right){cos}^2\left(\frac{nx}{l}\right),\hfill & n\in \mathbb{N},\hfill \\ \tilde{\tau }\left(\begin{array}{c}{\chi }_{20}\\ {\varrho }_{20}\end{array}\right)-\frac{1}{2}(p_1\left(0\right)g_{20}+p_2\left(0\right){\overline{g}}_{02})·f_0,\hfill & n=0.\hfill \end{array}\right.\end{array}$$

(48)

$$\begin{array}{c}\hfill H_{11}\left(0\right)=\left\{\begin{array}{cc}\tilde{\tau }\left(\begin{array}{c}{\chi }_{11}\\ {\varrho }_{11}\end{array}\right){cos}^2\left(\frac{nx}{l}\right),\hfill & n\in \mathbb{N},\hfill \\ \tilde{\tau }\left(\begin{array}{c}{\chi }_{11}\\ {\varrho }_{11}\end{array}\right)-\frac{1}{2}(p_1\left(0\right)g_{11}+p_2\left(0\right){\overline{g}}_{11})·f_0,\hfill & n=0.\hfill \end{array}\right.\end{array}$$

(49)

Then, we can obtain

$${\dot{W}}_{20}=A_{\tilde{\tau }}{W}_{20}=2i{\omega }_n\tilde{\tau }{W}_{20}+\frac{1}{2}(p_1\left(\theta \right)g_{20}+p_2\left(\theta \right){\overline{g}}_{02})·f_n,-1\le \theta <0.$$

That is,

$$W_{20}\left(\theta \right)=\frac{i}{2i{\omega }_n\tilde{\tau }}(g_{20}{p}_1\left(\theta \right)+\frac{{\overline{g}}_{02}}{3}{p}_2\left(\theta \right))·f_n+E_1{e}^{2i{\omega }_n\tilde{\tau }\theta },$$

where

$$\begin{array}{c}\hfill E_1=\left\{\begin{array}{cc}{W}_{20}\left(0\right)\hfill & n\in \mathbb{N},\hfill \\ W_{20}\left(0\right)-\frac{i}{2i{\omega }_n\tilde{\tau }}(g_{20}{p}_1\left(\theta \right)+\frac{{\overline{g}}_{02}}{3}{p}_2\left(\theta \right))·f_0\hfill & n=0.\hfill \end{array}\right.\end{array}$$

Using the definition of $A_{\tilde{\tau }}$ and (46), we have that for $-1\le \theta <0$

$$\begin{array}{cc}\hfill & -(g_{20}{p}_1\left(0\right)+\frac{{\overline{g}}_{02}}{3}{p}_2\left(0\right))·f_0+2i{\omega }_n\tilde{\tau }{E}_1-A_{\tilde{\tau }}(\frac{i}{2{\omega }_n\tilde{\tau }}(g_{20}{p}_1\left(0\right)+\frac{{\overline{g}}_{02}}{3}{p}_2\left(0\right))·f_0)\hfill \\ \hfill & -A_{\tilde{\tau }}{E}_1-L_{\tilde{\tau }}(\frac{i}{2{\omega }_n\tilde{\tau }}(g_{20}{p}_1\left(0\right)+\frac{{\overline{g}}_{02}}{3}{p}_2\left(0\right))·f_n+E_1{e}^{2i{\omega }_n\tilde{\tau }\theta })\hfill \\ \hfill & =\tilde{\tau }\left(\begin{array}{c}{\chi }_{20}\\ {\varrho }_{20}\end{array}\right)-\frac{1}{2}(p_1\left(0\right)g_{20}+p_2\left(0\right){\overline{g}}_{02})·f_0.\hfill \end{array}$$

As

$$A_{\tilde{\tau }}{p}_1\left(0\right)+L_{\tilde{\tau }}(p_1·f_0)=i{\omega }_0{p}_1\left(0\right)·f_0,$$

and

$$A_{\tilde{\tau }}{p}_2\left(0\right)+L_{\tilde{\tau }}(p_2·f_0)=-i{\omega }_0{p}_2\left(0\right)·f_0,$$

we have

$$\begin{array}{c}\hfill 2i{\omega }_n{E}_1-A_{\tilde{\tau }}{E}_1-L_{\tilde{\tau }}{E}_1{e}^{2i{\omega }_n}=\tilde{\tau }\left(\begin{array}{c}{\chi }_{20}\\ {\varrho }_{20}\end{array}\right){cos}^2\left(\frac{nx}{l}\right),n\in {\mathbb{N}}_0.\end{array}$$

That is,

$$E_1=\tilde{\tau }E\left(\begin{array}{c}{\chi }_{20}\\ {\varrho }_{20}\end{array}\right){cos}^2\left(\frac{nx}{l}\right)$$

where

$$E={\left(\begin{array}{cc}2i{\omega }_n\tilde{\tau }+d_1\frac{n^2}{l^2}+aA& aB\\ -C{e}^{-2i{\omega }_n\tilde{\tau }}& -D+2i{\omega }_n\tilde{\tau }+d_2\frac{n^2}{l^2}\end{array}\right)}^{-1}.$$

Similarly, we have

$$-{\dot{W}}_{11}=\frac{i}{2{\omega }_n\tilde{\tau }}(p_1\left(\theta \right)g_{11}+p_2\left(\theta \right){\overline{g}}_{11})·f_n,-1\le \theta <0.$$

That is,

$$W_{11}\left(\theta \right)=\frac{i}{2i{\omega }_n\tilde{\tau }}(p_1\left(\theta \right){\overline{g}}_{11}-p_1\left(\theta \right)g_{11})+E_2.$$

Then,

$$E_2=\tilde{\tau }{E}^{*}\left(\begin{array}{c}{\chi }_{11}\\ {\varrho }_{11}\end{array}\right){cos}^2\left(\frac{nx}{l}\right),$$

where

$$E^{*}={\left(\begin{array}{cc}{d}_1\frac{n^2}{l^2}+aA& aB\\ -C& -D+d_2\frac{n^2}{l^2}\end{array}\right)}^{-1}.$$

Therefore, we have

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}{c}_1\left(0\right)=\frac{i}{2{\omega }_n\tilde{\tau }}(g_{20}{g}_{11}-2{\left|g_{11}\right|}^2-\frac{|g_{02}{|}^2}{3})+\frac{1}{2}{g}_{21},\hfill & {\mu }_2=-\frac{Re\left(c_1\left(0\right)\right)}{Re\left({\lambda }^{\prime }\left({\tau }_n^j\right)\right)},\hfill \\ T_2=-\frac{1}{{\omega }_n\tilde{\tau }}[Im\left(c_1\left(0\right)\right)+{\mu }_{2}Im\left({\lambda }^{\prime }\left({\tau }_n^j\right)\right)],\hfill & {\beta }_2=2Re\left(c_1\left(0\right)\right).\hfill \end{array}\right.\end{array}$$

(50)

By [19], we have the following theorem.

Theorem 4.

For any critical value ${\tau }_n^j$, we have the Hopf bifurcation is forward (${\mu }_2>0$) or backward (${\mu }_2<0$). The bifurcating periodic solutions are orbitally asymptotically stable (${\beta }_2<0$) or unstable (${\beta }_2>0$). The period increases ($T_2>0$) or decreases ($T_2<0$).

4. Numerical Simulations

In order to verify the previous conclusion, we provide some numerical simulations by Matlab. In particular, the numerical simulation of the systems with $\tau =0$ is implemented by the pdepe function in Matlab, and $\tau >0$ is implemented by the finite-difference methods. Choose the following parameters.

$$d_1=2,h=1.667,s=0.267,\alpha =0.1,a=1,\beta =0.48,\mu =0.18,c=0.5,l=4.$$

(51)

By direct computation, we have $(N_{*},P_{*})\approx (1.2130,0.3344)$ is the unique positive equilibrium, and $A\approx 0.8244$, $B\approx 2.0220$, $C\approx 0.1672$, $D\approx 0.3064$, $BC-AD\approx 0.0854>0$, and $a_0\approx 0.3717$. Hence, hypothesis (11) holds. Now, we give the curves of $a_{-}$ and $\frac{d_1}{d_2}{a}_0$ with the predator’s diffusion coefficient $d_2$ (Figure 1). We can see that $a_0<a<a_{-}$ holds when $d_2<d_2^{*}$, then the Turing instability of $(N_{*},P_{*})$ may occur. When $d_2>d_2^{*}$, then $a>a_{-}$ holds, which implies $(N_{*},P_{*})$ is locally asymptotically stable. Choose $d_2=0.1$, we have $a_{-}=2.4603$, ${\sigma }_{-}=0.1723$, ${\sigma }_{+}=2.4800$, and $k\in \{2,3,4,5,6\}$ such that $\frac{k^2}{l^2}\in ({\sigma }_{-},{\sigma }_{+})$. Then $(N_{*},P_{*})$ is Turing unstable (Figure 2).

We choose $d_2=0.4$, and change the parameter $\beta$, which represents the release rate of toxic chemicals by the TPP population. The bifurcation diagram of system $\left(4\right)$ with parameter $\beta$ is given in Figure 3. We can see that the increasing parameter $\beta$ is not beneficial to the stability of $(N_{*},P_{*})$ initially. However, when $\beta$ crosses some critical value, increasing parameter $\beta$ is of benefit to the stability of $(N_{*},P_{*})$. In particular, when the parameter $\beta$ is sufficiently large, $(N_{*},P_{*})$ will always be stable.

If we choose $\beta =0.48$, we have ${\tau }_{*}={\tau }_0^0\approx 1.5710$ and ${\omega }_0\approx 0.2460$. Then $(P_{*},N_{*})$ is locally asymptotically stable for $\tau \in [0,{\tau }_{*})$ (Figure 4), and Hopf bifurcation occurs when $\tau ={\tau }_{*}$. We obtain

$${\mu }_2\approx 0.5391>0,{\beta }_2\approx -0.1217<0,and{T}_2\approx 12.3699>0.$$

Hence, the stable bifurcating periodic solutions exist for $\tau >{\tau }_{*}$ (Figure 5). However, if we choose $\beta =0.6$ and $\tau =2.3$, $(P_{*},N_{*})$ is locally asymptotically stable (Figure 6).

5. Conclusions

Diffusion and time delay was incorporated into a nutrient–phytoplankton model. The instability and Hopf bifurcation induced by the time delay was studied. Through the central manifold theory and normal form method, some parameters were given to determine the property of bifurcating periodic solutions. The results indicate diffusion may induce Turing unstable. The release rate $\beta$ of toxic chemicals by the TPP population has a stabilizing and destabilizing effect on the stability of the positive equilibrium. In addition, the time delay can also affect the stability of the positive equilibrium, and it can induce periodic oscillation of prey and predator population density.