10
$\begingroup$

Bug introduced in 10.3 or earlier and persisting in 11.0
10.0 is not affected


I'm fairly sure the following is a bug, and I would normally just report it to WRI. However, this is related to my answer to When using NDsolve, how to determine the positions of steady states? So I thought I would post it as a question, for two reasons. One is to let the community verify that it is in fact a bug. The other is to give a reference to the issue for the other question.

Here is a simplified example. With the WhenEvent rule x[t] -> 1., NDSolve fails to integrate the ODE and gives no message.

{sol} = NDSolve[
  {x'[t] == -0.08 x[t], x[0] == 1.,
   WhenEvent[Norm[{x'[t]}] < 0.0001, {x[t] -> 1.}]},
  {x}, {t, 0, 200}]

Mathematica graphics

One the other hand, using an empty action, the system is integrated:

{sol} = NDSolve[
  {x'[t] == -0.08 x[t], x[0] == 1.,
   WhenEvent[Norm[{x'[t]}] < 0.0001, {}]},
  {x}, {t, 0, 200}]

Mathematica graphics

Problem also seems connected to the use of the derivative x'[t] in the event condition. It works fine with x[t].

Mathematica graphics

Am I doing something wrong, or is it a bug?

[Mathematica 10.4.1, OSX 10.11.4.]

$\endgroup$
  • 1
    $\begingroup$ Same behaviour on 10.4.0 for Microsoft Windows (64-bit) . Any deviation from your first example seems to work, e.g. WhenEvent[Norm[{x'[t]}] < 0.0001, {Print@t; x[t] -> 1.}], which indicates that the default step size and resolution is adequate to locate events. Seems like a bug to me. $\endgroup$ – István Zachar Jun 13 '16 at 17:57
  • $\begingroup$ Works fine in v9.0.1: i.stack.imgur.com/IpevD.png I think it's safe to say it's a bug. $\endgroup$ – xzczd Jun 16 '16 at 3:13
  • $\begingroup$ @IstvánZachar Change Print@t to any other expression also works (e.g. an empty list {}). Pretty strange. $\endgroup$ – luyuwuli Sep 15 '16 at 6:53
7
$\begingroup$

Here is a workaround suggested by the response I received from WRI:

{sol} = NDSolve[{x'[t] == -0.08 x[t], x[0] == 1., 
   WhenEvent[Norm[{x'[t]}] < 0.0001, {x[t] -> 1.}]},
  {x}, {t, 0, 200},
  Method -> {"EquationSimplification" -> "Residual"}
  ]

Plot[x[t] /. sol, {t, 0, 200}]

Mathematica graphics

Warning: This option works by converting the system to a DAE, for which only machine-real code is available. That is probably sufficient for most cases, but it does limit one's options. For example, you cannot change WorkingPrecision, and it will not solve a BVP.

The response from WRI suggested using NDSolve[.., SolveDelayed -> True], which indeed works. Its use with NDSolve is pretty much undocumented, and it shows up in red in the front end. One can find references to it in the documentation for the messages NDSolve::ntdv and NDSolve::ndnum. While the text for NDSolve::ntdv in the documentation advises trying SolveDelayed -> True, the text in my current system (V10.4.1) reads like this:

NDSolve::ntdv
(*
"Cannot solve to find an explicit formula for the derivatives. \
Consider using the option \
Method->{\"EquationSimplification\"->\"Residual\"}."
*)

It's hard to know which option should be considered the preferred solution. After a little testing, they seem to produce identical solutions, so I feel they might be invoking the same internal code. I opted for "EquationSimplification" because it is documented for DAEs.

Also, the response from WRI did not indicate that they thought this was a bug (or "issue"); it did not indicate otherwise either.

$\endgroup$
  • 3
    $\begingroup$ I only saw this now: "EquationSimplification" -> "Residual" is the same as the old deprecated SolveDelayed -> True $\endgroup$ – user21 Jul 20 '16 at 17:22
  • $\begingroup$ @user21 Thanks. I thought that was probably the case. $\endgroup$ – Michael E2 Jul 20 '16 at 18:44
  • $\begingroup$ Its seems that the first action in WhenEvent will be ignored. Another workaround may be add a harmless action to feed it; (e.g. WhenEvent[Norm[{x'[t]}] < 0.0001, {}; x[t] -> 1.]), which is also pointed out by István Zachar in the comment under the post. $\endgroup$ – luyuwuli Sep 15 '16 at 7:00
  • $\begingroup$ Wow, you cracked the 100K. Well done! $\endgroup$ – user21 Oct 23 '16 at 14:53
  • $\begingroup$ @user21 Thanks! $\endgroup$ – Michael E2 Oct 23 '16 at 18:01
3
$\begingroup$

Try:

{sol} = NDSolve[{x'[t] == -0.08 x[t], x[0] == 1., 
   WhenEvent[Norm[{x'[t]}] < 0.0001, x[t] -> 1.; x[t] -> 1]}, {x}, {t,
    0, 200}]

enter image description here

$\endgroup$
  • $\begingroup$ But x'[t] does seem to be available in this similar code: NDSolve[{x'[t] == Sin[t], x[0] == 1, WhenEvent[x'[t] < 0, Null; x[t] -> 3 - x[t]]}, x, {t, 0, 20}]. Similarly it is used here, albeit in the pre-WhenEvent days, to locate all local extrema. -- During integration, the value of x'[t] is one of the "SolutionDataComponents" and can be used in the event (but not the action, I think) of a WhenEvent call. See Components and Data Structures for more. $\endgroup$ – Michael E2 Jun 13 '16 at 16:30
  • $\begingroup$ @MichaelE2 And apparently similarly {sol} = NDSolve[{x'[t] == -0.08 x[t], x[0] == 1., WhenEvent[Norm[{x'[t]}] < 0.0001, x[t] -> 1.; x[t] -> 1]}, {x}, {t, 0, 200}] Works fine. Now also allowing plotting of velocity. $\endgroup$ – Feyre Jun 13 '16 at 16:39

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.