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 $$ 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).
For more details see: Andronov-Hopf bifurcation.