Using WhenEvent to Bound Variables

Question

I would like to use the WhenEvent function (link here) to set a bound on a dynamic variable in a set of ODEs. That is, when the variables k21[t] and k12[t] reach above specific value a, I would like to set them back at that value. A simplified version of the problem looks like

w1 = 6/24.5;
w2 = 6/23.5;
a = 0.1;

Eqs = {
   x1'[t] == w1 + (k21[t]/2)*Sin[x2[t] - x1[t]],
   x2'[t] == w2 + (k12[t]/2)*Sin[x1[t] - x2[t]],
   k21'[t] == a*(Cos[x2[t] - x1[t] + \[Pi]] + 1),
   k12'[t] == a*(Cos[x1[t] - x2[t] + \[Pi]] + 1)};

ICs = {x1[0] == 3/2, x2[0] == 3/4, k21[0] == 0.0001, k12[0] == 0.0001};

events = {WhenEvent[Abs[k21[t]] > a, k21[t] -> a], WhenEvent[Abs[k12[t]] > a, k12[t] -> a]};

EqsICs = Join[Eqs, ICs, events];

SolutionValue[t_] = NDSolveValue[EqsICs, {x1[t], x2[t], k21[t], k12[t]}, {t, 0, 10^6}];

And I plot the solutions as:

Show[
Plot[SolutionValue[t][[3]], {t, 0, tmax}, PlotRange -> {{0, tmax}, {-0.1, 0.1}}, AxesOrigin -> {0, 0}], 
Plot[SolutionValue[t][[4]], {t, 0, tmax}, PlotRange -> {{0, tmax}, {-0.1, 0.1}}, AxesOrigin -> {0, 0}]
]

However, my WhenEventdoesn't work. The plot shows that the values of kij keep growing beyond a. Is my syntax wrong? Thanks :)

Answer 1

The syntax of WhenEvent can be a little misleading. A WhenEvent is triggered by an event that happens at a particular time value, not a condition. In your case, Abs[k21[t]] > a is only triggered at the moment Abs[k21[t]] passes through a from below, but it is not detected afterwards. Therefore this doesn't keep k21 in bounds.

The approach that works for me (example here) is to multiply each equation by a DiscreteVariable that indicates whether that variable is in-bounds (1) or not (0). One event is when the variable reaches the boundary, set the indicator to 0. The other event you'll need is when the variable changes to being pushed back in from the boundary, set the indicator back to 1. This prevents getting stuck on the boundary forever.

In your case there are four possible events (k12 and k21 could each become > a or < -a). The following seems to work (note I switched to NDSolve, replaced 1 by a constant c to test other cases, and reduced tmax to 100):

w1 = 6/24.5;
w2 = 6/23.5;
a = 0.1;

c = 1;
tmax = 100;

Eqs = {
   x1'[t] == w1 + (k21[t]/2)*Sin[x2[t] - x1[t]],
   x2'[t] == w2 + (k12[t]/2)*Sin[x1[t] - x2[t]],
   k21'[t] == in[k21][t] (a*(Cos[x2[t] - x1[t] + π] + c)),
   k12'[t] == in[k12][t] (a*(Cos[x1[t] - x2[t] + π] + c))
   };

ICs = {x1[0] == 3/2, x2[0] == 3/4, k21[0] == 0.0001, k12[0] == 0.0001,
    in[k12][0] == 1, in[k21][0] == 1};

events = {
  WhenEvent[k21[t] > a, in[k21][t] -> 0],
  WhenEvent[a*(Cos[x2[t] - x1[t] + π] + c) < 0 && in[k21][t] == 0, in[k21][t] -> 1],
  WhenEvent[k21[t] < -a, in[k21][t] -> 0],
  WhenEvent[a*(Cos[x2[t] - x1[t] + π] + c) > 0 && in[k21][t] == 0, in[k21][t] -> 1],
  WhenEvent[k12[t] > a, in[k12][t] -> 0],
  WhenEvent[a*(Cos[x1[t] - x2[t] + π] + c) < 0 && in[k12][t] == 0, in[k12][t] -> 1],
  WhenEvent[k12[t] < -a, in[k12][t] -> 0],
  WhenEvent[a*(Cos[x1[t] - x2[t] + π] + c) > 0 && in[k12][t] == 0, in[k12][t] -> 1]
};

EqsICs = Join[Eqs, ICs, events];

sol = NDSolve[
   EqsICs, {x1, x2, k21, k12, in[k12], in[k21]}, {t, 0, tmax}, 
   DiscreteVariables -> {in[k12], in[k21]}];

Plot[Evaluate[{k21[t], k12[t]} /. sol], {t, 0, tmax}, PlotRange -> All]

Evidently you do get stuck on the top boundary. If we reduce c to 0.95, we can see that you can leave the boundary: