Are there practical Mathematica tools for detection of limit cycles in two dimensional dynamical systems?

Question

I have a dynamical system with one boundary saddle point, and one unstable interior point, and I would like to detect existence of cycles. This is of course a hard problem if the cycle is very small, and I can only hope for partial answers. My first attempt is a program with two inputs:the dynamical system, and the time T to iterate--- see below. I I use a very crude idea: starting paths near all the fixed points, and enlarging the time hoping to observe a cycle. If nothing else, this will detect separatices -- see the posts Plotting separatrices for nonlinear system , How to plot the stable and unstable manifolds of a hyperbolic fixed point of a nonlinear system of differential equations? The example below revealed a big cycle, but I have others from the literature where this crude program finds nothing (and contradicts the literature). So, my question was: if you needed to check quickly many numerics offered in the literature without accompanying notebooks (the typical case, unfortunately), what macro would you write?

T = 380; dyn = {16 - s/10 - (i s)/(100 (1 + i/1000)), -((21 i)/50) - (
       13 i)/(2 (13/2 + i)) + (i s)/(100 (1 + i/1000))}; vars = {s, i};
    equi = Solve[Thread[dyn == 0 && s >= 0 && i >= 0], vars] // N; le = 
     Length[equi];
    Jacob = Grad[dyn, vars];
    varst = Through[vars[t]];
    diff = D[varst, t] \[Minus] (dyn /. Thread[vars -> varst]);
    Print["Fixed points:"]
    Do[p[j] = vars /. equi[[j]];
     Jac[j] = Jacob /. Thread[vars -> p[j]]; 
     Print[p[j], " with eigenvalues", Eigenvalues[Jac[j]]];
     e[j] = Text["E[j]", Offset[{10, 11}, p[j]]]; 
     st[j] = {PointSize[Large], Style[Point[p[j]], Green]}; 
     b[j] = p[j] - {.1, 0}; a[j] = p[j] + {.1, 0}; 
     eqb[j] = Join[Thread[diff == 0], Thread[{s[0], i[0]} == b[j]]]; 
     eqa[j] = Join[Thread[diff == 0], Thread[{s[0], i[0]} == a[j]]];
     eqb[j] = Join[Thread[diff == 0], Thread[{s[0], i[0]} == b[j]]]; 
     eqa[j] = Join[Thread[diff == 0], Thread[{s[0], i[0]} == a[j]]]; 
     sob[j] = NDSolve[eqb[j], vars, {t, 0, T}];
     pb[j] = 
      ParametricPlot[{s[t], i[t]} /. sob[j], {t, 0, T}, 
        PlotStyle -> Green] /. 
       Line[x_] :> {Arrowheads[Table[.03, {5}]], Arrow[x]};
     soa[j] = NDSolve[eqa[j], vars, {t, 0, T}];
     pa[j] = 
      ParametricPlot[{s[t], i[t]} /. soa[j], {t, 0, T}, 
        PlotStyle -> Red] /. 
       Line[x_] :> {Arrowheads[Table[.05, {5}]], Arrow[x]}, {j, le}]
    X = Table[p[j][[1]], {j, le}]; Y = Table[p[j][[2]], {j, le}]; xm = 
     3 Min[X]/4; xM = 99 Max[X]/98; ym = 0; yM = 4 Max[Y]/3;
    epi = Flatten[Table[{e[j], st[j]}, {j, le}]];
    sp = StreamPlot[{dyn}, {s, xm, xM}, {i, ym, yM}, StreamPoints -> 500, 
       Epilog -> epi, Frame -> True,
        FrameLabel -> {"s", "i"}, 
       LabelStyle -> Directive[Black, Medium]];
    Show[sp, Table[pb[j], {j, le}], Table[pa[j], {j, le}]]

But, I accept Chris's implied answer that this problem is too difficult for one macro, and the quickest way is to use his package, which seems to provide even the period, thanks @Chris K

Answer 1

In general, I'm not sure there's a good algorithm to find all limit cycles for a given set of equations. But if there's one in particular you want to find, then it's not too hard with a decent initial guess. This answer will be similar to this one, using my EcoEvo package.

First, you need to install the package (once):

PacletInstall["EcoEvo", "Site" -> "http://raw.githubusercontent.com/cklausme/EcoEvo/master"]

Then, load the package and define your model:

<< EcoEvo`;

SetModel[{
   Pop[s] -> {Equation :> 16 - s/10 - (i s)/(100 (1 + i/1000))}, 
   Pop[i] -> {Equation :> -((21 i)/50) - (13 i)/(2 (13/2 + i)) + (i s)/(100 (1 + i/1000))}
}];

(In general I'd recommend not hard-wiring in parameter values, but for sake of simplicity, I followed your code).

Find valid equilibria (Pop's are assumed to be non-negative):

eq = SelectValid[NSolveEcoEq[]]
{{s -> 160., i -> 0}, {s -> 84.878, i -> 8.92962}}

and check their linear stability:

EcoEigenvalues[eq]
{{0.18, -0.1}, {0.0239248 + 0.169326 I, 0.0239248 - 0.169326 I}}

Both are unstable, as you found. Now make the phase plane and indicate equilibria:

pp = Show[
  PlotEcoPhasePlane[{s, 0, 200}, {i, 0, 40}, PlotPoints -> 200],
  RuleListPlot[eq, PlotMarkers -> EcoStableQ[eq]]
]

Simulate the dynamics to get close to the limit cycle:

sol = EcoSim[{s -> 85, i -> 8}, 1000];
PlotDynamics[sol]

Finally, use the final time slice of sol as an initial guess to FindEcoCycle to find your limit cycle:

lc = FindEcoCycle[FinalSlice[sol]]
Show[pp, RuleListPlot[lc, {s, i}, PlotStyle -> Pink]]

The period is included in lc:

FinalTime[lc]
(* 54.8005 *)

and you can check its dynamical stability by calculating Floquet exponents:

EcoEigenvalues[lc]
(* {-4.74648*10^-8, -0.0496394} *)

As for methods and options, you can look up FindEcoCycle in Mathematica's built-in documentation once you've installed the package, otherwise read online here. There are two methods, one based on NDSolve and WhenEvents, the other based on FindRoot (which requires an initial guess of the period):

lc = FindEcoCycle[FinalSlice[sol], Method -> "FindRoot", Period -> 50];

Also, you can download the package's source code here if you want to look inside.