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:
- Download the .doc file but do not open it in Word or any other word processor.
- Open the file in a text editor (or even in the MatLab editor).
- Copy the code completely from the text editor.
- Open a new script file in MatLab.
- 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.
% LOGISTIC MODEL TEACHING MODULE
% 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
clf
clear
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);
end
%———————————————————————-
% 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(x(1:(maxT-1)),x(2:maxT),’or’)
% 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’)
end
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
return
end
% Now we know there's at least one root
rset=[];
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];
end
end