**Hopf bifurcation analysis**

   The differential system:

    f1[x_,y_]:=a x (1 - x/k) - b x y; 
    f2[x_,y_]:=-c y + d x y; 
    F[{x_,y_},{a_,b_,c_,d_,k_}]:=Evaluate@{f1[x,y],f2[x,y]};
    X={x, y};
    μ={a,b,c,d,k};
$$
\begin{align}
 &\dot{x}=a x\left(1-\frac{x}{k} \right)- bxy\\
 &\dot{y}= dxy - cy
\end{align}
$$

The Jacobian matrix:

    J[{x_,y_},{a_,b_,c_,d_,k_}]:=Evaluate@D[F[X,μ],{X}]

   The non-trivial equilibrium point:

    X0=Normal[Simplify[SolveValues[F[X,μ]==0&&Variables[F[X,μ]]>0,X]]][[1]];
    MatrixForm@X0
$$
\begin{align}
 P_{0}(x,y)=\left(\frac{c}{d} ,\frac{a}{b}\left(1-\frac{c}{dk} \right)\right)
\end{align}
$$
   The linear approximation at $P_{0}$ (coexistence equilibrium point):

    J0=Simplify@J[X0,μ];
    MatrixForm@J0
$$
\begin{align}
 J(P_{0})=\left(
\begin{array}{cc}
\hspace{-0.25cm}-\displaystyle\frac{a c}{d k} & -\displaystyle\frac{b c}{d}\hspace{0.3cm} \\\\
 \hspace{0.2cm}\displaystyle\frac{a (d k-c)}{b k} &\hspace{0.2cm} 0 \\
\end{array}
\right)
\end{align}
$$
Under the Hopf bifurcation conditions, $\text{tr}(J(P_{0},\mu_{0}))=0$ and $\text{det}(J(P_{0},\mu_{0}))>0$, where $\mu_{0}$ is the critical bifurcation value for some parameter of our system. In our case, the parameters are strictly positive and $\text{tr}(J(P_{0}))$ cannot be zero. Therefore, Hopf bifurcation not take place at $P_{0}$. The non-trivial equilibrium $P_{0}$ is always locally stable and the only condition that must be fulfilled is given by the following inequality
$$
\frac{c}{d k}<1
$$
Code for time series and phase portrait:

Time series

    s = ParametricNDSolve[{x'[t] == a x[t] (1 - x[t]/k) - bx[t]*y[t], 
    y'[t] == -c y[t] + d x[t]*y[t], x[0] == 11/5, y[0] == 4/5}, {x, y}, {t, 0, 1000}, {a, b, c, d, k}];
    Plot[Evaluate[x[1/4, 1/3, 1/2, 1/4, 10][t] /. s], {t, 0, 300}, PlotRange -> All, PlotPoints -> 500, 
    PlotStyle -> {Blue, Thickness[0.003]},AxesStyle -> Directive[Black, Small], Background -> Lighter[Gray, 0.95]]

[![Time series][1]][1]

Phase portrait:

    ParametricPlot[Evaluate[{x[1/4, 1/3, 1/2, 1/4, 10][t], y[1/4, 1/3, 1/2, 1/4, 10][t]} /. s], 
    {t, 0, 300}, PlotRange -> All, PlotPoints -> 500, PlotStyle -> {Blue, Thickness[0.003]}, 
    AxesStyle -> Directive[Black, Small],Background -> Lighter[Gray, 0.95]]

[![Phase portrait][2]][2]

**Example: Hopf bifurcation in the Brusselator system**

`Calculation of the first Lyapunov coefficient`

The Brusselator system is given by:
$$
\begin{align}
 &\dot{x}=A-(B+1)x + x^2 y\\
 &\dot{y}= B x - x^2 y
\end{align}
$$

Assuming $A> 0$ fixed and taking $B$ as a bifurcation parameter, we show that at $B = 1 + A^2$ the system exhibits a supercritical Hopf bifurcation.

The Brusselator system code:

    f1[x_, y_] := A - (B + 1) x + x^2 y;
    f2[x_, y_] := B x - x^2 y;
    F[{x_, y_}, {A_, B_}] := Evaluate@{f1[x, y], f2[x, y]};
    X = {x, y};
    μ = {A, B};

The Jacobian matrix and its transpose:

    J[{x_, y_}, {A_, B_}] = D[F[X, μ], {X}];
    Jt[{x_, y_}, {A_, B_}] = Transpose[J[X, μ]];
    MatrixForm[J[X, μ]]
    MatrixForm[Jt[X, μ]]

Stability analysis (Routh-Hurwitz criterion):

    X0[{A_, B_}] = SolveValues[F[X, μ] == 0, X][[1]]
    polJX0 = Collect[CharacteristicPolynomial[J[X0[μ], μ], λ], λ,Simplify];
    a0 = CoefficientList[polJX0, λ][[3]];
    a1 = CoefficientList[polJX0, λ][[2]];
    a2 = CoefficientList[polJX0, λ][[1]];
    Reduce[a1 > 0 && a2 > 0 && A > 0 && B > 0, B]
    (*A > 0 && 0 < B < 1 + A^2*)

Note that $a_{1}=0$ if and only if  $B=B_{0}$, where $B_{0}=1+A^2$. Then, the Brusselator is locally asymptotically stable at $X_{0}(\mu)$ for $B<B_{0}$ and locally asymptotically unstable for $B>B_{0}$ (appears a stable limit cycle surrounded the unstable equilibrium point). We verify the previous conclusion (transversality condition) with the sign of the following derivative:

    D[-a1, B]
    (*1*)

The analysis at the critical bifurcation value $B_{0}$:

Solve the following system of equations:
$$
\left\{\begin{align}
F\left((x,y),(A,B)\right) &=0, \\
\operatorname{tr}(J((x,y),(A,B))) &=0,
\end{align}\right.
$$
for $(x,y,B)$ and we must check that det $\operatorname{det}J((x,y),(A,B))>0$ when $B = B_{0}$ for the solution found, where $B_{0}$ is the Hopf critical bifurcation value.

The code for the above system of equations:

    X0μ0 = Delete[Part[SolveValues[F[X, μ] == 0 && Tr[J[X, μ]] == 0, {x, y, B}], 1], {3}]
    μ0 = Prepend[Delete[Part[SolveValues[F[X, μ] == 0 && Tr[J[X, μ]] == 0, {x, y, B}], 1], {{1}, {2}}], A]
    Det[J[X0μ0, μ0]]
Here, the Hopf critical bifurcation value is $B_{0}=1+A^2$ and $\operatorname{det}J((x,y),(A,B_{0}))=A^2>0$. Thus, the Brusselator at $B_{0}=1+A^2$ has the equilibrium
$$
\begin{align}
X_{0}(\mu_{0})=\left(A, \displaystyle\frac{1+A^2}{A} \right)
\end{align}
$$
and the linear approximation at $X_{0}(\mu_{0})$ has purely imaginary eigenvalues $\lambda_{1,2}=\pm \omega i$, $\omega=A$.

The code for the linear approximation and its transpose at $X_{0}(\mu_{0})$:

    A= ω;
    JX0μ0 = Simplify@J[X0μ0, μ0];
    JtX0μ0 = Simplify@Transpose@JX0μ0;
    MatrixForm@JX0μ0
    MatrixForm@JtX0μ0
    Eigenvalues[JX0μ0]

The next step is to translate the equilibrium $X_{0}(\mu_{0})$ to the origin of coordinates:

    bb = {0, 0};
    F0[{x_, y_}, {A_, B_}] = Collect[Expand@F[X + X0μ0, μ0], {x, x^2, y, y^2, x y, x^2 y},Factor]
    MatrixForm@F0[bb,μ]

Now, to  obtain the normal form of the Hopf bifurcation, we need the Taylor expansion of the third order for $F_{0}((x,y),(A,B))$:

    (*Rank 3 tensor*)
    D2[{x_, y_}, {A_, B_}] = Simplify@D[F0[X, μ], {X, 2}]
    D2X0μ0 = Simplify@D2[bb, μ0]
    (*Rank 4 tensor*)
    D3[{x_, y_}, {A_, B_}] = Simplify@D[F0[X, μ], {X, 3}]
    D3X0μ0 = Simplify@D3[bb, μ0]

We verify that the first three terms of the Taylor series expansion of $F_{0}((x,y),(A,B))$ are correct:

    FullSimplify[F0[X, μ]-(JX0μ0.X + 1/2!D2X0μ0.X.X + 1/3! D3X0μ0.X.X.X) /. {x -> t x, y -> t y}]
    (*{0,0}*)

Now, we compute the critical eigenvectors of $J((0,0),\mu_{0})$ and its transpose:

    (*Eigenvectors for J[X0,μ0]*)
    vp = ComplexExpand@Eigenvectors[JX0μ0]
    q = vp[[2]];
    qc = vp[[1]];
    MatrixForm@q
    MatrixForm@qc
    MatrixForm@Simplify[JX0μ0.q - I ω q]
    (*Eigenvectors for Transpose[J[X0,μ0]]*)
    vpt = ComplexExpand[Eigenvectors[JtX0μ0]]
    (*Normalization constant*)
    cn = ComplexExpand[Conjugate[vp[[1]] . vpt[[1]]]]
    p = Expand@Simplify[vpt[[2]]/cn];
    pc = ComplexExpand[Conjugate[p]];
    MatrixForm@p
    MatrixForm@pc
    Simplify[JtX0μ0.p-I ω p]

We verify the normalization condition $\langle p,q\rangle=1$

    Simplify@(p.q)
    (*1*)
Finally, we compute the first Lyapunov coefficient:
$$
\begin{align}
l_1(0,\mu_{0})= &\frac{1}{2\omega_0} {\rm Re}\left[\langle p,C(q,q,\bar{q}) \rangle - 2 \langle p, B(q,A_0^{-1}B(q,\bar{q}))\rangle 
+\\\hspace{0.5cm} \langle p, B(\bar{q},(2i\omega_0 I_n-A_0)^{-1}B(q,q))\rangle \right]
\end{align}
$$

Before to calculate $l_1(0,\mu_{0})$, we clean $A$

    Clear[A]
    ω0=A;

The code for the first Lyapunov coefficient:

    Factor@ComplexExpand[Re[1/(2 ω) (p.D3X0μ0.q.q.qc- 2(p.D2X0μ0.q.Inverse[JX0μ0].D2X0μ0.q. qc) +
    p.D2X0μ0.qc. Inverse[2 I ω*IdentityMatrix[2]-JX0μ0].D2X0μ0.q. q)] /. ω -> ω0]
    (*-((2 + A^2)/(2 A (1 + A^2)))*)

$$
\begin{align}
l_1(0,\mu_{0})=-\frac{A^2+2}{2 A \left(A^2+1\right)}
\end{align}
$$

The above expression is the result that Kuznetsov arrives at on page 105 of his book (see [Elements of Applied Bifurcation Theory][3]).

Limit cycle:

[![Limit cycle][4]][4]

For more details see: [Andronov-Hopf bifurcation][5].

###### In case anyone wishes to collaborate on these topics, please, communicate.


  [1]: https://i.sstatic.net/pRAqT.png
  [2]: https://i.sstatic.net/rQ6Nk.png
  [3]: https://doi.org/10.1007/978-1-4757-3978-7
  [4]: https://i.sstatic.net/SCEwC.png
  [5]: http://www.scholarpedia.org/article/Andronov-Hopf_bifurcation