Message Boards Message Boards

Why does NDSolve return a constant value when the derivative is positive?

Posted 5 years ago

In the code below, a system of two first order non-linear differential equations are integrated. Would like some help understanding the result. There are two variables, z[t] and p[t]. The results appear to settle out with p[t] about 168 and z[t] about 0.99, and then after some time the system goes unstable and p begins to increase and z begins to oscillate.

My primary concern is the time range when p[t[ and z[t] appear to be constant. Plugging the values of z = 0.99 and p = 168 into the expression for dp/dt produces a value of about 2.4, so it is not clear why p[t] is staying constant during this period.

Below is the code.

dpRHS[z_, p_] := 4. (0.07 z Sqrt[600. - p] - 0.005 Sqrt[p p - 100.])

dzRHS[z_, p_] := Module[{dfVal, dzVal},
  dfVal = 4. (0.07 z Sqrt[600. - p] - 0.005 Sqrt[p p - 100.]);
  dzVal = -0.3 dfVal + 0.4 (170. - p);
  Piecewise[{{Max[0, dzVal], z <= 0.01}, {Min[0, dzVal], z >= 0.99}}, 
   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}]
tmax = 20;
{res, steps} =  Reap[solODE = NDSolve[eqs, {p, z}, {t, 0, tmax},   StepMonitor :> Sow[{t, z[t], p[t]}]]];

Plot[{p[t] /. solODE, 170}, {t, 0, tmax}, PlotRange -> {All, All},  PlotLabel -> "P[t]", PlotStyle -> {Sequence, {Red, Dashing[0.01]} }]
Plot[z[t] /. solODE, {t, 0, tmax}, PlotRange -> {All, All},  PlotLabel -> "Z[t]"]

(*inspect details of the numeric soluion*)
tzVals =  steps[[1, All, {1, 2}]];
tpVals  =  steps[[1, All, {1, 3}]];
tdpVals = {#[[1]],  dpRHS[#[[2]], #[[3]]  ]} & /@ steps[[1, All ]];

ListPlot[ tzVals , PlotRange -> {{12, 15}, All}, PlotLabel -> "Z[t]" ]
ListPlot[ tpVals , PlotRange -> {{12, 15}, {165, 175}},  PlotLabel -> "P[t]", 
 Epilog -> {Red, Dashing[0.01], Line[{{0, 170}, {50, 170}}]}]
ListPlot[ tdpVals , PlotRange -> {{12, 15}, All},  PlotLabel -> "P'[t]   (computed from {z,p}" ]

Below are the plots of p[t] and z[t]. Which should p[t] increasing, but then stops increasing. Plots of Z and P

Note that for the time range from 10 to 15, the value of p appears to settle out near 168 and the value of z appears is first just below 1. Then at later times, the value of z becomes unstable. At the moment I am not too concerned about the instability (the coefficients can be changed to remove the instability.) At the moment I am concerned about the value of p for times between 10 and 15, when p appears to be constant; however, with a value of z around 0.99 and a value of p around 168, dp/dt should be around 2.4, so I would expect p to continue increasing. (Then when p > 170, p should begin to decrease. The system is intended to have a steady-state value of p = 170 and z = 0.58).

Additional information: this code was written to limit the value of z between 0 and 1. So the equation for the dz/dt is given as a piecewise function. When z gets close to 0, then dz/dt is only allowed to be greater than or equal to zero. When z gets close to 1, then dz/dt is only allowed to be less than or equal to zero.

Additional information, I used Reap/Sow, to investigate the individual steps that NDSolve was using. The results were similar to what one sees when inspecting the plot.

Attachments:
POSTED BY: Robert McHugh
2 Replies

Just an observation: your dzRHS is discontinuous, and it can take negative values:

Plot[dzRHS[z, 175], {z, -2, 2}, PlotRange -> All]
POSTED BY: Gianluca Gorni
Posted 5 years ago

That is interesting. Thanks for pointing it out. In the range of operation, e.g. 0 <= z <= 1 and 140 < p < 180 (pressure limits are arbitrary), the function dzRHS is continuous. Demonstrated by Plot3D[dzRHS[z, p], {z, 0, 1}, {p, 140, 180}].

It could be that z is getting outside of this region, in this example z might be getting slightly above one. The intent of the piecewise implementation of dzRHS was to keep z between 0 and 1. Maybe this approach isn't successful at limiting z.

Is it possible that with the WhenEvent implementation, NDSolve looks more carefully at the boundary (z near 1 in this example).

Your thoughts?

POSTED BY: Robert McHugh
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract