5

When I use the following system

model = {x'[t] == x[t] (1 - x[t]) - x[t] y[t], y'[t] == x[t] y[t] - y[t], x[0] == 0.5, y[0] == 0.5}

with the WhenEvent

perturb = WhenEvent[Mod[t, 1],
    {x[t] -> x[t], y[t] -> y[t]}
];

chosen to be trivial, I have what seems to be strange behavior to me, namely if I do the following integration

NDSolve[{model, perturb}, {x, y}, {t, 0, 10000}]

I get an error that the integration reached maximum steps at t==8670., even though nothing is changing. If I increase the MaxSteps I can get it to integrate. So my question is if this is expected? The longer you integrate a system when using a WhenEvent the larger the MaxSteps you will need? I can replicate this same behavior when I use initial conditions that are on a stable equilibrium value ... so truly nothing is changing.

J. M.'s missing motivation
  • 124,525
  • 11
  • 401
  • 574
Gabriel
  • 2,225
  • 18
  • 19

1 Answers1

10

The issue is that event detection itself adds plenty of overhead. We can see this by comparing this:

model = {x'[t] == x[t] (1 - x[t]) - x[t] y[t], 
  y'[t] == x[t] y[t] - y[t], x[0] == 0.5, y[0] == 0.5};
perturb = WhenEvent[Mod[t, 1], {x[t] -> x[t], y[t] -> y[t]}];
eventPoints = Reap[NDSolve[{model, perturb}, {x, y}, {t, 0, 10000},
  StepMonitor :> Sow[t]]][[2, 1]];
eventPoints // Length

(* Out: 18670 *)

To this:

nonEventPoints = Reap[NDSolve[model, {x, y}, {t, 0, 10000},
 StepMonitor :> Sow[t]]][[2, 1]];
nonEventPoints // Length

(* Out: 347 *)

We can even illustrate exactly why this would happen.

pts1 = {#, 0} & /@ Select[eventPoints, # < 6 &];
pts2  = {#, 1} & /@ Select[nonEventPoints, # < 6 &];
ListPlot[{pts1, pts2}, AspectRatio -> Automatic,
  Axes -> False]

enter image description here

We can see the points piled up at the events.

Mark McClure
  • 32,469
  • 3
  • 103
  • 161
  • I see ... so since I am doing a perturbation at each time step, then that is why I am having to grow MaxSteps for longer output ... as I using up my "bank" near each event. I was thinking more like this was analogous to a Fortran solver that does integration a "step" at a time ... not considering the how the adaptive step size was working. Thanks! – Gabriel Dec 12 '12 at 03:39
  • @NasserM.Abbasi I don't think that point of view is quite reasonable. The event is part of the specification of the problem, so you have to expect that changing it could change the computational complexity. Consider the event x[t]==0; there's no particular reason to expect that one of the interpolation points generated by the eventless equations will match that event exactly. If you want machine precision, you've got to back up and perform a FindRoot to figure out exactly when the event occurs. That's quite costly. – Mark McClure Dec 12 '12 at 19:03
  • @MarkMcClure true, but for cases that happen at regular intervals (ie Mod[t,1] etc) this is slightly counter intuitive until you understand your answer. In these cases the events are easy to detect ... but the way the solver is allocating effort is subtle since it is I imagine trying to ensure that possible local discontinuities are being handled correctly. – Gabriel Dec 12 '12 at 19:06
  • @Gabriel I certainly undertand the hope, even expectation, that event detection with Mod[t,1] might work better. I think the detection techniques are purely numerical, though, so while it's clear to us where they happen, they still generate a lot of extra work for the algorithm. – Mark McClure Dec 12 '12 at 19:11
  • @MarkMcClure don't get me wrong ... your point is well taken! It would just be nice to have an extra option that would tell the maximum steps taken between events ... that way the user doesn't need to have to do any hand calculations on the number of time steps, and the predicted number of events to figure out how many maximum steps need to be taken. This will make the calculations much slower ... but seems like it will be a common need for these types of problems. Again I think of it the way old Fortran codes would have worked ... each call you get the maximum steps so it naturally scales – Gabriel Dec 12 '12 at 19:16
  • @Gabriel Like much of V9, the event detection could certainly be refined - and probably will be! – Mark McClure Dec 12 '12 at 19:20
  • @MarkMcClure for sure. I love WhenEvent so much ... it could make me have to do a back flip each time I use it and it would still keep me coming back! – Gabriel Dec 12 '12 at 19:34
  • @NasserM.Abbasi weird ... guess it just keeps trying! – Gabriel Dec 12 '12 at 19:36
  • People seem to forget that event detection has to do rootfinding and interpolation in addition to the integration proper; it should be no surprise that detecting events interferes with the normal step size selection, and thus incurs a fair amount of overhead. – J. M.'s missing motivation Apr 21 '13 at 07:08