ORDINARY DIFFERENTIAL EQUATIONS POWERTOOL

Unit 25 -- Application: Population Models

Outline of Unit 25

Initialization

> restart;

> with( DEtools ):

> with( plots ):

> with( linalg ):

Warning, the name changecoords has been redefined

Warning, the name adjoint has been redefined

Warning, the protected names norm and trace have been redefined and unprotected

>

A population with constant birth and death rates can be modeled as

> ode1 := diff( p(t), t ) = k * p(t);

> ic := p(0) = P[0];

>

where denotes the size of the population at time , is the net growth rate, and is the initial size of the population.

The solution to any Malthusian model is an exponential function

> sol1 := dsolve( { ode1, ic }, p(t) );

>

If the net growth rate is positive, i.e. , > , then the population grows (exponentially) without bound. If the net growth rate is negative, then the population decays exponentially to zero. And, if the net growth rate is zero, the population remains constant for all time: .

While an exponential model can be useful for short-term projections, the unbounded growth is a major limitation to making long-term projections.

>

The logistic model addresses the unbounded growth of the Malthusian model by assuming that the growth rate starts at the net growth rate but decreases (linearly) as the population increases.

> ode2 := diff( p(t), t ) = k * ( 1 - p(t)/A ) * p(t);

> ic;

>

This autonomous ODE is separable; the solution is

> q1 := dsolve( { ode2, ic }, p(t) );

>

A more typical form for this solution is

> sol2 := subsop(2=map( collect, rhs(normal(q1)), exp ), q1);

>

There are two equilibria for the logistic model: and . The qualitative behavior can be determined from a sign chart, or from the direction field. Unfortunately, Maple is unable to produce a direction field for general parameters, but the essential features of the model are observable with any reasonable, i.e. , positive, values of and .

> param2 := { k=1, A=7 }:

> DEplot( eval( ode2, param2 ), p(t), t=0..10, p=0..10 );

>

The direction field is suggests that, when > 0, any solution that begins with a positive initial population ultimately tends to the value of . This shows how the logistic model overcomes the unbounded growth for exponential models. The parameter is frequently called the carrying capacity of the population. This conjecture is supported with the symbolic computation

> assume( k > 0 ):

> limit( sol2, t=infinity );

> unassign( 'k' );

>

To see what happens when k < 0, consider

> param2 := { k=-1, A=7 }:

> DEplot( eval( ode2, param2 ), p(t), t=0..10, p=0..10 );

>

Thus, when < 0, the carrying capacity is a bifurcation value for solutions. If < , solutions decay to zero and if > , solutions grow without bound. However, the symbolic computation of the limiting population yields

> assume( k < 0 ):

> limit( sol2, t=infinity );

> unassign( 'k' ):

>

The problem is that when the net growth rate is negative, the solution does not exist for all time. The solution exists only so long as the denominator is not zero.

> q2 := solve( denom(rhs(sol2))=0, t );

>

This expression is not easily simplified within Maple, but is easily seen to be equivalent to

> Tblowup := 1/k * ln( 1 - A/P[0] );

>

Since and are both positive, 0 < < 1 when > . Then, < 0 and the blow-up time is positive (and finite) when < 0. This finite time blowup is difficult to see in the direction field.

If < , the denominator is never zero and the solution exists for all time.

>

Of course, if = 0, then the population remains unchanged, .

>

Another interesting feature of the logistic equation is the change of concavity for some solutions with 0 < < . To investigate the location of this inflection point, zeroes of the second derivative of the solution must be identified:

> DDsol2 := diff( sol2, t$2 );

>

> q3 := solve( rhs(DDsol2)=0, t );

>

which can be simplified by hand to

> Tinflection := 1/k * ln( A/P[0] - 1 );

>

This time is the location of a possible inflection point. The complete verification that an inflection point does occur at this time is left as an exercise. (Note that the assumptions on the parameters ensure that this time is real and positive.) At this time, the population is

> simplify( eval( sol2, t=Tinflection ) );

>

That is, when 0 < < the solution changes concavity when the population is exactly half of the carrying capacity. The significance of this result is that it provides a means to estimate the carrying capacity for a population from data before the population actually reaches the carrying capacity. The only requirement is to be able to locate the inflection point from the data.

>

Another attempt to eliminate the unbounded growth of solutions to the Malthusian model is the Gompertz model.

> ode3 := diff( p(t), t ) = k * ( 1 - ln(p(t))/a ) * p(t);

> ic;

The main difference is the appearance of in the factor that modifies the net growth rate. With this change, and assuming > 0, the effective growth rate is when , then decreases (logarithmically) as increases. (For 0 < < 1, the effective growth rate exceeds the net growth rate.)

The solution to this IVP is

> q4 := simplify( dsolve( { ode3, ic }, p(t) ) );

>

which simplifies (manually) to

> sol3 := p(t) = exp( a ) * exp( (ln(P[0])-a)*exp(-k*t/a) );

>

To begin the qualitative analysis of solutions to the Gompertz model, observe that when and are nonzero and have the same sign

> assume( k/a > 0 ):

> limit( sol3, t=infinity );

> unassign( 'k', 'a' ):

>

A representative direction field for this situation is

> param3 := { k=1, a=ln(4.) }:

> DEplot( eval( ode3, param3 ), p(t), t=1..10, p=1..10 );

>

When and have different signs,

> assume( k/a < 0 ):

> limit( sol3, t=infinity );

> unassign( 'k', 'a' ):

>

The reason for Maple not simplifying this limit is that the value depends on the sign of . A representative direction field for this situation is

> param3 := { k=-1, a=ln(4.) }:

> DEplot( eval( ode3, param3 ), p(t), t=1..10, p=1..10 );

>

Note that solutions with < decay to zero while > gives rise to solutions that are unbounded as time increases without bound. A key difference from the logistic equation is that these unbounded solutions exist for all time.

The direction fields suggest that some solutions have a change in concavity. The second derivative of the solution is

> DDsol3 := diff( sol3, t$2 );

>

which is zero when

> Tinflect3 := solve( rhs(DDsol3)=0, t );

>

After it has been confirmed that this time actually is an inflection point for the solution, it is interesting to note that the inflection point occurs when the population is

> simplify( eval( sol3, t=Tinflect3 ) );

>

Thus, the limiting population for the Gompertz model occurs when the population is a factor larger than the population at the inflection point. (Recall that the logistic model's carrying capacity is twice the population at the change in concavity.)

>