Constraining function found by NDSolve to stay positive

Question

I am trying to constrain my variables in NDSolve so that they stay positive for the integration interval. Is is possible to do that? Using Assumptions as in Simplify does not work. WhenEvent also does not help, as it just finds the event once.

Here is a code snippet of a minimal example:

a = 10;
sol = NDSolve[{x'[t] == -10, x[0] == 10}, x, {t, 0, 10}];
Plot[Evaluate[x[t] /. sol], {t, 0, 10}, PlotRange -> All]

What I would like to get out is a function, that linearly decreases from 10 to 0 and stays on 0 -- instead of turning negative.

Edit

To clarify the question: I am using Mathematica for fast prototyping. In a C environment I can modify the time integration with an additional step, like if(x<0){ x=1e-10; } to ensure my variable does not turn negative. So I am looking for a way to introduce a lower bound into my system.

My full system consists of a system of 1st order nonlinear ODEs, which represent positive real quantities. I am aware that this modification will change the resulting solutions. The above just provides a minimal example of what I am looking for.

Changing the rhs as suggested below will unfortunately not do the trick, as it only adjusts the rate of change but does not introduce a lower bound for the solution. (Although, strictly speaking this procedure will result in a solution I was describing above ;)

Answer 1

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.

---

@DanielLichtblau Ok. Here are 4 tests based on the same system of equations. The code for equations (ndSolveLst) was actually generated, so I had to copy/paste it a few times and do some replacements to make it more readable but I doubt that that help much. My apologies.

The equations are a mix of chemical reactions (which are proportional to some integer powers of concentrations - in this particular example the possible powers are 1 and 2, if I am not mistaken) and crystallization/dissolution processes. It is the dissolution process, which is causing the trouble. Basically the dissolution process (at small amount of solid sediment) has a constant rate until all the sediment dissolves and then it must stop, of course! That's exactly the example above.

The test solves the same system of equations using 4 types of rules:

1. ruleBase - "discontinuity avoidance rules". That's the one I described above but with slight edits to simplify the tests.
2. ruleIf - Using If to enforce that all the terms, which must be identically zero if some x[t]<0, stay at zero.
3. rulePiecewise - Using Piecewise to enforce that all the terms, which must be identically zero if some x[t]<0, stay at zero.
4. ruleIfNone - Using If to enforce that discontinuous terms in the derivative function, which must be identically zero if some x[t]<0, stay at zero, while letting the values of x[t] in all other places take their values as they are.

tMaxValue is the maximum value of time for which to solve the equations. I put the value of tMaxValue = 10^8 because this is approximately a relaxation time of this system. However, all test except the first one never came back on my machine for that value of tMaxValue. The maximum time up to which they can work is slightly above 3*10^3. For that value the 4-th test produces meaningless results as some of the variables become substantially negative.

My point is that for complicated system it is better to avoid discontinuities rather than to handle them. If you have too many discontinuities (which, for example, this "simple" system exhibits) then NDSolve just dies trying to handle them. In the particular case of this system the physical essence of the model allows introducing smooth cut offs without changing the essence of the model. In some other cases that may not work, of course.

ClearAll["Global`*"];

tMaxValue = 10^8;

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

 NnnXIf[x_] := If[x > 0, x, 0];
 NnnXPiecewise[x_] := Piecewise[{{x, x > 0}}, 0];
 NnnXNone[x_] := x;

UstEIf[x_, norm_] := If[x > 0, 1, 0];
UstEPiecewise[x_, norm_] := Piecewise[{{1, x > 0}}, 0];

ruleBase = {NnnX -> NnnXBase, UstE -> UstEBase};
ruleIf = {NnnX -> NnnXIf, UstE -> UstEIf};
rulePiecewise = {NnnX -> NnnXPiecewise, UstE -> UstEPiecewise};
ruleIfNone = {NnnX -> NnnXNone, UstE -> UstEIf};

RunNDSolve[rules_, tMax_?NumericQ] := 
  Module[{ndSolveLst, sol, plotOptions, time},
   ndSolveLst = {{Derivative[1][rZ][t] == 0., 
       Derivative[1][rY][t] == 
        1.0009999999999999*^-7*NnnX[ra[t]] + 9.99*^-8*NnnX[rA[t]] - 
         2.0000000000000003*^-6*NnnX[rY[t]] + 
         5.776792245814018*^-6*NnnX[r$$$a$a[t]] + 
             5.519929566477231*^-6*NnnX[r$$$A$a[t]] + 
         5.776792245814018*^-6*NnnX[r$$$A$A[t]], 
           rX[t] == 
            2.*NnnX[ra$a[t]] + 2.*NnnX[rA$a[t]] + 2.*NnnX[rA$A[t]], 
           Derivative[1][rA][t] == -9.99*^-8*NnnX[rA[t]] - 
             0.2371801575241698*NnnX[ra[t]]*NnnX[rA[t]] - 
             0.13823802227503887*NnnX[rA[t]]^2 + 
             0.022937500731028893*NnnX[rA$a[t]] + 
             0.0975399548966282*NnnX[rA$A[t]] + 
             1.0000000000000002*^-6*NnnX[rY[t]], 
           Derivative[1][ra][t] == -1.0009999999999999*^-7*NnnX[ra[t]] - 
             0.13823802227503887*NnnX[ra[t]]^2 - 
             0.2371801575241698*NnnX[ra[t]]*NnnX[rA[t]] + 
             0.0975399548966282*NnnX[ra$a[t]] + 
             0.022937500731028893*NnnX[rA$a[t]] + 
             1.0000000000000002*^-6*NnnX[rY[t]], 
           Derivative[1][rA$A][t] == 
            0.06911901113751943*NnnX[rA[t]]^2 - 
             0.0487699774483141*NnnX[rA$A[t]] + 
             0.0002528839419740976*(0.13441061034286977 + 
                Sqrt[(0.13441061034286977 - 0.5*rX[t])^2] - 
                0.5*rX[t])*(r$$$A$A[t] + Sqrt[r$$$A$A[t]^2]) - 
             0.00498252784982204*(-0.13441061034286977 + 
                Sqrt[(-0.13441061034286977 + 0.5*rX[t])^2] + 0.5*rX[t])*
              UstE[rA$A[t], 0.13441061034286977], 
           Derivative[1][r$$$A$A][
         t] == -2.888396122907009*^-6*NnnX[r$$$A$A[t]] - 
             0.0002528839419740976*(0.13441061034286977 + 
                Sqrt[(0.13441061034286977 - 0.5*rX[t])^2] - 
                0.5*rX[t])*(r$$$A$A[t] + Sqrt[r$$$A$A[t]^2]) + 
             0.00498252784982204*(-0.13441061034286977 + 
                Sqrt[(-0.13441061034286977 + 0.5*rX[t])^2] + 0.5*rX[t])*
              UstE[rA$A[t], 0.13441061034286977], 
           Derivative[1][ra$a][t] == 
            0.06911901113751943*NnnX[ra[t]]^2 - 
             0.0487699774483141*NnnX[ra$a[t]] + 
             0.0002528839419740976*(0.13441061034286977 + 
                Sqrt[(0.13441061034286977 - 0.5*rX[t])^2] - 
                0.5*rX[t])*(r$$$a$a[t] + Sqrt[r$$$a$a[t]^2]) - 
             0.00498252784982204*(-0.13441061034286977 + 
                Sqrt[(-0.13441061034286977 + 0.5*rX[t])^2] + 0.5*rX[t])*
              UstE[ra$a[t], 0.13441061034286977], 
           Derivative[1][r$$$a$a][
         t] == -2.888396122907009*^-6*NnnX[r$$$a$a[t]] - 
             0.0002528839419740976*(0.13441061034286977 + 
                Sqrt[(0.13441061034286977 - 0.5*rX[t])^2] - 
                0.5*rX[t])*(r$$$a$a[t] + Sqrt[r$$$a$a[t]^2]) + 
             0.00498252784982204*(-0.13441061034286977 + 
                Sqrt[(-0.13441061034286977 + 0.5*rX[t])^2] + 0.5*rX[t])*
              UstE[ra$a[t], 0.13441061034286977], 
           Derivative[1][rA$a][t] == 
            0.2371801575241698*NnnX[ra[t]]*NnnX[rA[t]] - 
             0.022937500731028893*NnnX[rA$a[t]] + 
             0.0002492109966820526*(0.08882322843930857 + 
                Sqrt[(0.08882322843930857 - 0.5*rX[t])^2] - 
                0.5*rX[t])*(r$$$A$a[t] + Sqrt[r$$$A$a[t]^2]) - 
             0.005023486011561131*(-0.08882322843930857 + 
                Sqrt[(-0.08882322843930857 + 0.5*rX[t])^2] + 0.5*rX[t])*
              UstE[rA$a[t], 0.08882322843930857], 
           Derivative[1][r$$$A$a][
         t] == -2.7599647832386156*^-6*NnnX[r$$$A$a[t]] - 
             0.0002492109966820526*(0.08882322843930857 + 
                Sqrt[(0.08882322843930857 - 0.5*rX[t])^2] - 
                0.5*rX[t])*(r$$$A$a[t] + Sqrt[r$$$A$a[t]^2]) + 
             0.005023486011561131*(-0.08882322843930857 + 
                Sqrt[(-0.08882322843930857 + 0.5*rX[t])^2] + 0.5*rX[t])*
              UstE[rA$a[t], 0.08882322843930857], rZ[0.] == 0., 
           rY[0.] == 0., rA[0.] == 1., ra[0.] == 0., rA$A[0.] == 0., 
           r$$$A$A[0.] == 0., ra$a[0.] == 0., r$$$a$a[0.] == 0., 
           rA$a[0.] == 0., r$$$A$a[0.] == 0.}, {rZ, rY, rX, rA, ra, rA$A, 
           r$$$A$A, ra$a, r$$$a$a, rA$a, r$$$A$a}, {t, 0, tMax}, 
      MaxSteps -> Infinity, 
      Method -> {"EquationSimplification" -> "Residual"}} /. rules;

   Print["rules = ", rules];
   Print["ndSolveLst = ", ndSolveLst];

   time = Timing[sol = Apply[NDSolve, ndSolveLst]][[1]];
   Print["NDSolve run time = ", time];

   plotOptions = {PlotRange -> All, Frame -> True, 
     GridLines -> Automatic, PlotStyle -> Thick, ImageSize -> 500};

   Print["Substances: A, a, and (A+a)"];
   Print[Plot[{Evaluate[rA[t] /. sol], Evaluate[ra[t] /. sol], 
      Evaluate[rA[t] /. sol] + Evaluate[ra[t] /. sol]}, {t, 0, tMax}, 
     Evaluate[plotOptions]]];

   Print["Pairs: A$A, A$a, a$a, and (A$A+A$a+a$a)"];
   Print[Plot[{Evaluate[rA$A[t] /. sol], Evaluate[rA$a[t] /. sol], 
      Evaluate[
       ra$a[t] /. sol], (Evaluate[rA$A[t] /. sol] + 
        Evaluate[rA$a[t] /. sol] + Evaluate[ra$a[t] /. sol])}, {t, 0, 
      tMax}, Evaluate[plotOptions]]];

   Print["Solids: $$$A$A, $$$A$a, $$$a$a, and ($$$A$A+$$$A$a+$$$a$a)"];
   Print[Plot[{Evaluate[r$$$A$A[t] /. sol], 
          Evaluate[r$$$A$a[t] /. sol], 
      Evaluate[
       r$$$a$a[t] /. sol], (Evaluate[r$$$A$A[t] /. sol] + 
        Evaluate[r$$$A$a[t] /. sol] + 
            Evaluate[r$$$a$a[t] /. sol])}, {t, 0, tMax}, 
     Evaluate[plotOptions]]];
   Print["=================="];
   ];

RunNDSolve[ruleBase, tMaxValue];
RunNDSolve[ruleIf, tMaxValue];
RunNDSolve[rulePiecewise, tMaxValue];
RunNDSolve[ruleIfNone, tMaxValue];

Answer 2

Here are couple of ways to do what the OP want. As was pointed out by @belisarius, your differential equation is inconsistent with x[t] being constant. I can think of two ways to deal with this, either stop integration or change the differential equation when the condition is reached.

Method 1

Here the derivative parameter is a[t] and changes when the event x[t] == 0 is reached.

Clear[a];
sol = NDSolve[{x'[t] == -a[t], x[0] == 10, a'[t] == 0, a[0] == 10, 
    WhenEvent[x[t] == 0, a[t] -> 0]}, x, {t, 0, 10}];

Plot[Evaluate[x[t] /. sol], {t, 0, 10}, PlotRange -> All]

Method 2

In this method, we stop integration at the event and use Piecewise to complete the solution. This method may not work as is in the OP's use case. But instead of the default value of 0 used by Piecewise, another function could be supplied, such as the solution to another differential equation.

solval = Piecewise[{{#[t], #[[1, 1]] <= t <= #[[1, 2]] &@#["Domain"]}}] &@
  NDSolveValue[{x'[t] == -10, x[0] == 10, 
    WhenEvent[x[t] < 0, "StopIntegration"]}, x, {t, 0, 10}];

Plot[solval, {t, 0, 10}, PlotRange -> All]

Update: Method 2'

Similar to Method 2, but using "ExtrapolationHandler" instead of Piecewise:

solval = NDSolveValue[{x'[t] == -10, x[0] == 10, 
   WhenEvent[x[t] < 0, "StopIntegration"]}, x, {t, 0, 10}, 
  "ExtrapolationHandler" -> {0 &, "WarningMessage" -> False}]

In all cases the graphics output is the same:

Answer 3

Of course belisarius is right with his comment that the ODE in your example doesn't have a positive solution. I guess that you want to actually solve a different equation, and this should do what I think you try to achieve:

a = 10;
sol = NDSolve[{x'[t] == If[x[t] >= 0, -10, 0], x[0] == 10}, x, {t, 0, 10}];
Plot[Evaluate[x[t] /. sol], {t, 0, 10}, PlotRange -> All]

I also don't think that what you say about WhenEvent is true: it very well can find events more than once. I also guess that your example is most probably simplified too much and for a more complicated example WhenEvent certainly would be the way to go. What is the code you did try with WhenEvent?