**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};

[![enter image description here][1]][1]

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

[![enter image description here][2]][2]

   The linear approximation at $P_{0}$ (coexistence equilibrium point):

    J0=Simplify@J[X0,μ];
    MatrixForm@J0
   
[![enter image description here][3]][3]

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}$.

**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, μ]]

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})$:

    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[[2]] . vpt[[2]]]]
    p = Expand@Simplify[vpt[[2]]/cn];
    pc = ComplexExpand[Conjugate[p]];
    MatrixForm@p
    MatrixForm@pc
    Simplify[JtX0μ0.p-I ω p]
    

Example in development.

  [1]: https://i.sstatic.net/We5lI.png
  [2]: https://i.sstatic.net/a45Eb.png
  [3]: https://i.sstatic.net/YqJdt.png