FE2D is a collection of MATLAB routines using the finite element method for simulating the dynamics of predator-prey interactions in two space dimensions and time. The codes generalize my earlier predator-prey 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, 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 link 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
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 many simulations we found it acceptable to use no "restarts" of the iterative method, or, preconditioners, and a tolerance for the relative error of 1E-6. 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, time-step, initial data functions, and boundary data functions.
The reaction-diffusion system models spatially extended predator-prey interactions of
Rosenzweig-MacArthur form, namely
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 clock-wise or counter-clockwise, 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 = [ 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 ]; |
For convenience we provide an example of a 'front-end' 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 grid-nodes 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 front-end that can be copied and adapted are given below:
The codes are mostly self-explanatory, 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 line-numbered codes are listed at the bottom of the page in PDF format.
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*((x-50)^2+(y-50)^2<25)
The boundary of the whole domain is denoted Gamma. In the mixed boundary condition cases we assume
Gamma is the union of two pieces Gamma 1 and Gamma2. The code allows for the following boundary
conditions:
>> Enter the Dirichlet b.c. g1u(x,y,t) for u 0.2*sin(x-y)+t/10+0.2
>> Enter the Dirichlet b.c. g1v(x,y,t) for v 0.3*cos(x-y)+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
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 run-time 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 time-step 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 ill-conditioned 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).
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 LUINC. 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 Gauss-Seidel algorithms are widely used, although the run-times 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 aii of the coefficient matrix A are non-zero, and they will not always converge even then.
Files you may copy include:
Files you may copy include the MATLAB M-files and PDF versions with line numbering:
Files you may copy include:
FE2D is distributed under the GNU GPL; see the License and Copyright notice for more information.