# Differential equation of a physical property which can not be below zero

A simple question but maybe also interesting for others. I have the simple differential equation:

a'[t] == -r


The quantity a has a physical meaning and can not be below 0. The ODE is of course simple to solve but I don't know how to tell DSolve to restet a to 0 if it is below zero.

My approach now is to simple manipulate the result:

sol = DSolve[{a'[t] == -r, a[0] == 1}, a[t], t]
solspecial = sol /. r -> 0.1
Table[{t, Max[(a[t] /. solspecial), 0]}, {t, 0, 15}]
ListLinePlot[%, PlotRange -> All, Frame -> True,
FrameLabel -> {"time", "a"}]


Is there a better way to do this and already state that a can not be below 0 in DSolve or NDSolve?

edit: to clarify only negative values are reset to 0 and if r=Sin[t] the value should "come back".

-
You could solve the equation, then define aPhys[t_]=Max[0, a[t]]. – b.gatessucks Jan 24 '14 at 11:52
@b.gatessucks , yes thanks, but thats basically the same as in my code example only in a nicer formatted way. But can I already demand that a can not be below 0 inside Dsolve? – gogoolplex Jan 24 '14 at 11:57
Isn't your ODE incompatible with that boundary condition? a'[t]=-r cannot be true for all t. Shouldn't this be a'[t]=-r UnitStep[a[t]] or equivalent? – Ymareth Jan 24 '14 at 12:20

The answer to this question depends first of all on the physical mechanism that keeps $a$ from becoming negative. Your equation does not account for this mechanism, otherwise your question would not arise. In principle the choice of the way of how to make the cut-off must take into account this mechanism. If the mechanism is of no importance, then it is also not important which way to choose among those already proposed. Then it is the matter of convenience only.

If, on the contrary, it counts, then you need to think first of all what happens at $a\rightarrow 0$ and derive the method from this.

Just to complete the methods already proposed I would define $a(t)$ as a non-linear function of a new variable, say, $f(t)$, such that this function cannot become negative. I would then solve for $f(t)$ instead of $a(t)$. There are lots of such functions, the choice evidently depends upon the mechanism. Thus, we return to the first question.

To give and example let us take $a(t)=e^{f(t)}$. Equation then has the form $f'(t)=r \left(-e^{-f(t)}\right)$ and the initial condition $f(0)=0$. Note that this approach works good in application to objects of statistical physics, such as e.g. concentration.

    Needs["DifferentialEquationsNDSolveProblems"];
sl = NDSolve[{f'[t] == -0.1*Exp[-f[t]], f[0] == 0}, f, {t, 0, 16},
Method -> "StiffnessSwitching"];
pl = Plot[Evaluate[Exp[f[t]] /. sl], {t, 0, 16},
Epilog ->
Inset[Plot[Evaluate[Exp[f[t]] /. sl], {t, 10, 16},
PlotRange -> {0, 0.001}], Scaled[{0.7, 0.7}]]]


The system exhibits a stiffness at a close to zero, which is expected, of course, since the mechanism to avoid the negative value is forced artificially.

This returns the expected result:

The plot shows the $a(t)=e^{f(t)}$ behavior. The inset in the image shows more closely, what happens in the vicinity of zero.

-

If you don't mind using WhenEvent combined with NDSolve you can try the following. The parameter crossing gives you the control where to set that limitting value. In your case it should be zero.

crossing = 0.;
indicator[x_?NumberQ] := If[x > crossing, 1, 0];
sol = With[{r = 0.1},
NDSolve[{a'[t] == -r indicator[a[t]], a[0] == 1,
WhenEvent[a[t] == crossing, "CrossSlidingDiscontinuity"]},a, {t, 0, 20}]
];
Plot[a[t] /. sol, {t, 0, 20}, PlotRange -> All, Frame -> True,FrameLabel -> {"time", "a"}]


Unfortunately there is still no built in event trigger for the mighty DSolve.

BR

-
Ok, interesting, but quite complicated, I was hoping for an more simple solution and with using DSolve. Also if this method is used with say r=Sin[t] then, the a value sticks to 0 and does not come back. So maybe my approach with manipulating the DSolve result is not so bad. – gogoolplex Jan 24 '14 at 14:06
@gogoolplex Please check the WhenEvent doc properly. I am sure you can make it "come back" even with a r=Sin[t]. There are such coming back examples in the doc! – PlatoManiac Jan 24 '14 at 14:21

Also;

r = .1;
sol = NDSolve[{
a'[t] == -r u[t],
a[0]  == 1,
u[0]  == 1,
WhenEvent[a[t] < 0, u[t] -> 0]},
{a, u}, {t, 0, 15}, DiscreteVariables -> {u}]
Plot[a[t] /. sol, {t, 0, 15}, Evaluated -> True, PlotRange -> All, Frame -> True]


Edit:

r[t_] := Sin[t];
sol = NDSolve[
{a'[t] == -r[t] u[t], a[0] == 1, u[0] == 1,
WhenEvent[r[t] < 0, u[t] -> 0],
WhenEvent[r[t] > 0, u[t] -> 1]}, {a, u}, {t, 0, 15},
DiscreteVariables -> {u}]
Plot[a[t] /. sol, {t, 0, 15}, Evaluated -> True, PlotRange -> All, Frame -> True]


Or equivalently:

r[t_] := Max[0, Sin[t]];
sol = NDSolve[{a'[t] == -r[t] , a[0] == 1}, {a}, {t, 0, 15}]
Plot[a[t] /. sol, {t, 0, 15}, Evaluated -> True, PlotRange -> All,  Frame -> True]

-
I prefer this over PlatoManiacs solution, but also here a stays at 0 and does not "come back" when r=Sin[t], please see may edit in the Question. I guess I was not clear about that. – gogoolplex Jan 24 '14 at 14:19