# How to use WhenEvent with a vector ODE in NDSolve

I have an ODE system I'd like to specify as a vector equation in NDSolve. I'm not clear on how to use WhenEvent for a system specified in this way. Ultimately I'd like WhenEvent to change the value of one of the state variables in the vector equation when it reaches a threshold.

For example, the following code doesn't work because WhenEvent[y[t]==0, ...] is applying a conditional test to a list of values. Replacing this with y[t][[1]]==0 also doesn't work.

sol = NDSolve[{y'[t] == {{.1, -.2}, {-.1, .2}}.y[t], y[0] == {1, 1},
WhenEvent[y[t] == 0, y[t] -> 1]}, y, {t, 0, 10}]


Ideas anyone?

s = NDSolve[{y'[t] == {{.1, -.2}, {-.1, .2}}.y[t], y[0] == {1, 1},
WhenEvent[Norm[y[t] - {0.9460552574072016, 1.053944742592798}] <= .01, y[t] -> {1, 1}]}
, y[t], {t, 0, 1}]
Plot[y[t] /. s[[1]] /. t -> u, {u, 0, 1}]


• Although this works, the Norm[] is an attribute of the whole system of state variables. This did give me the idea to try Part[y[t],1]<=0.1, which seems to work. Not sure why Part[] works when y[t][[1]] doesn't... – dvasseur Oct 31 '13 at 18:58
• @dvasseur y[t][[1]] works OK. Your error was using == instead of <= – Dr. belisarius Oct 31 '13 at 19:04

The proper way to generate events is to symbolically refer to vector elements, that is: using Indexed (since v10). Here I assume that the vector length n can change (in that case one must adjust the matrix used in y'[t] also), and I also use a programmatic way to set up events for all vector elements (see this answer).

ClearAll[y, t];

n = 2; (* vector length *)
x = 0; (* event threshold *)
events = WhenEvent[Indexed[y[t], #] == x,
y[t] -> ReplacePart[y[t], # -> 1]] & /@ Range@n;
sol = NDSolve[{
y'[t] == {{.1, -.2}, {-.1, .2}}.y[t],
y[0] == Table[1, {n}],
events}, y, {t, 0, 10}];

Plot[y[t] /. sol, {t, 0, 10}]


With Part instead of Indexed, premature evaluation of symbolic expressions (here Part[y[t], #] == x) during NDSolve fails. One must either use Indexed or hardcode the threshold value x into the WhenEvent.

ClearAll[y, t];

n = 2;
x = 0;
events = WhenEvent[Part[y[t], #] == x,
y[t] -> ReplacePart[y[t], # -> 1]] & /@ Range@n;
sol = NDSolve[{y'[t] == {{.1, -.2}, {-.1, .2}}.y[t],
y[0] == Table[1, {n}], events}, y, {t, 0, 10}];

Plot[Evaluate[y[t] /. sol], {t, 0, 10}]


Here is my solution, just change the WhenEvent part to WhenEvent[First@y[t] == 0, y[t] -> {1, Last@y[t]}]

sol = NDSolve[{y'[t] == {{.1, -.2}, {-.1, .2}}.y[t], y[0] == {1, 1},
WhenEvent[First@y[t] == 0, y[t] -> {1, Last@y[t]}]}, y, {t, 0, 10}]
Plot[Evaluate[y[t] /. sol], {t, 0, 10}]


Mathematica gives

• +1. I always feel a sense of relief when the result confirms that First@y[t] is evaluated the way you hoped it would be. :) – Michael E2 Jan 8 '15 at 18:28
• @MichaelE2 Thanks. I like something straightforward, and Mathematica rarely let me down:) – luyuwuli Jan 9 '15 at 1:28
• @MichaelE2 For the record, Part (First, Last, etc.) does not always work as one would expect, especially when it comes to symbolic evaluation, see my answer. – István Zachar Oct 11 '18 at 20:18
• @MichaelE2 Thanks for telling me this. However, I think in the case of fixed dimension, HoldAll attribute of WhenEvent is exactly the key to success. With this attribute, MMTC would first get the solution of explicit dimension, then Part or First without any problem. I admit this cannot handle the arbitrary dimension and it prevent us from encapsulating the whole events before the show up of NDSolve. – luyuwuli Oct 15 '18 at 8:59