3
$\begingroup$

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.

$\endgroup$
  • 1
    $\begingroup$ 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! $\endgroup$ – Chris K Jan 18 '16 at 23:30
  • $\begingroup$ Sorry for the mistake before! Yes, I would like to find a limit cycle with period 1. $\endgroup$ – X Bruno Jan 19 '16 at 1:12
1
$\begingroup$

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.

|improve this answer|||||
$\endgroup$
0
$\begingroup$

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}]

enter image description here

|improve this answer|||||
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.