Detecting existence of period orbit

Question

I wish to detect the existence of a periodic orbit (limit cycle). For example, the test should return FALSE (i.e. it should fail) for the left (m = 0.35) and return TRUE for the right (m = 0.3) in the following toy model.

sys[m_] := {
  x'[t] == x[t] - 0.5 x[t]^2 - x[t] y[t] (x[t] - 1)  / (x[t]^2 - 1),
  y'[t] == -m y[t] +  x[t] y[t] (x[t] - 1) / (x[t]^2 - 1),
  x[0] == 0.1,
  y[0] == 0.1}
T = 500;
m = 0.35;
sol = First[NDSolve[{sys[m]}, {x, y}, {t, 0, T}]];
p = Plot[Evaluate[{x[t], y[t]} /. sol], {t, 0, T}];
m = 0.3;
sol2 = First[NDSolve[{sys[m]}, {x, y}, {t, 0, T}]];
p2 = Plot[Evaluate[{x[t], y[t]} /. sol2], {t, 0, T}];
GraphicsRow[{p, p2}]

Note that I am not interested in performing local stability analysis on the system (with which the Hopf bifurcation of the above toy model could be discerned), but rather wish to assess the existence of an orbit from the numerical solution itself.

It seems that using a Poincaré section and WhenEvent is the way to go, but I have not figured out a way to implement it that is robust to short transients (when x'[t]==0 is true for all t in the evaluated time period, an empty solution is returned) or very long transients, as demonstrated in the following by changing the period of time (Tstart to Tend) that is assessed:

test[m_] :=
 Module[{res},
  res = Reap[nd[{m}] =
      NDSolve[{sys[m], 
        WhenEvent[x'[t] == 0 && t > Tstart, 
         Sow[{m, Round[y[t], 10^-5]}]]},
       {x[t], y[t]},
       {t, 0, Tend},
       MaxSteps -> \[Infinity]]][[-1, 1]];
  CountDistinct[res] == 2
  ]

(* "Correct" inference *)
Tstart = 4500;
Tend = 5000;
test[0.35]
test[0.3]
Out[2916]= False
Out[2917]= True

(* False negative for m=0.3; still in transient phase *)

Tstart = 450;
Tend = 500;
test[0.35]
test[0.3]

Out[2920]= False

Out[2921]= False

(* non-existent part for m=0.35 when past transient phase *)

Tstart = 45000;
Tend = 50000;
test[0.35]
test[0.3]

During evaluation of In[2922]:= Part::partw: Part 1 of {} does not exist.

During evaluation of In[2922]:= Part::partw: Part 1 of {} does not exist.

Out[2924]= False

Out[2925]= True

I have not managed to make sense of How to locate the position of a periodic orbit and Locating periodic orbits with Mathematica.

Answer 1

Try

{X, Y} = NDSolveValue[{sys[m]}, {x, y}, {t, 0, T}]

Assuming periodoc solution near t~T the period tp follows with

gz=NMinimize[{  (X[T] - X[T-tp ])^2 + (Y[T] - Y[T-tp ])^2, 1 < tp < 20  }, tp ]
(* {7.54119*10^-12, {tp -> 16.2574}}*)
ParametricPlot[ {X[t], Y[t]}  , {t, T - tp /. gz[[2]], T}]

Hope it helps!

Answer 2

A few initial thoughts:

1. Since limit cycles are defined for $t\to\infty$, you'll never be quite sure with finite time. This will be especially problematic near a bifurcation point, where transients last longer.
2. You can simplify (x[t] - 1) / (x[t]^2 - 1) to 1 / (x[t] + 1).
3. With that simplification, this becomes the well-known Rosenzweig-MacArthur predator-prey model.

With those points out of the way, I can suggest that this will be easier if if you're willing to use my EcoEvo package. First, install it with:

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

Now load the package and set the model for analysis:

<< EcoEvo`
SetModel[{
  Aux[x] -> {Equation :> x - 0.5 x^2 - x y /(x + 1)},
  Aux[y] -> {Equation :> -m y + x y/(x + 1)},
  Parameters -> {m}
}]

Then you can find a limit cycle with FindEcoLimitCycle, and check its final derivatives to make sure it's not an equilibrium.

A limit cycle:

m = 0.3;
PlotEcoPhasePlane[{x, 0, 2}, {y, 0, 2}]
ec = FindEcoCycle[{x -> 0.1, y -> 0.1}];
PlotDynamics[ec]
FinalDerivatives[ec]

(* {x' -> -0.00828357, y' -> 0.179955} *)

An equilibrium:

m = 0.35;
PlotEcoPhasePlane[{x, 0, 2}, {y, 0, 2}]
ec = FindEcoCycle[{x -> 0.1, y -> 0.1}];
PlotDynamics[ec]
FinalDerivatives[ec]

(* {x' -> -4.6931*10^-7, y' -> 2.49523*10^-7} *)

I'm working on some new functionality to track limit cycles vs a parameter, but that's not quite ready.

Answer 3

First try to find the steady solution

dgl[m_] := {x'[t] ==x[t] - 1/2 x[t]^2 - x[t] y[t] (x[t] - 1)/(x[t]^2 - 1),y'[t] == -m y[t] + x[t] y[t] (x[t] - 1)/(x[t]^2 - 1) }
T = 500;
m = 0.35;
{X, Y} = NDSolveValue[{dgl[m], x[0] == 0.1, y[0] == 0.1 }, {x, y}, {t, 0, T}];
{X, Y} = NDSolveValue[{dgl[m], x[0] == X[T], y[0] == Y[T] }, {x, y}, {t, 0, T}];  (* steady solution*)

equilibrium

equi =  First@Solve[{0 == x[t] - 1/2 x[t]^2 - x[t] y[t] (x[t] - 1)/(x[t]^2 - 1),0 == -m y[t] + x[t] y[t] (x[t] - 1)/(x[t]^2 - 1) }, {x[t], y[t]}] /. f_[t] -> f; 
{x\[Infinity], y\[Infinity]} = {{x /. equi}, {y /. equi}}
(*{{0.538462}, {1.12426}}*)
Plot[ Norm[ {X[t] - x[Infinity], Y[t] - y[Infinity]} ]  , {t, 0, T},PlotLabel->"m="<>ToString[m]]

Result m==.3

m = 0.3 ;
{X, Y} = NDSolveValue[{dgl[m], x[0] == 0.1, y[0] == 0.1 }, {x, y}, {t, 0,T}];
{X, Y} = NDSolveValue[{dgl[m], x[0] == X[T], y[0] == Y[T] }, {x, y}, {t, 0,T}];  (* steady solution*)
equi =  First@
Solve[{0 == x[t] - 1/2 x[t]^2 - x[t] y[t] (x[t] - 1)/(x[t]^2 - 1),0 == -m y[t] + x[t] y[t] (x[t] - 1)/(x[t]^2 - 1) }, {x[t], y[t]}] /. f_[t] -> f; 
{x\[Infinity], y\[Infinity]} = {{x /. equi}, {y /. equi}}
(*{{0.428571}, {1.12245}}*)
Plot[ Norm[ {X[t] - x\[Infinity], Y[t] - y\[Infinity]} ]  , {t, 0,T},PlotLabel->"m="<>ToString[m]]

Answer 4

Although the suggestion made by Ulrich Neumann works for detecting orbits in general, it is not robust to "vanishing orbits", i.e. damped oscillations, as his subsequent suggestion points out. (My interest is in the detection of stable limit cycles only.)

Finite-time challenges aside, I tried some variations on his method of estimating the period (incl. allowing for multiple tp solutions and using FFT). However, in the end, modifying my originally-posted method (increasing T, testing for convergence to a fixed point using a lagged time comparison, increasing Accuracy and Precision goals) proved to get me closest to what I desire.

test2[m_] :=
 Module[{X, Y, res, T, Tlag, Diff, Tol},
  T = 150000;
  Tlag = 100;
  Tol = 10^-5;
{X, Y} = 
   NDSolveValue[{sys[m], x[0] == 0.1, y[0] == 0.1}, {x, y}, {t, 0, 
     T}, AccuracyGoal -> 10, PrecisionGoal -> 10];
Diff = Abs[X[T] - X[T - Tlag]];
  If[Diff < Tol,
   res = {1},
   res = 
    Reap[nd[{m}] = NDSolve[{{sys[m], x[0] == X[T], y[0] == Y[T]},
         WhenEvent[x'[t] == 0, Sow[{m, Round[y[t], Tol]}]]}, {x[t], 
         y[t]}, {t, 0, T}, AccuracyGoal -> 10, PrecisionGoal -> 10,
        MaxSteps -> [Infinity]]][[-1, 1]]
   ];
  CountDistinct[res] == 2
  ]