Population Dynamics

Mathematical Modeling

Richard Hitt
Department of Mathematics and Statistics
University of South Alabama
Mobile, AL 36688

Fibonacci Growth

Perhaps one of the first mathematical models of a population growth was due to Leonardo Pisano (Fibonacci) in 1202.  In this setting, we assume we have a newborn pair of rabbits who will eventually mate.  We assume that rabbits are born in oppositely sexed pairs which become sexually mature at age 1 month and that the gestational term is 1 month.  

Let F[n] denote the number of pairs of rabbits at the end of each month.  Lets look at different ways we can compute the number of rabbits at the end of the nth month.

Simple but Inneficient Mathematica Code for the Fibonacci Recurrence

The most natural was to describe the Fibonacci relation to Mathematica is as follows.

F[1] = F[2] = 1;
F[n_] := F[n-1] + F[n-2]

Note that if we want to make a table of Fibonacci numbers, this formula for F[n] is horribly inefficient.  Explain why.  You can use the Timing function in Mathematica to measure how long a calculation takes (in CPU seconds, not wall clock time).  Of course, the value of the Timing output will vary from machine to machine.

Timing[Table[F[n],{n,1,15}]]

Improved Code

Mathematica can be convinced to cache its previously computed information about a sequence as follows.

Clear[F2];
F2[1] = F2[2] = 1;
F2[n_] := F2[n] = F2[n-1] + F2[n-2]
Timing[Table[F2[n], {n,1,100}]]

This approach is great if n is not too large.  However, it is easy to exceed Mathematica's recursion depth if n is large.  You can experiment with this approach, but you should clear Mathematica's memory of the cached F2 values if you run Timing tests so you are not using previously stored information.

Timing[Table[F2[500]]]

Another Approach

The disadvantage in the above improvement is that in order to compute a Fibonacci number, all the intermediate numbers are stored as rules.  This will eventually require a large amount of storage for a term in the sequence if the index n is very large.  In fact, for sufficiently large n, the process will fail due to exceeded capacity.  Also, the lookup procedure for rules is slow.  These problems can be overcome with an iterative procedure.

F3[n_] :=
        Block[{f = 1, f1 = 0},
          Do[ {f, f1} = {f + f1, f}, {n - 1} ];
          f
        ]

Timing[Table[F3[n],{n,1,100}]][[1]]

Timing[F3[1000]]

This procedure really flies.  It does not have the recursion depth restriction of the previous procedure.

Graphs

To determine the growth rate of the Fibonacci sequence, we can plot the points.

ListPlot[Table[F3[n],{n,1,10}],AxesOrigin->{0,0}]

Is this exponential?  To see, we can plot the log of the data.  In fact, we can let the data go further out.

ListPlot[Log[Table[F3[n],{n,1,100}]]]

Since the log plot looks linear, the data are probably exponential in growth, at least over the range we are examining.  This observation is made more precise in the next subsection.

A Closed Form Method

In this subsection, we will derive the standard closed form formula for the Fibonacci sequence to show that the terms grow exponentially.

The method used to compute F3 can be rewritten in matirx form.

A = {{1,1},{1,0}};

Now the Fibonacci sequence can be described by a matrix equation with appropriate initial conditions.

X[n_] := A . X[n-1]; X[1] = {{1},{0}};

Since the eigenvectors of the matrix A are distinct, the basis can be changed to an eigenbasis (a basis consisting of distinct eigenvectors) in order to solve the system.

T = Transpose[Eigenvectors[A]];

Under this new basis, the coupled recurrence equations become uncoupled since the coefficient matrix A is transformed into a disgonal matrix of eigenvalues.

MatrixForm[Simplify[Inverse[T] . A . T]]

This shows that the formula, X[n], for the nth term of the Fibonacci sequence is a linear combination of powers, say nth, of the two eigenvalues.  A formula for X[n] can now be determined as follows.

Clear[c1, c2, e1, e2];
{e1,e2} = Eigenvalues[A];
sol=
  Solve[ 
    {c1 + c2 == 0, c1 e1 + c2 e2 == 1},
    {c1, c2} 
  ]

{c1, c2} = {c1, c2} /. sol[[1]]

This shows that a closed form formula for the nth term of the Fibonacci sequence is given by

OutputForm[c1 e1^n + c2 e2^n]

Since e1 is less than one in absolute value, the first term will go to zero as n gets large.  Thus for large n, the Fibonacci sequence is approximately c2*e2^n -- exponential!

To implement this closed form description as a Mathematica definition is easy.

F4[n_] := Expand[c1 e1^n + c2 e2^n]

Timing[F4[100]]

Timing[F4[1000]]

Notice that there is a considerable amount of algebra being done by Mathematica when using this formula to compute large terms in the sequence.

Malthusian Growth

Thomas Malthus (1766-1834) is generally credited with the idea that populations tend to grow exponentially.  Exponential growth (or decay) occurs whenever the rate of change over time of a variable is directly proportional to the value of the variable.  To see this in Mathematica, we use will K as the constant of proportionality.

solution1=
  DSolve[{P'[t] == K P[t], P[0] == a}, P[t], t]

We can plot several solution curves at once to see the behaviour of the solutions.

K = 1;
Plot[
  Evaluate[
    Table[
      P[t] /. solution1,
      {a,.4,4,.2}
    ]
  ],
  {t,0,4},
  PlotRange -> {0,20}
]

Logistic Growth

According to Malthus, although populations tend to increase exponentially, resources which support populations (e.g. food) tend to grow only linearly.  As the exponential growth of the population outstrips the linear growth of the resources, the resulting constraints would damp the population growth.

We can build in this dampening effect of excessive growth into the differential equation of the previous section.  Instead of P' being linear in P, we will make it quadratic in P.  This approach is called the logistic equation.  It was introduced by the Belgian statistician Pierre-Francois Verhulst (1804-1849).

P'[t_] = (1 - P[t]/K)*P[t]

With this model, the rate of change of P will be 0 when P=0 and when P=K.  Before we actually solve this DE, we will analyze it graphically using vector fields.

<<Graphics`PlotField`

Clear[K,f,P];
K=1;
f[t_,P_] := (1-P/K)*P

PlotVectorField[
  {1,f[t,P]},
  {t,0,4},
  {P,0,2},
  Axes->True,
  AxesLabel->{"t","P"}
]

Note that P' does not depend at all on t.  It only depends on P.  So all the columns of vectors in the above vector field are identical.

Clear[K,a,P]
K=1;

solution = 
  Table[
    DSolve[
      { P'[t] == (1 - P[t]/K)*P[t], P[0] == a }, 
      P[t],
      t
    ],
    {a,.1,1.7,.2}
  ]

We have to get rid of the outer set of braces to provide suitable input for the ParametricPlot function.  We do this with the Flatten command shown below.

Plot[
  Evaluate[P[t] /. Flatten[solution,1]],
  {t,0,4}
]

If, instead of having P' quadratic in P, we make P' cubic in P.  We can build a model which has three stationary points:  P = 0, T, K  (T<K).  The equilibrium point P=T represents the threshhold of a population - a number below which the population becomes extinct.  For 0<P<T, we want P'<0.  For T<P<K, we want P'>0.  And for P>K, we want P'<0.  A cubic function with these properties is easily constructed.

K=5; T=1;
g[P_] := (P/T - 1)*(1 - P/K)*P
Plot[g[P],{P,0,6},AxesLabel->{"P","P'"}]

We use this cubic function to assign values to P' in a differential equation.  The we plot solution curves for various values of the intial condition.

solution = 
  Table[
    NDSolve[
      { P'[t] == (P[t]/T - 1)*(1 - P[t]/K)*P[t], 
        P[0] == a 
      }, 
      P[t],
      {t,0,2}
    ],
    {a,.3,7,.3}
  ];

Plot[
  Evaluate[P[t] /. Flatten[solution,1]],
  {t,0,2},
  PlotRange->{0,7}
]

Lotka-Volterra Predator-Prey Populations

Vito Volterra was an Italian mathematician who lived from 1860-1940.  His son-in-law, Humberto D'Ancona was an Italian biologist who, in 1926, completed a statistical study of fish populations in the Adriadic Sea.  D'Ancona asked Volterra if there was a mathematical model which could explain the increase in predator fish (and corresponding decrease in prey fish) which he observed during the World War I period.  Within a couple of months, Volterra produced a series of models for the interaction of two or more species.  Alfred J. Lotka was an American biologist and actuary who independently produced many of the same models.

The simplest of their models can be described by letting 

  x(t) = predator population at time t;
  y(t) = prey population at time t;
  a     =  the natural rate of decay of the predator population if there are no prey;
  b     =  the natural rate of growth of the prey population if there are no prey;
  c     =  the efficiency with which predators convert prey encounters into offspring,
  d    =  the rate of decrease of prey due to encounter with predators.

The simplest Lotka-Volterra system for modeling predator-prey interactions is

x'[t] = -a x[t] + c x[t] y[t]
  y'[t] =  b y[t] - d x[t] y[t]

For a specific example, we consider the following.

eqn={x'[t]== -x[t] + x[t] y[t],
     y'[t]== 3y[t] - x[t] y[t],
     x[0]==x0,
     y[0]==1} ;

In order to graph several solution curves, we make a table of solutions.

tbl =
  Table[
    NDSolve[
      eqn, 
      {x,y},
      {t,0,6}
    ],
    {x0,.2,3,.3}
  ] ;

Then we can plot the solutions together.

ParametricPlot[
  Evaluate[
    { x[t],y[t] } /. Flatten[tbl,1]
  ],
  {t,0,6}
]

There are many variations on the Lotka-Volterra equations.  Here is one example which introduces an additional quadratic term in the rate of change of y.

eqn={x'[t]== -x[t] + x[t] y[t],
     y'[t]== 3y[t] - x[t] y[t] - y[t]^2,
     x[0]==x0,
     y[0]==1} ;

tbl =
  Table[
    NDSolve[
      eqn, 
      {x,y},
      {t,-3,3}
    ],
    {x0,.2,1.7,.3}
  ] ;

ParametricPlot[
  Evaluate[
    { x[t],y[t] } /. Flatten[tbl,1]
  ],
  {t,-3,3}
]

Competitive Exclusion

Instead of a predator-prey relationship, two species may compete with each other for the same niche in an ecosystem.  How can this relationship be modeled?  As you would expect, there are many ways to do this.  Perhaps the most direct way is to modify the predator-prey system as follows:

x'[t] =  a x[t] - b x[t] y[t]
  y'[t] =  c y[t] - d x[t] y[t]

You get an opportunity to explain the differential equations in the exercises below!

First, we get a feel for the solutions to the system by plotting the vector field.

<< Graphics`PlotField`

\(PlotVectorField[{0.1\ x - 0.01\ x*y, 
        0.2\ y - 0.05\ x*y}, 
\ \ \ \ \ \ \ \ \ \ \ \ \ {x, \ 0, \ 
        12}, \ {y, \ 0, \ 15}, \ Frame -> True, \ 
      PlotPoints -> 20]\)

For a specific example, we consider the following.

eqn={x'[t]== 0.1 x[t] - 0.01 x[t] y[t],
     y'[t]== 0.2 y[t] - 0.05 x[t] y[t],
     x[0]==x0,
     y[0]==y0} ;

In order to graph several solution curves, we make a table of solutions.

tbl =
  Table[
    NDSolve[
      eqn, 
      {x,y},
      {t,0,20}
    ],
    {x0,1,9,1},
    {y0,1,19,2}
  ] ;

Then we can plot the solutions together.

ParametricPlot[
  Evaluate[
    { x[t],y[t] } /. Flatten[tbl,2]
  ],
  {t,0,20},PlotRange->{{0,10},{0,20}}
]

Exercise 1

Exercise 2

years   = Table[1790 + 10n, {n,0,21}];
pop     = {3.929214,5.308483,7.239881,9.638453,
			12.860702,17.063353,23.191876,31.443321,
			38.558371,50.189209,62.979766,76.212168,
			92.228496,106.021537,123.202624,132.164569,
			151.325798,179.323175,203.302031,226.542199,
            248.709873,281.421906};
data    = Table[{years[[i]], pop[[i]]}, {i,1,22
}];

ListPlot[
  data,
  AxesOrigin -> {1780,0},
  Prolog -> PointSize[.016]
]

Part (a)

DSolve[{P'[t] == r*P[t], P[1790] == 3.9}, P[t], t] // 
      Simplify

So we have one parameter, r, to vary in order to find the best fit to the dataset.

<< Statistics`NonlinearFit`

fitdata1 = 
      NonlinearFit[data, 3.9*E^(r*(t - 1790)), t, {r, 0}, 
        ShowProgress -> True]

You can also leave the initial population unspecified and let the power of the NonlinearFit function use that extra variability to improve the exponential fit to the data.

DSolve[{P'[t] == r*P[t], P[1790] == a}, P[t], t] // 
      Simplify

fitdata1 = 
      NonlinearFit[data, a*E^(r*(t - 1790)), t, {{a, 3}, {r, 0}}, 
        ShowProgress -> True]

Part (b)

Exercise 3

Exercise 4

Exercise 5

The equation dP/dt = r*P*(1-P/K) - H gives a logistic type model of a population with natural rate r, carrying capacity K and constant harvesting rate H (r, K, H positive).

Part A

Part B

Exercise 6

Part A

Part B

Exercise 7