(AM)^2 REU


Below are exercises and MatLab scripts that support the mathematical biology module of the joint ASU-SCC REU.

MatLab Scripts (programs)

The following file (LogisticMap) is a MatLab .m file that runs the logistic map from the introductory exercise. It also makes a cobweb diagram. Due to security constraints on this server, it had to be saved as a .doc file instead of MatLab's standard .m. Therefore, to run the script you'll need to do the following:

  1. Download the .doc file but do not open it in Word or any other word processor.
  2. Open the file in a text editor (or even in the MatLab editor).
  3. Copy the code completely from the text editor.
  4. Open a new script file in MatLab.
  5. Paste the code into the new file and save it under any name with a .m extension.

MatLab logistic model script: LogisticMap The following example scripts can be copied and pasted into the MatLab editor and saved with any file name followed by a .m extension. Beware, though, that you may find formatting problems if you just cut-and-paste, especially the quotation marks. Experienced programmers may note that these are not the most efficient algorithms, but they have the advantage of transparency, making it easier for someone who is not a MatLab jockey to follow the logic.

% This script runs the discrete-time logistic map. 

% Written by John Nagy 
% June 2013 (revised May 2016)

% Initialize and set parameters 

% Clear graphical and work spaces 
hold off

% Set demographic parameters (birth, death, competition rates)

lambda = 2; 
mu = 0.6; % Must be between 0 and 1, inclusive 
c = 0.1;

% Compute the derived parameters

r = lambda – mu; 
K = r/c;

% Set the time horizon

maxT = 50;

% Preallocate the main data vector and set the initial number of animals

x = zeros(maxT,1); 
x(1) = 5; % In hundreds.

% Run the model 

for i=2:maxT 
x(i) = x(i-1) + r*x(i-1)*(1-x(i-1)/K); 

% Plot 1: Time series 

plot(1:maxT, x,’-o’) % Plots the basic time series as open points and lines 
ax = gca; % Gain control of the current plot axes 
line([0 maxT], [K K], ‘color’, ‘red’) % Plots carrying capacity in red 
ax.XAxis.FontSize = 12; % Adjust font size of axis marks (numbers) 
ax.YAxis.FontSize = 12; 
xlabel(‘Time (Number of generations)’, ‘interpreter’, ‘latex’, ‘FontSize’, 14) 
ylabel(‘Number of moose (Hundreds)’, ‘interpreter’, ‘latex’, ‘FontSize’, 14)

% Plot 2: Cobweb plot 

figure % Create a new plot window

% Calculate the positive portion of the r.h.s of the model (return map) 
% [Next 3 lines]

xbaseMax = K*(1+r)/r; 
xbase = 0:0.01:xbaseMax; 
f = xbase + r*xbase.*(1-xbase/K);

% Plot the return map and the identity line 
% [Next 3 lines]

plot(xbase, f, ‘-b’) 
ax=gca; % see above 
ax.XAxis.FontSize = 12; % see above 
ax.YAxis.FontSize = 12; 
xlabel(‘Number of moose in generation $n$ (hundreds)’, ‘interpreter’, ‘latex’, ‘FontSize’, 14) 
ylabel(‘Number of moose in generation $n+1$ (hundreds)’, ‘interpreter’, ‘latex’, ‘FontSize’, 14) 
line([0 ax.XAxis.Limits(2)], [0 ax.XAxis.Limits(2)], ‘color’, ‘black’)

% Overlay the solution x onto the return map 
% [Next 2 lines]

hold on 

% Plot the cobweb

line([x(1) x(1)], [0 x(2)], ‘Color’, ‘r’) 
line([x(1) x(2)], [x(2) x(2)], ‘Color’, ‘r’) 
for i = 2:(maxT-1) 
line([x(i) x(i)], [x(i), x(i+1)], ‘Color’, ‘r’) 
line([x(i) x(i+1)], [x(i+1) x(i+1)], ‘Color’, ‘r’) 


gen_roots function

The following script requires the function gen_roots, which must be copied verbatim with the filename "gen_roots.m" (including underscore) and saved in the same file as this one. The gen_roots function is listed after this script.


% MCTP example Nullcline Plot 

% See Linear Statbility Theory handout 

% Written by John Nagy % June 2013

% Parameterize

K = 1; 
r = 0.2; 
a = 2;

Max_u = r/a;

% Define the horizontal axis (predators) 
u = eps:0.001*Max_u:Max_u;

v_null = K*(1 - u*a/r);

u_null = u./(1-exp(-a*u));

plot(u,v_null, u,u_null,'-r')

xlabel('Number of predators (\itu)') 
ylabel('Number of prey (\itv)')

% Find the actual fixed point 
% First, subtract the two nullclines:

Null = v_null - u_null;

%Then fine the root:

uFP = gen_roots(u, Null); 
vFP = K*(1 - uFP*a/r);

% Plot the fixed pointsline([uFP uFP], [0 K], 'color', 'black', 'LineStyle', '--') 
line([0, Max_u], [vFP vFP], 'color', 'black', 'LineStyle', '--')

% In the plot, write the values of the fixed points: 
% First, convert the values to strings and write the string to be placed 
% in the plot

uFP = num2str(uFP); 
vFP = num2str(vFP);

%The next construction concatinates the strings into a single string. 
uStr = strcat('u fixed point = ', uFP); 
vStr = strcat('v fixed point = ', vFP);

%Now, write the text on the plot 
text(0.06, 0.8, uStr) 
text(0.06, 0.76, vStr)


function rset = gen_roots(d,f) 
% Returns a column vector, rset, that represents the roots of a general 
% function, f, obtained by simple linear interpolation It takes two 
% vectors, d (domain, x-values) and f (function, y-values). If no root 
% (zero) exists, then rset = NaN. 

% Written December 2009 
% Written by John D. Nagy

len = length(d); 
s = sign(f);

if abs(sum(s)) == len rset = NaN; 
% the abs val of sum of sign function = len iff no roots 

% Now we know there's at least one root


for i=1:len-1 
if f(i) == 0 
rset=[rset; d(i)]; 
elseif s(i)+s(i+1) == 0 
interp = (f(i)*d(i+1)-f(i+1)*d(i))/(f(i)-f(i+1)); 
rset=[rset; interp];