Switching Differential Equation in NDSolve

Question

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

$\dot{z}_t=.5(1-z_t)$

$\dot{y}_t=.05y_t+z_t-x_t$

subject to the following constraint:

$-y_t-z_t\le0$

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:

zdot=.5*(1-z[t]);
ydot=.05*y[t]+z[t]-x[t];
xdotbind=D[Solve[-ydot-zdot==0,x[t]][[1,1,2]],t];
xdot=Piecewise[{{0,bind==0},{xdotbind,bind==1}}];
bind=0;
sol=NDSolve[{x'[t]==xdot,y'[t]==ydot,z'[t]==zdot,x[0]==.75,y[0]==0,z[0]==0.1,WhenEvent[-y[t]-z[t]==0,{x[t]->.05*y[t]+z[t],bind=1,"RemoveEvent"}]},{x,y,z},{t,0,25}];
Grid[{{Plot[Evaluate[x[t]/.sol],{t,0,2.5},PlotRange->All,AxesLabel->{"t","x"}],Plot[{Evaluate[y[t]/.sol],-Evaluate[z[t]/.sol]},{t,0,2.5},PlotRange->All,AxesLabel->{"t","y"},PlotStyle->{Automatic,{Gray,Dashed}}],Plot[Evaluate[z[t]/.sol],{t,0,25},PlotRange->All,AxesLabel->{"t","z"}]}}]

which generates the following:

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?

Answer 1

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

Clear[bind];
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, 
       "RemoveEvent"}]},
    {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"}]}}]

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.

Answer 2

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:

zdot=.5*(1-z[t]);
ydot=.05*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, y[t] > -z[t]}}, xdotbind];
sol=NDSolve[{x'[t]==xdot,y'[t]==ydot,z'[t]==zdot,x[0]==.75,y[0]==0,z[0]==0.1,WhenEvent[y[t]==-z[t],{x[t]->0.05*y[t]+0.5*(1-z[t])+z[t],"RemoveEvent"}]},{x,y,z},{t,0,10},Method->{"DiscontinuityProcessing"->False}];
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"}]}}]

Which gives the following output:

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:

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

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