Message Boards Message Boards

0
|
4862 Views
|
5 Replies
|
1 Total Likes
View groups...
Share
Share this post:

NDSolve misbehaviour

Marc Widdowson
Posted 3 years ago

I'm using the latest version, Mathematica 12.1.1.0 on Windows.

I've been doing some work with NDSolve involving coupled differential equations, and it hasn't behaved as I would have expected. I've got x'=f(x,y) and y'=g(x,y), and I'm finding that the variation in y is not feeding back into x' as I think it ought to, i.e. x' seems to be insensitive to the changes in y even though it's supposed to depend on y. In the course of investigating this, I have got it down to a very simple situation where I have got one variable, x, and I make x' depend on t. This also does not behave as I would expect. Specifically, I have

eqns = {
   x'[t] == If[OddQ[Floor[t/10]], 1, -1],
   x[0] == 0
   };
soln = NDSolve[eqns, {x[t]}, {t, 0, 100}];

I would expect this to generate a sawtooth. x should increase with slope 1 when time is an odd multiple of 10 and decrease with slope -1 when time is an even multiple of 10.

But what I get is that x just decreases with constant slope -1. The first plot below is of 

Plot[Evaluate[x[t] /. First[soln]], {t, 0, 100}]

which does not display the expected sawtooth, and the second plot below is of

OddQ[Floor[t/10]]

which does switch between 1 and -1 as expected.

enter image description here

If I use a one-off switch in the value of x', so that it changes from 1 to -1 as t goes from below 10 to above 10, i.e.

eqns = {
   x'[t] == If[Floor[t/10] < 1, 1, -1],
   x[0] == 0
   };
soln = NDSolve[eqns, {x[t]}, {t, 0, 100}];
Plot[Evaluate[x[t] /. First[soln]], {t, 0, 100}]

Then I get the behaviour I'd expect...

enter image description here

Similarly, if I make the driving of x' sinusoidal, I also get the behaviour I'd expect...

eqns = {
   x'[t] == Sin[t],
   x[0] == 0
   };
soln = NDSolve[eqns, {x[t]}, {t, 0, 100}];
Plot[Evaluate[x[t] /. First[soln]], {t, 0, 100}]

results in...

enter image description here

So why doesn't the initial square wave driving and expected sawtooth for x work? Is this a bug, or an expected limitation of NDSolve? Or is there something I am missing?

This is disturbing because I don't know whether the unexpected behaviour of what I'm actually trying to model is a problem with my model or a problem with NDSolve.

POSTED BY: Marc Widdowson
5 Replies
Sort By:

I don't know what happens here, but this is a workaround

eqns = {x'[t] == (-1)^(IntegerPart[t/10]), x[0] == 0};
soln = NDSolve[eqns, x, {t, 0, 100}];
Plot[(x /. First[soln])[t], {t, 0, 100}]
POSTED BY: Hans Dolhaine
Marc Widdowson
Marc Widdowson
Posted 3 years ago

Here is some further weirdness.

I define a function to return +1 or -1 according to whether the input is an odd or even multiple of 10. This function tests the parity in two ways, one by seeing whether the number is 1 modulo 2, and one using the OddQ function.

func[t_] := Module[{},
  {If[Mod[Floor[t/10], 2] == 1, 1, -1], If[OddQ[Floor[t/10]], 1, -1]}
  ]

Obviously, both return values should be the same, and indeed this is the case. E.g.

{func[11], func[23], func[35]}

gives

{{1, 1}, {-1, -1}, {1, 1}}

as expected.

I now use this in my differential equation. Firstly, I take the first part of the return value of func, i.e.

soln = NDSolve[{x'[t] == func[t][[1]], x[0] == 0}, x, {t, 0, 100}];
Plot[x[t] /. soln, {t, 0, 100}]

This gives a sawtooth as expected:

enter image description here

Next, I take the second part of the return value of func, i.e.

soln = NDSolve[{x'[t] == func[t][[2]], x[0] == 0}, x, {t, 0, 100}];
Plot[x[t] /. soln, {t, 0, 100}]

This gives a monotonic decline as previously:

enter image description here

This makes no sense. The first and second parts of the return value from func are exactly the same, as shown above. How come Mathematica knows that the second part was calculated from OddQ and treats it differently, somehow ignoring the actual value and substituting a constant value?
POSTED BY: Marc Widdowson
Marc Widdowson
Marc Widdowson
Posted 3 years ago

Thank you very much. I also found that testing the remainder mod 2 works, see my other post, and it's weird that Mathematica can distinguish between two things that ought to be a pure number--both came from If statements but one of them used OddQ...and Mathematica somehow doesn't like that one. What's more, if you plot the number directly, it behaves normally, but if you use it in a differential equation, Mathematica ignores the fact it is changing.

Of course, this isn't really what I was interested in or trying to do. I was trying to understand why, when a variable changed value in my set of coupled differential equations, it didn't affect the behaviour of the other variable even though the other variable's differential equation depended on it. So I hard-wired in this switching behaviour and found that that had no effect either. Now I find that it does have an effect provided I use the test Mod[u,2]==1 and not OddQ[u], even though those should be equivalent.
POSTED BY: Marc Widdowson

This could be a matter of when the function is evaluated. If OddQ is evaluated on symbolic as opposed to numeric input it will return False right away. Similar with the If as opposed to a more "mathematical" formulation.

General advice: When possible, use mathematical rather than procedural code inside functions like NDSolve. For example, use Piecewise rather than If. When not possible, force the function to remain "inert" unless given explicit numeric input. This is done by restricting input patterns, with constructs like f[x_?NumberQ]:=.... See NDSolve documentation for further details about this.
POSTED BY: Daniel Lichtblau

Ok, but

g[x_?NumberQ] := Piecewise[
  {
   {1, OddQ[Floor[x/10]]},
   {-1, EvenQ[Floor[x/10]]}
   }
  ]

eqns = {x'[t] == g[t], x[0] == 0};
soln = NDSolve[eqns, x, {t, 0, 100}];
Plot[(x /. First[soln])[t], {t, 0, 100}]

gives just something coinciding with the x-axes. Why that????
POSTED BY: Hans Dolhaine
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Tag limit exceeded
Note: Only the first five people you tag will receive an email notification; the other tagged names will appear as links to their profiles.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract