# Switching Differential Equation in NDSolve

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==.75,y==0,z==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?

• 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. – Michael E2 Apr 26 '16 at 0:46
• 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. – clr66 Apr 26 '16 at 17:45
• 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. – 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

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, x == 3/4, y == 0, z == 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.

• This was tremendously informative and beneficial. Thanks so much, @Michael – 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:

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==.75,y==0,z==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.