11
$\begingroup$

Bug introduced after 10.0, in 10.3 or earlier and persisting in 13.2.1

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

Improve this question
$\endgroup$
3
  • 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, 2016 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, 2016 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, 2016 at 6:53

2 Answers 2

Reset to default
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.

Improve this answer
$\endgroup$
5
  • 3
    $\begingroup$ I only saw this now: "EquationSimplification" -> "Residual" is the same as the old deprecated SolveDelayed -> True $\endgroup$
    – user21
    Jul 20, 2016 at 17:22
  • $\begingroup$ @user21 Thanks. I thought that was probably the case. $\endgroup$
    – Michael E2
    Jul 20, 2016 at 18:44
  • 1
    $\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, 2016 at 7:00
  • 1
    $\begingroup$ Wow, you cracked the 100K. Well done! $\endgroup$
    – user21
    Oct 23, 2016 at 14:53
  • $\begingroup$ @user21 Thanks! $\endgroup$
    – Michael E2
    Oct 23, 2016 at 18:01
4
$\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

Improve this answer
$\endgroup$
2
  • $\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, 2016 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, 2016 at 16:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.