# Bifurcation Analysis and Lyapunov Exponents Visualization [closed]

(*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

• 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}. Mar 25 at 0:16
• 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} Mar 25 at 8:41
• 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) Mar 25 at 16:30
• @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. Mar 25 at 22:30
• @AthanasiosParaskevopoulos Interesting, could you give a reference? Mar 26 at 17:07

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