# How to increase the speed of NDsolve when using WhenEvent

I am trying to solve a set of ODE's with switching. I implemented using WhenEvents.

xd = {1/√6, √5/√6};
x0 = {1/√2, 1/√2};
min = -0.0001;
max = 0.0001;
tmax = 10;
sol = NDSolve[{
x1'[t] == -x2[t]*u1[t],
x2'[t] == x1[t]*u1[t],
WhenEvent[xd[]* x2[t] - xd[] x1[t] > max, u1[t] -> -1],
WhenEvent[xd[]* x2[t] - xd[] x1[t] < min, u1[t] -> 1],
x1 == x0[], x2 == x0[], u1 == -1},
{x1, x2, x3, u1}, {t, 0, tmax}, DiscreteVariables -> {u1}]
Plot[{xd[], xd[], x1[t] /. sol, x2[t] /. sol}, {t, 0, tmax},
PlotPoints -> 10000]
Plot[u1[t] /. sol, {t, 0, tmax}, PlotPoints -> 10000] How ever, I need to another event WhenEvent[ xd[]* x2[t] - xd[] x1[t] < max && xd[]* x2[t] - xd[] x1[t] > min, u1[t] -> 0], it takes hours and does not produce results.

I think, the last event is raising too many events, so it's taking a lot of time. Is there a work around for this?

• I think it's just the new event; try it as the only event. Possibly related: mathematica.stackexchange.com/a/39393/4999 Sep 20, 2020 at 21:41
• Perhaps WhenEvent[Null; <ineq.>, u1[t] -> 0] from my answer in the link does what you want? Sep 20, 2020 at 21:55
• @MichaelE2 No, it did not work. For f(x)>eps, f(x)<-eps it works. I want to include the event -eps<f(x)<eps.
– kosa
Sep 20, 2020 at 23:30
• @MichaelE2 I checked it. I think it's not related to my issue.
– kosa
Sep 20, 2020 at 23:31
• Thanks for the accept! This might be worth reporting to Wolfram. I might not be considered a bug, but I think they would want the behavior to be a bit better here. I suspect the algorithm could be improved. Sep 21, 2020 at 5:27

It seems to be related to the event location method. I used foo to indicate when an event was detected, and murf to show when a step is taken. With the default "LocationMethod", the integration gets stuck on the first event. This happens even when the other two events are removed. It also happens only if the event action is u[t] -> 0; change it to another nonzero value, and everything works fine. (I guess this makes sense. NDSolve is trying to find where the event occurs, but each time it steps across the event, the system stops moving, since u[t] == 0 sets all derivatives equal to zero. I guess that's confusing the root-finding algorithm.)

foo = murf = 0.;
PrintTemporary@
Dynamic@{foo, Style[murf, PrintPrecision -> 17], Clock@Infinity};
xd = {1/√6, √5/√6};
x0 = {1/√2, 1/√2};
min = -0.0001;
max = 0.0001;
tmax = 10;
sol = NDSolve[{x1'[t] == -x2[t]*u1[t], x2'[t] == x1[t]*u1[t],
WhenEvent[
xd[]*x2[t] - xd[] x1[t] > min &&
xd[]*x2[t] - xd[] x1[t] < max,
foo = t; u1[t] -> 0, "LocationMethod" -> "LinearInterpolation"],
WhenEvent[xd[]*x2[t] - xd[] x1[t] > max,(*foo=t;*)
u1[t] -> -1],
WhenEvent[xd[]*x2[t] - xd[] x1[t] < min,(*foo=t;*)u1[t] -> 1],
x1 == x0[], x2 == x0[], u1 == -1},
{x1, x2, x3, u1}, {t, 0, tmax},
DiscreteVariables -> {u1 ∈ {-1, 0, 1}},
StepMonitor :> (murf = t), MaxSteps -> 1000]


Of course it's stops pretty soon:

Plot[{xd[], xd[], x1[t] /. sol, x2[t] /. sol}, {t, 0, tmax}]
Plot[u1[t] /. sol, {t, 0, tmax}, PlotPoints -> 100] 