Application Center - Maplesoft

App Preview:

Lesson 3: Application: Exponential and Logistic Growth

You can switch back to the summary page by clicking here.

Learn about Maple
Download Application


 

Lesson03.mw

ORDINARY DIFFERENTIAL EQUATIONS POWERTOOL

Lesson 3 -- Application: Exponential and Logistic Growth

Prof. Douglas B. Meade

Industrial Mathematics Institute

Department of Mathematics

University of South Carolina

Columbia, SC 29208

URL:   http://www.math.sc.edu/~meade/

E-mail: meade@math.sc.edu

Copyright  2001  by Douglas B. Meade

All rights reserved

-------------------------------------------------------------------

>

Outline for Lesson 3

3.A Exponential Growth and Decay

                3.A-1 General Solution

                3.A-2 Doubling Time

                3.A-3 Half-Life

3.B Logistic Equation

                3.B-1  General Solution

                3.B-2 Carrying Capacity

>

Initialization

> restart;

> with( DEtools ):

> with( plots ):

Warning, the name changecoords has been redefined

>

3.A Exponential Growth and Decay

3.A-1 General Solution

The general model for exponential growth and decay with rate constant k , is

> ode := diff( x(t), t ) = k * x(t);

ode := diff(x(t), t) = k*x(t)

>

and initial condition

> ic := x(0) = A;

ic := x(0) = A

>

If k is positive, the equation models exponential growth, whereas if k is negative, it models exponential decay.

This equation is separable, and Maple agrees with this classification, as seen with

> odeadvisor( ode, [separable] );

[_separable]

>

To find the general solution from first principles, divide through by x(t) to separate variables, and then integrate both sides of the resulting equation to obtain

> impl_soln := subs( _t=t, map( int, ode/x(t), t=0.._t,continuous ) );

impl_soln := -ln(x(0))+ln(x(t)) = k*t

>

Substitution of the initial condition leads to

> impl_part_soln := subs( ic, impl_soln );

impl_part_soln := -ln(A)+ln(x(t)) = k*t

>

and solving for x(t) gives

> expl_part_soln := simplify(op(solve( impl_part_soln, {x(t)} )));

expl_part_soln := x(t) = exp(k*t)*A

>

Of course, the same result could be obtained from

> infolevel[dsolve] := 3:

> soln2 := dsolve( {ode, ic}, x(t), [separable] );

> infolevel[dsolve] := 0:

`Classification methods on request`
`Methods to be used are: [separable]`

`----------------------------`

`* Tackling ODE using method: separable`

`--- Trying classification methods ---`

`trying separable`

`<- separable successful`

soln2 := x(t) = exp(k*t)*A

>

and the final simplification

> simplify( soln2 );

x(t) = exp(k*t)*A

>

3.A-2 Doubling Time

The doubling time for a process growing exponentially is the time t[d] needed for the quantity to double from its original size.  Thus, the doubling time is found from the relationship

x(t[d]) = 2*x(0) .

Implementing that for the solution found above leads to the equation

> double_eqn := subs( t=t[d], x(t[d])=2*x(0), ic, expl_part_soln );

double_eqn := 2*A = exp(k*t[d])*A

>

whose solution is

> double_time := solve( double_eqn, {t[d]} );

double_time := {t[d] = ln(2)/k}

>

One of the important characteristics of the doubling time is that not only is it the time needed for the initial size to double, it is the time needed for the size to double at any point in an exponential process.  The ratio of the solution at time T+t[d] and at time t[d] is given by

> q1 := x(T+t[d])/x(T)

>    = subs( t=T+t[d],

>            rhs(expl_part_soln) ) / subs( t=T, rhs(expl_part_soln) );

q1 := x(T+t[d])/x(T) = exp(k*(T+t[d]))/exp(k*T)

>

If t[d] is the true doubling time, then this ratio simplifies to 2, as is seen via

> q1;

> `` = simplify(eval( rhs(q1), double_time ));

x(T+t[d])/x(T) = exp(k*(T+t[d]))/exp(k*T)

`` = 2

>

3.A-3 Half-Life

The half-life for a quantity that is decaying according to an exponential model is the time after which exactly half the original amount remains.  The half life t[h] must obey the equation

> half_eqn := subs( t=t[h], x(t[h])=1/2*x(0), ic, expl_part_soln );

half_eqn := 1/2*A = exp(k*t[h])*A

>

Solving for t[h] leads to

> half_time := solve( half_eqn, {t[h]} );

half_time := {t[h] = -ln(2)/k}

>

Note that k < 0 holds for a decaying process. Thus, the half-life for an exponential model with "growth" rate -k is the same as the doubling time for an exponential model with growth rate k .

>

3.B Logistic Equation

3.B-1 General Solution

The general logistic equation is a modification of the exponential model in which the growth is tempered by the factor (K-x(t) ).  Therefore, the model consists of the ODE

> logistic_ode := diff( x(t), t ) = A * x(t) * ( K - x(t) );

logistic_ode := diff(x(t), t) = A*x(t)*(K-x(t))

>

with, for example, an initial condition of the form

> logistic_ic := x(0)=X[0];

logistic_ic := x(0) = X[0]

>

Initially, when x(t) is small, the factor K-x(t) is essentially constant, so the population obeys an exponential growth law.  As x(t) increases and approaches K , the factor K-x(t) tends towards zero, so the rate of change, dx/dt , approaches zero also.  Thus, x(t) approaches a constant, and the growth is said to be self-limiting.

This ODE is easily separated by the obvious division, resulting in

> sep_log_ode := logistic_ode / (x(t)*(K-x(t)));

sep_log_ode := (diff(x(t), t))/(x(t)*(K-x(t))) = A

>

Integrating from the initial time, 0, to any other time t , yields

> part_log_soln := subs( _t=t, logistic_ic, map(int,sep_log_ode,t=0.._t, continuous) );

part_log_soln := -(-ln(-K+X[0])+ln(X[0])+ln(-K+x(t))-ln(x(t)))/K = A*t

>

and solving for x(t) to obtain the explicit solution, gives

> part_expl_soln := collect(op( solve( part_log_soln, {x(t)} ) ),exp);

part_expl_soln := x(t) = X[0]*K*exp(A*t*K)/(K-X[0]+X[0]*exp(A*t*K))

>

A quick check that this is, in fact, a solution to the original ODE shows

> odetest( part_expl_soln, logistic_ode );

0

>

and, for the initial condition,

> eval( part_expl_soln, t=0 );

x(0) = X[0]

>

3.B-2 Carrying Capacity

Before discussing the generic properties of the logistic model, it is instructive to use graphical methods to examine a specific example.  For parameter values, choose

> param := { A = 1/2, K=3 };

param := {A = 1/2, K = 3}

>

so the logistic equation becomes

> l_ode := subs(param,logistic_ode);

l_ode := diff(x(t), t) = 1/2*x(t)*(3-x(t))

>

Then, the direction field for this model is seen in Figure 3.1.

> DEplot( l_ode, x(t), t=0..10, x=0..10, arrows=SMALL, title="Figure 3.1" );

[Plot]

>

An equilibrium solution of the differential equation F(x, `x'`) = 0 is any solution x(t) for which `x'`(t) = 0 , identically.  Hence, an equilibrium solution is a constant solution.  Such constant solutions show up on the direction field as horizontal lines, corresponding to the constant value of x(t) .

Figure 3.1 suggests that x = 3 might be an equilibrium solution.  Closer inspection of Figure 3.1 suggests that x = 0 might also be an equilibrium solution.  These potential equilibrium solutions can be found analytically by solving the equation `x'`(t) = f(x) = 0 for its roots.  The roots of the equation

> EQUILIB_EQN := rhs(l_ode) = 0;

EQUILIB_EQN := 1/2*x(t)*(3-x(t)) = 0

>

are found to be

> solve(EQUILIB_EQN, x(t));

0, 3

>

Initial conditions that will produce the equilibrium solutions are

> equil_ic := [ [x(0)=0], [x(0)=3] ];

equil_ic := [[x(0) = 0], [x(0) = 3]]

>

Use of these initial conditions in the DEplot command yields Figure 3.2.

> equil_plot := DEplot( l_ode, x(t), t=0..5, equil_ic, x=0..10,

>                      arrows=SMALL, linecolor=CYAN, title="Figure 3.2" ):

> equil_plot;

[Plot]

>

Figure 3.3 shows a sample of solutions with initial conditions between the two equilibria.

> ic1 := [ [x(0)=i/2] $ i=1..5 ]:

> soln_plot1 := DEplot( l_ode, x(t), t=0..5, ic1, x=0..10, arrows=SMALL,

>                      linecolor=GREEN, title="Figure 3.3" ):

> soln_plot1;

[Plot]

>

Note that all of these solutions are increasing and appear to approach x = 3 as t continues to increase.

Figure 3.4 shows a sample of solutions with initial conditions above the positive equilibrium.

> ic2 := [ [x(0)=2*i] $ i=2..5 ]:

> soln_plot2 := DEplot( l_ode, x(t), t=0..5, ic2, x=0..10, arrows=SMALL,

>                      linecolor=BLUE, stepsize=0.1, title="Figure 3.4" ):

> soln_plot2;

[Plot]

>

These solutions are all decreasing and also appear to approach the x = 3 as t increases.

The composite plot is appears in Figure 3.5.

> display( [equil_plot, soln_plot1, soln_plot2], title="Figure 3.5" );

[Plot]

>

To investigate some of the general properties of solutions to the logistic equation, recall that the equilibrium solutions are

> equil_sol := solve( rhs(logistic_ode)=0, {x(t)} );

equil_sol := {x(t) = 0}, {x(t) = K}

>

A quick inspection of the ODE shows that diff(x(t), t) > 0 when 0 < x(t) < K and diff(x(t), t) < 0 when x(t) > K . The long-term behavior of the solutions can be determined from the explicit solution to the IVP by evaluating the following limit.

> limit_size := Limit( part_expl_soln, t=infinity );

limit_size := Limit(x(t) = X[0]*K*exp(A*t*K)/(K-X[0]+X[0]*exp(A*t*K)), t = infinity)

>

The limit exists as a finite number only if both A and K are nonnegative.  With these physically reasonable assumptions, the limit is

> value( limit_size ) assuming A>=0, K>=0;

limit(x(t), t = infinity) = K

>

Because all solutions with x(0)*`>`*K decrease to x = K and all solutions with 0 < x(0) < K increase to x = K , the parameter K is known as the carrying capacity of the logistic equation.

The equilibrium solution x = K is called a stable equilibrium because, as t increases, nearby solutions, those close and on either side of it, tend towards the equilibrium.  This is shown very well in Figure 3.5.

The equilibrium solution x = 0 is called an unstable equilibrium because, as t increases, nearby solutions, those close and on either side of it, move away from the equilibrium.  Figure 3.6, exhibiting solutions starting both above and below the equilibrium solution x = 0 , shows that this equilibrium is unstable.

> DEplot( l_ode, x(t), t=0..1, [[x(0)=1],[x(0)=0],[x(0)=-1]], x=-2..2, arrows=SMALL,

>                      linecolor=[green,red,green], stepsize=0.1, title="Figure 3.6" );

[Plot]

>

An equilibrium solution can also be called semi-stable if, as t increases, nearby solutions on one side approach the equilibrium, but nearby solutions on the other move away from the equilibrium.  Adifferential equation with a semi-stable equilibrium would be the equation

> SS_ode := diff(x(t),t) = x(t)*(1-x(t))^2;

SS_ode := diff(x(t), t) = x(t)*(1-x(t))^2

>

The equilibrium solutions are x = 0 and x = 1 , with the latter being the semi-stable equilibrium, as shown in Figure 3.7.

> DEplot(SS_ode, x(t),t=0..1, [[x(0)=1],[x(0)=.5],[x(0)=1.5]], x=0..2, arrows=medium, dirgrid=[10,10], linecolor=[red,green,green], stepsize=0.1, title="Figure 3.7" );

[Plot]

>

Solutions starting above x = 1 move away from this equilibrium, but solutions starting below x = 1 move towards it.  Solutions on one side of the equilibrium are attracted to the equilibrium, but solutions on the other side are repelled.

>

[Back to ODE Powertool Table of Contents]

>

>