4
$\begingroup$

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}?

$\endgroup$
12
  • 1
    $\begingroup$ The proper event syntax according to the docs is WhenEvent[t > 10^(-6), "StopIntegration"]. $\endgroup$ – Michael E2 Mar 25 at 13:38
  • $\begingroup$ @MichaelE2 I agree about the "StopIntegration", my mistake, but it still doesn't work. The condition you put would stop integration whenever t > 10^(-6), but I am interested if the derivative at this point is larger than zero. Again, the example I gave is just to illustrate the actual problem which is more complicated. $\endgroup$ – Jovan Markov Mar 25 at 13:43
  • $\begingroup$ Your event is still ill-defined. WhenEvent[event && condition,..] is triggered when the event "happens" (normally, flips from False to True but there are special cases) provided condition is true; otherwise, the event is ignored. In your example, which is the event and which is the condition, t > 10^(-6) / y'[t] > 0? Or is the event something else? $\endgroup$ – Michael E2 Mar 25 at 14:28
  • $\begingroup$ I thought NDSolve[{y'[t] == y[t], y[0] == 1, WhenEvent[t > 10^-6 && y'[t] > 0, "StopIntegration"]}, y, {t, 0, 10}] would work like NDSolve[{y'[t] == y[t], y[0] == 1, WhenEvent[t > 10^-6 && y[t] > 0, "StopIntegration"]}, y, {t, 0, 10}] does (no prime in event), but it doesn't. Strange... $\endgroup$ – Chris K Mar 25 at 14:33
  • 1
    $\begingroup$ @ChrisK I was going to deal with that after knowing what the event is. It's because the highest order derivative is not treated as a state variable in an ODE. So only t and y[t] are available in WhenEvent. All variables are available in a DAE, so retry your code with the option Method -> {"EquationSimplification" -> "Residual"}. $\endgroup$ – Michael E2 Mar 25 at 15:06
4
$\begingroup$

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.

$\endgroup$

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.