2
$\begingroup$
(*Define the model equations*)
modelEquations[{x_, y_}, {a_, k_, b_, d_, p_, c_, e_, 
   f_}] := {x Exp[
    a (1 - x/k) - (b (1 - p) y)/(1 + c (1 - p) x + d y)], 
  y Exp[((e b (1 - p) x)/(1 + c (1 - p) x + d y)) - f]}

(*Compute the Jacobian matrix*)
jacobian[{x_, y_}, params_] := 
 D[modelEquations[{x, y}, params], {{x, y}}]

(*Initialize parameters and ranges*)
params = {a, k, b, d, p, c, e, f};
initialConditions = {2.6, 2.7};
paramValues = {3.2, 1.8, 0.2, 0.2, 0.9, 2.9};  (*b,d,p,c,e,f*)
aRange = Range[0.1, 10, 0.1];
kRange = Range[0.1, 10, 0.1];

(*Simulation and bifurcation analysis*)
bifurcationAnalysis = 
  Flatten[Table[{a, k, 
     LyapunovExponents = 
      Re[Mean[Log[
         Abs[Eigenvalues[
           jacobian[initialConditions, {a, k}~Join~paramValues]]]]]], 
     BifurcationType = 
      Which[And @@ (LyapunovExponents < 0), "Stable Node", 
       And @@ (LyapunovExponents > 0), "Unstable Node", True, 
       "Saddle Point"]}, {a, aRange}, {k, kRange}], 1];

(*Visualization*)
ListContourPlot[bifurcationAnalysis[[All, {1, 2, 3}]], 
 FrameLabel -> {"a", "k"}, PlotLegends -> Automatic, Contours -> 20, 
 ColorFunction -> "TemperatureMap", PlotRange -> All, 
 PlotLabel -> "Bifurcation Diagram with Lyapunov Exponents"]

I don't understand how to fix this error
enter image description here

$\endgroup$
5
  • 4
    $\begingroup$ The error is from the D in D[modelEquations[{x, y}, params], {{x, y}}]. {x,y} is passed in as initialConditions, which is {2.6, 2.7}. $\endgroup$
    – MelaGo
    Mar 25 at 0:16
  • $\begingroup$ Change your definition from: D[modelEquations[{x0, y0}, params], {{x0, y0}}] /. {x0 -> x, y0 -> y} to D[modelEquations[{x0, y0}, params], {{x0, y0}}] /. {x0 -> x, y0 -> y} $\endgroup$ Mar 25 at 8:41
  • $\begingroup$ Not necessarily a Mathematica comment, but could you explain what you're calculating? Usually calculating Lyapunov exponents requires iterating the model on the attractor, but I don't see that here. (could be I just don't follow the code well) $\endgroup$
    – Chris K
    Mar 25 at 16:30
  • $\begingroup$ @ChrisK the method used in the code is a shortcut that leverages the Jacobian matrix to approximate Lyapunov exponents locally rather than iterating the full dynamical system over time. This approach can be useful for a bifurcation analysis but does not replace the full Lyapunov exponent calculation for understanding the global dynamical behavior. $\endgroup$ Mar 25 at 22:30
  • $\begingroup$ @AthanasiosParaskevopoulos Interesting, could you give a reference? $\endgroup$
    – Chris K
    Mar 26 at 17:07

1 Answer 1

6
$\begingroup$

To solve this problem it could be better to define system of equations and Jacobian as expressions in a form of

(*Define the model equations*)
modelEquations = {x Exp[
     a (1 - x/k) - (b (1 - p) y)/(1 + c (1 - p) x + d y)], 
   y Exp[((e b (1 - p) x)/(1 + c (1 - p) x + d y)) - f]};

(*Compute the Jacobian matrix*)
jacobian = D[modelEquations, {{x, y}}];

With these definitions we have

(*Initialize parameters and ranges*)
params = {a, k, b, d, p, c, e, f};
initialConditions = {2.6, 2.7};
paramValues = {3.2, 1.8, 0.2, 0.2, 0.9, 2.9};  (*b,d,p,c,e,f*)
aRange = Range[0.1, 10, 0.1];
kRange = Range[0.1, 10, 0.1];

(*Simulation and bifurcation analysis*)
bifurcationAnalysis = 
  Flatten[Table[{a1, k1, 
     LyapunovExponents = 
      Re[Mean[Log[
         Abs[Eigenvalues[
           jacobian /. 
             Thread[params -> Join[{a1, k1}, paramValues]] /. 
            Thread[{x, y} -> initialConditions]]]]]], 
     BifurcationType = 
      Which[And @@ (LyapunovExponents < 0), "Stable Node", 
       And @@ (LyapunovExponents > 0), "Unstable Node", True, 
       "Saddle Point"]}, {a1, .05, 10, 0.05}, {k1, 0.05, 10, .05}], 
   1];

Visualization

ListContourPlot[bifurcationAnalysis[[All, 1 ;; 3]], 
 FrameLabel -> {"a", "k"}, PlotLegends -> Automatic, Contours -> 20, 
 ColorFunction -> "TemperatureMap", PlotRange -> Automatic, 
 PlotLabel -> "Bifurcation Diagram with Lyapunov Exponents"]

Figure 1

$\endgroup$

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