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.

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?

share|improve this question
    
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
show 1 more comment

2 Answers

up vote 3 down vote accepted

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]]];
share|improve this answer
add comment

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

share|improve this answer
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
add comment

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.