Detecting nearly simultaneous WhenEvents in NDSolve

I am trying to solve a system of (many) coupled nonlinear ODEs. I need to decouple some of the equations (i.e. set the time derivatives of some of the dependent variables to zero) at various points in the time integration when particular criteria are reached. To do this, I am using NDSolve with equations that contain some terms multiplied by DiscreteVariables whose values (1 or 0) are controlled by WhenEvents. The problem is, some of the WhenEvents tend to occur at nearly simultaneous times, and NDSolve appears to have difficult detecting multiple simultaneous events without the use of a ridiculously small MaxStepSize.

I know that the DiscreteVariables functionality for NDSolve is new, as of v9, and I have found similar question at mathematica.stackexchange.com/questions/76777, although in a system without DiscreteVariables.

Here is a minimal working example (per Michael E2's first answer, note that we can use Sow and Reap to check that only one event is being triggered) that fails to evaluate both WhenEvents:

With[
{
ti = Log@100,
tf = Log@(10^9),
a0 = 3.05917*^7,
b0 = 3.05242*^7
},
{{sol}, {evts}} = Reap@NDSolve[
{
a'[t] == (a[t] - b[t] s2[t]) s1[t],
b'[t] == (-a[t] s1[t] + b[t]) s2[t],
a[tf] == a0^2,
b[tf] == b0^2,
s1[tf] == 1,
s2[tf] == 1,
WhenEvent[Evaluate[a[t] > Exp[2 t]], {Sow["a" -> t], s1[t] -> 0}, "DetectionMethod" -> "Sign"],
WhenEvent[Evaluate[b[t] > Exp[2 t]], {Sow["b" -> t], s2[t] -> 0}, "DetectionMethod" -> "Sign"]
},
{a, b, s1, s2},
{t, ti, tf},
DiscreteVariables -> {Element[s1, {0, 1}], Element[s2, {0, 1}]},
Method -> Automatic,
WorkingPrecision -> Automatic
MaxStepSize -> Automatic,
MaxSteps -> Infinity
]
]

(*
{{{a -> InterpolatingFunction[{{4.60517, 20.7233}}, <>],
b -> InterpolatingFunction[{{4.60517, 20.7233}}, <>],
s1 -> InterpolatingFunction[{{4.60517, 20.7233}}, <>],
s2 -> InterpolatingFunction[{{4.60517, 20.7233}}, <>]}},
{{"b" -> 17.2351}}}
*)


The scales (10^7, 10^9) are taken from my work, but the equations are much simplified. Both s1 and s2 should step from 1 to 0 as the system integrates down from tf to ti; however, usually only the "b" event triggering s2, the first one to occur, is recognized.

According to the documentation, the "DetectionMethod" -> "Interpolation" option for WhenEvent should be more "robust" against multiple events per step than the default option, but I actually generally get better behavior with the default "DetectionMethod" -> "Sign".

Using more precision and a small step size does get the job done:

a0 = SetPrecision[3.05917*^7, 2 MachinePrecision],
b0 = SetPrecision[3.05242*^7, 2 MachinePrecision],
Method -> "ExplicitRungeKutta",
MaxStepSize -> 0.0001


However, this small of a step size (seems to require the step to be smaller than the time separation between events by about an order of magnitude) takes too much time for my big system of coupled equations and is not always successful.

I can also get both events to be evaluated if I set "LocationMethod" -> "StepEnd" in each WhenEvent, but that tends to produce moderate discontinuities in the solutions when I run it with my larger system (probably an interesting problem in itself). The fact that it works makes me suspect that the trouble has something to do with the order in which the WhenEvents are being evaluated (e.g. the event for s2 is processed first, and the result of that evaluation causes the system not to recognize that the event for s1 also occurs during that timestep). Giving each WhenEvent a different "Priority" doesn't help.

My question is: is there a robust way to ensure that all the WhenEvents are recognized (such that the decoupling variables s1 and s2 both go to zero) without resorting to a small timestep or pushing the effects of the decoupling to the end of the step via "LocationMethod" -> "StepEnd"?

Update: This appears to be a bug. I sent this in to Wolfram, and they acknowledged that "NDSolve is not behaving properly in this case", and forwarded it to their developers.

• Welcome to Mathematica.SE! I suggest that: 1) You take the introductory Tour now! 2) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! 3) As you receive help, try to give it too, by answering questions in your area of expertise. – bbgodfrey Apr 12 '15 at 21:54

Here's a workaround. I'm not sure why the variables s1[t], s2[t] are not reset in my first answer (see edit history). We can take care of things manually by making s1 and s2 numerical functions.

Block[{ti = Log@100, tf = Log@(10^9), a0 = 3.05917*^7,
b0 = 3.05242*^7, s1, s2, s10 = 1, s20 = 1},
s1[t_?NumericQ] := s10;
s2[t_?NumericQ] := s20;
{{sol}, {evts}} =
Reap@NDSolve[{a'[t] == (a[t] - b[t] s2[t]) s1[t],
b'[t] == (-a[t] s1[t] + b[t]) s2[t], a[tf] == a0^2, b[tf] == b0^2,
WhenEvent[a[t] > Exp[2 t], Sow["a" -> t]; s10 = 0],
WhenEvent[b[t] > Exp[2 t], Sow["b" -> t]; s20 = 0]},
{a, b}, {t, ti, tf}]

(*
{{{a -> InterpolatingFunction[{{4.60517, 20.7233}}, <>],
b -> InterpolatingFunction[{{4.60517, 20.7233}}, <>]}},
{{"b" -> 17.2351, "a" -> 17.2351}}}
*)

• That is interesting. In your example it looks like both WhenEvents are being evaluated, which is technically what I was asking for. However, if you look at the interpolation functions for s1 and s2 in the solution both are constant at 1, rather than going to 0 as they should (i.e only the first Sow action occurs?), which is what I need. – Virgil Apr 12 '15 at 23:54
• @Virgil I think you're right. I didn't notice since graphs of a and b flatten out at about the right time. NDSolve does not like me today. This is the second time things do not work as advertised. – Michael E2 Apr 13 '15 at 0:05
• Good use of Reap and Sow, though. I'm not sure why NDSolve sees both events if there is a Sow statement first in the WhenEvent action. The answer to the question I referenced with this question seems to indicate that WhenEvent has different behavior for simple vs. compound action statements. I might be running in to that. – Virgil Apr 13 '15 at 0:15
• I have, of course, linked you to your own answer. – Virgil Apr 13 '15 at 0:40
• @Virgil, yes, I realized that. So far everything else I've tried (other than the current aswer) either misses an event or fails to reset the variables. – Michael E2 Apr 13 '15 at 1:49