1
$\begingroup$

I want to change part of a set of equations by using WhenEvent in NDSolve. I have made a much less complicated minimum working example as follows:

sol = x /. First@NDSolve[{
      x''[t] +
        If[x[t] < 0,
         50 x[t] ,
         16 50 x[t] ] == 0,
      x'[0] == 0,
      x[0] == 1
      },
     {x}, {t, 0, 2}];
Plot[sol[t], {t, 0, 1}]

enter image description here

This does not use WhenEvent but does switch a term in the equation. Now I try and introduce a WhenEvent. I also introduce a module that contains the two versions of the term in the equation.

ClearAll[ff];
ff[status_, t_?NumericQ, x_?NumericQ] := 
  If[x < 0,
   Sow[{t, x, -1, status}, tag1]; 50 x ,
   Sow[{t, x,   1, status}, tag1]; 16 50 x ];

status = "More";
a = Reap[
   res = Reap[
      sol = x /. First@NDSolve[{
            x''[t] +
              ff[status, t, x[t]] == 0,
            x'[0] == 0, x[0] == 1,

            WhenEvent[x[t] == 0, 
             If[x[t] < 0, status = "Less", status = "More"]; 
             Sow[{t, x[t], status, "WhenEvent"}, tag2]]},
           {x}, {t, 0, 2}];,
      tag1];,
   tag2];
Plot[sol[t], {t, 0, 1}, 
 Epilog -> {Line[res[[2, 1, All, {1, 3}]]], Red, PointSize[0.01], 
   Point[a[[2, 1, All, {1, 2}]]]}]

enter image description here

This has also worked and the two Sow functions collect various values which shows everything is working. I have introduced a status value which shows which of the two states I am in. I would like to use the status value for controlling the switching of the term in the differential equation. So I modify the function ff so that the status is examined rather than the value of the variable. This does not work.

ClearAll[ff];
ff[status_, t_?NumericQ, x_?NumericQ] := 
  If[status == "Less",
   Sow[{t, x, -1, status}, tag1]; 50 x ,
   Sow[{t, x,   1, status}, tag1]; 16 50 x ];

status = "More";
a = Reap[
   res = Reap[
      sol = x /. First@NDSolve[{
            x''[t] +
              ff[status, t, x[t]] == 0,
            x'[0] == 0, x[0] == 1,

            WhenEvent[x[t] == 0, 
             If[x[t] < 0, status = "Less", status = "More"]; 
             Sow[{t, x[t], status, "WhenEvent"}, tag2]]},
           {x}, {t, 0, 2}];,
      tag1];,
   tag2];
Plot[sol[t], {t, 0, 1}, 
 Epilog -> {Line[res[[2, 1, All, {1, 3}]]], Red, PointSize[0.01], 
   Point[a[[2, 1, All, {1, 2}]]]}]

enter image description here

Why does this not work? How can I fix it?

$\endgroup$
2
$\begingroup$

I think you need to make status a DiscreteVariable in order to change it within NDSolve. In that case, you'll also need to use numeric values instead of "More" and "Less".

Does this work for you?

ClearAll[ff];
ff[status_, t_?NumericQ, x_?NumericQ] := 
  If[status == -1, Sow[{t, x, -1, status}, tag1]; 50 x, 
   Sow[{t, x, 1, status}, tag1]; 16 50 x];

Clear[status];
a = Reap[res = 
     Reap[sol = 
        x /. First@
          NDSolve[{x''[t] + ff[status[t], t, x[t]] == 0, x'[0] == 0, 
            x[0] == 1, status[0] == 1,            
            WhenEvent[
             x[t] < 0, {status[t] -> -1, 
              Sow[{t, x[t], status[t], "WhenEvent"}, tag2]}],
            WhenEvent[
             x[t] > 0, {status[t] -> 1, 
              Sow[{t, x[t], status[t], "WhenEvent"}, tag2]}]}, {x, 
            status}, {t, 0, 2}, DiscreteVariables -> {status}];, 
      tag1];, tag2];

Plot[sol[t], {t, 0, 1}, 
 Epilog -> {Line[res[[2, 1, All, {1, 3}]]], Red, PointSize[0.01], 
   Point[a[[2, 1, All, {1, 2}]]]}]

enter image description here

| improve this answer | |
$\endgroup$
  • $\begingroup$ Thank you. Yes its all there, with examples, if you look up DiscreteVariables.in help. Pity I had not come across it before. Hope this helps someone else who needs to do something similar. $\endgroup$ – Hugh Apr 4 at 18:39

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.