5
$\begingroup$

I'm having some trouble stopping NDSolve when the solutions gets outside of a given domain. There are five unknowns, and NDSolve should stop when the solution becomes not canonical, where canonical is defined by the following function:

canonical[list_List] := Last@Sort@ Flatten[NestList[Reverse /@ # &, 
                                 NestList[RotateLeft, list, Length[list] - 1], 1], 1]

Understanding canonical is not required for the question, but let me explain in few words: the problem whose solution are given by NDSolve is invariant by circular permutations and symmetry; so {2,3,4,5,1} does not bring any additional information if I already have {5,4,3,2,1}. That's why I want to stop NDSolve when it gets out of my "canonical domain".

The following function checks if a list is canonical or not:

nonCanonicalQ[vals_] := (vals != canonical[vals])

For example:

nonCanonicalQ[{5,4,3,2,1}] (* False: means that it's a canonical list *)
nonCanonicalQ[{3,4,5,1,2}] (* True: means it's not a canonical list *)

Now, let's use it in WhenEvent using NDSolve, with a very basic example (all solutions are constant, except the last one which is linear):

vars = Array[s, 5];
varsOt = Through[vars[t]];

sysDAE0 = {{s[1]'[t], s[2]'[t], s[3]'[t], s[4]'[t], s[5]'[t]} == {0,  0, 0, 0, 1}, 
           {s[1][0], s[2][0], s[3][0], s[4][0], s[5][0]} == {4, 3, 2, 1, 0}}

{solDAE} = NDSolve[{sysDAE0}~Join~{WhenEvent[Evaluate@nonCanonicalQ[varsOt], 
  "StopIntegration"]}, vars, {t, -1000, 1000}];
solDAE = vars /. solDAE

This outputs a solution valid on domain $[-1000,1000]$. However, it should have triggered the WhenEvent because

nonCanonicalQ[Through[solDAE[0]]] (* False *)
nonCanonicalQ[Through[solDAE[100]]] (* True *)

The integration should have stopped somewhere between $t=0$ and $t=100$.

$\endgroup$
4
  • $\begingroup$ Did you define vars and varsOt? $\endgroup$
    – Cassini
    Commented Dec 9, 2016 at 19:40
  • $\begingroup$ @Cassini Thank you, I forgot to include their definition in the OP. Edited. $\endgroup$
    – anderstood
    Commented Dec 9, 2016 at 19:42
  • $\begingroup$ It seems to have something to do with the evaluation of the predicate because WhenEvent[{s[1][t], s[2][t], s[3][t], s[4][t], s[5][t]} != canonical[{s[1][t], s[2][t], s[3][t], s[4][t], s[5][t]}], "StopIntegration"] stops at t=3. $\endgroup$
    – Cassini
    Commented Dec 9, 2016 at 19:50
  • 1
    $\begingroup$ The third component in each side of nonCanonicalQ[varsOt] is the same, s[3][t], and the relation between them NotEqual should always come up False. This doesn't happen with the solution because nonCanonicalQ is evaluated on the solution values, not on the symbolic expressions as in the NDSolve code. $\endgroup$
    – Michael E2
    Commented Dec 10, 2016 at 1:52

2 Answers 2

6
$\begingroup$

Cassinis solution is one possibility to solve your problem and he is right in assuming that it has to do with evaluation order. The following is a summary of techniques which I found helpful when investigating situations when NDSolve doesn't behave as one would expect.

WhenEvent does have the attribute HoldAll and I guess that is the reason why you did put the Evaluate within it. The Evaluate seems to work but does evaluate with symbolic arguments before NDSolve starts to solve. It will evaluate to:

 {s[1][t],s[2][t],s[3][t],s[4][t],s[5][t]} != {s[5][t],s[4][t],s[3][t],s[2][t],s[1][t]}

which then will be used in every step with numerical values for the s[n][t] in the event check while NDSolve is running and of course will always evaluate to false and thus never trigger the event (strictly speaking it could even evaluate to true for cases where it shouldn't, but that doesn't happen for your example). Here is a simple change to nonCanonicalQ which shows what happens:

nonCanonicalQ[vals_] := Module[{result},
  result = vals != canonical[vals];
  Print["nonCanonicalQ"[vals] -> result];
  result
]

rerunning the call to NDSolve after making that change will show that nonCanonicalQ will be evaluated only once and give the mentioned result.

Cassinis answer shows how to solve the problem by explicitly writing out the variables and avoid the Evaluate. Using the above definition for nonCanonicalQ will show that nonCanonicalQ then is also evaluated only once at the beginning and return:

{s[1][t],s[2][t],s[3][t],s[4][t],s[5][t]} !=
  canonical[{s[1][t],s[2][t],s[3] [t],s[4][t],s[5][t]}]

the canonical in that expression will only be evaluated at every time step with numeric arguments and make WhenEvent trigger as intended (you could check that with a similar change to canonical as I have shown for nonCanonicalQ).

Of course there are several other possibilities to achieve the same thing and without writing out the variables to make that work. One very simple solution would be to change either canonical or nonCanonicalQ so that they only evaluate for numeric arguments, e.g.:

ClearAll@canonical;
canonical[list : {__?NumericQ}] := (
  Last@Sort@Flatten[NestList[
   Reverse /@ # &, NestList[RotateLeft, list, Length[list] - 1], 1], 1]
)

Using that with the OP's rest of the code should work as intended, which can be checked by using my printing version of nonCanonicalQ.

Note that it is important to ClearAll before making the new definition, otherwise the old one will still evaluate early. If you want to go back to the original definition to try the following solution you should also use a ClearAll to really verify that it really does work.

Instead of evaluating early and avoiding evaluation with symbolic arguments you can insert the variables into the WhenEvent with something like this:

With[{vars = varsOt},
  solDAE = NDSolveValue[
    {sysDAE0}~Join~{WhenEvent[nonCanonicalQ[vars], "StopIntegration"]}, 
    vars, {t, -1000, 1000}
  ]
]

which will work even without the Evaluate and without the restriction of canonical to numeric arguments. I would also suggest to adopt NDSolveValue which is avoiding some extra steps to extract the result.

$\endgroup$
5
$\begingroup$

This works:

  {solDAE} = 
    NDSolve[
      {sysDAE0} 
       ~ Join ~ 
      {WhenEvent[
         nonCanonicalQ[
           {s[1][t], s[2][t], s[3][t], s[4][t], s[5][t]}], 
           "StopIntegration"]}, 
      vars, {t, -1000, 1000}]
$\endgroup$
2
  • $\begingroup$ Indeed... Any idea why the other solution does not? $\endgroup$
    – anderstood
    Commented Dec 9, 2016 at 19:55
  • 1
    $\begingroup$ I have a vague idea that it has to do with Hold..., but I'm sure one of the MMA wizards here can tell you. $\endgroup$
    – Cassini
    Commented Dec 9, 2016 at 19:58

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.