6

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.

J. M.'s missing motivation
  • 124,525
  • 11
  • 401
  • 574
polkabug
  • 61
  • 2
  • 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. – march Nov 16 '15 at 21:41
  • 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. – polkabug Nov 17 '15 at 05:49
  • 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. – march Nov 17 '15 at 06:07

2 Answers2

4

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.

march
  • 23,399
  • 2
  • 44
  • 100
  • 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. – polkabug Nov 17 '15 at 07:52
  • @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. – march Nov 17 '15 at 17:45
1

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??

Michael E2
  • 235,386
  • 17
  • 334
  • 747
  • 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. – march Nov 17 '15 at 15:35