# 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 at 17:22