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:
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).
It's only the first segment of the solution from t=0 to t=0.5 that varies from run to run.
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?
=================
I get the same strange behavior.
– chris
Mar 13 '13 at 21:03
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.
– Sjoerd C. de Vries
Mar 13 '13 at 21:54
I filed this as a bug. Thanks.
– user21
Mar 14 '13 at 10:19
1
@ruebenko: Have you seen ChrisK's finding about the "DiscontinuityProcessing" method option? It probably might help to track down the bug somewhat faster...
– Albert Retey
Mar 14 '13 at 14:04
@AlbertRetey, yes this seems to be an issue with the new in V9 discontinuity processing. Thanks for letting me know.
– user21
Mar 14 '13 at 19:37
=================
2 Answers
2
=================
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]]];

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…

3

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.

– Chris K

Mar 14 ’13 at 13:48

1

@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…

– Albert Retey

Mar 14 ’13 at 14:03