6
$\begingroup$

I was wondering whether there is a way to use WhenEvent in a system of delay differential equations, for example:

z[t_] = Piecewise[{{1, t <= 5}, {3/2, 5 < t <= 10}}, 2];
NDSolve[{x'[t] == 1, y'[t] == x'[t] - x'[t - 1], 
 WhenEvent[x[t - 1] == z[t], x[t] -> 0], x[t /; t <= 0] == 0, y[t /; t <= 0] == 0 }, 
 {x[t], y[t]}, {t, 0, 20}];

The issue seems to be that the WhenEvent doesn't recognize x[t-1] as a function and tries to evaluate the condition by throwing in the most recent x and t values. I'm not sure whether this is because WhenEvent can't access the value at x[t - 1] or because I need to input the condition in a way that explicitly tells it to go back and find this value.

A more complicated example to illustrate one of the systems I am working with:

n = 3; m = 2;
staff = {{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 3, 3, 3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0,0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
P = {{{0, 0}, {1, 0}, {0, 1}}, {{0, 0}, {0, 0}, {7/10, 0}}, {{0, 0}, {0, 4/5}, {0, 0}}};
PT = {{75, 360, 720, 0.3}, {75, 420, 930, 0.4}};
PA = {{1, 7/10}, {0, 0}, {0, 3/10}};
v = {{7/2, 5}, {7/2, 9/2}, {15, 20}};

a[t_, i_, k_] := PT[[k, 1]]*PA[[i, k]]*PDF[
    TriangularDistribution[{PT[[k, 2]], PT[[k, 3]]}, ((1 - PT[[k, 4]])*PT[[k, 2]] + PT[[k, 4]]*PT[[k, 3]])]
  , t];

Ns[t_, i_] := Piecewise[Table[{staff[[i, k]], (k - 1)*30 <= t < k*30}, {k, 1, 48}]];

evts = Flatten[Table[With[{i = i}, {
   WhenEvent[Evaluate[Subscript[V, i][t] == 0] && Evaluate[Subscript[S, i][t] <= Ns[t, i]], Subscript[dv, i][t] -> 0]
   , WhenEvent[{Evaluate[Subscript[V, i][t] > 0], Evaluate[Subscript[S, i][t] > Ns[t, i]]}, Subscript[dv, i][t] -> 1]
   , Subscript[dv, i][t /; t <= 0] == 1
  }], {i, 1, n}]];

eqns = Flatten[Table[With[{i = i}, With[{p = p}, {Subscript[U, i, p]'[ t] == (Subscript[dv, i][t])*(Subscript[Q, i, p][t])*Ns[t, i]/(Subscript[V, i][t] + .001) + (1 - Subscript[dv, i][t])*Subscript[L, i, p]'[t - v[[i, p]]], Subscript[U, i, p][t /; t <= 0] == 0}]], {i, 1, n}, {p, 1, m}]];

rls = Flatten[Table[With[{i = i}
   , {Subscript[S, i][t_] -> Sum[With[{p = p}, Subscript[L, i, p]'[t - v[[i, p]]]*v[[i, p]]], {p, 1, m}]
    , Subscript[V, i][t_] -> Sum[With[{p = p}, Subscript[Q, i, p][t - v[[i, p]]]*v[[i, p]]], {p, 1, m}]
    , Table[With[{p = p}
       , {Subscript[Q, i, p][t_] -> Subscript[L, i, p][t] - Subscript[U, i, p][t]
        , Subscript[L, i, p][t_] -> Integrate[a[t, i, p], t] + Sum[With[{jj = jj}, P[[jj, i, p]]*Subscript[U, jj, p][t]], {jj, 1, n}], Subscript[L, i, p]'[t_] -> a[i, t, p] + Sum[With[{jj = jj}, P[[jj, i, p]]*Subscript[U, jj, p]'[t]], {jj, 1, n}]
      }]
      , {p, 1, m}]
     }]
   , {i, 1, n}]];

NDSolve[
  Flatten[{evts, eqns} //. rls]
  , Flatten[Table[ With[{i = i}, {Subscript[dv, i][t], Table[With[{p = p}, {Subscript[U, i, p][t]}], {p, 1, m}]}], {i, 1, n}]]
  , {t, 0, 1000}
  , DiscreteVariables -> Flatten[Table[With[{i = i}, Subscript[dv, i][t]], {i, 1, n}]]];

This example gives fairly decent results without the time delays, but bad things start happening when any element in v gets too large.

Improve this question
$\endgroup$
3
  • $\begingroup$ Is this too simple an example? It seems to me that since x'[t] == 1, y'[t] will always be zero, and so y is constant, unless it's supposed to know about the discontinuity in x[t] that arises as a consequence of the reset x[t] -> 0, but in that case, you actually have an impulse for y (because of the infinite derivative at those points), and I'm not sure that NDSolve can deal with that. $\endgroup$
    – march
    Nov 16, 2015 at 21:41
  • $\begingroup$ Sorry - I should have been clear that the example was just to illustrate what I'm trying to do with the WhenEvent and that the equations (especially y[t]) are not relevant. I am only interested in whether there is a way to get the WhenEvent to recognize x[t-1] as a function instead of directly substituting the most recent values. $\endgroup$
    – polkabug
    Nov 17, 2015 at 5:49
  • $\begingroup$ My first thought was to define an auxiliary function a[t] and include a[t] == x[t - 1] and a[t/;t<=0]==0 in the list of equations. I tried NDSolve[{x'[t] == t, a[t] == x[t - 1], WhenEvent[a[t] == 1, x[t] -> 0], x[t /; t <= 0] == 0, a[t /; t < 0] == 0}, {x[t], a[t]}, {t, 0, 20}]; as a simple example, but NDSolve had problems once it was time for a to reset (one time unit after x resets). Perhaps we could make this work, but I'm not sure how to fix the resetting problem. $\endgroup$
    – march
    Nov 17, 2015 at 6:07

2 Answers 2

Reset to default
4
$\begingroup$

Here's something I came up with to get you started. Consider this an extended comment until we can figure out how to make this work..

The idea is to use an auxiliary variable a[t] that is a DiscreteVariable, and use it to find the times at which x crosses z. We set the value of a[t] to be that particular time and then use WhenEvent to reset x one time unit later.

I will use a simplified version of your code that doesn't actually involve the delay piece involving y.

timeLimit = 10;
z[t_] = Piecewise[{{1, t <= 2}, {3/2, 2 < t <= 5}}, 2];
NDSolve[{x'[t] == 1/1.5 t
   , WhenEvent[x[t] == z[t], a[t] -> t]
   , WhenEvent[t == a[t] + 1, {x[t] -> 0, a[t] -> 0}]
   , a[t /; t <= 0] == timeLimit
   , x[t /; t <= 0] == 0}
  , {x[t], a[t]}
  , {t, 0, timeLimit}
  , DiscreteVariables -> {a}];
Show[
 Plot[Evaluate[{x[t], z[t], a[t]} /. %], {t, 0, timeLimit}]
 , Plot[{t, t - 1}, {t, 0, timeLimit}, 
  PlotStyle -> {{Black, Dashed, Opacity[0.5]}}]
 ]

enter image description here

The blue curve is the solution for x[t], the orange curve is z[t], the green curve is a[t], and the dashed lines are the functions t and t - 1 that are acting as aids for the eye to show that a[t] gets set to the current value of t and is then reset to 0 after 1 time unit has elapsed.

There are some important observations.

  • We can see that a[t] gets set to t when x[t] crosses z[t]: at these points, x, z, and a all cross.
  • We can see that x[t] gets reset to 0 when one time unit has elapsed since a[t] gets set to t.
  • Most importantly, we can see that the solution is not perfect: between the second and third resets, x crosses z twice, and a gets set only once. The first crossing is not detected for some reason.

This last point illustrates a second possible weakness of this method: if there is a second crossing of x with z within one unit time step (i.e. before x gets reset), then a might get over-written with the new time, and so x won't reset until one unit time step after the new time.

To illustrate this last problem, change the first two lines in the code above to

timeLimit = 4;
z[t_] = Piecewise[{{1, t <= 2}, {2, 2 < t <= 5}}, 2];

and run the rest of it. This results in:

enter image description here

It's worse than I feared! a[t] gets set by WhenEvent to the current time t when x crosses z at about t == 1.8. a[t] gets set again at t == 2 when z jumps from having a value of 1 to having a value of 2; in other words, this is interpreted as a crossing between x and z even though z exhibits a discontinuity here! Finally, a[t] gets reset one more time when x crosses the new value of z, which occurs before t == 3, and since a[t] was set to 2 at the second "crossing", x won't get reset before t == 3. Finally, x resets at about t == 3.5.

The upshot: this is a good idea, but I don't know exactly how to make it robust.

Improve this answer
$\endgroup$
2
  • $\begingroup$ Thanks for the suggestion - unfortunately, this doesn't extend easily to something like WhenEvent[x1[t-1]+x2[t-3]==z[t],{...}]. I will try to post a more complicated example to demonstrate why I'm so fixated on putting time delays in the WhenEvent. $\endgroup$
    – polkabug
    Nov 17, 2015 at 7:52
  • $\begingroup$ @polkabug. Yeah, it seems like this method wouldn't generalize to that use case. Thus, the best bet is demonstrated by Michael E2, but as I noted, it seems to break when the new variables try to reset. I don't yet know how to get around this. $\endgroup$
    – march
    Nov 17, 2015 at 17:45
1
$\begingroup$

Almost a solution: I feel that one ought to be able to integrate the delayed values of the variables via equations such as

x1[t] == x[t - 1], y2[t] == y[t - 2]

and use the variables in WhenEvent. It works in the following:

z[t_] = Piecewise[{{1/22, t <= 1}}, 2];
foo = NDSolve[{x'[t] == 2 t, y'[t] == x'[t] - x'[t - 1],
    x1[t] == x[t - 1], y2[t] == y[t - 2],
    WhenEvent[x1[t] + y2[t] > z[t], {x[t] -> 0, y[t] -> 0}], 
    x[t /; t <= 0] == 0, y[t /; t <= 0] == 1/2,
    x1[t /; t <= 0] == 0, y2[t /; t <= 0] == 1/2},
   {x, y, x1, y2}, {t, 0, 3}];

ListLinePlot[{x, y, x1, y2} /. First[foo]]

Mathematica graphics

But that seems to be an accident. Extend the integration to {t, 0, 4}, and you'll get

NDSolve::nderr: Error test failure at t == 3.20710674525188`; unable to continue. >>

Integrate x'[t] == 2 Sqrt[t] instead of x'[t] == 2 t and I get a kernel crash!! Maybe it's a bug??

Improve this answer
$\endgroup$
1
  • $\begingroup$ I noticed this too, which is why I went the DiscreteVariable route. It crashes when the new variables are supposed to reset based on when x was reset. $\endgroup$
    – march
    Nov 17, 2015 at 15:35

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.