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.
x'[t] == 1
,y'[t]
will always be zero, and soy
is constant, unless it's supposed to know about the discontinuity inx[t]
that arises as a consequence of the resetx[t] -> 0
, but in that case, you actually have an impulse fory
(because of the infinite derivative at those points), and I'm not sure thatNDSolve
can deal with that. $\endgroup$WhenEvent
and that the equations (especiallyy[t]
) are not relevant. I am only interested in whether there is a way to get theWhenEvent
to recognizex[t-1]
as a function instead of directly substituting the most recent values. $\endgroup$a[t]
and includea[t] == x[t - 1]
anda[t/;t<=0]==0
in the list of equations. I triedNDSolve[{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, butNDSolve
had problems once it was time fora
to reset (one time unit afterx
resets). Perhaps we could make this work, but I'm not sure how to fix the resetting problem. $\endgroup$