To find a limit cycle with period 1 in NDSolve

I am trying to find a limit cycle with period 1 of an ode system. The idea is to compare the values of the system during two periods, for example, x[t]=x[t-1], y[t]=y[t-1]. I am using Whenevent to do this, but I find that the event action can not be activated. Below is the code for this:

f[x_] := 0.1 + 0.05*Sin[2*Pi*(x + 1/3)];
xxold = yold = 0;
tmax = 200;
ressol = NDSolve[{xx'[t] == 0.2*xx[t]*y[t] - f[t]*xx[t],
y'[t] == (1 - y[t]) - 0.5*xx[t]*y[t], xx[0] == 0.5, y[0] == 0.5,
WhenEvent[
Mod[t, 1] == 0, {{xxold, yold} = {xx[t], y[t]},
If[Abs[xx[t] - xxold] + Abs[y[t] - yold] <
Power[10, -5], {{xxeq, yeq, teq} = {xx[t], y[t], t},
"StopIntegration"}]}]}, {xx, y}, {t, 0, tmax}]


Can anybody have a look at it for me? Thank you so much!

1 Answer

This works with a couple of subtle changes. Notable WhenEvent needs to return just the string "StopIntegration", not a list with that among other things. From the looks of it you are maybe trying to use {} for grouping things, where the curly brackets are only for lists in mathematica.

f[x_] := 0.1 + 0.05*Sin[2*Pi*(x + 1/3)];
xxold = yold = 0;
tmax = 400;
ressol = NDSolve[{xx'[t] == 0.2*xx[t]*y[t] - f[t]*xx[t],
y'[t] == (1 - y[t]) - 0.5*xx[t]*y[t], xx[0] == 0.5, y[0] == 0.5,
WhenEvent[Mod[t, 1] == 0,
If[Abs[xx[t] - xxold] + Abs[y[t] - yold] < Power[10, -5],
{xxeq, yeq, teq} = {xx[t], y[t], t};
"StopIntegration",
{xxold, yold} = {xx[t], y[t]}]]}, {xx, y}, {t,0, tmax}]

ParametricPlot[{(xx /. ressol[[1]])[t], (y /. ressol[[1]])[t]}, {t,
teq-1, teq}, AspectRatio -> 1/GoldenRatio]


• cool! Thank you so much! Feb 9, 2016 at 17:54