3

I have defined a periodic function $f(t)$ like this (Mathematica 11):

fTmp[t_ /; 0 <= t <= 2] := \[Piecewise] {
   {1, t <= 1},
   {0, t <= 2}
  }
f[t_] := fTmp[Mod[t, 2]]

This works fine, as can be seen from

Plot[f[t], {t, 0, 5}]

Plot of the function...

But when I try to use this function in NDSolve, I get a weird error:

NDSolve[{a'[t] == f[t]*a[t], a[0] == 1}, a, {t, 0, 5}]
... NDSolve::nlnum: The function value {2.71828 fTmp[2.18952]} is not a list of numbers with dimensions {1} at {t,a[t],NDSolve`s$49306[t]} = {2.18952,2.71828,1}.

As you can see, $fTmp(t)$ is called with an argument $t>2$, which of course never should happen.
Any ideas how to fix this problem? Thanks in advance.

Noiralef
  • 379
  • 1
  • 10

2 Answers2

5

Two alternatives:

  • Remove the argument restriction so that ff[t] evaluates to Piecewise[..], which allows the discontinuity processing phase to see Piecewise and set up the integration scheme properly.
  • Use DiscreteVariables to code the function f.

Codes (OP's answer, Piecewise, DiscreteVariables):

Clear[f, fTmp];
fTmp[t_ /; 0 <= t <= 2] := \[Piecewise]{{1, t <= 1}, {0, t <= 2}};
f[t_] := fTmp[Mod[t, 2]];
{sol} = NDSolve[{a'[t] == f[t]*a[t], a[0] == 1}, a, {t, 0, 5}, 
   Method -> {"DiscontinuityProcessing" -> False}];

Clear[ff, ffTmp];
ffTmp[t_ (*/;0 <= t <= 2*)] := \[Piecewise]{{1, t <= 1}, {0, t <= 2}};
ff[t_] := ffTmp[Mod[t, 2]];
{sol2} = NDSolve[{a'[t] == ff[t]*a[t], a[0] == 1}, a, {t, 0, 5}];

{sol3} = NDSolve[{a'[t] == fff[t]*a[t], a[0] == 1,
    fff[0] == 1, WhenEvent[Mod[t, 1] == 0, fff[t] -> 1 - fff[t]]},
   a, {t, 0, 5}, DiscreteVariables -> {fff}];

Comparison (time vs. steps and size of solution): The first alternative is better than the OP's, and the last is better than the first.

ListPlot[Flatten[a["Grid"] /. sol], GridLines -> {None, Automatic}, 
 PlotLabel -> Row[{ByteCount@sol, "bytes"}, " "],
 AxesLabel -> {"steps", "time"}
 ]

Mathematica graphics

ListPlot[Flatten[a["Grid"] /. sol2], GridLines -> {None, Automatic}, 
 PlotLabel -> Row[{ByteCount@sol2, "bytes"}, " "],
 AxesLabel -> {"steps", "time"}
 ]

Mathematica graphics

ListPlot[Flatten[a["Grid"] /. sol3], GridLines -> {None, Automatic}, 
 PlotLabel -> Row[{ByteCount@sol3, "bytes"}, " "], 
 AxesLabel -> {"steps", "time"}
 ]

Mathematica graphics

Another way for the given example is to use DSolve with the first alternative:

{sol4} = DSolve[{a'[t] == ff[t]*a[t], a[0] == 1}, a, {t, 0, 5}]

Mathematica graphics

They produce overlapping plots:

Plot[a[t] /. {sol, sol2, sol3, sol4} // Evaluate, {t, 0, 5}]

Mathematica graphics

Michael E2
  • 235,386
  • 17
  • 334
  • 747
3

I found the answer myself in the answer to this question: Weird NDSolve behavior with Piecewise (MMA9)

NDSolve[{a'[t] == f[t]*a[t], a[0] == 1}, a, {t, 0, 5}, Method -> {"DiscontinuityProcessing" -> False}]

does the trick.

Noiralef
  • 379
  • 1
  • 10