7
$\begingroup$

I have to numerically integrate an equation system and monitor the accumulating datapoints. For example, I fit a line to a subsample of the points and terminate via "StopIntegration" if fitted line has ~0 slope. However, when I define a regular sampling time for WhenEvent (Mod[t, 1]), the collected log values have huge errors.

ode = {
   x'[t] == -.2 x[t]^2 + 2 y[t],
   y'[t] == x[t] + .1 x[t]^2 - 1.5 y[t],
   x[0] == y[0] == 1,
   WhenEvent[Mod[t, 1] == 0, 
    If[140 < t < 200, AppendTo[events, {t, Log@x@t}]]]
   };

events = steps = {};
if = NDSolveValue[ode, {x, y}, {t, 0, 200}, StepMonitor :> AppendTo[steps, {t, Log@x@t}]];

LogPlot[Through@if@t, {t, 0, 200}, Epilog -> {
   {Blue, AbsolutePointSize@6, Point@steps},
   {Red, AbsolutePointSize@3, Point@events}
   }, ImageSize -> 500]

Mathematica graphics

I can see two ways to overcome this. 1) Set a method for NDSolve with uniform known stepsize and make the WhenEvent test at the same times. For certain reasons I want to avoid this. For one, it's nice to rely on the automatic adaptive stepsize algorithm. 2) Force WhenEvent to only evaluate when an internal step is taken by the integrator without modifying the step size manually. I failed to achieve this reliably. How to do this?

As you can see, it's perfectly unnecessary to increase accuracy/precision, as at steps the solution values are correct and I don't need to have extra values in between them. Note, that I cannot use StepMonitor solely, as I have to stop integration when a condition is met (based on the collected points).

Details

The problem is in the way Mathematica calculates solution values but not via the actual integrator method: when a value is calculated, WhenEvent interface does not know where the next point will be, so it has to extrapolate. It just does it in a horrible way. Somewhere in the process the actual derivative is lost, or even worse, a constant derivative is assumed, resulting in wild errors. See the following, even simpler example:

Block[{events = {}, steps = {}, if, x, t},
   if = NDSolveValue[{x'[t] == -1/10 x[t], x[0] == 1, 
      WhenEvent[Mod[t, #] == 0., AppendTo[events, {t, x[t]}]]}, 
     x, {t, 0, 100},
     StepMonitor :> {AppendTo[steps, t]}, 
     Method -> "ExplicitRungeKutta"(*,WorkingPrecision -> 50*)];
   Plot[if@t, {t, 0, 100}, 
    Epilog -> {Green, Point@events}, 
    PlotRange -> {All, {1, -1}}, GridLines -> {steps, {0}}, 
    PlotLabel -> Row@{"d = ", #}]
   ] & /@ {10, 1, 1/10, 1/100}

Mathematica graphics

Vertical gridlines indicate internal steps, green points are the detected events. If you uncomment the option WorkingPrecision -> 50, the result is what one would expect - at the cost of at least double time required and I had to convert all reals to exact numbers.

Update

According to TechSupport, there is no way at the moment to evaluate a WhenEvent test only when the integrator takes a step. Maybe in a future version.

$\endgroup$
  • $\begingroup$ Just a note, the Log doesn't seem to matter - same problem without it. $\endgroup$ – Chris K Mar 15 '16 at 0:40
  • $\begingroup$ Is there a reason you want to fit a line to points to estimate the slope rather than using something like WhenEvent[Abs[x'[t]] < 10^-4, "StopIntegration"]? $\endgroup$ – Chris K Mar 15 '16 at 0:44
  • 1
    $\begingroup$ @ChrisK Sure: you can fit a line to a large sample of points, but the derivative is checked at a point. For certain cases, I want to check whether the last 10000 datapoints are equilibrated (slope ~ 0) or in exponential decay (slope is of log.values is negative). Furthermore, my solution could heavily oscillate, while still could be a limit cycle. $\endgroup$ – István Zachar Mar 15 '16 at 8:17
  • $\begingroup$ I see, thanks! That sounds generally useful, I'd be happy if you posted the final product. $\endgroup$ – Chris K Mar 15 '16 at 11:42
  • $\begingroup$ @MichaelE2 I am aware of this method but try to collect more than one point this way. The problem is that whenever the test of the WhenEvent becomes true, it won't check that test any more, even if you reset the flag. $\endgroup$ – István Zachar Mar 16 '16 at 9:22
5
$\begingroup$

This is interesting even if it's not a perfect answer. It makes me suspect a bug somewhere.

I thought, for reasons that elude me -- it just occurred to me to try -- of enforcing a certain checking of x[t] by adding a "trivial" discrete variable and updating it from time to time. It is trivial in the sense that it never changes value and does not affect the OP's system. I guess in part I thought the variable would force NDSolve to track it, and perhaps x[t] would be calculated instead of extrapolated (for it seemed that extrapolation might be the explanation of the bad values for x[t] in the OP's code). So I changed the first DE to

x'[t] == -.2 x[t]^2 + 2 y[t] + a[t]

where a[t] == 0. I was able to get what I wanted by changing the event action

WhenEvent[Mod[t, 1] == 0, 
 If[140 < t < 200, AppendTo[events, {t, Log@x@t}], a[t] -> 0]]

But that was a typo. I meant to try the following which did not work:

WhenEvent[Mod[t, 1] == 0, 
 If[140 < t < 200, AppendTo[events, {t, Log@x@t}]; a[t] -> 0]]

I couldn't (and still cannot) explain why one worked and the other failed. I got the following, respectively:

Mathematica graphics

Then I thought, what if the discrete variable a[t] was update periodically, but not very often, so that the OP's goal of letting the step size grow somewhat large might might be realized. This is, admittedly, similar to @user21's approach with MaxStepSize, but that approach still had "extrapolation" errors, albeit much smaller. This is unlike what a[t] -> 0 produced above (left image):

Mathematica graphics

According to my speculative theory, an event of the form WhenEvent[Mod[t, dt] == 0, a[t] -> 0], might cause x[t] to be assigned the correct value. (I was beginning to think that extrapolation is not the reason for the wild values in the OP's code.) With dt = 10 or 20 everything seems to work fine. The value chosen for dt seemed to affect the step size (because NDSolve is worried that some discontinuous change to the ODE is going to occur at such an event). So I thought, how large can I make dt? I found that around 85, strange things happen.

ode2 = {x'[t] == -.2 x[t]^2 + 2 y[t] + a[t], 
   y'[t] == x[t] + .1 x[t]^2 - 1.5 y[t],
   a[0] == 0, x[0] == y[0] == 1, 
   WhenEvent[Mod[t, 1] == 0, 
    If[140 < t < 200, AppendTo[events, {t, Log@x@t}](*,a[t]->0*)]],
   WhenEvent[Mod[t, dt] == 0, a[t] -> 0]};

Table[
   events = steps = {};
   if = NDSolveValue[ode2, {x, y}, {t, 0, 200}, 
     StepMonitor :> AppendTo[steps, {t, Log@x@t}], 
     DiscreteVariables -> {a}];
   LogPlot[Through@if@t, {t, 0, 200}, 
    Epilog -> {{Blue, AbsolutePointSize@6, Point@steps}, {Red, 
       AbsolutePointSize@3, Point@events}}, ImageSize -> 250, 
    PlotLabel -> HoldForm[dt] == dt],
   {dt, 85, 88}
   ] // Partition[#, 2] & // GraphicsGrid

Mathematica graphics

Three of the plots are have no error message; the upper right one has this error:

Coordinate {171., Complex[0.38451996298641833`, 3.141592653589793]} should be a pair of numbers, or a Scaled or Offset form.

The error just means that the value of some x[t] are negative when dt = 86, which in not inconsistent with the graph of the solution. For dt = 85 or 88, the values of event seem accurate, but in the case of dt = 88, there's a little glitch in the solution. As I said, it appears to be buggy (unless someone can explain why this is to be expected).

In any case, adding a[t] -> 0 in

WhenEvent[Mod[t, 1] == 0, 
 If[140 < t < 200, AppendTo[events, {t, Log@x@t}], a[t] -> 0]]

seems to be a workaround, even though I cannot explain why it works. None of these events, it seems to me, should affect the computed solutions to x[t] and y[t]. But they do.

$\endgroup$
  • $\begingroup$ Interesting findings, thanks for the analysis! If we can close i on the issue, I'll file it for TechSupport. $\endgroup$ – István Zachar Mar 15 '16 at 8:28
  • $\begingroup$ You'd have to use MaxStepSize->1; Mod[t,dt] in a WhenEvent is documented to sample at dt. Maybe I am missing the point? $\endgroup$ – user21 Mar 15 '16 at 10:26
  • $\begingroup$ @user21 Doesn't Mod[t,dt] add a step size every dt? I don't know how a Mod event affects the adaptive step size. I do notice that when the event action is not null, NDSolve does take a step at the event time; however, it takes other steps in between, if the time between events is long enough. I know I'm confused about something. I'm surprised that changing the timing of an event that seems to do nothing affects the quality of the solution. I assume the event action a[t] -> 0 is treated as a discontinuity and that in turn interferes with the next step in an obscure way. $\endgroup$ – Michael E2 Mar 15 '16 at 12:32
  • $\begingroup$ That is the point: specifying a high-resolution event-condition does not (by default) increases the resolution of the adaptive stepsize algorithm. It either should, or there should be a way to only test event conditions when a step is actually taken. $\endgroup$ – István Zachar Mar 15 '16 at 18:12
  • $\begingroup$ @IstvánZachar You can force NDSolve to take a step at each event by having the WhenEvent action return "RestartIntegration" (but I thought you wanted to avoid the extra sampling). $\endgroup$ – Michael E2 Mar 15 '16 at 18:19
3
$\begingroup$

One way to do this is not to restrict NDSolve to use a fixed size method but to restrict the MaxStepSize

if = NDSolveValue[ode, {x, y}, {t, 0, 200}, 
   StepMonitor :> AppendTo[steps, {t, Log@x@t}], MaxStepSize -> 2];

So if you use WhenEvent[Mod[t, dt] == 0.... I'd think you want MaxStepSize -> dt

From the documentation of WhenEvent: Mod[t,dt]==0 sample at regular intervals dt in the time variable t. But perhaps I miss the point?

$\endgroup$
  • $\begingroup$ I appreciate the effort, but MaxStepSize explicitly discards the adaptive step selection method and unnecessarily multiplies step points. Something I rather avoid as I wrote. Emphasized it in my edit. $\endgroup$ – István Zachar Mar 14 '16 at 23:43
  • 2
    $\begingroup$ @IstvánZachar, it does not discard the adaptive step size method, it limits it. That's something different. $\endgroup$ – user21 Mar 14 '16 at 23:55
  • $\begingroup$ You are right, but what I want is not to restrict the default step method in any way, if possible and not to collect too many datapoints. $\endgroup$ – István Zachar Mar 15 '16 at 8:19
3
$\begingroup$

Here is a potential workaround, that works on the OP's MWE. We'll restart integration at the beginning of logging of events. Since I'm not sure why the OP's code fails, I'm not sure why this fixes it. It potentially adds a step at t == 140, and it causes NDSolve to calculate a new starting step size. It's hard to know whether this would slow things down in all cases, but it has a negligible effect on the number of steps in the OP's MWE. I'll log the steps in steps2 and events2 to compare with the OP.

ode = {x'[t] == -.2 x[t]^2 + 2 y[t], 
   y'[t] == x[t] + .1 x[t]^2 - 1.5 y[t], x[0] == y[0] == 1, 
   WhenEvent[Mod[t, 1] == 0, If[140 < t < 200, AppendTo[events2, {t, Log@x@t}]]], 
   WhenEvent[t == 140, "RestartIntegration"]};

events2 = steps2 = {};
if = NDSolveValue[ode, {x, y}, {t, 0, 200}, 
   StepMonitor :> AppendTo[steps2, {t, Log@x@t}]];

LogPlot[Through@if@t, {t, 0, 200}, 
 Epilog -> {{Blue, AbsolutePointSize@6, Point@steps2}, {Red, 
    AbsolutePointSize@3, Point@events2}}]

Mathematica graphics

The steps are about the same, with an adjustment at the restart of integration.

NumberLinePlot[{steps[[All, 1]], steps2[[All, 1]]}, 
 GridLines -> {{140}, None}]

Mathematica graphics

There is a small error between the computed solution and the event values:

With[{xif = First@if},
 maxerr = Max@Abs[Differences /@ MapAt[Log @* xif, events, {All, 1}]]
 ]
(*  3.81291*10^-10  *)

We can see that this does not really fix the error observed in the OP:

Show[
 plot = LogPlot[First[if][t], {t, 139, 200}], 
 Graphics[{{Blue, AbsolutePointSize@6, Point@steps2}, {Red, 
    AbsolutePointSize@3, Point@events2}}],
 PlotRange -> {{139, 200}, Last@PlotRange@plot + maxerr {-1, 1}}]

Mathematica graphics

Note that the initial slopes of the event arcs are close to the derivative of x[t] at t == 140, where integration (re)started. A similar pattern may be observed in the OP. This is not strictly true, nor does it seem that each arc starts with the same rate of change.

SplitBy[Differences@Exp@events2[[All, 2]], Positive][[1 ;; ;; 2, 1]]
First[if]'[140]
(*
  {6.90459*10^-10, 6.54047*10^-10, 6.05088*10^-10, 6.11475*10^-10}
  7.31003*10^-10
*)

If the arcs are connected to x'[140], then this workaround may only work accidentally, as it were, in the OP's MWE.

$\endgroup$
  • $\begingroup$ Unfortunately, restarting only works if the function does not change much afterwards. I've quite a few cases where restart does not work - especially, if I include the whole range for the WhenEvent and not just after t=140. $\endgroup$ – István Zachar Mar 15 '16 at 17:24
  • $\begingroup$ @IstvánZachar I was afraid of that, given my last figure. I think that NDSolve is interpolating between steps to get x[t] at the events, but it appears that it has botched the interpolation formula (a somewhat wild guess). $\endgroup$ – Michael E2 Mar 15 '16 at 17:51

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.