2
$\begingroup$

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

z(t) p(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}

Results for all three methods are similar. enter image description here

|improve this question
$\endgroup$
  • $\begingroup$ Your syntax for the action of WhenEvent is off - that should be WhenEvent[…,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
  • $\begingroup$ "Eventually, the system goes unstable -- which is actually the point of this exercise" -- is the point to demonstrate the instability or to eliminate it? $\endgroup$ – Michael E2 Jan 6 at 20:41
  • $\begingroup$ this exercise was to demonstrate the instability, subsequent exercise will be to eliminate it. I observed the instability in another context (outside of Mathematica) and was trying to reproduce it in Mathematica. What was new to me, was that without any limits on z, (e.g. if z can go above 1 and if z can move very fast) then the system appears to be stable. However, the physical system does limit the value of z between 0 and 1 and it does limit how fast z can move. This question was about how to add the limit on the value of z. $\endgroup$ – user6546 Jan 6 at 22:05
  • $\begingroup$ I don’t know, my answer doesn’t show any instability. It does keep z[t]<1 though. $\endgroup$ – Chris K Jan 6 at 22:17
  • $\begingroup$ thanks for pointing this out. investigating the reason for this. $\endgroup$ – user6546 Jan 7 at 16:54
4
$\begingroup$

To constrain NDSolve to keep one of the variables in bounds, I add an indicator variable in[t] that changes from 1 to 0 at the boundary and back to 1 when z'[t] becomes negative again, using two WhenEvents.

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};
tmax = 20;
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}];
solODE = NDSolve[eqs, {p, z, in}, {t, 0, tmax}, DiscreteVariables -> {in}];
Plot[p[t] /. solODE, {t, 0, tmax}, PlotLabel -> "P[t]"]
Plot[z[t] /. solODE, {t, 0, tmax}, PlotLabel -> "Z[t]"]
Plot[in[t] /. solODE, {t, 0, tmax}, PlotLabel -> "in[t]"]

Mathematica graphics Mathematica graphics Mathematica graphics

Here are the dynamics superimposed on the phase plane (including isoclines), with the boundary at z=1 indicated:

Show[
  myStreamPlot[{dz, df}, {z[t], 0.5, 1.05}, {p[t], 140, 175}, StreamPoints -> Fine, StreamStyle -> Gray],
  ContourPlot[{dz == 0, df == 0}, {z[t], 0.5, 1.05}, {p[t], 140, 175}],
  Graphics[Line[{{1, 140}, {1, 175}}]],
  ParametricPlot[{z[t], p[t]} /. solODE, {t, 0, tmax}, PlotStyle -> Pink],
  FrameLabel -> {z, p}
]

Mathematica graphics

myStreamPlot is from here, originally by @Rahul.

As to the weird way to define the WhenEvents, this answer by Mark McClure notes that WhenEvent has the attribute HoldAll. This answer by Michael E2 explains why the order dz < 0 && z[t] == 1 is required in the second WhenEvent.

|improve this answer
$\endgroup$
  • $\begingroup$ Thanks! this provides some direction and tools for going forward. $\endgroup$ – user6546 Jan 6 at 20:32
2
$\begingroup$

It can be easier

df = 4.*(0.07 z[t] Sqrt[600. - p[t]] - 0.005 Sqrt[p[t]^2 - 100.]);
F = If[z[t] <= 
    1, -0.3*4.*(0.07 z[t] Sqrt[600. - p[t]] - 
       0.005 Sqrt[p[t]^2 - 100.]) + 0.4 (170. - p[t]), 
   Min[0, -0.3*4.*(0.07 z[t] Sqrt[600. - p[t]] - 
        0.005 Sqrt[p[t]^2 - 100.]) + 0.4*(170. - p[t])]];
de = {p'[t] == df, z'[t] == F};
ic = {p[0] == 140., z[0] == 0.5};
eqs = Flatten[{de, ic}];
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]"], 
 Plot[z[t] /. solODE, {t, 0, .1}, PlotRange -> {All, All}, 
  PlotLabel -> "Z[t]"]}

fig1 As pointed out by Chris, the solution depends on the method

solODE = NDSolve[eqs, {p, z}, {t, 0, 20}, 
  Method -> {"DiscontinuityProcessing" -> False}]

fig2

|improve this answer
$\endgroup$
  • $\begingroup$ any thoughts about why z[t] did not settle out to the same value as the original posting? $\endgroup$ – user6546 Jan 5 at 22:15
  • $\begingroup$ Why should it be like in your mistaken solution of equations? $\endgroup$ – Alex Trounev Jan 5 at 22:46
  • $\begingroup$ because when p[t] is 170, for p'[t] to be zero, then z[t] should be 0.58. Put another way, p'[t] = df, and NSolve[ 0 == df /. p[t] -> 170, z[t]], yields {{z[t] -> 0.584567}} $\endgroup$ – user6546 Jan 6 at 0:26
  • $\begingroup$ the system eventually goes unstable, which I suspected, but wasn't seeing when z was increasing beyond 1. There is a delay for the instability to occur, probably because the NDSolve is taking increasingly smaller step sizes. The instability occurs much earlier when a simple solutions method, like ExplicitEuler, is used. $\endgroup$ – user6546 Jan 6 at 1:49
  • 1
    $\begingroup$ Interestingly, if you add Method -> {"DiscontinuityProcessing" -> False} to your NDSolve it matches my solution. $\endgroup$ – Chris K Jan 7 at 0:24

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.