NDSolve with Piecewise gives the incorrect answer randomly

Question

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?

Answer 1

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]]];

Answer 2

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...

Answer 3

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

Answer 4

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]

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.