5
$\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

2 Answers 2

4
$\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
5
$\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 at 17:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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