# Showing that a Hopf bifurcation exists?

I have the following system of ODEs:

dx/dt = a/(1 + z) - Q*z

dy/dt = Qx - qy

dz/dt = qy - cz/(K + z).

Assuming K = 1, Q = q < c, and a = c*(Sqrt[c/Q] - 1), is there a way to use Mathematica to test whether a Hopf bifurcation exists?

Any help would be greatly appreciated!

Thank you

• To be clear, there are two parameters left (q and c)? Oct 27, 2021 at 2:00
• @ChrisK Yes, we are missing q and c. I found that a bifurcation happens when q = c/4, but I have no idea how to test if this is a Hopf bifurcation. Oct 27, 2021 at 2:16

We can use my EcoEvo package, despite the fact that this doesn't seem like an ecological model.

First, install the package (only needs to be done once):

PacletInstall["EcoEvo", "Site" -> "http://raw.githubusercontent.com/cklausme/EcoEvo/master"]


<< EcoEvo

SetModel[{
Aux[x] -> {Equation :> a/(1 + z) - Q z},
Aux[y] -> {Equation :> Q x - q y},
Aux[z] -> {Equation :> q y - c z/(k + z)}
}]

k = 1;
Q := q;
a := c (Sqrt[c/Q] - 1);


We can find equilibria with SolveEcoEq:

eq = SolveEcoEq[]


There are two. To get an overview of where bifurcations happen in the c--q plane, plot where the maximum real part of the eigenvalues at each equilibrium equal zero:

λ1[c_?NumericQ, q_?NumericQ] := Max[Re[EcoEigenvalues[eq[[1]]]]];
λ2[c_?NumericQ, q_?NumericQ] := Max[Re[EcoEigenvalues[eq[[2]]]]];

ContourPlot[{λ1[c, q] == 0, λ2[c, q] == 0}, {c, -2, 2}, {q, -2, 2},
FrameLabel -> Automatic, MaxRecursion -> 3]


To investigate more closely, let's take an arbitrary slice at q = 0.2.

q = 0.2;
Plot[λ1[c, q], {c, -1, 2}]
Plot[λ2[c, q], {c, -1, 2}]


Seems like each equilibrium has two bifurcations in this range. We can use FindRoot to solve for them, eyeballing the graphs for initial guesses:

bif1a = FindRoot[λ1[c, q], {c, 0.19}]
bif1b = FindRoot[λ1[c, q], {c, 1.1}]
bif2a = FindRoot[λ2[c, q], {c, 0.01}]
bif2b = FindRoot[λ2[c, q], {c, 0.5}]

(* {c -> 0.2} *)
(* {c -> 1.03904} *)
(* {c -> 0.070988} *)
(* {c -> 0.466855} *)


Now we can evaluate the eigenvalues of the jacobian matrix at the bifurcation points to see if they're Hopf bifurcations (two complex conjugate eigenvalues with zero real part, non-zero imaginary part).

c = c /. bif1a;
EcoEigenvalues[eq[[1]]]
(* {-284390., -0.199832, -0.000168003} *)

c = c /. bif1b;
EcoEigenvalues[eq[[1]]]
(* {-0.42987, 0. + 0.214415 I, 0. - 0.214415 I} *)

c = c /. bif2a;
EcoEigenvalues[eq[[2]]]
(* {-0.303951, 0. + 0.144188 I, 0. - 0.144188 I} *)

c = c /. bif2b;
EcoEigenvalues[eq[[2]]]
(* {-0.358281, 0. + 0.177922 I, 0. - 0.177922 I} *)


Nope, yep, yep, yep!

We can simulate the dynamics for select values of c to make sure.

bif2a:

c = 0.05;
EcoEigenvalues[eq[[2]]]
sol = EcoSim[RuleListAdd[eq[[2]], {x -> 0.01}], 1000];
PlotDynamics[sol]
(* {-0.297406, 0.0143882 + 0.148584 I, 0.0143882 - 0.148584 I} *)


c = 0.1;
EcoEigenvalues[eq[[2]]]
sol = EcoSim[RuleListAdd[eq[[2]], {x -> 0.01}], 1000];
PlotDynamics[sol]
(* {-0.316999, -0.0155364 + 0.139724 I, -0.0155364 - 0.139724 I} *)


bif2b:

c = 0.45;
EcoEigenvalues[eq[[2]]]
sol = EcoSim[RuleListAdd[eq[[2]], {x -> 0.01}], 1000];
PlotDynamics[sol]
(* {-0.358413, -0.00121928 + 0.176904 I, -0.00121928 - 0.176904 I} *)


c = 0.48;
EcoEigenvalues[eq[[2]]]
sol = EcoSim[RuleListAdd[eq[[2]], {x -> 0.01}], 10000];
PlotDynamics[sol]
(* {-0.358171, 0.000917014 + 0.178671 I, 0.000917014 - 0.178671 I} *)


bif1b:

c = 1.0;
EcoEigenvalues[eq[[1]]]
sol = EcoSim[RuleListAdd[eq[[1]], {x -> 0.1}], 1000];
PlotDynamics[sol]
(* {-0.435545, -0.00286826 + 0.213892 I, -0.00286826 - 0.213892 I} *)


c = 1.1;
EcoEigenvalues[eq[[1]]]
sol = EcoSim[RuleListAdd[eq[[1]], {x -> 0.1}], 10000];
PlotDynamics[sol]
(* {-0.422319, 0.00401044 + 0.215014 I, 0.00401044 - 0.215014 I} *)


Interesting to note that at c = 1.1 both equilibria are unstable. If we start next to eq[[2]] we actually end up on a different limit cycle:

sol2 = EcoSim[RuleListAdd[eq[[2]], {x -> 0.1}], 10000];
PlotDynamics[sol2]


RuleListPlot[{FinalSlice[sol, 100], FinalSlice[sol2, 100]}]


So, there's a lot going on here. BTW, I don't see anything special about q = c/4.

• Thank you so much!!!! Oct 29, 2021 at 8:18

You will have a Hopf bifurcation if two Eigenvalues of the Jacobian at a fixpoint have pure imaginary conjugate complex values.

Therefore, we first need the equilibrium or fix points. With the parametervalues from your post we get the fixpoints:

rhs= {a/(1+z)- Q z, Q x- q y, q y - c z/(K+z)} //. {K->1,Q->q,a->c(Sqrt[c/Q]-1),q->c/4}

fixpt = {x, y, z} /. Solve[rhs == 0, {x, y, z}]


Now we can linearize the right hand sides around the fixpoints. For the first fixpoint we get:

lin = Series[rhs,
Sequence @@ Transpose[{{x, y, z}, fixpt[[1]], {1, 1, 1}}]] //
Normal;


The Jacobi matrix (all first derivatives) can be get by:

jm = D[lin, {{x, y, z}}];


And the Eigenvalues:

Eigenvalues[jm] // N


You see that we can not get 2 pure imaginary conjugate complex Eigenvalues.

For the second Eigenvalue we get:

Again 2 pure imaginary conjugate complex Eigenvalues are not possible.

• So because two of the Eigenvalues have conjugate complex values, we know that the bifurcation is a Hopf bifurcation? Thank you! Oct 27, 2021 at 22:09
• Sorry, I remembered this wrong. The eigenvalues must be conjugate and purely imaginary. I will look at it if I have time. Oct 28, 2021 at 7:38
• I corrected my answer. Oct 28, 2021 at 16:18
• So, to be clear, there is no Hopf when q=c/4, but three Hopf when q = 0.2 -- see first answer Feb 26, 2022 at 17:22

I will make a guess that this system comes from an application involving positive variables and parameters. This "essential nonnegativity" is only possible if the first equation is wrong, and Q z is in fact Q x. The fixed point becomes then unique. I will ignore the assumption on a, since Mathematica can handle problems with three parameters. The essential equation for a "candidate of Hopf bif. in three dimensions" is c1 c2 ==c3, where ci are the coefficients of the Char. Pol. of the Jacobian -- see for example Routh–hurwitz-liénard–chipart criteria SL Wiggers, P Pedersen A short program will produce a condition on a (not the one given), from which we can get examples and test them graphically

xv = a/(1 + z) - q*x;
yv = q*x - q*y;
zv = q*y - c*z/(1 + z); X = {x, y, z}; par = {q, c, a}; cp =
dyn = {xv, yv, zv}; fp = Solve[Thread[dyn == 0], X][[1]]
jac = D[dyn, {X}]; jac // MatrixForm
tr = Tr@jac /. fp // FullSimplify
det = Det[jac] /. fp // FullSimplify
pol = -CharacteristicPolynomial[jac, s] /. fp // Apart
cof = CoefficientList[pol, s]
re = Reduce[Join[{cof[[1]] *cof[[2]] - cof[[3]] == 0 && q < c}, cp]]
fi = FindInstance[re, par, 3]; fi // N


The original question could be slightly harder, due to the two fixed points.

The following answer is just a guide to detect the Hopf bifurcation in your system:

Considering $$a=\frac{3}{8}$$, $$c=\frac{4}{3}$$, $$k=\frac{3}{2}$$, $$Q=\frac{1}{2}$$ and $$q$$ as a bifurcation parameter, your system becomes as follows:

F[{x_, y_, z_}][q_] := {-(z/2) + 3/(8 (1 + z)), 1/2 (x - 2 q y), q y - (8 z)/(9 + 6 z)}
X = {x, y, z};


The nontrivial equilibrium point:

X0[q_] := First@SolveValues[F[X][q] == 0, X]
X0[q]
(*{2/3, 1/(3 q), 1/2}*)


Jacobian matrix:

J[{x_, y_, z_}][q_] := Evaluate[D[F[X][q], {X}]]
J[X][q]
(*{{0, 0, -(1/2) - 3/(8 (1 + z)^2)}, {1/2, -q, 0}, {0, q, (48 z)/(9 + 6 z)^2 - 8/(9 + 6 z)}}*)


The linear approximations at $$X_{0}(q)$$:

J0[q_] := J[X0[q]][q]
J0[q]
(*{{0, 0, -(2/3)}, {1/2, -q, 0}, {0, q, -(1/2)}}*)


Local stability using BialternateSum:

Reduce[Det[BialternateSum[J0[q]]] < 0 && q > 0, q]
(*q > 1/6*)


Looking for the Hopf bifurcation value using BialternateSum:

q0 = First@SolveValues[Det[BialternateSum[J0[q]]] == 0 && q > 0, q]
(*1/6*)


The eigenvalues at $$q=q_0$$:

Eigenvalues[J0[q0]]
(*{-(2/3), I/(2 Sqrt[3]), -(I/(2 Sqrt[3]))}*)
`

The limit cycle appears for $$q:

There is a condition that makes the two non-trivial equilibria collapse into one, so it is possible to study a Bogdanov-Takens bifurcation in your system.

I hope you enjoy it!