MAT 581 PROJECT 2 TOPICS
Introduction: Projects are mostly computational (programming) assignments that
are larger in scope than a homework problem. Your work should consist of a Matlab
file with a program and a separate file that explains your method. A project can be
done either individually or by a group of 2 students (in the latter case, the names of
both students must be included).
You should choose one of the following topics for your project. But be sure to read
all of the descriptions, because the description of one project may provide tools for
another one. The common theme of project 2 topics is differential equations.
1
2 MAT 581 PROJECT 2 TOPICS
Topic 1. The period and average population size in the Lotka-Volterra system.
Description. Consider the Lotka-Volterra system
R′ = aR− bRF
F′ = −cR + dRF(1)
with two positive values as initial conditions (R0, F0). The solution of this system is
known to be periodic. Compute:
• The least period T of the solution (smallest number T such that y(t+ T) = y(t)
for all t)
• The average values Ra and Fa of the two populations (average is taken over a
period).
The script should also plot the functions R, F on the interval [0, T].
Based on your experiment, answer the question: does the period T depend only on
the coefficients a, b, c, d or also on the initial condition?
You can use ode45 to solve the system. It can be assumed that the period is less
than 1000. A natural starting point is to obtain a numerical solution y with equally
spaced t-values on the interval [0, 1000]. Then, after a period is found, the values of y
over one period can be averaged, for example using mean(y(m:n, :), 1) where m, n
are the beginning and end of one period. The tricky part is finding the period.
Test cases (note that the output values may vary slightly based on your method):
• Data: a = 0.1, b = 0.004, c = 0.2, d = 0.001, R0 = 100, F0 = 30. Output:
T ≈ 47, Ra ≈ 200, Fa ≈ 25.
• Data: a = 0.2, b = 0.01, c = 0.25, d = 0.007, R0 = 70, F0 = 40. Output: T ≈ 31,
Ra ≈ 36, Fa ≈ 20.
Potentially useful Matlab commands:
• [vals, locs] = findpeaks(z) finds the local maxima within the vector z.
• find(z) returns all indices k such that z(k) is not zero. It is often used with
logical expressions, like find(z > 0)
• locs = find(z(2:end-1) > z(1:end-2) & z(2:end-1) > z(3:end)) + 1; is
an alternative to findpeaks.
• diff(z) takes differences of consecutive elements: z(2)-z(1), z(3)-z(2), ...
• mean(A, 1) is the average of all rows of matrix A, while mean(A, 2) is the
average of all columns.
MAT 581 PROJECT 2 TOPICS 3
Topic 2. Resource extraction in the Lotka-Volterra model. Suppose that the equation for R′
has an additional term −E, representing the rate at which the resource is extracted from the
system, probably by humans (hunted, fished, gathered, etc). Find the value at E at which the
population R goes extinct at given time.
Description. In the modified Lotka-Volterra model
R′ = aR− bRF− E
F′ = −cR + dRF(2)
the following parameters are considered fixed: a = 0.2, b = 0.005, c = 0.3, d = 0.003,
R0 = 100, F0 = 30. Your goal is to determine the parameter E so that R(200) = 0 (the
rabbits go extinct at time T = 200). The script should show the critical value of E and
the population graph for this value.
When working on this project, it may be helpful to plot the solution for a few values
of E to get an idea of the range in which the solution may be found.
Once you have a reasonable bracketing interval, either the bisection method or
fzero can be used to solve for E from the condition R(200) = 0. Note that one can
use fzero either with anonymous function or with a named function, for example
E = fzero(@f, [A B]);
...
function finalR = f(E)
... (solve the system and return R(200))
end
Since each evaluation of the function involves solving a system of differential equa-
tions, this can take a long time if the bracketing interval [A, B] is too wide. In addition,
large values of E can drive the solution of the system deep into negative numbers (not
meaningful) where it will rapidly grow in size, posing a challenge for ode45 solver.
Remember that Matlab execution can be interrupted with Ctrl-C or its equivalent,
when the Command Window is in focus.
Test case: if T was 100 instead of 200, the critical value of E would be about 6.4.
4 MAT 581 PROJECT 2 TOPICS
Topic 3. Implement SIR model with an arbitrary number of groups.
Description. Given a positive integer n (provided by a user with input) find a numeric
solution of the SIR model with n ≥ 2 different groups with the following parameters.
• Individuals of group k can be infected either by those from within group k, or
by those from its neighboring groups k− 1 and k + 1. (Group 1 and group n
have only one neighboring group each.)
• In-group infection rate is β = 0.2, meaning equations will contain the term
βSk Ik.
• Between-groups infection rate is γ = 0.03, so the terms γSk Ik−1 and γSk Ik+1
will appear in the equations.
• The recovery rate is γ = 0.1 for all groups, meaning there will be terms γIk.
• The initial conditions are: S1 = 0.99, the rest of Sk are 1; I1 = 0.01, the rest of
Ik are 0.
• To simplify the model, do not include “R” group: recovered individuals are
simply subtracted from Ik without being added anywhere.
• Solve the system on the interval [0, 100].
You will need a named function, like function yp = f(t, y) to implement this
system. The ode45 solver can be used with a named function, like
[t, y] = ode45(@f, [0, 100], y0). It is best to arrange the vector y so that its com-
ponents represent S1, . . . , Sn, I1, . . . , In in this order.
The output must contain:
• The plot of solution, plot(t, y)
• The peak infection rate in each group; that is, the maximum values of I1, . . . , In.
Test case: when n = 4, the peak infection rates are 0.2011, 0.2367, 0.2438, 0.2261.
Potentially useful Matlab commands:
• max(A, [], 1) finds the maximum value in each column of A. (The empty
argument [] is necessary because max(A, 1) is understood as taking the max-
imum of A and 1.)
• max(A, [], 2) finds the maximum value in each row of A.







































































































学霸联盟