Analysis of a toxin producing phytoplankton–zooplankton interaction with Holling IV type scheme and time delay

Abstract

The present paper aims to study the interaction of toxin producing phytoplankton (TPP)–zooplankton (a prey–predator interaction) and its role in plankton ecology. The delay in the zooplankton predation is considered and its effect on the overall dynamic of phytoplankton–zooplankton interaction is studied. Moreover, Holling IV type response function is used for zooplankton grazing to account for the effect of toxication by the TPP population. It is shown that time delay can destabilise the given system and induce oscillation in population due to Hopf-bifurcation. Further stability of the bifurcating periodic solution is determined by using normal form theory and centre manifold arguments. Some numerical simulations are executed to validate the analytical findings.

Appendix

In order to compute the properties of the Hopf-bifurcation, we use the following transformation, $u_{1}=P-P_*$, $u_{2}=Z-Z_*$ and using the Taylor series expansion about $E_*=(P_*,Z_*)$, the given system (2.1) can be rewritten as,

$$\begin{aligned} {\left\{ \begin{array}{ll} \frac{{\hbox {d}}u_1}{{\hbox {d}}t}=a_{10}u_{1}(t)+a_{01}u_{2}(t-\tau )\\ \qquad \quad +{\displaystyle \sum _{i+j\ge 2}}a_{ij}{u^{i}}_{1}(t){u^{j}}_{2}(t-\tau )=F_{1}(u_1,u_2)\\ \frac{{\hbox {d}}u_2}{{\hbox {d}}t}=b_{10}u_{1}(t)+b_{01}u_{2}(t)\\ \qquad \quad +{\displaystyle \sum _{i+j\ge 2}}b_{ij}{u^{i}}_{1}(t){u^{j}}_{2}(t)=F_{2}(u_1,u_2) \end{array}\right. } \end{aligned}$$

(8.1)

where

$$\begin{aligned} a_{ij}={\displaystyle \frac{1}{!i !j}}{\displaystyle \frac{{\partial }^{i+j}F_{1}}{\partial P^{i}\partial Z^{j}(t-\tau )}}, \quad b_{ij}={\displaystyle \frac{1}{!i !j}}{\displaystyle \frac{{\partial }^{i+j}F_{2}}{\partial P^{i}\partial Z^{j}}}, \end{aligned}$$

The coefficients of the nonlinear terms are given by

$$\begin{aligned} a_{20}&= -2\delta _1-\beta Z_*\frac{\left( -\frac{6P^3_*}{i^2}-\frac{4P^2_*}{i}+\frac{2P_*b}{i}+2b\right) }{(\frac{P^2_*}{i}+P_*+b)^3} ,\\ a_{11}&= - \frac{\beta \left( -\frac{P^2_*}{i}+b\right) }{\left( \frac{P^2_*}{i}+P_*+b\right) ^2},\\ b_{20}\!&= \!\frac{\beta _1 Z\left[ \left( \!-\!2\frac{P_*}{i}\left( \frac{P^2_*}{i}\!+\!P_*\!+\!b\right) \right) \!-\!\left( \frac{P^2_*}{i}\!+\!b\right) ^2\left( \frac{2P^*}{i}\!+\!1\right) \right] }{\left( \frac{P^2_*}{i}\!+\!P_*\!+\!b\right) ^3},\\ b_{11}&= \frac{\beta _1\left( -\frac{P^2_*}{i}+b\right) }{\left( \frac{P^2_*}{i}+P_*+b\right) ^2} \end{aligned}$$

Let $\tau =\tau _{k}+\mu$ , $\bar{u}_{i}(t)=u_{i}(\tau t)$ and $u_{t}(\theta )=u(t+\theta )$ for $\theta \in [-1,0]$ and dropping the bars for simplification of notations, system (8.1) becomes a functional differential equation in $C=C([-1,0],\mathfrak {R}^2)$ as

$$\begin{aligned} \dot{u}(t)=L_{\mu }(u_{t})+f(\mu ,u_{t}) \end{aligned}$$

(8.2)

where $u(t)=(u_1(t),u_2(t))^{T}\in \mathfrak {R}^{2}$ and $L_{\mu }:C\rightarrow \mathfrak {R}^{3}, f:\mathfrak {R}\times C\rightarrow \mathfrak {R}^{2}$ are given, respectively, by

$$\begin{aligned} L_{\mu }(\phi )=(\tau _{k}+\mu )[A_{11}\phi (0)+A_{22}\phi (-1)] \end{aligned}$$

(8.3)

$$\begin{aligned} A_{11}=\left[ \begin{array}{ccc} a_{10} &{}\quad 0 \\ b_{10} &{}\quad b_{01} \end{array} \right] , A_{22}=\left[ \begin{array}{ccc} 0 &{}\quad a_{01}\\ 0 &{}\quad 0 \end{array} \right] \end{aligned}$$

and

$$\begin{aligned} f(\mu ,\phi )=(\tau _{k}+\mu )\left[ \begin{array}{cc} a_{20}\phi ^2_{1}(0)+a_{11}\phi _(0)\phi _{2}(-1)\\ b_{20}\phi ^2_{1}(0)+b_{11}\phi _{1}(0)\phi _{2}(0) \end{array}\right] \nonumber \\ \end{aligned}$$

(8.4)

By Riesz representation theorem, there exists a function $\phi (\theta ,\mu )$ of bounded variation for $\theta \in [-1,0]$ such that

$$\begin{aligned} L_{\mu }(\phi )=\int _{-1}^0{\hbox {d}}\eta (\theta ,\mu )\phi (\theta )\quad for \quad \phi \in C \end{aligned}$$

(8.5)

In fact, we can choose

$$\begin{aligned} \phi (\theta ,\mu )=(\tau _{k}+\mu )[A_{11}\delta (\theta )-A_{22}\delta (\theta +1)] \end{aligned}$$

(8.6)

where $\delta$ denote the Dirac delta function.

For $\phi \in C([-1,0],\mathfrak {R}^2)$, define

$$\begin{aligned} A(\mu )\phi (\theta )= {\left\{ \begin{array}{ll} {\displaystyle \frac{{\hbox {d}}\phi (\theta )}{{\hbox {d}}\theta }} &{} \theta \in [-1,0), \\ \int _{-1}^0{\hbox {d}}\varsigma (s,\mu )\phi (s) &{} \theta =0 \end{array}\right. } \end{aligned}$$

and

$$\begin{aligned} R(\mu )\phi (\theta )={\left\{ \begin{array}{ll} 0 &{} \theta \in [-1,0), \\ f(\mu ,\phi ) &{} \theta =0 \end{array}\right. } \end{aligned}$$

Then system (8.2) is equivalent to

$$\begin{aligned} \dot{u}(t)=A(\mu )u_{t}+R(\mu )u_{t}, \end{aligned}$$

(8.7)

For $\psi \in C^{1}([0,1],(\mathfrak {R}^{2})^*)$, define

$$\begin{aligned} A^{*}\psi (s) = {\left\{ \begin{array}{ll} -{\displaystyle \frac{{\hbox {d}}\psi (s)}{{\hbox {d}}s}} &{}\quad s\in (0,1], \\ \int _{-1}^0{\hbox {d}}\varsigma (t,0)\psi (-t) &{}\quad s=0 \end{array}\right. } \end{aligned}$$

and a bilinear inner product

$$\begin{aligned}&<\psi (s),\phi (\theta )>=\bar{\psi }(0)\phi (0)\nonumber \\&-\!\int _{-1}^0\int _{\xi =0}^\theta \bar{\psi }(\xi \!-\!\theta ){\hbox {d}}\eta (\theta )\phi (\xi ){\hbox {d}}\xi \nonumber \\ \end{aligned}$$

(8.8)

where $\varsigma (\theta )=\varsigma (\theta ,0)$. Then $A(0)$ and $A^{*}$ are adjoint operators. From the results of last section, we know that $\pm \iota \omega _{0}\tau _{k}$ are the eigenvalues of $A(0)$ and $\mp \iota \omega _{0}\tau _{k}$ are the eigenvalues of $A^*$.

Let $q(\theta )=(1,\alpha )^Te^{\iota \omega _{0}\tau _{k}\theta }$ be the eigenvector of $A(0)$ corresponding to the eigenvalue $\iota \omega _{0}\tau _k$ and $q^{*}(s)=D(1,\beta )e^{\iota \omega _{0}\tau _{k}s}$ be the eigenvector of $A^{*}$ corresponding to the eigenvalue $-\iota \omega _{0}\tau _k$, then

$$\begin{aligned} A(0)q(\theta )=\iota \omega \tau _kq(\theta ) \end{aligned}$$

It follows from the definition of $A(0)$, $L_{\mu }(\phi )$ and $\eta (\theta , \mu )$ that

$$\begin{aligned} \left[ (A_{11}+A_{22}e^{-\iota \omega \tau _0})-\iota \omega I\right] q(0)=0\\ \end{aligned}$$

where I is identity matrix of order 2, that is,

$$\begin{aligned} \left( \begin{array}{cc} a_{10}-\iota \omega &{}\quad a_{01}e^{-\iota \omega \tau _0} \\ b_{10} &{}\quad b_{01}-\iota \omega _0 \\ \end{array} \right) \left( \begin{array}{c} 1 \\ \alpha \\ \end{array} \right) =\left( \begin{array}{c} 0 \\ 0 \\ \end{array} \right) \end{aligned}$$

Solving the above equations, we get

$$\begin{aligned} q(0)=\left( \begin{array}{c} 1 \\ \alpha \\ \end{array} \right) =\left( \begin{array}{c} 1 \\ {\displaystyle \frac{b_{10}}{\iota \omega _{0}-b_{01}}}\\ \end{array} \right) \end{aligned}$$

From the definition of $A^*$, we obtain

$$\begin{aligned} A^*(\mu )\psi (s)&= \int _{-1}^0{\hbox {d}}\eta (t,0)\psi (-t)\\&= A_{11}^T\psi (0)+A_{22}^T\psi (-1)\\ \end{aligned}$$

Since $q(\theta )$ is the eigenvector of $A^*$ corresponding to -$\iota \omega _0\tau _k$, then we have

$$\begin{aligned}&A^{*}q^{*}(0)=-\iota \omega _{0}\tau _{kq}^{*}(\theta )\\&{\hbox {or}} \\&\left[ ({A_{11}^{T}}+{A_{22}^{T}} {e}^{\iota \omega \tau _{0}})+\iota \omega I\right] {\left( q^{*}(0)\right) }^T=0 \end{aligned}$$

where I is identity matrix of order 2, that is,

$$\begin{aligned} \left( \begin{array}{cc} a_{10}+\iota \omega &{}\quad b_{10} \\ a_{01}e^{-\iota \omega \tau _0} &{}\quad b_{01}+\iota \omega _0 \\ \end{array} \right) \left( \begin{array}{c} 1 \\ \beta \\ \end{array} \right) =\left( \begin{array}{c} 0 \\ 0 \\ \end{array} \right) \end{aligned}$$

Solving the above equations, we get

$$\begin{aligned} q(0)=\left( \begin{array}{c} 1 \\ \beta \\ \end{array} \right) =\left( \begin{array}{c} 1 \\ -{\displaystyle \frac{\iota \omega _0+a_{10}}{b_{10}}}\\ \end{array} \right) \end{aligned}$$

Now,

$$\begin{aligned}&(q^*, \bar{q})=\bar{D}\left[ (1,\bar{\beta })(1,\alpha )-\int _{-1}^0\int _{\xi =0}^\theta (1,\bar{\beta })\right. \nonumber \\&\quad \quad \left. e^{-\iota (\xi -\theta )\omega _0\tau _k}){\hbox {d}}\eta (\theta )(1, \alpha )^{T}e^{\iota \xi \omega _0\tau _k}\right] \\&\quad =\bar{D}\left[ 1+\alpha \bar{\beta }-\int _{-1}^0(1,\bar{\beta })\theta e^{\iota \theta \omega _0\tau _k}){\hbox {d}}\eta (\theta )(1, \alpha )^T\right] \\&\quad =\bar{D}\left[ 1+\alpha \bar{\beta }+\tau _k\alpha e^{-\iota \omega _0\tau _k}\right] \\ \end{aligned}$$

In order to have $(q^*, q)=1$, we can choose $D={\displaystyle \frac{1}{1+\bar{\alpha }\beta +\bar{\alpha }a_{01}\tau _{k}e^{\iota \omega _{0}\tau _{k}}}}$. such that $(\psi , A\phi )=(A^*\psi , \phi )$, or $-\iota \omega _0(q^*, \bar{q})=(q^*, A\bar{q})=(A^*q^*, \bar{q})=(-\iota \omega _0q^*, \bar{q})=\iota \omega _0(q^*, \bar{q})$ i.e. $(q^*, \bar{q})=0$.

Next we will compute the coordinates to describe the centre manifold $C_{0}$ at $\mu =0$. Let $u_{t}$ be the solution of (8.7) when $\mu =0$.

Define

$$\begin{aligned}&z(t)=<q^*,u_{t}>,\nonumber \\&\quad W(t,\theta )=u_{t}(\theta )-2Re\{z(t)q(\theta )\} \end{aligned}$$

(8.9)

On the centre manifold $C_{0}$, we have

$$\begin{aligned} W(t,\theta )=W(z(t),\bar{z}(t),\theta ) \end{aligned}$$

where

$$\begin{aligned} W(z(t),\bar{z}(t),\theta )&= W_{20}(\theta )\frac{z^2}{2}+W_{11}(\theta )z\bar{z}\nonumber \\&+\,W_{02}(\theta )\frac{{\bar{z}}^2}{2}+........ \end{aligned}$$

(8.10)

z and $\bar{z}$ are local coordinates for centre manifold $C_{0}$ in the direction of $q$ and $q^*$ respectively. Here we consider only real solutions as W is real if $u_{t}$ is real. For the solution $u_{t}\in C_{0}$ of (8.7), since $\mu =0$, we have

$$\begin{aligned} \dot{z}(t)&= \iota \omega _{0}\tau _{k}z+\bar{q}^*(0)f(0, W(z,\bar{z},0)+2Re(zq(\theta )))\\&= \iota \omega _{0}\tau _{k}z+\bar{q}^*(0)f_{0}(z,\bar{z}) \end{aligned}$$

This equation can be rewritten as

$$\begin{aligned} \dot{z}(t)=\iota \omega _0\tau _k z+g(z,\bar{z}) \end{aligned}$$

where

$$\begin{aligned} g(z,\bar{z})&= \bar{q}^*(0)f_{0}(z,\bar{z})\nonumber \\&= g_{20}(\theta )\frac{z^2}{2}+g_{11}(\theta )z\bar{z}+g_{02}(\theta )\frac{{\bar{z}}^2}{2}\nonumber \\&+\,g_{21}(\theta )\frac{z^2\bar{z}}{2}+........ \end{aligned}$$

(8.11)

From (8.9) and (8.10), we have

$$\begin{aligned} u_{t}(\theta )&= W(t,\theta )+2Re\{z(t)q(\theta )\} \nonumber \\&= W_{20}(\theta )\frac{z^2}{2}+W_{11}(\theta )z\bar{z}+W_{02}(\theta )\frac{{\bar{z}}^2}{2}\nonumber \\&+\,(1,q_1,q_2)^{T}e^{\iota \omega _0\tau _{k}\theta }z\nonumber \\&+\,(1,\bar{q_1},\bar{q_2})^{T}e^{-\iota \omega _0\tau _{k}\theta }\bar{z}+.... \end{aligned}$$

(8.12)

Now from (8.4) and (8.11), it follows that

$$\begin{aligned}&g(z,\bar{z})=\tau _k\bar{D}(1,\bar{\beta }) \left[ \begin{array}{cc} &{}a_{20}u^2_{1t}(0)+a_{11}u_{1t}(0)u_{2t}(-1) \nonumber \\ &{}b_{20}u^2_{1t}(0)+b_{11}u_{1t}(0)u_{2t}(0) \end{array}\right] \\&\quad =\tau _k\bar{D}(1,\bar{\beta }) \left[ \begin{array}{c} p_1z^2+p_2z\bar{z}+p_3\bar{z}^2+p_4\frac{z^2}{2}\bar{z}+........ \\ q_1z^2+q_2z\bar{z}+q_3\bar{z}^2+q_4\frac{z^2}{2}\bar{z}+........ \end{array}\right] \nonumber \\ \end{aligned}$$

(8.13)

or

$$\begin{aligned} g(z,\bar{z})&= \tau _k\bar{D}\left[ \left( p_1+\bar{\beta }q_1\right) z^2\right. \nonumber \\&+\left( p_2+\bar{\beta }q_2\right) z\bar{z}+\left( p_3+\bar{\beta }q_3\right) \bar{z}^2\nonumber \\&\left. +\left( p_4+\bar{\beta }q_4\right) \frac{z^2}{2}\bar{z}+............\right] \end{aligned}$$

(8.14)

where

$$\begin{aligned} p_1&= a_{20}+a_{11}\alpha e^{-\iota \omega \tau _0}, p_2=2a_{20}+2Re\left( \alpha e^{\iota \omega \tau _0}\right) a_{11},\\ p_3&= a_{20}+a_{11}\bar{\alpha }e^{\iota \omega \tau _0},\\ p_4&= \left( 2W^1_{20}(0)+4W^1_{11}(0)\right) a_{20}+a_{11}\left( \bar{\alpha }W^1_{20}(0)e^{\iota \omega \tau _0}\right. \\&\quad \left. +\,2\alpha W^2_{11}(0)e^{-\iota \omega \tau _0}+W^2_{20}(-1)\right) +2W^2_{11}(-1)\\ q_1&= b_{20}+b_{11}\alpha , q_2=2b_{20}+2Re(\alpha )b_{11},\\&q_3=b_{20}+b_{11}\bar{\alpha },\\ q_4&= \left( 2W^1_{20}(0)+4W^1_{11}(0)\right) b_{20}\nonumber \\&+\,b_{11}\left( \bar{\alpha }W^1_{20}(0)+2\alpha W^2_{11}(0)+W^2_{20}(0)\right) \\&+\,2W^2_{11}(0)\\ \end{aligned}$$

Comparing coefficients of (8.11) and (8.14), we find

$$\begin{aligned} {\left\{ \begin{array}{ll} g_{20}=2\bar{D}\tau _{k}\left( p_1+q_1\bar{\beta }\right) \\ g_{11}=\bar{D}\tau _{k}\left( p_2+q_2\bar{\beta }\right) \\ g_{02}=2\bar{D}\tau _{k}\left( p_3+q_3\bar{\beta }\right) \\ g_{21}=\bar{D}\tau _{k}\left( p_4+q_4\bar{\beta }\right) \\ \end{array}\right. } \end{aligned}$$

(8.15)

Since $W_{20}$ and $W_{11}$ are in $g_{21}$, we still to compute them.

Now from (8.7) and (8.9), we have

$$\begin{aligned} \dot{W}&= \dot{u}(t)-\dot{z}q-\bar{\dot{z}}\bar{q} \nonumber \\&= {\left\{ \begin{array}{ll} A(0)W-2Re\left\{ \bar{q^*}(0)f_{0}q(\theta )\right\} &{} \theta \in [-1,0) \\ A(0)W-2Re\left\{ \bar{q^*}(0)f_{0}q(0)\right\} \\ \quad +f_0(z,\bar{z}) &{} \theta =0 \end{array}\right. }\nonumber \\&\triangleq A(0)W+H(z,\bar{z},\theta ) \end{aligned}$$

(8.16)

where

$$\begin{aligned} H(z,\bar{z},\theta )&= H_{20}(\theta )\frac{z^2}{2}+H_{11}(\theta )z\bar{z}\nonumber \\&+\,H_{02}(\theta )\frac{{\bar{z}}^2}{2}+........ \end{aligned}$$

(8.17)

Substituting (8.17) into (8.16) and comparing the coefficients, we get

$$\begin{aligned}&(A(0)-2\iota \omega _{0}\tau _{k}I)W_{20}(\theta )=-H_{20}(\theta ),\nonumber \\&A(0)W_{11}(\theta )=-H_{11}(\theta ) \end{aligned}$$

(8.18)

From (8.16) and for $\theta \in [-1,0)$

$$\begin{aligned} H(z,\bar{z},\theta )&= -\bar{q}^*(0)f_0q(\theta )-q^*(0)\bar{f_0}\bar{q}(\theta )\nonumber \\&= -g(z,\bar{z})q(\theta )-\bar{g}(z,\bar{z})\bar{q}(\theta ) \end{aligned}$$

(8.19)

Using (8.11) in (8.19) and comparing coefficients with (8.17), we can obtain

$$\begin{aligned} H_{20}(\theta )=-g_{20}q(\theta )-\bar{g}_{02}\bar{q}(\theta ) \end{aligned}$$

(8.20)

and

$$\begin{aligned} H_{11}(\theta )=-g_{11}q(\theta )-\bar{g}_{11}\bar{q}(\theta ) \end{aligned}$$

(8.21)

From the definition of $A(0)$, (8.18) and (8.20), we obtain

$$\begin{aligned} \dot{W}_{20}(\theta )&= 2\iota \omega _0\tau _k W_{20}(\theta )+g_{20}q(\theta )+\bar{g}_{02}\bar{q}(\theta ) \end{aligned}$$

Solving it and for $q(\theta )=(1,q_1,q_2)^Te^{\iota \omega _0\tau _k\theta }$, we have

$$\begin{aligned} W_{20}(\theta )&= \frac{\iota g_{20}}{\omega _0\tau _k}q(0)e^{\iota \omega _0\tau _k\theta }+\frac{\iota \bar{g}_{02}}{3\omega _0\tau _k}\bar{q}(0)e^{-\iota \omega _0\tau _k\theta }\nonumber \\&+\,E_{1}e^{2\iota \omega _0\tau _k\theta } \end{aligned}$$

(8.22)

Similarly, from (8.18) and (8.21) it follows that,

$$\begin{aligned} W_{11}(\theta )&= -\frac{\iota g_{11}}{\omega _0\tau _k}q(0)e^{\iota \omega _0\tau _k\theta }\nonumber \\&+\frac{\iota \bar{g}_{11}}{\omega _0\tau _k}\bar{q}(0)e^{-\iota \omega _0\tau _k\theta }+E_{2} \end{aligned}$$

(8.23)

where $E_{1}=({E_{1}}^{(1)},{E_{1}}^{(2)})^{T}$ and $E_{2}=({E_{2}}^{(1)},{E_{2}}^{(2)})^{T}$ are three dimensional constant vectors, and can be determined by setting $\theta =0$ in $H(z,\bar{z},\theta )$.

Again from the definition of $A(0)$ and (8.18), we have

$$\begin{aligned} \int _{-1}^0{\hbox {d}}\varsigma (\theta )W_{20}(\theta )=2\iota \omega _{0}\tau _{k}W_{20}(0)-H_{20}(0) \end{aligned}$$

(8.24)

and

$$\begin{aligned} \int _{-1}^0{\hbox {d}}\varsigma (\theta )W_{11}(\theta )=-H_{11}(0) \end{aligned}$$

(8.25)

where $\varsigma (\theta )=\varsigma (0,\theta )$.

From (8.16), we know when $\theta =0$,

$$\begin{aligned} H(z,\bar{z},0)&= -2Re(\bar{q}^*(0)f_{0}q(0))+f_{0}(z,\bar{z}) \nonumber \\&= -\bar{q}^*(0)f_{0}q(0)-q^*(0)\bar{f}_{0}\bar{q}(0)\!+\!f_{0}(z,\bar{z}) \end{aligned}$$

That is,

$$\begin{aligned}&H_{20}(\theta )\frac{z^2}{2}+H_{11}(\theta )z\bar{z}+H_{02}(\theta )\frac{{\bar{z}}^2}{2}+.......\nonumber \\&\quad =-q(0)\left\{ g_{20}\frac{z^2}{2}+g_{11}z\bar{z}+g_{02}\frac{{\bar{z}}^2}{2}+......\right\} \nonumber \\&\quad \quad \bar{q}(0)\left\{ \bar{g}_{20}\frac{{\bar{z}}^2}{2}+\bar{g}_{11}z\bar{z}+\bar{g}_{02}\frac{z^2}{2}+....\right\} +f_{0}(z,\bar{z})\nonumber \\ \end{aligned}$$

(8.26)

From (8.9), we have

$$\begin{aligned} u_{t}(\theta )&= W(t,\theta )+2Re\{z(t)q(\theta )\}\\&= W(t,\theta )+z(t)q(\theta )+\bar{z}(t)\bar{q}(t)\\&= W_{20}(\theta )\frac{z^2}{2}+W_{11}(\theta )z\bar{z}+W_{02}(\theta )\frac{{\bar{z}}^2}{2}+.... \end{aligned}$$

Thus we can obtain,

$$\begin{aligned} f_{0}=2\tau _{k}\left[ \begin{array}{cc} &{}p_1\\ &{}q_1 \end{array}\right] {\displaystyle \frac{z^2}{2}}+\tau _k\left[ \begin{array}{cc} &{}p_2\\ &{}q_2 \end{array}\right] z\bar{z}+......... \end{aligned}$$

(8.27)

Comparing the coefficients in (8.26) and using (8.27), we get

$$\begin{aligned}&H_{20}(0)=-g_{20}q(0)-\bar{g}_{02}\bar{q}(0)+2\tau _k\left[ \begin{array}{cc} &{}p_1\\ &{}q_1 \end{array}\right] \end{aligned}$$

(8.28)

$$\begin{aligned}&H_{11}(0)=-g_{11}q(0)-\bar{g}_{11}\bar{q}(0)+\tau _k\left[ \begin{array}{cc}&{}p_2\\ &{}q_2 \end{array}\right] \end{aligned}$$

(8.29)

Since $\iota \omega _{0}\tau _{k}$ is the eigenvalue of $A(0)$ corresponding to eigenvector $q(0)$, then

$$\begin{aligned}&\{\iota \omega _{0}\tau _{k}I- \int _{-1}^0e^{\iota \omega _{0}\tau _{k}\theta }{\hbox {d}}\varsigma (\theta )\}q(0)=0 \quad {\hbox {and}} \\&\{-\iota \omega _{0}\tau _{k}I- \int _{-1}^0e^{-\iota \omega _{0}\tau _{k}\theta }{\hbox {d}}\varsigma (\theta )\}\bar{q}(0)=0. \end{aligned}$$

Substituting (8.22) and (8.28) into (8.24), we find

$$\begin{aligned} \{2\iota \omega _{0}\tau _{k}I-\int _{-1}^0e^{2\iota \omega _{0}\tau _{k}\theta }{\hbox {d}}\varsigma (\theta )\}E_{1} =2\tau _k\left[ \begin{array}{cc} &{}p_{1}\\ &{}q_{1} \end{array}\right] \end{aligned}$$

or

$$\begin{aligned} \left[ \begin{array}{cccc} &{}2\iota \omega _{0}-a_{10} &{} -a_{01}e^{-2\iota \omega \tau _0}\\ &{}-b_{10} &{} 2\iota \omega _{0}-b_{01} \end{array}\right] E_{1} =\left[ \begin{array}{cc} &{}2p_{1}\\ &{}2q_{1} \end{array}\right] \end{aligned}$$

or $E_1=2\left[ \begin{array}{cc} &{}p_1 \\ &{}q_1 \end{array}\right] \left[ \begin{array}{cccc} &{}2\iota \omega _{0}-a_{10} &{} -a_{01}e^{-2\iota \omega \tau _0}\\ &{}-b_{10} &{} 2\iota \omega _{0}-b_{01} \end{array}\right] ^{-1}$

Similarly substituting (8.23) and (8.29) into (8.25), we obtain

$$\begin{aligned} \left[ \begin{array}{cccc} &{} -a_{10} &{} -a_{01}\\ &{}-b_{10} &{} -b_{01} \end{array}\right] E_{2}=\left[ \begin{array}{cc}&{}p_2\\ &{}q_2 \end{array}\right] \end{aligned}$$

or

$$\begin{aligned} E_2=\left[ \begin{array}{cc}&{}p_2\\ &{}q_2 \end{array}\right] \left[ \begin{array}{cccc} &{} -a_{10} &{} -a_{01}\\ &{}-b_{10} &{} -b_{01} \end{array}\right] ^{-1} \end{aligned}$$

Thus, we can determine $W_{20}(\theta )$, $W_{11}(\theta )$ from (8.22), (8.23) and $g_{21}$ can be computed from (8.15).

Finally, we can compute the following quantities:

$$\begin{aligned} \begin{aligned}&c_{1}(0)={\displaystyle \frac{\iota \left\{ g_{20}g_{11}-2|g_ {11}|^2-\frac{|g_{02}|^2}{3}\right\} }{2\omega _{0}\tau _{k}}+\frac{g_{21}}{2}} \\&\mu _{2}=-{\displaystyle \frac{Re\left\{ c_{1}(0)\right\} }{Re\left\{ \frac{{\hbox {d}}{\uplambda }(\tau _k)}{{\hbox {d}}\tau }\right\} }}, \\&\beta _2=2Re\left\{ c_1(0)\right\} , \\&T_{2}=-{\displaystyle \frac{\mathfrak {I}\{ c_1(0)\}+\mu _{2}\mathfrak {I}\left\{ \frac{{\hbox {d}}{\uplambda }(\tau _k)}{{\hbox {d}}\tau }\right\} }{\omega _{0}\tau _{k}}}, \quad k=0,1,2,..... \end{aligned} \end{aligned}$$