I am trying to solve the following system of differential equations using NDSolve:



subject to the following constraint:


where $z_t$ essentially follows a predetermined path, and $x_t=0.75$ until the constraint becomes binding, then dynamically adjusts to keep $y_t$ at the level constrained by $z_t$.

I am using the following code:


which generates the following:

enter image description here

Basically, I am using a WhenEvent trigger to change the value of $x_t$ once $y_t$ hits the constraint (represented by the dotted line). This works fine, but then I try to force NDSolve to switch the expression used for $\dot{x}_t$ to the one that would make $y_t$ follow the constraint path going forward by using a Piecewise function, but this doesn't seem to work for some reason.

Alternatively, I've been able to achieve what I want using "StopIntegration" in a WhenEvent trigger to shut down NDSolve when it hits the constraint, then using the values of the state variables at this point as the initial values in a second NDSolve with the appropriate constraint expression for $\dot{x}_t$. However, I then have to stitch together the two sets of results from the two NDSolve commands, which is a bit cumbersome, so I'm hoping there's a way to do this within a single NDSolve.

Any ideas?

  • $\begingroup$ With bind = 0 before xdot, xdot is defined to be 0. If bind is defined after xdot, then xdot is defined in terms of xdot, which causes an infinite recursion. $\endgroup$
    – Michael E2
    Apr 26 '16 at 0:46
  • $\begingroup$ Thanks for pointing that out; I made an error when converting from a more complicated version of the problem. I believe I have correctly defined xdot now (and define bind=0 at the appropriate time). However, the central issue remains unresolved. $\endgroup$
    – clr66
    Apr 26 '16 at 17:45
  • $\begingroup$ In other words, it doesn't matter what I put in as the expression for xdot when bind==1 in the Piecewise function, NDsolve doesn't generate any dynamics in x[t] after the WhenEvent triggers. $\endgroup$
    – clr66
    Apr 26 '16 at 17:54

Changing parameter values during integration works better with DiscreteVariables. But I think the problem with OP's code, in the question and the OP's answer, has more to do with Mathematica numerics.

My solution

zdot = 1/2 (1 - z[t]);
ydot = 1/20*y[t] + z[t] - x[t];
xdotbind = 
  D[Solve[-ydot - zdot == 0, x[t]][[1, 1, 2]], t] /. {y'[t] -> ydot, 
    z'[t] -> zdot};
xdot = Piecewise[{{0, bind[t] == 0}, {xdotbind, bind[t] == 1}}];
{sol, {data}} = Reap@NDSolve[{
     x'[t] == xdot, y'[t] == ydot, z'[t] == zdot,
     bind[0] == 0, x[0] == 3/4, y[0] == 0, z[0] == 1/10,
     WhenEvent[-y[t] - z[t] == 0,
      {x[t] -> 1/20*y[t] + z[t] + 1/2*(1 - z[t]), bind[t] -> 1, 
    {x, y, z}, {t, 0, 5},
    DiscreteVariables -> {bind \[Element] {0, 1}},
    StepMonitor :> 
     Sow[{bind[t], t, {x[t], y[t], z[t]}, {xdot, ydot, -zdot}}]];

Grid[{{Plot[Evaluate[x[t] /. sol], {t, 0, 5}, PlotRange -> All, 
    AxesLabel -> {"t", "x"}], 
   Plot[{Evaluate[y[t] /. sol], Evaluate[-z[t] /. sol]}, {t, 0, 5}, 
    PlotRange -> All, AxesLabel -> {"t", "y"}, 
    PlotStyle -> {Automatic, {Gray, Dashed}}], 
   Plot[Evaluate[z[t] /. sol], {t, 0, 5}, PlotRange -> All, 
    AxesLabel -> {"t", "z"}]}}]

Mathematica graphics

Numerical issues

The OP puts a finger on the problem with the observation in the OP's answer:

if I change the region condition in the definition of xdot from y[t] > -z[t] to y[t]+z[t] > 0 (which seems equivalent to me...)

The inequalities are mathematically equivalent but not numerically equivalent. This can be demonstrated by the following:

1. + $MachineEpsilon > 1.
(*  False  *)

(1. + $MachineEpsilon) - 1. > 0
(*  True  *)

This is because there is a tolerance in comparing numbers. If they are approximately equal, as defined by Internal`$EqualTolerance, then strict comparison will return False.

Block[{Internal`$EqualTolerance = 0.},
 1. + $MachineEpsilon > 1.
(*  True  *)

The documentation page for Less points out another problem, rounding error:

0.00001 < 2.00006 - 2.00005
(*  True  *)

0.00001 + 2.00005 == 2.00006
(*  True  *)

The problem is almost bound to arise in the OP's answer with the variant

xdot = Piecewise[{{0, y[t] + z[t] > 0}}, xdotbind]

because NDSolve is trying to tiptoe along the path where y[t] == z[t]. If we look at the two forms of the condition in xdot for the OP's first solution (the correct one), we can see that the comparison is an issue:

grid = x["Grid"] /. sol // Flatten; (* the times of the steps *)

x'["ValuesOnGrid"] /. sol // First; (* values of the derivative   x'  *)
Split[Flatten@Position[%, 0.],      (* positions where  x'[t] == 0  *)
  Subtract[##] == -1 &] /. {a_Integer, ___, b_Integer} :> a ;; b
(*  {1 ;; 38}  *)

y[grid[[38 ;; -1]]] > -z[grid[[38 ;; -1]]] /. sol // First // Thread
y[grid[[38 ;; -1]]] + z[grid[[38 ;; -1]]] > 0 /. sol // First // Thread

(*  y[t] > -z[t]  from the correct solution
{True, False, False, False, False, False, False, False, False, False, 
... 30 False's ...
False, False, False, False, False, False, False}

    y[t] + z[t] > 0  from the alternative
{True, False, False, False, False, False, False, False, False, True, 
... 33 True's ...
True, True, True, True}

The first True to False change is where the jump in x[t] occurs around t == 0.73. In the correct solution, the condition is False at all subsequent steps, which means the x'[t] == xdotbind and the solution has y track z. The alternative condition is False for a few steps and then is True. This would set x'[t] == 0 and y would stop tracking z. This is what happens when NDSolve is run with the alternative condition defining xdot.

  • $\begingroup$ This was tremendously informative and beneficial. Thanks so much, @Michael $\endgroup$
    – clr66
    Apr 27 '16 at 13:49

So, after playing around with this problem all day, I finally stumbled upon a version of the code that gives me my desired result, although I'm not at all clear on why it works where other versions fail.

First, here is the working version:

xdot = Piecewise[{{0, y[t] > -z[t]}}, xdotbind];

Which gives the following output:

enter image description here

The two important changes seem to have been the following:

1) Abandoning the use of the bind variable as a method of switching between the Piecewise definitions of xdot, and instead defining the change in xdot explicitly on the values of the other state variables, and

2) Including the Method->{"DiscontinuityProcessing"->False} option in NDSolve.

However, the solution seems incredibly sensitive to the manner in which the binding condition is defined, i.e. if I change the region condition in the definition of xdot from y[t]>-z[t] to y[t]+z[t]>0 (which seems equivalent to me...), I get the following incorrect output:

enter image description here

Why is it being so finicky? Am I overlooking something obvious?

In any case, I hope someone else finds this useful.


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.