# track equilibrium of periodic ode system

I am trying to track equilibium in a periodid ode system. In such systems the equiblirium is defined as x[t]=x[t-1] and y[t]=y[t-1]. So it requires me to evaluate x[t-1] and y[t-1] when EventCondition takes effect within the NDSolve process.

For example, for such a system:

f[x_]:=0.1+0.05*Sin[2*Pi*(x+1/3)];
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}, {xx, y}, {t, 0, 200},
Method -> {"EventLocator", "Event" -> Abs[xx[t]-xx[t-1]]+Abs[y[t]-y[t-1]]<Power[10,-5],
"EventCondition" -> (t > 10),
"EventAction" :> {Print[ t, " ", xx[t], y[t]],
"StopIntegration"}}];
Plot[Evaluate[{xx[t], y[t]} /. ressol], {t, 0, 200}, Frame -> True,
PlotStyle -> {Red, Blue}]


But I get the errors: NDSolve::evre: The value of the event function at t = 9.777182689739677 was not a real number. The event will be ignored in steps where it does not evaluate to real numbers at both ends.The reason is xx[t-1] and y[t-1] can not be evaluated within the NDSolve process. Can anyone help me for this issue? Thank you.

• I'm a bit confused what you want to do -- what does "track equilibrium" mean? From your definition it sounds like you want to find a limit cycle with period 1, but your forcing function doesn't have period 1 (f[x] divides by 2pi instead of multiplying). Anyhow, from numerically solving it, it looks like xx[t] goes to zero rather than a limit cycle. Looks like a predator-prey system with periodic predator mortality. If you want periodic behavior out, you need to let the predator xx live first! – Chris K Jan 18 '16 at 23:30
• Sorry for the mistake before! Yes, I would like to find a limit cycle with period 1. – X Bruno Jan 19 '16 at 1:12

With your new f[x], I made the event Mod[t,1], with the action to compare the current state against the previous one saved in xxold and yold.

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}, {xx, y}, {t, 0, tmax},
Method -> {"EventLocator", "Event" -> Mod[t, 1],
"EventAction" :> {If[Abs[xx[t] - xxold] + Abs[y[t] - yold] < 10^-5,
Throw[{xxeq, yeq, teq} = {xx[t], y[t], t},
"StopIntegration"]]; xxold = xx[t]; yold = y[t]}}];

{xxeq, yeq, teq}

(* {1.99184, 0.499614, 192.} *)


You might want a more stringent condition than <10^-5, because it still seems to be changing a bit.

If you have Mathematica 9 or later, I think WhenEvent is easier to read.

I'm not certain this is correct but it seems promising. The idea is to placate NDSolve with a delay differential equation.

f[x_] := 1 + Sin[x/(2*Pi) + 1/3] + 1/2;
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], xxd'[t] == xx'[t - 1],
xxd[t /; t <= 0] == 0, yd'[t] == y'[t - 1], yd[t /; t <= 0] == 0,
xx[0] == 0.5, y[0] == 0.5}, {xx, y, xxd, yd}, {t, 0, 200},
Method -> {"EventLocator",
"Event" ->
Abs[xx[t] - xxd[t]] + Abs[y[t] - yd[t]] < Power[10, -5],
"EventCondition" -> (t > 10),
"EventAction" :> (Print[t, " ", xx[t], y[t]];
"StopIntegration")}];


We get this message:

NDSolve::ihist: Conditions given at t = 0. will be interpreted as initial history functions for t/;t<=0.. >>


Now plot the result.

Plot[Evaluate[{xx[t], y[t]} /. ressol], {t, 0, 200}, Frame -> True,
PlotStyle -> {Red, Blue}]
`