3
$\begingroup$

Perhaps it's something obvious, but why does this simple piecewise function (as well as other conditionals) not evaluate in NDSolve:

fcnP[t_] := Piecewise[{{1, 0.4 < FractionalPart[t] < 0.6},
   {0, FractionalPart[t] >= 0.6}, {0, FractionalPart[t] <= 0.4}}]

temp = NDSolve[{y'[t] == fcnP[t], y[0] == 0}, y, {t, 0, 5}, 
   MaxStepSize -> 0.1][[1]]

If I fit an interpolating function to fcnP, I get the correct solution.

$\endgroup$
  • $\begingroup$ Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! $\endgroup$ – Michael E2 Jun 21 at 2:12
  • $\begingroup$ You can format inline code and code blocks by selecting the code and clicking the {} button above the edit window. The edit window help button ? is useful for learning how to format your questions and answers. You may also find this meta Q&A helpful $\endgroup$ – Michael E2 Jun 21 at 2:12
  • 1
    $\begingroup$ I think this is a bug, and you should report it to Wolfram Research. It turns out that NDSolve does process FractionalPart, but it seems to do it incorrectly. $\endgroup$ – Michael E2 Jun 21 at 4:40
4
$\begingroup$

Perhaps NDSolve dos not know how to expand FractionalPart. PiecewiseExpand benefits from the addition of an assumption giving the interval of integration:

temp = NDSolve[{y'[t] == PiecewiseExpand[fcnP[t], 0 <= t <= 5], 
   y[0] == 0}, y, {t, 0, 5}, MaxStepSize -> 0.1]

ListLinePlot[y /. temp]

enter image description here

Here's another reformulation that works (no need for MaxStepSize either):

fcnP[t_] := Floor[t + 0.6] - Floor[t + 0.4]

temp = NDSolve[{y'[t] == fcnP[t], y[0] == 0}, y, {t, 0, 5}]

Addendum: Evidence of a bug

ClearAll[fcnP];
fcnP[t_] := Piecewise[{{1, 0.4 < FractionalPart[t] < 0.6}}];
{state} = 
  NDSolve`ProcessEquations[{y'[t] == fcnP[t], y[0] == 0}, 
   y, {t, 0, 5}, MaxStepSize -> 0.1, EvaluationMonitor :> Sow[foo]];

nf = state@"NumericalFunction"

Mathematica graphics

In the numerical function, we see the FractionalPart in the right-hand side has been replaced by an expression $(1+s)/2$, where $s$ is NDSolve`s$nnn and nnn is a module number (NDSolve`s$1299443 in my run). The value of $(1+s)/2$ as a function of t is shown in the blue graph; the gold graph shows the value of the right-hand side of the ode (which should be equivalent to fcnP[t]). The only time the RHS is 1 is at t == 0; however, it should be zero then. It is also clear that $(1+s)/2$ is not periodic like FractionalPart. It seems to me that there is a bug with how this ODE is handled.

Table[
   NDSolve`Iterate[state, t];
   With[{
     s = (1 + 
         First@NDSolve`SolutionDataComponent[
           state@"SolutionData"["Forward"], "ID"])/2,
     rhs = First@NDSolve`EvaluateWithSolutionData[nf, 
        state@"SolutionData"["Forward"]]
     },
    {{t, s}, {t, rhs}}
    ],
   {t, 0, 5, 0.01}] // Transpose // ListLinePlot

enter image description here

$\endgroup$
  • $\begingroup$ Another possibility is to define fcnP[t_?NumericQ] := ... $\endgroup$ – Carl Woll Jun 21 at 3:24
  • $\begingroup$ @CarlWoll Yes, but the discontinuity is not handled efficiently. This is better: fp[t_?NumericQ] := FractionalPart[t]; fcnP[t_] := Piecewise[{{1, 0.4 < fp[t] < 0.6}}], as good as PiecewiseExpand, but the Floor version takes the fewest steps. $\endgroup$ – Michael E2 Jun 21 at 4:38
  • $\begingroup$ Thanks @CarlWoll, this is what I also just found worked. The above problem is just a toy example of a more complex function that is difficult to use piecewise on. I've used ?NumericQ approach for Minimize and other optimization functions. It appears that NDSolve has change from many years ago, and now places greater importance on symbolic interrogation of the ODE. It seems that discontinuous functions like Max and FractionalPart require more explicit handling if the ODE is left in symbolic form. Thanks also MichaelE2. $\endgroup$ – Joe Vallino Jun 21 at 14:41

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.