fd2d_predator_prey, an Octave code which uses the finite difference method (FDM) to simulate a predator-prey system in two space dimensions and time, by Marcus Garvie.
The MATLAB code is mostly self explanatory, with the names of variables and parameters corresponding to the symbols used in the finite difference methods described in the paper.
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 "for-loops" and conditional "if-then-else" statements, which again helps speed up the computations.
The linear systems are solved using MATLAB's built in function gmres(). The gmres() algorithm in MATLAB requires a number of input arguments. For our simulations we found it acceptable to use the default settings, e.g., no "restarts" of the iterative method or use of 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 and space-steps, and initial data. Due to a limitation in MATLAB, vector indices cannot be equal to zero; thus the nodal indices 0,...,J are shifted up one unit to give 1,...,(J+1) so xi=(i-1)*h + a.
It is easy to have MATLAB draw a contour plot of either the predator or prey density at any time. By making a plot at every timestep, an animation can be constructed, which gives an intriguing experience of how the dynamics of the system can support stable spiral patterns which eventually dissolve into a sort of chaos.
The initial data is again entered by the user as a string, but now the functions depend on a grid of x and y values denoted by X and Y (as a common MATLAB convention, a variable name using capital letters indicates that the variable is a matrix). Functions are entered using the same element-by-element rules for the arithmetic operation as in the 1-D code.
An example with an acceptable format is the following:
>> Enter initial prey 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=[-500,500]2, in order to choose an initial
prey density of 0.2 within the circle with radius 10 and center
(-50,200), and a density of 0.01 elsewhere on Omega
we input the following:
>> Enter initial predator function V0(X,Y)
0.2*(((X+50).^2+(Y-200).^2)<100)+0.01*(((X+50).^2+(Y-200).^2)>=100)
The computer code and data files described and made available on this web page are distributed under the MIT license
fd2d_predator_prey is available in a Fortran90 version and a MATLAB version and an Octave version.
fd_predator_prey, an Octave code which solves a pair of predator prey ODE's using a finite difference approximation.
fd1d_predator_prey, an Octave code which uses finite differences to solve a 1D predator prey problem.
Marcus Garvie