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};
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
The linear approximation at $P_{0}$ (coexistence equilibrium point):
J0=Simplify@J[X0,μ];
MatrixForm@J0
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.