FE2D:
Predatorprey simulations in arbitrary shaped 2D habitats
by Marcus R. Garvie (mgarvie@uoguelph.ca)
Back
Introduction
FE2D is a collection of MATLAB routines using the finite
element method for simulating the dynamics of predatorprey interactions
in two space dimensions and time. The codes generalize my earlier predatorprey codes
solved on the square (see FD2D)
to the situation of arbitrary shaped domains with general boundary
conditions (functions depending on space and time of Neumann, Dirichlet,
Robin, Periodic, or mixed type).
The MATLAB code is mostly self explanatory, with the names of
variables and parameters corresponding to the symbols used in
the finite element method described in the paper referenced below.
Copies of the MATLAB codes are freely available via the links below.
The code employs the sparse matrix facilities of MATLAB when solving
the linear systems, which provides advantages in both matrix storage
and computation time. The code is vectorized to minimize the
number of "forloops" and conditional "ifthenelse" statements,
which again helps speed up the computations.
The linear systems are solved using MATLAB's built in function
gmres.m. The gmres.m algorithm in
MATLAB requires a number of input arguments. For some simulations
we found it acceptable to use no
"restarts" of the iterative method, or, preconditioners, and a
tolerance for the relative error of 1E6. In practise, the
user will need to experiment with the restart value, tolerance, and
"maxit" (maximum number of outer iterations) in order to achieve
satisfactory rates of convergence of the gmres.m function.
For definitions of these input arguments (and others) see
the description in MATLAB's Help pages. We remark that a pure C or
FORTRAN code is likely to be faster than our codes, but with the
disadvantage of much greater complexity and length.
The user is prompted for all the necessary parameters, timestep, initial data
functions, and boundary data functions. However, the grid must be supplied by the user
(further details given below).
The reactiondiffusion system models spatially extended predatorprey interactions of
RosenzweigMacArthur form, namely:
Kinetics 1:
where u(x,t) and v(x,t) are the prey and predator densities at time t
and vector position x, and Delta is the usual Laplacian operator, and the parameters
alpha, beta, gamma, and delta are strictly positive.
Grid generation
Before the finite element codes can be run in MATLAB it is necessary to
triangulate the region using an appropriate unstructured grid generator in 2D.
As a starting point, we assume that the user has provided a list of node coordinates,
ordered clockwise or counterclockwise, corresponding to the boundary of the
computational domain. For simple geometries this is easily worked out by hand, but
for more complicated domains the use of a drawing package may be advisable.
The grid generator then needs to output two arrays that define the mesh, namely

t, a list of node indices corresponding to M triangles, given by an 3byM array;

p, a list of N node coordinates, given by an 2byN array;
We found it convenient to use an unstructured mesh generator (Mesh2d v24) in MATLAB freely available from
http://www.mathworks.com/matlabcentral/fileexchange/25555mesh2dautomaticmeshgeneration". To help understand how the arrays "p" and
"t" are used to define a mesh we provide a simple example below. For example, the 2nd triangle (colured blue) in the list for array
"t" has vertices labelled 2,4,3. And from array "p" we can see for example, that node 4 has coordinates (4,5).

t = [ 1 2 1 2 5; ...
2 4 5 5 6; ...
3 3 2 4 4 ];
p = [ 1 3 2 4 4 5; ...
2 2 4 5 1 2 ];

Using the mesh generator Mesh2d
For convenience we provide an example of a 'frontend' in MATLAB for constructing an unstructured grid using using Mesh2d
(see above), for the case of a hypothetical lake with an island. Initially a list of nodes defining the boundary of the lake
and the island are supplied. Then after calling Mesh2d to construct the mesh, the arrays p and t that define
the mesh are exported as files 'p_coord.dat' and 't_triang.dat' respectively. We also export the list of gridnodes on the edge
of the lake and on the island as files 'bn2_nodes.dat' and 'bn1_nodes.dat' respectively. An example mesh and the MATLAB frontend
that can be copied and adapted are given below:
mesh2d_example.m,
Frontend code for meshing a hypothetical lake with an island.
lake_mesh.jpg
JPEG image of a coarse grid generated from this example.
Linebyline description of the code for fe2d_nd.m
The codes are mostly selfexplanatory, with comments to explain what each section of the code does. For example, the
structure of fe2d_nd.m is outlined below. The structure of the other codes are similar. The linenumbered codes are listed
at the bottom of the page in PDF format.

Lines 23  46: External data files for the mesh are loaded.

Lines 50  78: User prompted for the parameter values, initial data functions,
Dirichlet data functions, and Neumann data functions. The functions are entered as strings
(allowable formats discussed below) and then converted to anonymous functions.

Lines 82  140: The Stiffness Matrix K and the Finite Element Matrix L are assembled.

Lines 147  197: The linear system is solved repeatedly until time T .

Lines 149  157: The coefficient matrices A2, B2 and C2 are updated with the
current solutions, as is the righthandside of the linear system.

Lines 158  173: The Neumann boundary conditions are imposed for each
edge of Gamma2.

Lines 174186: The Dirichlet boundary conditions are imposed for each node
on Gamma1.

Lines 187  189: The Incomplete LU factorizations of the coefficient matrices
A2 and C2 are computed to provide preconditioners for the GMRES iterative solver.

Lines 190196: The linear system is solved using the GMRES iterative solver.

Lines 201  214: The finite element solutions for u and v are plotted using a triangular
surface plot at time T (a vertical 'colorbar' provides a scale).
The initial data
The user is prompted for the initial data, which can be functions depending on space, or just constants.
Unlike with the simpler PRED_PREY_SIM codes no special format is required for entering the initial
data functions.
An example with an acceptable format is the following:
>> Enter initial data function u0(x,y) 0.2*exp((x^2+y^2))
>> Enter initial predator function v0(x,y) 1.0
We can also define functions in a piecewise fashion. For example, with
Omega=[0,100]^{2}, in order to choose an initial
predator density of 0.2 within a circle with radius 5 and center
(50,50), and a density of 0 elsewhere on Omega
we input the following:
>> Enter initial data function u0(x,y) 0.2*((x50)^2+(y50)^2<25)
The boundary conditions
The boundary of the whole domain is denoted Gamma. In the mixed boundary condition cases we assume
Gamma is the union of two pieces Gamma1 and Gamma2. The code allows for the following boundary
conditions:
 Pure Neumann boundary conditions on Gamma.
 Pure Dirichlet boundary conditions on Gamma.
 Pure Robin boundary conditions on Gamma.
 Mixed boundary conditions with a Dirichlet condition on Gamma1 and a Neumann condition on Gamma2.
 Mixed boundary conditions with a Robin condition on Gamma1 and a Neumann condition on Gamma2.
 Periodic boundary conditions on Gamma.
The Neumann, Dirichlet, and Robin boundary conditions have the following form:
where n denotes the unit outward normal to the domain Gamma and k_{1} and
k_{2} are real numbers (positive, negative, or zero). The Neumann and Dirichlet boundary conditions can
be functions of space (x and y) and time t.The periodic boundary conditions are used for
the problem defined over a square of arbitrary side length L (the code is self contained and thus a grid does not need
to be imported). If we label the sides of a square bn1, bn2, bn3 and bn4
then at each time step the boundary values on sides bn2 and bn1 are used to overwrite the boundary values on sides bn4 and bn3 respectively,
thus enforcing the periodicity of the boundary conditions (see the above figure).
There is no special format for entering the boundary functions, for example when running fe2dx_nd.m we might have:
>> Enter the Dirichlet b.c. g1u(x,y,t) for u 0.2*sin(xy)+t/10+0.2
>> Enter the Dirichlet b.c. g1v(x,y,t) for v 0.3*cos(xy)+t/20+0.3
>> Enter the Neumann b.c. g2v(x,y,t) for v 0.0
>> Enter the Neumann b.c. g2v(x,y,t) for v 0.6
Some practical issues
Firstly, bear in mind that if you run a simulation with a large domain size
and large final time T, coupled with small temporal and spatial discretization
parameters, then the runtime in MATLAB can be prohibative.
Another point concerns the choice of parameters alpha, beta, and gamma
used to run the code. In order for the local kinetics of the systems
to be biologically meaningful, there are restrictions on the parameters that need to be satisfied
(for further details see the references below).
Another issue concerns the stability and accuracy of the numerical solutions. The user needs to choose the time
step sufficiently small to obtain converged solutions. This can easily be checked by reducing the timestep until
the solutions corresponding to scheme 1 and scheme 2 (see the references below) are the same. Experience has shown
that with many examples the time step needs to be of the order of 1/384.
One final point concerns the code with the Robin boundary conditions. If you run the code with k1 or k2
large and positive on part of the boundary then this part of the boundary acts as a source for the domain.
Thus solvers may fail due to growing solutions, which leads to illconditioned coefficient matrices.
No such problems occur if we choose k1 or k2 large and negative (although solutions may go negative if the
time step is too large).
Nonconvergence of GMRES
The GMRES algorithm does not always converge, or it may converge very slowly. One approach that we have
incorporated into the codes to help overcome this problem is to use MATLAB's 'incomplete lu factorization' algorithm ILU
(which replaces the previously used LUINC function). This
provides a preconditioner for GMRES. See MATLAB's Help page under GMRES for usage.
If this still doesn't help you can try other preconditioners, or use a less sophisticated iterative
method for solving the linear systems. For example, the Jacobi and GaussSeidel algorithms are widely used,
although the runtimes will likely be considerably longer than for GMRES. As these methods are not part of
the standard suite of MATLAB functions they are provided below. The algorithms require that all diagonal
elements a_{ii} of the coefficient matrix A are nonzero, and they will not always converge even then.
Files you may copy include:

jacobi.m,
MATLAB code for iteratively solving the (square) system of linear equations Ax = b.

gauss_seidel.m,
A variation (usually an improvement) of JACOBI.
References

Garvie M.R.
, "Finite Difference Schemes for ReactionDiffusion Equations
Modeling PredatorPrey Interactions in MATLAB," Bulletin of Mathematical
Biology (2007) Vol.69, No.3., pp. 931956

Garvie M.R.
, Burkardt J., Morgan J. : "Simple Finite Element Methods for
Approximating PredatorPrey Dynamics in Two Dimensions using MATLAB," Bulletin of Mathematical
Biology (2015) Vol.77, No.3, pp. 548578.
Download the 'slow' codes for FE2D
These are the codes described in paper no. 2 referenced above. For more computationally intensive problems the 'fast' versions of these codes
are recommended (see below).
Files you may copy include the MATLAB Mfiles and PDF versions with line numbering:

fe2d_n.m, fe2d_n.pdf
Scheme 2 applied to Kinetics 1 with pure Neumann boundary conditions.

fe2dx_n.m, fe2dx_n.pdf
Scheme 1 applied to Kinetics 1 with pure Neumann boundary conditions.

fe2d_d.m, fe2d_d.pdf
Scheme 2 applied to kinetics 1 with pure Dirichlet boundary conditions.

fe2dx_d.m, fe2dx_d.pdf
Scheme 1 applied to kinetics 1 with pure Dirichlet boundary conditions.

fe2d_r.m, fe2d_r.pdf
Scheme 2 applied to kinetics 1 with pure Robin boundary conditions.

fe2dx_r.m, fe2dx_r.pdf
Scheme 1 applied to kinetics 1 with pure Robin boundary conditions.

fe2d_p.m, fe2d_p.pdf
Scheme 2 applied to kinetics 1 with periodic boundary conditions over the square.

fe2dx_p.m, fe2dx_p.pdf
Scheme 1 applied to kinetics 1 with periodic boundary conditions over the square.

fe2d_nd.m, fe2d_nd.pdf
Scheme 2 applied to kinetics 1 with mixed Neumann/Dirichlet boundary conditions.

fe2dx_nd.m, fe2dx_nd.pdf
Scheme 1 applied to kinetics 1 with mixed Neumann/Dirichlet boundary conditions.

fe2d_nr.m, fe2d_nr.pdf
Scheme 2 applied to kinetics 1 with mixed Neumann/Robin boundary conditions.

fe2dx_nr.m, fe2dx_nr.pdf
Scheme 1 applied to kinetics 1 with mixed Neumann/Robin boundary conditions.

fe2dx_nr_alt.m, fe2dx_nr_alt.pdf
Same as fe2dx_nr, but using an implicit approximation of the boundary terms

mycopyright.txt,
FE2D copyright details.
Download the 'fast' codes for FE2D
These are the same as the codes listed above, but optimized for speed.
Files you may copy include the MATLAB Mfiles, the PDF versions with line numbering, and front ends
with test data for running the fast codes. See the comments in the fast Mfiles for further implementation details:

fe2d_n_fast.m, fe2d_n_fast.pdf
Scheme 2 applied to Kinetics 1 with pure Neumann boundary conditions.

fe2d_n_fast_test.m, fe2d_n_fast_test.pdf
A front end with some test data that can be used to run fe2d_n_fast.m.

fe2dx_n_fast.m, fe2dx_n_fast.pdf
Scheme 1 applied to Kinetics 1 with pure Neumann boundary conditions.

fe2dx_n_fast_test.m, fe2dx_n_fast_test.pdf
A front end with some test data that can be used to run fe2dx_n_fast.m.

fe2d_d_fast.m, fe2d_d_fast.pdf
Scheme 2 applied to kinetics 1 with pure Dirichlet boundary conditions.

fe2d_d_fast_test.m, fe2d_d_fast_test.pdf
A front end with some test data that can be used to run fe2d_d_fast.m.

fe2dx_d_fast.m, fe2dx_d_fast.pdf
Scheme 1 applied to kinetics 1 with pure Dirichlet boundary conditions.

fe2dx_d_fast_test.m, fe2dx_d_fast_test.pdf
A front end with some test data that can be used to run fe2dx_d_fast.m.

fe2d_r_fast.m, fe2d_r_fast.pdf
Scheme 2 applied to kinetics 1 with pure Robin boundary conditions.

fe2d_r_fast_test.m, fe2d_r_fast_test.pdf
A front end with some test data that can be used to run fe2d_r_fast.m.

fe2dx_r_fast.m, fe2dx_r_fast.pdf
Scheme 1 applied to kinetics 1 with pure Robin boundary conditions.

fe2dx_r_fast_test.m, fe2dx_r_fast_test.pdf
A front end with some test data that can be used to run fe2dx_r_fast.m.

fe2d_p_fast.m, fe2d_p_fast.pdf
Scheme 2 applied to kinetics 1 with Periodic boundary conditions.

fe2d_p_fast_test.m, fe2d_p_fast_test.pdf
A front end with some test data that can be used to run fe2d_p_fast.m.

fe2dx_p_fast.m, fe2dx_p_fast.pdf
Scheme 1 applied to kinetics 1 with Periodic boundary conditions.

fe2dx_p_fast_test.m, fe2dx_p_fast_test.pdf
A front end with some test data that can be used to run fe2dx_p_fast.m.

fe2d_nd_fast.m, fe2d_nd_fast.pdf
Scheme 2 applied to kinetics 1 with mixed Neumann/Dirichlet boundary conditions.

fe2d_nd_fast_test.m, fe2d_nd_fast_test.pdf
A front end with some test data that can be used to run fe2d_nd_fast.m.

fe2dx_nd_fast.m, fe2dx_nd_fast.pdf
Scheme 1 applied to kinetics 1 with mixed Neumann/Dirichlet boundary conditions.

fe2dx_nd_fast_test.m, fe2dx_nd_fast_test.pdf
A front end with some test data that can be used to run fe2dx_nd_fast.m.

fe2d_nr_fast.m, fe2d_nr_fast.pdf
Scheme 2 applied to kinetics 1 with mixed Neumann/Robin boundary conditions.

fe2d_nr_fast_test.m, fe2d_nr_fast_test.pdf
A front end with some test data that can be used to run fe2d_nr_fast.m.

fe2dx_nr_fast.m, fe2dx_nr_fast.pdf
Scheme 1 applied to kinetics 1 with mixed Neumann/Robin boundary conditions.

fe2dx_nr_fast_test.m, fe2dx_nr_fast_test.pdf
A front end with some test data that can be used to run fe2dx_nr_fast.m.

fe2dx_nr_alt_fast.m, fe2dx_nr_alt_fast.pdf
Same as fe2dx_nr_fast, but using an implicit approximation of the boundary terms

fe2dx_nr_alt_fast_test.m, fe2dx_nr_alt_fast_test.pdf
A front end with some test data that can be used to run fe2dx_nr_alt_fast.m.

mycopyright.txt,
FE2D copyright details.
Associated codes you will need
Files you may copy include:

boundedges.m,
Finds all the boundary edges in a triangular mesh.

subsetconnectivity.m
Finds the boundary edges in a triangular mesh for a portion of the boundary.

timestamp.m,
Prints the current YMDHMS date as a timestamp.
Boundary node lists for Experiments 2, 3, 4 and 6.
(To be uploaded soon  see reference 2. listed above.)
FE2D is distributed under the GNU GPL;
see the mycopyright.txt notice for more information.
Last revised on June 24, 2014.