# How to explain NDSolve that functions must be positive only?

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 == x0,
y == 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 == x0,
y == 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 rvencio@usp.br

• According to your differential equations, if x[t] and y[t] stay positive, x'[t] approaches -1 or less as t -> Infinity. (I get x'[t] == -(1/(1 + E^(-8 (-5 + t)))) - x[t]/y[t].) Hence x would decrease to negative infinity over time. (Of course when the sign of x changes, it alters things.) If x should be positive, perhaps you need to check your diff. eq. – Michael E2 Sep 20 '14 at 11:43
• @MichaelE2 I am a little confused also: I assume that this functionality (if it really exists) is only useful for very singular ODE - otherwise, by the existence and uniqueness, you don't get to indicate the sign of the solution. – Igor Rivin Sep 20 '14 at 16:48
• @IgorRivin The functionality can be used to model discrete events such as switches or other controls, or, say, when you drop an egg on the floor. From the description in the question, that does not seem to be the case here. – Michael E2 Sep 20 '14 at 20:07
• I would like to thank all those who took the time to answer me and try to help, not only on the Mathematica questions itself but also on the mathematical question. Both x[t] and y[t] are quantities that only make sense if positives, that is why I was trying to impose that somehow. If the system solution goes negative, makes no physical sense anymore even if it is allowed mathematically. Maybe the best way to go would be just stop everything when reach x[t] < 0 since the model itself does not apply anymore. Thank you all. – Vêncio Sep 21 '14 at 14:14
• Out of curiosity, what is the point of the coefficients (1) and (0) in front of the derivatives? They make things seemingly unnecessarily complicated. Are they just particular values for general parameters in your problem? I would have left them out of my answer, but I wanted to show how to deal with (0)*y'[t] in the case the coefficient is nonzero. – Michael E2 Oct 15 '14 at 13:40

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 == x0, y == 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]}] 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[x][t] == UsbE[x[t], SomeNorm]*(-b NnnX[x[t]] + g[t] NnnX[y[t]])/NnnX[y[t]],
Derivative[y][t] == UsbE[y[t], SomeNorm]*(a NnnX[x[t]] + f[t] NnnX[y[t]])/NnnX[y[t]],
x == x0,
y == 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)];


 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.

• @MichaelE2 Sorry I did not get the question. Are you asking if looked through other solution in this thread and other threads? Yes, I did. None worked. Neither WhenEvent, nor If[...], nor anything else. The only possibly working solution (in my case) is the example above. For some other models there might be other solutions. – Konstantin Konstantinov Oct 14 '14 at 21:05
• @MichaelE2 Sorry I scrolled too fast and missed the name of the person, who proposed using If. I updated both posts to remove the name(s) and generalize the problem. – Konstantin Konstantinov Oct 14 '14 at 21:45
• I thought that might be it. :) (I'll delete my comments now.) – Michael E2 Oct 14 '14 at 21:47