How to explain NDSolve that functions must be positive only?

Question

I have a simple system of differential equations to solve but both functions should be positive: x[t] > 0 and y[t] > 0 for all t. I was not able to explain NDSolve this fact so it goes x[t] < 0 for large t's.

I read some related post here on how to impose somehow that NDSolve solutions stay positive (Constraining function found by NDSolve to stay positive and Define a function with the condition that it must be positive), but didn't work so I'm here asking for help.

My code is:

a = 1;
b = 1;
p = 8;
q = 5;
r = 1;
f[t_] := a*r*1/(1 + E^(-p*(t - q)))
g[t_] := -b*r*1/(1 + E^(-p*(t - q)))
x0 = 2;
y0 = 2;
Plot[f[t], {t, -10, 20}]
system = {
           (1)*x'[t] + (0)*y'[t] + (b/y[t]*x[t]) == g[t], 
           (0)*x'[t] + (1)*y'[t] + (-a/y[t]*x[t]) == f[t], 
           x[0] == x0, 
           y[0] == y0
          };
s = NDSolve[system, {x, y}, {t, 0, 10}];
Plot[Evaluate[y[t] /. s], {t, 0, 10}, PlotRange -> All]
Plot[Evaluate[x[t] /. s], {t, 0, 10}, PlotRange -> All]

After aprox t=5, x[t] becomes negative, which makes no sense physically. I tried "WhenEvent", with no success:

system = {
           (1)*x'[t] + (0)*y'[t] + (b/y[t]*x[t]) == g[t], 
           (0)*x'[t] + (1)*y'[t] + (-a/y[t]*x[t]) == f[t],
           x[0] == x0,
           y[0] == y0,
           WhenEvent[x[t] < 0, x[t] -> 0] 
};

I would like to understand how can I impose that my unknown functions have positive codomain. I didn't find a way to use "Assuming" or similar resources to impose codomain or image sets to yet unknown functions.

I would appreciate feedback from experienced Mathematica users.

Thank you a lot in advance

Ricardo Vêncio [email protected]

Answer 1

Concerning the follow-up question in a comment:

...just for curiosity, if someone knows how to impose a positive codomain to not yet known function, I would very much like to know.

One way is to define the system so that x'[t] is positive whenever x[t] < 0. Then when x[t] reaches 0, it will be stuck at 0, unless x'[t] becomes positive for x[t] > 0 in a neighborhood of x == 0. See the tutorial Events and Discontinuities in Differential Equations. Below I use Piecewise to define the OP's system for x[t] >= 0 and set x'[t] to the arbitrarily chosen positive value 1 for x[t] < 0.

a = 1;
b = 1;
p = 8;
q = 5;
r = 1;
f[t_] := a*r*1/(1 + E^(-p*(t - q)))
g[t_] := -b*r*1/(1 + E^(-p*(t - q)))
x0 = 2;
y0 = 2;

system = {
   (1)*x'[t] == Piecewise[{{g[t] - ((0)*y'[t] + (b/y[t]*x[t])), x[t] >= 0}}, 1],
   (0)*x'[t] + (1)*y'[t] + (-a/y[t]*x[t]) == f[t],
   x[0] == x0, y[0] == y0};

s = NDSolve[system, {x, y}, {t, 0, 10}];
GraphicsRow[{
  Plot[Evaluate[y[t] /. s], {t, 0, 10}, PlotRange -> All],
  Plot[Evaluate[x[t] /. s], {t, 0, 10}, PlotRange -> All]}]

Answer 2

I posted the solution I found in another thread (Constraining function found by NDSolve to stay positive) and I am reposting it here for convenience. However, in your particular case the model is not complete and you need to ensure that x'[t] is identically zero when x[t]<0. Here is your model with the necessary completion. The description is below.

ClearAll["Global`*"];
NnnX[x_?NumericQ] := ((x + Sqrt[x^2])/2);
UsbEEpsilonValue = 10^-3;
UnitStepBase[x_?NumericQ] := 
  Piecewise[{{0, x < 0}, {(2*x^2 - x^4), 0 <= x < 1}, {1, x >= 1}}];
UsbE[x_?NumericQ, norm_?NumericQ] := UsbE[x, norm, UsbEEpsilonValue];
UsbE[x_?NumericQ, norm_?NumericQ, eps_?NumericQ] := 
  UnitStepBase[x/(norm*eps)];
SomeNorm = 10;

f[t_] := a*r*1/(1 + E^(-p*(t - q)))
g[t_] := -b*r*1/(1 + E^(-p*(t - q)))
a = 1;
b = 1;
p = 8;
q = 5;
r = 1;
x0 = 2;
y0 = 2;
Plot[f[t], {t, -10, 20}]

system = {
  Derivative[1][x][t] == UsbE[x[t], SomeNorm]*(-b NnnX[x[t]] + g[t] NnnX[y[t]])/NnnX[y[t]],
  Derivative[1][y][t] == UsbE[y[t], SomeNorm]*(a NnnX[x[t]] + f[t] NnnX[y[t]])/NnnX[y[t]], 
  x[0] == x0, 
  y[0] == y0
};
s = NDSolve[system, {x, y}, {t, 0, 10}];
Plot[Evaluate[y[t] /. s], {t, 0, 10}, PlotRange -> All, Frame -> True, GridLines -> Automatic]
Plot[Evaluate[x[t] /. s], {t, 0, 10}, PlotRange -> All, Frame -> True, GridLines -> Automatic]

I had a very similar problem for a system of very large number of ODEs and none of the methods proposed above worked. The main reason is that the coefficient

If[x[t] >= 0, -10, 0] 

is a stiff function. Several bad things happen because of that.

1. For stiff system the step size decreases substantially and it may take forever to integrate. The simplest model I run has 11 variables and that's exactly what I observed: the integration never finished!

2. NDSolve cannot know that the coefficient will be exactly zero for x[t]<0. What happens is that as the complexity of the system increases there is a high chance that NDSolve will overshoot zero and make some variables negative. What happens next depends on your model. In my case the model usually "exploded" and the variables started to run to plus and minus infinities.

3. In many cases it is necessary to just impose positivity constraint and do not stop integration. In such case using WhenEvent and then stopping integration does not help as we need to continue integrating (while keeping a function at zero). At some point in time the conditions may change and the function may become positive again.

The solution I found involved several steps:

Any and all stiff functions should be replaced by "numerical" smoothed out functions. To do so I use a smoothed out UnitStep with some scale (norm) and stiffness (eps - default value =10^-3):

UsbEEpsilonValue=10^-3;
UnitStepBase[x_?NumericQ] := 
  Piecewise[{{0, x < 0}, {(2*x^2 - x^4), 0 <= x < 1}, {1, x >= 1}}];
UsbE[x_?NumericQ, norm_?NumericQ] := UsbE[x, norm, UsbEEpsilonValue];
UsbE[x_?NumericQ, norm_?NumericQ, eps_?NumericQ] := UnitStepBase[x/(norm*eps)];

So instead of

 If[x[t] >= 0, -10, 0] 

I would write:

-10*UsbE[x[t],SomeNorm]

with, let's say

SomeNorm = 10;

The second step involves making sure that negative values of x[t] do not affect any other calculations. To do so I wrap all x[t] inside the derivative function

x'[t]==f[x[t]]

into another "numerical" function, which ensures that for negative x[t] we get exactly zero:

NnnX[x_?NumericQ] := ((x + Sqrt[x^2])/2);

and then eventually use:

x'[t]==UsbE[x[t],SomeNorm]*f[NnnX[x[t]]]

or

x'[t]==f[NnnX[x[t]]]

depending on the behavior of the term f. A small remaining problem is that if NDSolve does overshoot and makes x[t]<0 then there is no way for it to come back (in the first case): it will effectively stay at zero. This is not good if you have some terms, which may at some point increase x[t]. To deal with that it is important to identify the stiff terms, which may result in negative x[t] and deal with them in the described above way. You then just need to leave the positive terms with only NnnX[x[t]] wrappers.

You might wonder why using numerical functions? Well, I don't know. What I noticed is that NDSolve tries to evaluate all symbolic functions (except the ones, which can only be evaluated with numerical values). And when NDSolve encounters a lot of conditional functions (like If or Piecewise) then it starts to behave unpredictably.