# Piecewise with differential equations

I have the following differential system. The calibration I use is as follows

paramFinal = {rho -> 0.025, R -> 0.25, alpha -> 0, k -> -0.55, g -> 1.2, C0 -> 3, C1 -> 0.1, sbar -> 10, \[Eta] -> 8.5, hbar -> 0.5, Hbar -> 23.5};


The differential system is

dec = W'[t] == -( rho + k*rho*C1 *H[t]/W[t] - k C1 R/W[t] - rho (g + k C0)/W[t]);

des = H'[t] ==  -(R - (1 - alpha) W[t]);


With NDSolve, I can solve the problem without any problem but I would like to restrict the variable $$H[t]$$ such that it does not go below $$Hbar$$ and once it hits $$Hbar$$, it stays there. I try to use Piecewise function as follows

pdec = Piecewise[{W'[t] == -(rho + k*rho*C1 *H[t]/W[t] - k C1 R/W[t] - rho (g + k C0)/W[t]), Hbar < H[t], {Hbar, Hbar > H[t]}}];

pdes = Piecewise[{H'[t] ==  -(R - (1 - alpha) W[t]), Hbar < H[t], {Hbar, Hbar > H[t]}}];

nsba3aa = NDSolve[{pdec /. paramFinal, pdes /. paramFinal, W[0] == 0.275, H[0] == 25, {W[t], H[t]}, {t, 0, tmax}]


I have the following error message and I cannot figure out where is the problem

Piecewise::pairs: The first argument {(W^[Prime])[t]==-rho+(C1 k R)/W[t]+((g+C0 k) rho)/W[t]-(C1 k rho H[t])/W[t],HbarH[t]}} of Piecewise is not a list of pairs. >>

• H[t] appears to be an increasing function, and H[0]>Hbar, so you don't need to force H[t]>Hbar: this is implemented automatically. Jul 1, 2019 at 11:57
• You're using Piecewise blindly. Please start from checking document of WhenEvent. Jul 1, 2019 at 11:59

There are a lot of little errors in your code:

1. As the error message suggests, Piecewise takes a list of pairs as an argument. The structure in your code is Piecewise[{eqn1, cond1, {eqn2, cond2}}]; it needs to be Piecewise[{{eqn1, cond1}, {eqn2, cond2}}] instead.

2. The Piecewise should be on the right hand side, since the eqn2 isn't a differential equation by itself. Use W'[t] == Piecewise[...] instead.

3. I think you want eqn2 to be 0 not Hbar to keep H[t]>=Hbar.

4. You need to define tmax.

5. The initial conditions need to be in the same list as the differential equations in NDSolve.

6. I prefer to solve for {W, H} not {W[t], H[t]}.

Together:

tmax = 100;
pdec = W'[t] == Piecewise[{
{-(rho + k*rho*C1*H[t]/W[t] - k C1 R/W[t] - rho (g + k C0)/W[t]), Hbar < H[t]},
{0, Hbar > H[t]}}];

pdes = H'[t] == Piecewise[{
{-(R - (1 - alpha) W[t]), Hbar < H[t]}, {0, Hbar > H[t]}}];

nsba3aa = NDSolve[{pdec /. paramFinal, pdes /. paramFinal,
W[0] == 0.275, H[0] == 25}, {W, H}, {t, 0, tmax}][[1]];


We can't see that it works, since H[t]>Hbar in this example:

Plot[Evaluate[H[t] /. nsba3aa], {t, 0, tmax}]


but if you flip the sign of H'[t] you will find that it does what you ask: