Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question
add comment

1 Answer 1

up vote 9 down vote accepted

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.

share|improve this answer
    
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 3: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
show 4 more comments

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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