problem when using WhenEvent "StopIntegration" in NDSolve

Question

I am trying to solve some system of nonlinear differential equations by using NDSolve. I am looking for a particular shape of the solution which definitely excludes functions that at the initial point have a positive derivative. I want to use WhenEvent to stop the integration if the first derivative of the function is positive at the initial point (or actually a bit after it). However, when I try to do that, WhenEvent is never triggered. Here is a simplified problem where the same error occurs:

NDSolve[{y'[t] == y[t], y[0] == 1, WhenEvent[y'[10^{-6}] > 0, {"Stop Integration"}]}, y[t], {t, 0, 10}]

I would expect that the integration would stop at t=10^{-6} because the derivative is definitely positive there, but actually, NDSolve continues to integrate up until t=10. What is the problem? Why can't WhenEvent register that the derivative is larger than zero and stop the integration at t=10^{-6}?

Answer 1

The limitation that an event in WhenEvent cannot happen at the initial condition has been criticized on this site before. Sometimes a user wants the event to happen on the initial step and sometimes not. Whether or not this conflict was considered, the design choice sticks to the idea that an event is a change in state between steps that crosses a user-defined threshold. So an event can happen between step 0 and step 1, but not at step 0.

The OP would like an event that ends integration as soon as the derivative y'[t] is positive, including at the beginning. It may appear the way to do this is to check the system at the initial condition before calling NDSolve. One might think an event like t > 10^(-6) && y'[t] > 0 might be a way to accomplish it. However, an event of the form And[event, condition] has a special semantics, which has also been criticized. It seems to be the case that the highest-order derivative y'[t] may appear in event, but errors ensue if it appears in condition. However, there is a hack to get around this, but I think the most satisfying approach is to use the NDSolve components. Whatever the limitations of WhenEvent, one can get under the hood of NDSolve and the principal components are documented. With a little programming, one can usually get the NDSolve suite to do what one wants.

There's one more restriction to deal with. The form y'[t] > 0 has the special semantics of triggering an event if y'[t] changes from negative to positive. It does not trigger an event if y'[t] changes from zero to positive, at least from step 0 to step 1. I'm not sure what happens in the middle of integration, since it's virtually impossible to get a continuously changing quantity to land exactly at the value 0.. To get the desired event we use the logically equivalent Not[y'[t] >= 0], which is treated strictly as a logical, boolean condition with no special interpretation. ("Aaargh," I hear you say.)

Here is a modification of the OP's example that has four solutions, two which start with a positive y'[t]:

states = NDSolve`ProcessEquations[{
    y'[t]^2 == y[t]^2, y[0]^2 == 1,
    WhenEvent[Not[y'[t] <= 0], "StopIntegration"]
    }, y, {t, 0, 10}];
Table[
 (* advance integration if y' is negative or zero *)
 If[NonPositive@First@NDSolve`SolutionDataComponent[
     state@"SolutionData"["Active"], "X'"],
  NDSolve`Iterate[state, First@state@"VariableRanges"]
  ];
 NDSolve`ProcessSolutions[state],
 {state, states}]

Here's the hacky way to get WhenEvent to stop on the first step (replace $MinMachineNumber by 10^-155 to avoid underflow warning):

NDSolve[{
  y'[t] == y[t], y[0] == 1,
  s[0] == 1,   (* user determines initial sign *)
  WhenEvent[y'[t] > 0, "StopIntegration"],
  WhenEvent[t > $MinMachineNumber && s[t] > 0,
   "StopIntegration",
   (* not necessary but reflects the intention: *)
   "LocationMethod" -> "StepBegin"]
  },
 y[t], {t, 0, 10},
 DiscreteVariables -> {s \[Element] {-1, 1}}]

Extension of the first method. For higher dimensional systems higher than a single first-order ODE, the derivative data, the "X'" component, will be a vector longer than 1. One needs a way to find the particular derivative of interest, say y'. The value of state@"Variables" indicates what the components of state@"SolutionData"["Active"] represent.

states = NDSolve`ProcessEquations[{
    y''[t]^2 == y[t]^2, y'[0]^2 == 1, y[0] == 1,
    WhenEvent[Not[y'[t] <= 0], "StopIntegration"]
    }, y, {t, 0, 10}];
Table[
 (* advance integration if y' is negative or zero *)
 If[NonPositive@First@Extract[
     NDSolve`SolutionDataComponent[
      state@"SolutionData"["Active"], "X'"],
     Position[
      NDSolve`SolutionDataComponent[
       state@"Variables", "X'"],
      y'
      ]],
  NDSolve`Iterate[state, First@state@"VariableRanges"]
  ];
 NDSolve`ProcessSolutions[state],
 {state, states}]

Here is an example where the distinction between y'[t] > 0 and Not[y'[t] >= 0] makes a difference:

states = NDSolve`ProcessEquations[{
    y''[t]^2 == y[t]^2, y'[0] == 0, y[0]^2 == 1,
    WhenEvent[Not[y'[t] <= 0], "StopIntegration"]
    }, y, {t, 0, 10}];
Table[
 If[NonPositive@First@Extract[
     NDSolve`SolutionDataComponent[
      state@"SolutionData"["Active"], "X'"],
     Position[
      NDSolve`SolutionDataComponent[
       state@"Variables", "X'"],
      y'
      ]],
  NDSolve`Iterate[state, First@state@"VariableRanges"]
  ];
 NDSolve`ProcessSolutions[state],
 {state, states}]
(* underflow warnings omitted *)

If the event is changed to y'[t] > 0, the first and last solutions run through a long interval in which y'[t] > 0.