In the system below, would like to keep z[t] between 0 and 1. The intent of the code below is to use WhenEvent to detect when z[t] reaches the limit of 1, and then only allow negative values of the derivative, so that z[t] could decrease but not increase. However, in the plots below, z[t] continues to increase. Perhaps some user error, but I haven't been able to identify it.
df = 4. (0.07 z[t] Sqrt[600. - p[t]] - 0.005 Sqrt[p[t] p[t] - 100.] );
de = {p'[t] == df, z'[t] == -0.3 df + 0.4 (170. - p[t]) };
ic = {p[0] == 140., z[0] == 0.5};
events = {WhenEvent[z[t] > 1,
z'[t] == Min[0, -0.3 df + 0.4 (170. - p[t])] ]};
eqs = Flatten[{de, ic, events}];
solODE = NDSolve[eqs, {p, z}, {t, 0, 20}];
Plot[p[t] /. solODE, {t, 0, 20}, PlotRange -> {All, All}, PlotLabel -> "P[t]" ]
Plot[z[t] /. solODE, {t, 0, 20}, PlotRange -> {All, All}, PlotLabel -> "Z[t]" ]
Follow-up
3 methods that produce the intended results have been identified.
tmax = 20;
(*WhenEvent and appropriate expressions*)
df = 4. (0.07 z[t] Sqrt[600. - p[t]] - 0.005 Sqrt[p[t] p[t] - 100.]);
dz = (-0.3 df + 0.4 (170. - p[t]));
de = {p'[t] == df, z'[t] == in[t] dz};
ic = {p[0] == 140., z[0] == 0.5, in[0] == 1};
events = WhenEvent[event, action] /. {
{event -> z[t] > 1, action -> in[t] -> 0},
{event -> dz < 0 && z[t] == 1, action -> in[t] -> 1}};
eqs = Flatten[{de, ic, events}];
{pFuncA, zFuncA} = {p[t], z[t]} /. First@NDSolve[eqs, {p, z, in}, {t, 0, tmax}, DiscreteVariables -> {in}];
(*Piecewise and appropriate method for NDSolve*)
dpRHS[z_, p_] := 4. (0.07 z Sqrt[600. - p] - 0.005 Sqrt[p p - 100.])
dzRHS[z_, p_] := Module[{dzVal},
dzVal = -0.3 dpRHS[z, p] + 0.4 (170. - p);
Piecewise[{{Max[0, dzVal], z <= 0}, {Min[0, dzVal], z >= 1}}, dzVal]]
de = {p'[t] == dpRHS[z[t], p[t]], z'[t] == dzRHS[z[t], p[t]]};
ic = {p[0] == 140., z[0] == 0.5};
eqs = Flatten[{de, ic}];
{pFuncB, zFuncB} = {p[t], z[t]} /. First@NDSolve[eqs, {p, z}, {t, 0, 20}, Method -> {"DiscontinuityProcessing" -> False}];
(*Piecewise and appropriate form for test conditions*)
dpRHS[z_, p_] := 4. (0.07 z Sqrt[600. - p] - 0.005 Sqrt[p p - 100.])
dzRHS[z_, p_] := Module[{dzVal},
dzVal = -0.3 dpRHS[z, p] + 0.4 (170. - p);
Piecewise[{{0, (z <= 0 && dzVal < 0) || (z >= 1 && dzVal > 0)}},
dzVal]];
de = {p'[t] == dpRHS[z[t], p[t]], z'[t] == dzRHS[z[t], p[t]]};
ic = {p[0] == 140., z[0] == 0.5};
eqs = Flatten[{de, ic}];
{pFuncC, zFuncC} = {p[t], z[t]} /. First@NDSolve[eqs, {p, z}, {t, 0, tmax}];
(*results*)
plotP = Plot[{pFuncA, pFuncB, pFuncC, 170}, {t, 0, tmax}, PlotRange -> {All, All}, PlotLabel -> "P[t]",
PlotStyle -> {Sequence, Sequence, Sequence, {Red, Dashing[0.01]} }, ImageSize -> 400];
plotZ = Plot[{zFuncA, zFuncB, zFuncC}, {t, 0, tmax}, PlotRange -> {All, All}, PlotLabel -> "Z[t]", ImageSize -> 400];
{plotP, plotZ}
WhenEvent
is off - that should beWhenEvent[…,y[t]'->…]
.NDSolve
will complain however that you can't set the highest derivative. See @AlexTrounev's answer for a workaround $\endgroup$ – Lukas Lang Jan 5 at 22:10