How can I "stop" the graph which "hits the zero"?

Question

I have a simple model for the population:

a := 5; 
b := 1;
h := 25/4;, 
de := diff(P(t), t) = P(t)*(a-b*P(t))-h;
cond1 := P(0) = .5;
cond2 := P(0) = 2;
cond3 := P(0) = 5;
sol1 := dsolve({cond1, de}, P(t)); 
sol2 := dsolve({cond2, de}, P(t)); 
sol3 := dsolve({cond3, de}, P(t));
P1 := plot(rhs(sol1), t = 0 .. 5, color = blue);
P2 := plot(rhs(sol2), t = 0 .. 5, color = green);
P3 := plot(rhs(sol3), t = 0 .. 5, color = orange);
plots:-display(P1, P2, P3);

When the graph "hits the zero", population goes extinct, so I need to "stop" the graph there. I can use Optimization[Minimize](...); to minimize deviation from zero (and include this point, say, A, in t=0..A), but it's not always possible and kind of "brute-force". How can I "stop" the graph which "hits the zero"? Also, how can I find the point where non-extinct population becomes stable (say, 1% deviation from limit value)?

Answer 1

You can solve the differential equation, for an "exact" symbolic formula, without having to specify a particular numeric value for the initial condition.

restart;
a := 5:
b := 1:
h := 25/4:

de := diff(P(t), t) = P(t)*(a-b*P(t))-h:

solparexact := dsolve({P(0) = P0, de}, P(t));

                            10 P0 t + 4 P0 - 25 t
      solparexact := P(t) = ---------------------
                            2 (2 P0 t - 5 t + 2) 

Let's pick off the right hand side, for convenience. This is P(t).

solparexact := rhs(solparexact);

                         10 P0 t + 4 P0 - 25 t
          solparexact := ---------------------
                         2 (2 P0 t - 5 t + 2) 

We can solve that for a formula expressing where P(t) would become zero.

general_root := solve({t>0, solparexact}, t) assuming P0>0;

     general_root :=  / [ /         4 P0   \ ]           5 
                      | [{ t = - ---------- }]      P0 < - 
                      | [ \      10 P0 - 25/ ]           2 
                     <                                     
                      |                             5      
                      |           []                - <= P0
                      \                             2      

That prints ugly in character-mode, and nicer in the Maple GUI. It is a piecewise, showing that there is no root when P0 > 5/2. That is, when the initial population is greater than 5/2 the population doesn't ever become zero.

We can evaluate that result at various values of P0. The following does that.

It is known as two-argument eval, to evaluate an expression for specific values (substituions) of one of more of its unknowns. This is a very useful bit of Maple syntax, which I'll utilize several times in all the code below.

eval(general_root, P0=2);

                      [ /    8\ ]
                      [{ t = - }]
                      [ \    5/ ]

eval(general_root, P0=0.5);

                  [{t = 0.1000000000}]

eval(general_root, P0=5); # no root

                           []

The closer P0 gets to 5/2, from above, the longer it takes for the population to ever get to zero.

eval(general_root, P0=5/2 - 1e-9);

               [ /                  8\ ]
               [{ t = 9.999999996 10  }]
               [ \                   / ]

We could even extract the formula (for the time t where the root occurs and the population becomes zero) within the piecewise, valid for P0>0 and P0<5/2,

genroot := eval(t, general_root[1]) assuming P0>0, P0<5/2;

                                4 P0   
                genroot := - ----------
                             10 P0 - 25

We can find the limiting value for the population, as t goes to infinity.

L := limit(solparexact, t=infinity);

                              5
                         L := -
                              2

We can re-express the exact formula for P(t) at particular values of P0.

S[3.0] := eval(solparexact, P0 = 3.0);

                          5.0 t + 12.0 
                S[3.0] := -------------
                          2 (1.0 t + 2)

S[5.0] := eval(solparexact, par = 5.0);

                      10 P0 t + 4 P0 - 25 t
            S[5.0] := ---------------------
                      2 (2 P0 t - 5 t + 2) 

We can find the time at which the population gets close (relatively), to its limiting value. We can do that exactly, or in floating-point,

End[3.0] := fsolve( (S[3.0] - L)/L = 0.1, t = 0..infinity );

                End[3.0] := 2.000000000

solve( (eval(solparexact, P0 = 3) - L)/L = 1/10);

                           2

We can even find a general representation of that,

Endform := solve( {P0>0, t>0, ((solparexact - L)/L) = 1/10},
                  t) assuming P0>5/2;

    Endform :=  /                                11
                |         []               P0 <= --
                |                                4 
               <                                   
                | [ /    8 P0 - 22\ ]      11      
                | [{ t = --------- }]      -- < P0 
                \ [ \    2 P0 - 5 / ]       4      

That too can be evaluated at particular values of the initial population size. These agree with what was found above.

eval(Endform, P0=3), eval(Endform, P0=3.0);

             [{t = 2}], [{t = 2.000000000}]

We could pick of the time value, alone,

eval(t, op(eval(Endform, P0=3.0)));

                      2.000000000

Or we could pick off the formula, under the appropriate assumption.

genend := eval(t,Endform[1]) assuming P0>11/4;

                            8 P0 - 22
                  genend := ---------
                            2 P0 - 5 

And now for a plot, for a list of P0 values that attain zero population, and a list of values which do not.

Zerovals := [0.5, 1.5, 2.0, 2.1, 2.25, 2.3]:
P0vals := [3.0, 3.5, 5.0, 10.0]:

plots:-display(
  plots:-shadebetween(L, 1.1*L, t = 0 .. 4, linestyle=dot,
                      color=gray, transparency=0.5),
  seq( plot(eval(solparexact,P0=P0vals[i]),
            t = 0 .. eval(genend,P0=P0vals[i]),
            color=cat("Spring ",i), legend=P0vals[i]),
       i = nops(P0vals)..1, -1 ),
  plot(L, t = 0 .. 4, linestyle=dot, thickness = 3,
       legend=evalf[2](L)),
  seq( plot(eval(solparexact,P0=Zerovals[i]),
            t = 0 .. eval(genroot,P0=Zerovals[i]),
            color=ColorTools:-Color([1,1,1]-
                                    [i,i,i]/(nops(Zerovals)+1)),
            legend=Zerovals[i]),
       i = nops(Zerovals)..1, -1 ),
  axis[2]=[tickmarks=[op(Zerovals),L,op(P0vals)]],
  legendstyle=[location = right],
  size = [700, 500],
  view = [0 .. 4, 0 .. max(P0vals)]
);