Showing that a Hopf bifurcation exists?

Question

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

Answer 1

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"]

Now, load the package, set the model and add your assumptions on parameters:

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

Answer 2

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.

Answer 3

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 = 
 Thread[par > 0];
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.

Answer 4

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<q_0$:

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!