2
$\begingroup$

I'm interested in carrying out stability analysis on discrete dynamic systems, as was done in this paper. The Author used their function TwoDimensionalStability[x^2 + 3 x y, y^2 - 2 x y] to look for stability properties. The output was:

The fixed point {− 2/7 , 3/7} is saddle
The fixed point {0, 1}is source
The fixed point {1, 0} is oscillatory source
The fixed point{0, 0}is Stable sink

I tried to run their code in my Mathematica 12 but I only ever got this:

Out[109]=TwoDimensionalStability[2 x^2 + 3 xy, 0.315165 - 2 xy]

Here is the code from the paper I mentioned above:

TwoDimensionalStability[f , g ]:=Module[{cz, czm, F, G, J, A, l},cz = Solve[{f == x, g == y}, {x, y}]; czm = {x, y}/.cz; F [u , v ]:=f /.x → ucz = Solve[{f == x, g == y}, {x, y}]; czm = {x, y}/.cz; F [u , v ]:=f /.x → ucz = Solve[{f == x, g == y}, {x, y}]; czm = {x, y}/.cz; F [u , v ]:=f /.x → u/.y → v; G[u , v ]:=g/.x → u /.y → v;/.y → v; G[u , v ]:=g/.x → u /.y → v;/.y → v; G[u , v ]:=g/.x → u /.y → v;J[{u , v }] = D[{F [u, v], G[u, v]}, {{u, v}}]; l = Length[czm];J[{u , v }] = D[{F [u, v], G[u, v]}, {{u, v}}]; l = Length[czm];J[{u , v }] = D[{F [u, v], G[u, v]}, {{u, v}}]; l = Length[czm];If[l==0, Print[“No fixed point”], Do[If[Im[First[Eigenvalues[J[czm[[i]]]]]]==0,If[l==0, Print[“No fixed point”], Do[If[Im[First[Eigenvalues[J[czm[[i]]]]]]==0,If[l==0, Print[“No fixed point”], Do[If[Im[First[Eigenvalues[J[czm[[i]]]]]]==0,If[−1 < Max[Eigenvalues[J[czm[[i]]]]] < 1&& − 1 < Min[Eigenvalues[J[czm[[i]]]]] < 1,If[−1 < Max[Eigenvalues[J[czm[[i]]]]] < 1&& − 1 < Min[Eigenvalues[J[czm[[i]]]]] < 1,If[−1 < Max[Eigenvalues[J[czm[[i]]]]] < 1&& − 1 < Min[Eigenvalues[J[czm[[i]]]]] < 1,Print[“The fixed point” , czm[[i]], “is Stable sink”],Print[“The fixed point” , czm[[i]], “is Stable sink”],Print[“The fixed point” , czm[[i]], “is Stable sink”],Which[Max[Eigenvalues[J[czm[[i]]]]] > 1&&0 < Min[Eigenvalues[J[czm[[i]]]]] < 1,Which[Max[Eigenvalues[J[czm[[i]]]]] > 1&&0 < Min[Eigenvalues[J[czm[[i]]]]] < 1,Which[Max[Eigenvalues[J[czm[[i]]]]] > 1&&0 < Min[Eigenvalues[J[czm[[i]]]]] < 1,Print[“The fixed point”, czm[[i]], “ is saddle”],Print[“The fixed point”, czm[[i]], “ is saddle”],Print[“The fixed point”, czm[[i]], “ is saddle”],−1 < Max[Eigenvalues[J[czm[[i]]]]] < 0&&Min[Eigenvalues[J[czm[[i]]]]] < −1,−1 < Max[Eigenvalues[J[czm[[i]]]]] < 0&&Min[Eigenvalues[J[czm[[i]]]]] < −1,−1 < Max[Eigenvalues[J[czm[[i]]]]] < 0&&Min[Eigenvalues[J[czm[[i]]]]] < −1,Print[“The fixed point”, czm[[i]] , “is oscillatory saddle”],Print[“The fixed point”, czm[[i]] , “is oscillatory saddle”],Print[“The fixed point”, czm[[i]] , “is oscillatory saddle” ,
Min[Eigenvalues[J[czm[[i]]]]] > 1, Print[“The fixed point”, czm[[i]] , “is source”],Min[Eigenvalues[J[czm[[i]]]]] > 1, Print[“The fixed point”, czm[[i]] , “is source”],Min[Eigenvalues[J[czm[[i]]]]] > 1, Print[“The fixed point”, czm[[i]] , “is source”],Min[Eigenvalues[J[czm[[i]]]]]== − 1&&Max[Eigenvalues[J[czm[[i]]]]] < 1,Min[Eigenvalues[J[czm[[i]]]]]== − 1&&Max[Eigenvalues[J[czm[[i]]]]] < 1,Min[Eigenvalues[J[czm[[i]]]]]== − 1&&Max[Eigenvalues[J[czm[[i]]]]] < 1,Print[“Try to use Center manifold module for ”, czm[[i]]],Print[“Try to use Center manifold module for ”, czm[[i]]],Print[“Try to use Center manifold module for ”, czm[[i]]],Min[Eigenvalues[J[czm[[i]]]]]==1, Print[“The fixed point”, czm[[i]], “is unstable”],Min[Eigenvalues[J[czm[[i]]]]]==1, Print[“The fixed point”, czm[[i]], “is unstable”],Min[Eigenvalues[J[czm[[i]]]]]==1, Print[“The fixed point”, czm[[i]], “is unstable”],Max[Eigenvalues[J[czm[[i]]]]] == −1, Print[“The fixed point”, czm[[i]], “is unstable”],Max[Eigenvalues[J[czm[[i]]]]] == −1, Print[“The fixed point”, czm[[i]], “is unstable”],Max[Eigenvalues[J[czm[[i]]]]] == −1, Print[“The fixed point”, czm[[i]], “is unstable”],Max[Eigenvalues[J[czm[[i]]]]]==1&&Min[Eigenvalues[J[czm[[i]]]]] < −1,Max[Eigenvalues[J[czm[[i]]]]]==1&&Min[Eigenvalues[J[czm[[i]]]]] < −1,Max[Eigenvalues[J[czm[[i]]]]]==1&&Min[Eigenvalues[J[czm[[i]]]]] < −1,Print[“The fixed point”, czm[[i]], “is unstable”],Print[“The fixed point”, czm[[i]], “is unstable”],Print[“The fixed point”, czm[[i]], “is unstable”],Min[Eigenvalues[J[czm[[i]]]]] == −1&&Max[Eigenvalues[J[czm[[i]]]]] > 1,Min[Eigenvalues[J[czm[[i]]]]] == −1&&Max[Eigenvalues[J[czm[[i]]]]] > 1,Min[Eigenvalues[J[czm[[i]]]]] == −1&&Max[Eigenvalues[J[czm[[i]]]]] > 1,Print[“The fixed point”, czm[[i]], “is unstable”],Print[“The fixed point”, czm[[i]], “is unstable”],Print[“The fixed point”, czm[[i]], “is unstable”],Max[Eigenvalues[J[czm[[i]]]]]==1&&Min[Eigenvalues[J[czm[[i]]]]] > −1,Max[Eigenvalues[J[czm[[i]]]]]==1&&Min[Eigenvalues[J[czm[[i]]]]] > −1,Max[Eigenvalues[J[czm[[i]]]]]==1&&Min[Eigenvalues[J[czm[[i]]]]] > −1,Print[“Try to use Center manifold module for ”, czm[[i]]],Print[“Try to use Center manifold module for ”, czm[[i]]],Print[“Try to use Center manifold module for ”, czm[[i]]],Max[Eigenvalues[J[czm[[i]]]]] > 1&&Min[Eigenvalues[J[czm[[i]]]]] < −1,Max[Eigenvalues[J[czm[[i]]]]] > 1&&Min[Eigenvalues[J[czm[[i]]]]] < −1,Max[Eigenvalues[J[czm[[i]]]]] > 1&&Min[Eigenvalues[J[czm[[i]]]]] < −1,Print[“The fixed point”, czm[[i]] , “ is oscillatory source”]]],Print[“The fixed point”, czm[[i]] , “ is oscillatory source”]]],Print[“The fixed point”, czm[[i]] , “ is oscillatory source”]]],Which[Abs[First[Eigenvalues[J[czm[[i]]]]]] < 1, Print[“The fixed point”, czm[[i]],Which[Abs[First[Eigenvalues[J[czm[[i]]]]]] < 1, Print[“The fixed point”, czm[[i]],Which[Abs[First[Eigenvalues[J[czm[[i]]]]]] < 1, Print[“The fixed point”, czm[[i]],“ is spiral sink”], Abs[First[Eigenvalues[J[czm[[i]]]]]]==1, Print[“ is spiral sink”], Abs[First[Eigenvalues[J[czm[[i]]]]]]==1, Print[“ is spiral sink”], Abs[First[Eigenvalues[J[czm[[i]]]]]]==1, Print[“The fixed point” , czm[[i]] , “is Center”], Abs[First[Eigenvalues[J[czm[[i]]]]]] > 1,“The fixed point” , czm[[i]] , “is Center”], Abs[First[Eigenvalues[J[czm[[i]]]]]] > 1,“The fixed point” , czm[[i]] , “is Center”], Abs[First[Eigenvalues[J[czm[[i]]]]]] > 1,Print[“The fixed point”, czm[[i]], “ is spiral source”]]], {i, 1, l}]]]Print[“The fixed point”, czm[[i]], “ is spiral source”]]], {i, 1, l}]]]Print[“The fixed point”, czm[[i]], “ is spiral source”]]], {i, 1, l}]]]
Improve this question
$\endgroup$
8
  • 4
    $\begingroup$ I don't know why it doesnt work because there is no such command in Mathematica. May be the authors used a package or their own code? Try to look at book/paper website and see if they have the code there. I see now they have the code listed in the paper itself. May be you can copy it from the pdf? $\endgroup$
    – Nasser
    Jul 21, 2022 at 0:54
  • 4
    $\begingroup$ The code for TwoDimensionalStability is found in that paper. Please include it in your question. $\endgroup$
    – JimB
    Jul 21, 2022 at 3:09
  • 1
    $\begingroup$ Good. Please post that copy into your question. $\endgroup$
    – JimB
    Jul 21, 2022 at 13:55
  • 3
    $\begingroup$ The copied code has a lot of typos so it doesn't work, including using curly quotes “ instead of straight quotes ". But one thing I notice in your usage is that xy should be x*y or at least x y (with a space) -- otherwise Mathematica thinks xy is a single variable. $\endgroup$
    – Chris K
    Jul 21, 2022 at 15:30
  • 2
    $\begingroup$ The "typos" are from the way text is stored in a PDF document. There are curly quotes (as @ChrisK has pointed out), all of the underscores are missing, etc. But it is up to you to fix all of that. (Or write the authors for a clean copy of the code.) $\endgroup$
    – JimB
    Jul 21, 2022 at 15:55

1 Answer 1

Reset to default
11
$\begingroup$

The code in the article you cited was poorly formatted, poorly indented, and generally not at all readable. However, what it did was relatively simple, so I rewrote it to be more readable and perhaps more efficient, e.g. avoiding multiple calls to Eigenvalues.

I used a separate evaluator function to encapsulate the decision logic, which I find more readable and more maintainable. It looks to me that the branches of the decision tree are disjoint, so rather than the nested If functions, I opted for a Which expressions.

ClearAll[twoDimStab]
twoDimStab[f_, g_] := Module[
  {stablePoints, J, pointEvalPairs},

  stablePoints = {x, y} /. Solve[{f == x, g == y}, {x, y}];
  (*stablePoints=SolveValues[{f==x,g==y},{x,y}]*)
  
  J[{u_, v_}] := D[{f, g}, {{x, y}}] /. {x -> u, y -> v};

  pointEvalPairs = Transpose[{stablePoints, (SortBy[N]@*Eigenvalues@*J) /@ stablePoints}];
  iEvaluateStablePoint /@ pointEvalPairs
]

ClearAll[iEvaluateStablePoint]
iEvaluateStablePoint[{point_List, {lambda1_, lambda2_}}] :=
 Which[
   0 < lambda1 < 1 < lambda2, point -> "saddle",
   0 < lambda1 < lambda2 < 1, point -> "sink",
   1 < lambda1 < lambda2, point -> "source",
   lambda1 < -1 && -1 < lambda2 < 0, point -> "oscillatory saddle",
   lambda1 < -1 && lambda2 > 1, point -> "oscillatory source",
   -1 < lambda2 < 1 && -1 < lambda1 < 1, point -> "stable sink",
   AnyTrue[{lambda1, lambda2}, Im[#] != 0 && Abs[#] < 1 &], point -> "spiral sink",
   AnyTrue[{lambda1, lambda2}, Im[#] != 0 && Abs[#] > 1 &], point -> "spiral source",
   AnyTrue[{lambda1, lambda2}, Im[#] != 0 && Abs[#] == 1 &], point -> "center"
 ]

With that code in hand, I tried the example reported in the paper:

twoDimStab[x^2 + 3 y x, y^2 - 2 y x]

(* Out: 
{{-(2/7), 3/7} -> "saddle", {0, 1} -> "source", 
 {1, 0} -> "oscillatory source", {0, 0} -> "stable sink"}
*)

Note that the decision tree shown in the paper (see below) was missing the "stable sink" case, as kindly pointed out by MichaelM2, so their textual description is not completely representative of what the function really does in code.

classification conditions in plain text from paper

Finally, if you have access to recent versions of MMA, you can replace the stablePoints = {x, y} /. Solve[{f == x, g == y}, {x, y}]; with the equivalent and more succinct stablePoints = SolveValues[{f == x, g == y}, {x, y}] that is shown in comments in the code above. However, I left the traditional version because I think it ensures far more backwards compatibility for the code.

Improve this answer
$\endgroup$
8
  • $\begingroup$ You might want SortBy[N] or NumericalSort on the eigenvalues. (E.g., consider SortBy[{Sqrt[3], 2, Sqrt[5]}]). $\endgroup$
    – Michael E2
    Jul 21, 2022 at 17:43
  • 3
    $\begingroup$ And {0, 0} is missed because the text of the paper left out a case in the code. You which statement needs something like -1 < lambda2 < 1 && -1 < lambda1 < 1, point -> "stable sink", added to it. $\endgroup$
    – Michael E2
    Jul 21, 2022 at 17:49
  • $\begingroup$ @MichaelE2 Thank you for fishing out the missing case! I've added it to the Which statement. I also switched to SortBy[N], which is a very good point. In fact, I got bitten by that before, but I still forget about it. $\endgroup$
    – MarcoB
    Jul 21, 2022 at 20:19
  • 1
    $\begingroup$ @MarcoB: Why does it fail for cases like twoDimStab[y, 2 x + y]? SP shows StreamPlot[{y, 2 x + y}, {x, -4, 4}, {y, -4, 4}] $\endgroup$
    – Moo
    Jul 24, 2022 at 2:41
  • 2
    $\begingroup$ @Moo The eigenvalues in that case are {-1, 2}, which puts the equilibrium exactly on the border between "oscillatory saddle" and another missing case (-1 < lambda1 < 0 && lambda2 > 1). BTW, the StreamPlot is misleading, because twoDimStab is for discrete-time systems but StreamPlot is appropriate for continuous-time differential equations. $\endgroup$
    – Chris K
    Jul 24, 2022 at 17:27

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.