5
$\begingroup$

Bug introduced in 9.0 and fixed in 11.1


NDSolve in Mathematica 9.0.0 (MacOS) is behaving strangely with a piecewise right hand side. The following code (a simplified version of my real problem):

sol = NDSolve[{x'[t] == 
        Piecewise[{{2, 0 <= Mod[t, 1] < 0.5}, 
                   {-1, 0.5 <= Mod[t, 1] < 1}}
        ], x[0] == 0}, x, {t, 0, 1}];
Print[x[1] /. sol[[1]]];

gives the correct answer of 0.5 about 50% of the time, but often returns -0.5 and -1 instead. Rerunning it gives apparently random results. It always gives the correct result in Mathematica 8.

Here's what I've figured out so far:

  1. It apparently has something to do with the Mod[t,1], because it works fine with just "t" in the Piecewise. Unfortunately I'm looking at a piecewise periodic system (not just from t=0 to 1).
  2. It's only the first segment of the solution from t=0 to t=0.5 that varies from run to run.
  3. Using initial condition x[10^-100]==0 fixes the problem, but this is an ugly hack.

Can anyone replicate this strange behavior, know what's behind it, or have a better suggested fix?

$\endgroup$
  • $\begingroup$ I get the same strange behavior. $\endgroup$ – chris Mar 13 '13 at 21:03
  • $\begingroup$ Same here Win7-64 v9.0.1. In v8 no such thing indeed. I assume it's a bug and suggest you contact wolfram support. $\endgroup$ – Sjoerd C. de Vries Mar 13 '13 at 21:54
  • $\begingroup$ I filed this as a bug. Thanks. $\endgroup$ – user21 Mar 14 '13 at 10:19
  • 1
    $\begingroup$ @ruebenko: Have you seen ChrisK's finding about the "DiscontinuityProcessing" method option? It probably might help to track down the bug somewhat faster... $\endgroup$ – Albert Retey Mar 14 '13 at 14:04
  • 1
    $\begingroup$ @AlbertRetey, yes this seems to be an issue with the new in V9 discontinuity processing. Thanks for letting me know. $\endgroup$ – user21 Mar 14 '13 at 19:37
7
$\begingroup$

Following a lead from Albert Retey, I found an NDSolve option that fixes this problem:

sol = NDSolve[{x'[t] == 
    Piecewise[{{2, 0 <= Mod[t, 1] < 0.5}, 
               {-1, 0.5 <= Mod[t, 1] < 1}}
    ], x[0] == 0}, x, {t, 0, 1}, Method -> {"DiscontinuityProcessing" -> False}];
Print[x[1] /. sol[[1]]];
$\endgroup$
4
$\begingroup$

This is not really an answer (the answer is of course that this is a bug), but it is too long for a comment and probably gives a hint where the problem is and how one can avoid it in other cases. The workaround is to define the piecewise function as an external definition only for numeric arguments. It looks like otherwise some invalid optimization with the symbolic expression is done:

ClearAll@rhs
rhs[t_?NumericQ] := Piecewise[{
   {2, 0 <= Mod[t, 1] < 0.5},
   {-1, 0.5 <= Mod[t, 1] < 1}
   }]

Table[
 sol = NDSolve[{x'[t] == rhs[t], x[0] == 0},
   x, {t, 0, 1}
   ];
 Plot[x[t] /. sol[[1]], {t, 0, 1}, PlotLabel -> (x[1] /. sol[[1]]), 
  Frame -> True, PlotRange -> All],
 {10}
 ]

I have tested this with Mathematica 9.0.1 on Windows 7 64bit. Unlike for NIntegrate there seems to not be a Method option with which one could switch off the symbolic optimization, at least I didn't find one. It might well be there, even if not documented and of course might be a better workaround...

$\endgroup$
  • 3
    $\begingroup$ Thanks for the idea of looking for a Method option. On a hunch I tried adding Method -> {"DiscontinuityProcessing" -> False}, which seems to do the trick. $\endgroup$ – Chris K Mar 14 '13 at 13:48
  • 1
    $\begingroup$ @ChrisK: good find, and the name indicates a plausible source of errors for this case. I'm not sure whether you are familiar with the site: it is welcomed to give answers to your own questions and accept them. I would suggest you should do that in this case... $\endgroup$ – Albert Retey Mar 14 '13 at 14:03
3
$\begingroup$

Another possible fix to the bug is to use Simplify`PWToUnitStep to expand the Piecewise into a combination of UnitStep:

Table[NDSolveValue[{x'[t] == 
      Simplify`PWToUnitStep@
       Piecewise[{{2, 0 <= Mod[t, 1] < 0.5}, {-1, 0.5 <= Mod[t, 1] < 1}}], x[0] == 0}, 
    x, {t, 0, 1}][1], {100}] // Union

Mathematica graphics

$\endgroup$
3
$\begingroup$

Taking the OP's clue that starting at x[10^-100] == 0 solves the problem, we can try explicitly setting the derivative value at t == 0. The following produces a consistent and correct result:

solME2 = First@ NDSolve[{
   x'[t] == Piecewise[{
      {2, t == 0},
      {2, 0 <= Mod[t, 1] < 0.5},
      {-1, 0.5 <= Mod[t, 1] < 1}}],
   x[0] == 0}, x, {t, 0, 1}]

ListLinePlot[x /. solME2, Mesh -> All]

As an aside, discontinuities are processed as events, and generally an event at the initial condition is problematic. Even if this observation is related to the issue, it's is still far from explaining the lack of determinacy in the computation. It's certainly a bug, and a full explanation may be elusive.

Note also that turning off discontinuity processing, either with the option Method -> {"DiscontinuityProcessing" -> False} or by using ?NumericQ-protected function for the right-hand side, has its own issues:

sol10 = NDSolve[{x'[t] == Piecewise[{
       {2, 0 <= Mod[t, 1] < 0.5},
       {-1, 0.5 <= Mod[t, 1] < 1}}],
    x[0] == 0}, x, {t, 0, 10},
   Method -> {"DiscontinuityProcessing" -> False}];
ListLinePlot[x /. sol10[[1]], Mesh -> All]

Mathematica graphics

The maximum step size is too large and the discontinuities are jumped over. You should be prepared to manually intervene in such cases.

The OP's code also leaks variables of the form sNNN from Unique["s"], like another recent question, Generation of global variables when using NDSolveValue and Piecewise function.

$\endgroup$
  • $\begingroup$ Thanks for the thoughts. Another time I ran into the "events at the initial conditions are problematic" is here. Also, that ListLinePlot[x /. solME2, Mesh -> All] trick is slick! $\endgroup$ – Chris K Nov 2 '16 at 14:15

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

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