3

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!

X Bruno
  • 143
  • 4

1 Answers1

9

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]

enter image description here

george2079
  • 38,913
  • 1
  • 43
  • 110