4
$\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

$\endgroup$
5
  • $\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, 2019 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, 2019 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, 2019 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, 2019 at 22:17
  • $\begingroup$ thanks for pointing this out. investigating the reason for this. $\endgroup$
    – user6546
    Jan 7, 2019 at 16:54

2 Answers 2

5
$\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.

$\endgroup$
1
  • $\begingroup$ Thanks! this provides some direction and tools for going forward. $\endgroup$
    – user6546
    Jan 6, 2019 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

$\endgroup$
13
  • $\begingroup$ any thoughts about why z[t] did not settle out to the same value as the original posting? $\endgroup$
    – user6546
    Jan 5, 2019 at 22:15
  • $\begingroup$ Why should it be like in your mistaken solution of equations? $\endgroup$ Jan 5, 2019 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, 2019 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, 2019 at 1:49
  • 2
    $\begingroup$ Interestingly, if you add Method -> {"DiscontinuityProcessing" -> False} to your NDSolve it matches my solution. $\endgroup$
    – Chris K
    Jan 7, 2019 at 0:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.