1
$\begingroup$

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. >>

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – AccidentalFourierTransform Jul 1 at 11:57
  • $\begingroup$ You're using Piecewise blindly. Please start from checking document of WhenEvent. $\endgroup$ – xzczd Jul 1 at 11:59
3
$\begingroup$

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}]

Mathematica graphics

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

$\endgroup$

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.