11
$\begingroup$

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?

Improve this question
$\endgroup$

3 Answers 3

Reset to default
7
$\begingroup$ 
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}]

Mathematica graphics

Improve this answer
$\endgroup$
2
  • $\begingroup$ 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... $\endgroup$
    – dvasseur
    Oct 31, 2013 at 18:58
  • $\begingroup$ @dvasseur y[t][[1]] works OK. Your error was using == instead of <= $\endgroup$
    – Dr. belisarius
    Oct 31, 2013 at 19:04
6
$\begingroup$

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}]

Mathematica graphics

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}]

Mathematica graphics

Improve this answer
$\endgroup$
5
$\begingroup$

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

enter image description here

Improve this answer
$\endgroup$
4
  • $\begingroup$ +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. :) $\endgroup$
    – Michael E2
    Jan 8, 2015 at 18:28
  • $\begingroup$ @MichaelE2 Thanks. I like something straightforward, and Mathematica rarely let me down:) $\endgroup$
    – luyuwuli
    Jan 9, 2015 at 1:28
  • 1
    $\begingroup$ @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. $\endgroup$
    – István Zachar
    Oct 11, 2018 at 20:18
  • $\begingroup$ @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. $\endgroup$
    – luyuwuli
    Oct 15, 2018 at 8:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.