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!
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]
Comments
Post a Comment