5
$\begingroup$

Good day. I need help with the code in mathematica to plot the bifurcation diagram (e vs z*) or (e vs x*), for the system x'[t]=x(1-x)-k1xy,y'[t]=ey(1-y)+k2xy-k3yz,z'[t]=-k4z+ek3yz, where X*=(x*,y*,z*) is a fixed point of the system, with k4=0.1,k1=0.02,k2=0.1,k3=1.2, and ebelongs to [0,4].

Thank you

Improve this question
$\endgroup$
2
  • $\begingroup$ or the bifurcation diagram (e vs z), where z is a solution $\endgroup$
    – Patrick
    Sep 23, 2019 at 13:40
  • 1
    $\begingroup$ Welcome to Mathematica.SE, Patrick! I suggest the following: 1) Take the tour and check the faqs. 2) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! 3) As you receive help, try to give it too, by answering questions in your area of expertise. $\endgroup$
    – Chris K
    Sep 23, 2019 at 14:18

1 Answer 1

Reset to default
10
$\begingroup$

Looks like a little three-species food web model -- a perfect excuse to use my new EcoEvo package.

First, install the package (one-time only):

PacletInstall["https://github.com/cklausme/EcoEvo/releases/download/v1.0.1/EcoEvo-1.0.1.paclet"]

Then load the package and set the model for analysis with SetModel:

<< EcoEvo`

SetModel[{
  Pop[x] -> {Equation :> x[t] (1 - x[t]) - k1 x[t] y[t]},
  Pop[y] -> {Equation :> e y[t] (1 - y[t]) + k2 x[t] y[t] - k3 y[t] z[t]},
  Pop[z] -> {Equation :> -k4 z[t] + e k3 y[t] z[t]}
  }]

Set the parameters, then solve for the equilibria:

k4 = 0.1; k1 = 0.02; k2 = 0.1; k3 = 1.2;
eq = SolveEcoEq[]

enter image description here

Now you can plot the parts of the final, coexistence equilibrium vs e:

Plot[x /. eq[[-1]], {e, 0, 4}, PlotRange -> {0, All}]
Plot[z /. eq[[-1]], {e, 0, 4}, PlotRange -> {0, All}]

Mathematica graphics Mathematica graphics

Just for fun, you can simulate the model with EcoSim:

e = 1;
sol = EcoSim[{x -> 1, y -> 1, z -> 0.1}, 100];
PlotDynamics[sol]

Mathematica graphics

Improve this answer
$\endgroup$
10
  • 1
    $\begingroup$ There's a bug in your add-on: the documentation pages are asking for /Users/klaus/Library/Mathematica/SystemFiles/FrontEnd/StyleSheets/MyStyleSheet.nb ! $\endgroup$
    – murray
    Sep 23, 2019 at 15:44
  • $\begingroup$ @murray Ugh! I'll try to fix it tonight. Thanks for the note! $\endgroup$
    – Chris K
    Sep 23, 2019 at 16:43
  • 1
    $\begingroup$ @ChrisK: I just copied what you posted and was able to duplicate the results. I did get one error message: "Solve::ratnz: Solve was unable to solve the system with inexact coefficients. The answer was obtained by solving a corresponding exact system and numericizing the result." Thanks. $\endgroup$
    – Moo
    Sep 23, 2019 at 21:34
  • 1
    $\begingroup$ @moo That's more of a warning than an error. $\endgroup$
    – Chris K
    Sep 24, 2019 at 6:34
  • 1
    $\begingroup$ @Moo I put a few examples in the built-in docs - check out EcoEvo/tutorial/ExampleModels in the help browser. $\endgroup$
    – Chris K
    Sep 24, 2019 at 12:13

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.