Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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"}]

fig 1

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". enter image description here

share|improve this question
    
You could solve the equation, then define aPhys[t_]=Max[0, a[t]]. –  b.gatessucks Jan 24 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 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 at 12:20
add comment

3 Answers

up vote 3 down vote accepted

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["DifferentialEquations`NDSolveProblems`"];
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:

enter image description here

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.

share|improve this answer
add comment

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"}]

enter image description here

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

BR

share|improve this answer
    
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 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 at 14:21
add comment

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]

Mathematica graphics

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]

Mathematica graphics

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]
share|improve this answer
    
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 at 14:19
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.