%% Population Dynamics
% We are going to look at two models for population growth of a
% species.  The first is a standard exponential growth/decay model
% that describes quite well the population of a species becoming
% extinct, or the short-term behavior of a population growing in an
% unchecked fashion.  The second, more realistic, model describes the
% growth of a species subject to constraints of space, food supply,
% and competitors/predators. 

%% Exponential Growth/Decay
% We assume that the species starts with an initial population $P_0$.
% The population after $n$ time units is denoted $P_n$.  Suppose that,
% in each time interval, the population increases or decreases by a
% fixed proportion of its value at the beginning of the interval.
% Thus,
%
% $$ P_{n+1} = P_n + r P_n, \qquad n \geq 0.$$
%
% The constant $r$ represents the difference between the birth rate
% and the death rate.  The population increases if $r$ is positive,
% decreases if $r$ is negative, and remains fixed if $r = 0$.
%
% Here is a simple M-file that we will use to compute the population
% at stage $n$, given the population at the previous stage and the
% rate $r$:

type itseq

%%
% In fact, this is a simple program for computing iteratively the values
% of a sequence $x_{k+1} = f(x_k, r), n \geq 0$, given any function
% $f$, the value of its parameter $r$, and the initial value $x_0$ of
% the sequence.
%
% Now let's use the program to compute two populations at five-year
% intervals for $r = 0.1$ and then $r = -0.1$:

Xinit = 100; f = @(x, r) x*(1 + r);
X = itseq(f, Xinit, 100, 0.1);  
format long; X(1:5:101)  

%%

X = itseq(f, Xinit, 100, -0.1); X(1:5:101)  

%%
% In the first case, the population is growing rapidly; in the second,
% decaying rapidly.  In fact, it is clear from the model that, for any
% $n$, the quotient $P_n/P_{n+1} = (1 + r)$, and therefore it follows
% that $P_n = P_0(1 + r)^n, n \geq 0$.  This accounts for the
% expression ``exponential growth/decay.''  The model predicts a
% population growth without bound (for growing populations), and is
% therefore not realistic.  Our next model allows for a check on the
% population caused by limited space, limited food supply,
% competitors, and predators.

%% Logistic Growth
% The previous model assumes that the relative change in population is
% constant, that is
%
% $$(P_{n+1} - P_n)/P_n =  r.$$
%
% Now let's build in a term that holds down the growth, namely
%
% $$(P_{n+1} - P_n)/P_n =  r - uP_n.$$
%
% We shall simplify matters by assuming that $u = 1 + r$, so that our
% recursion relation becomes
%
% $$P_{n+1} = u P_n(1 - P_n),$$
%
% where $u$ is a positive constant.  In this model, the population $P$
% is constrained to lie between $0$ and $1$, and should be interpreted
% as a percentage of a maximum possible population in the environment
% in question.  So let us define the function we will use in the
% iterative procedure:

f = @(x, u) u*x*(1 - x);

%%
% Now let's compute a few examples, and use *|plot|* to display the
% results:

Xinit = 0.5; X = itseq(f, Xinit, 20, 0.5); plot(X)

%%

X = itseq(f, Xinit, 20, 1); plot(X) 

%%

X = itseq(f, Xinit, 20, 1.5); plot(X)

%%

X = itseq(f, Xinit, 20, 3.4); plot(X)  

%%
% In the first computation, we have used our iterative program to
% compute the population density for $20$ time intervals, assuming a
% logistic growth constant $u = 0.5$, and an initial population density
% of 50%.  The population seems to be dying out.  In the remaining
% examples, we kept the initial population density at 50%; the only
% thing we varied was the logistic growth constant.  In the second
% example, with a growth constant $u = 1$, once again the population
% is dying out -- although more slowly.  In the third example, with a
% growth constant of $1.5$, the population seems to be stabilizing at
% 33.3...%.  Finally, in the last example, with a constant of $3.4$, the
% population seems to oscillate between densities of approximately 45%
% and 84%.
%
% These examples illustrate the remarkable features of the logistic
% population dynamics model.  This model has been studied for more
% than $150$ years, its origins lying in an analysis by the Belgian
% mathematician Pierre-Francois Verhulst.  Here are some of the facts
% associated with this model.  We will corroborate some of them with
% MATLAB.   In particular, we shall use *|bar|* as well as *|plot|*
% to display some of the data.  
%
% * *The logistic constant cannot be larger than $4$.*
%
% In order for the model to work, the output at any point must be
% between $0$ and $1$.  But the parabola $ux(1 - x)$, for $0 \leq x
% \leq 1$, has its maximum height when $x = 1/2$, where its value is
% $u/4$.  To keep that number between $0$ and $1$, we must restrict $u
% \leq 4$.  Here is what happens if $u$ is greater than $4$: 

X = itseq(f, 0.5, 10, 4.5)

%%
% * *If $0 \leq u \leq 1$, the population density tends to zero for any
% initial value.*

X = itseq(f, 0.5, 100, 0.8); X(101)  

%%

X = itseq(f, 0.5, 20, 1); bar(X)  

%%
% * *If $1 < u \leq 3$, the population will stabilize at density $1 -
% 1/u$ for any initial density other than zero.*
%
% The third of the original four examples corroborates the assertion
% (with $u = 1.5$ and $1 - 1/u = 1/3$). In the following examples, we
% set $u = 2$, $2.5$, and $3$, respectively, so that $1 - 1/u$ equals
% $0.5$, $0.6$, and $0.666\ldots$, respectively.  The convergence in the
% last computation is rather slow (as might be expected from a boundary
% case -- or ``bifurcation point'').  

X = itseq(f, 0.25, 100, 2); X(101)  

%%

X = itseq(f, 0.75, 100, 2); X(101)  

%%

X = itseq(f, 0.5, 20, 2.5);  
plot(X)  

%%

X = itseq(f, 0.5, 100, 3);  
bar(X); axis([0 100 0 0.8])

%%
% * *If $3 < u < 3.56994\ldots$, then there is a periodic cycle.*
%
% The theory is quite subtle.  For a fuller explanation, the reader
% may consult _Encounters with Chaos_, by Denny Gulick, McGraw-Hill, 
% New York, 1992,
% Section 1.5.  In fact, there is a sequence
%
% $$u_0 = 3 < u_1 = 1 + \sqrt{6} < u_2 < u_3 < \cdots < 4,$$
%
% such that between $u_0$ and $u_1$ there is a cycle of period $2$;
% between $u_1$ and $u_2$ there is cycle of period $4$; and, in
% general, between $u_k$ and $u_{k+1}$ there is a cycle of period
% $2^{k+1}$.  In fact, at least for small $k$, we have
% the approximation $u_{k+1} \approx 1 + \sqrt{3+u_k}$.  So:

u1 = 1 + sqrt(6)

%%

u2approx = 1 + sqrt(3 + u1)

%%
% This explains the oscillatory behavior we saw in the last of the
% original four examples (with $u_0 < u = 3.4 < u_1$).  Here is the
% behavior for $u_1 < u = 3.5 < u_2$.  The command *|bar|* is
% particularly effective here for spotting the cycle of order 4.

X = itseq(f, 0.75, 100, 3.5);  
bar(X); axis([0 100 0 0.9])

%%
% * *There is a value $u < 4$ beyond which -- chaos!*
%
% It is possible to prove that the sequence $u_k$ tends to a limit
% $u_\infty$.  The value of $u_\infty$, sometimes called the
% ``Feigenbaum parameter,'' is approximately $3.56994\ldots$.  Let's see
% what happens if we use a value of $u$ between the Feigenbaum
% parameter and $4$. 

X = itseq(f, 0.75, 100, 3.7);
plot(X); axis([0 100 0 1])

%%
% This is an example of what mathematicians call a ``chaotic'' phenomenon!
% It is not random -- the sequence was generated by a precise, fixed
% mathematical procedure, but the results manifest no discernible
% pattern.  Chaotic phenomena are unpredictable, but with modern
% methods (including computer analysis), mathematicians have been able
% to identify certain patterns of behavior in chaotic phenomena.  For
% example, the last figure suggests the possibility of unstable
% periodic cycles and other recurring phenomena.  Indeed, a great deal
% of information is known.  The aforementioned book by Gulick is a
% fine reference, as well as a source of an excellent bibliography
% on the subject.

%% Re-running the Model with Simulink
%
% The logistic growth model that we have been exploring lends itself
% particularly well to simulation using Simulink.  Here is a simple
% Simulink model that corresponds to the above calculations: 

open_system popdy

%%
% Let's briefly explain how this works.  If you ignore the Pulse
% Generator block and the Sum block in the lower left for a moment,
% this model implements the equation
%
% $$x \mbox{ at next time} = u x (1 - x) \mbox{ at current time,}$$
%
% which is the equation for the logistic model.  The Scope block
% displays a plot of $x$ as a function of (discrete)
% time.  However, we need somehow to build in the initial condition
% for $x$.  The simplest way to do this is as illustrated here: we add
% to the right-hand side a discrete pulse that is the initial value of
% $x$ (here we use $0.5$) at time $t = 0$ and is $0$ thereafter.
% Since the model is 
% discrete, you can achieve this by setting the Pulse Generator block
% to ``Sample based'' mode, setting the period of the pulse to
% something longer than the length of the simulation, setting the
% width of the pulse to $1$, and setting the amplitude of the pulse to
% the initial value of $x$.  The outputs from the model in the two
% interesting cases of $u = 3.4$ and $u = 3.7$ are shown here:
set_param('popdy/Logistic Constant', 'Value', '3.4')
popout = sim('popdy', 'SaveTime','on','TimeSaveName','t', ...
   'SaveState','on','StateSaveName','x','StopTime', '120');
t = popout.find('t'); x = popout.find('x');
simplot(t, x); title('u = 3.4')
set(gcf, 'InvertHardcopy', 'off')
%%
% In the first case of $u = 3.4$, the periodic behavior is clearly 
% visible.

set_param('popdy/Logistic Constant', 'Value', '3.7')
popout = sim('popdy', 'SaveTime','on','TimeSaveName','t', ...
   'SaveState','on','StateSaveName','x','StopTime', '120');
t = popout.find('t'); x = popout.find('x');
simplot(t, x); title('u = 3.7')
set(gcf, 'InvertHardcopy', 'off')
%%
% On the other hand, when $u = 3.7$, we get chaotic behavior.