Krdytkyu

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)?

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)]

);