6
$\begingroup$

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

$\endgroup$
2
  • $\begingroup$ To be clear, there are two parameters left (q and c)? $\endgroup$
    – Chris K
    Oct 27, 2021 at 2:00
  • $\begingroup$ @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. $\endgroup$
    – tardigrade
    Oct 27, 2021 at 2:16

4 Answers 4

8
$\begingroup$

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

enter image description here

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]

enter image description here

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

enter image description here

enter image description here

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} *)

enter image description here

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} *)

enter image description here

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} *)

enter image description here

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} *)

enter image description here

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} *)

enter image description here

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} *)

enter image description here

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]

enter image description here

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

enter image description here

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

$\endgroup$
1
  • $\begingroup$ Thank you so much!!!! $\endgroup$
    – tardigrade
    Oct 29, 2021 at 8:18
7
$\begingroup$

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

enter image description here

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

enter image description here

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

For the second Eigenvalue we get:

enter image description here

Again 2 pure imaginary conjugate complex Eigenvalues are not possible.

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

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.

$\endgroup$
0
1
$\begingroup$

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

enter image description here

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!

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.