5
$\begingroup$

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.

Improve this question
$\endgroup$

1 Answer 1

Reset to default
10
$\begingroup$

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.

Improve this answer
$\endgroup$
9
  • $\begingroup$ 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! $\endgroup$ Commented Dec 12, 2012 at 3:39
  • $\begingroup$ @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. $\endgroup$ Commented Dec 12, 2012 at 19:03
  • $\begingroup$ @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. $\endgroup$ Commented Dec 12, 2012 at 19:06
  • $\begingroup$ @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. $\endgroup$ Commented Dec 12, 2012 at 19:11
  • $\begingroup$ @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 $\endgroup$ Commented Dec 12, 2012 at 19:16

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.