Figure 1: Snapshot at T = 150
Numerical solution of spiral waves for phytoplankton on [0,400]^{2} in a
phytoplanktonzooplankton reactiondiffusion model.
The standard Galerkin FEM using continuous piecewise linear basis
functions was employed with zero flux boundary conditions. Solved on a uniform
401by401 grid with a fixed timestep of 1/384.
Figure 2: Snapshot at T = 1000
Numerical solution of spatiotemporal 'chaos' for phytoplankton on [0,400]^{2}
in a phytoplanktonzooplankton reactiondiffusion model. The standard Galerkin FEM
using continuous piecewise linear basis functions was employed with zero flux boundary conditions.
Solved on a uniform 401by401 grid with a fixed timestep of 1/384.
Figure 3: Snapshot at T = 110
Numerical solution of an expanding wavefront for prey on [0,400]^{2} in a
predatorprey reactiondiffusion model
after the localized introduction of predators
into a homogeneous distribution of prey. The standard Galerkin FEM
using continuous piecewise linear basis functions was employed with zero flux boundary conditions.
Solved on a uniform 401by401 grid with a fixed timestep of 1/384.
Figure 4: Triangulation of a lake with an island
Unstructured grid generation on a complex geometry
using DISTMESH, a mesh generator in Matlab (see SIAM Review
Vol. 46, No. 2, pp.329345, 2004). 281 nodes were used and
312 triangles were created.
Figure 5: Snapshot at T = 200
Numerical solution of an expanding wavefront for phytoplankton in a
phytoplanktonzooplankton reactiondiffusion model with complex
lake geometry (bounding box is [0,400]^{2}), after the localized
introduction of predators at (200,380) (marked with an 'X') into
a homogeneous distribution of prey. Associated
unstructured mesh has 2645 nodes and 4928 triangles. The standard
Galerkin FEM using continuous piecewise linear basis
functions was employed with zero flux boundary conditions.
Figure 6: Snapshot at T = 200
Numerical solution of spatiotemporal 'chaos' for phytoplankton
in a phytoplanktonzooplankton reactiondiffusion model with
complex lake geometry (bounding box is [0,400]^{2}),
after the localized introduction of predators at (250,200)
(marked with an 'X'),
into a homogeneous distribution of prey. Associated
unstructured mesh has 2645 nodes and 4928 triangles. The standard
Galerkin FEM using continuous piecewise linear basis
functions was employed with zero flux boundary conditions.
Figure 7: Interactions incorporated in a fishnutrientsplankton model
Arrows indicate positive effects, circles indicate negative effects (redrawn
and modified from OIKOS Vol. 62, pp.271282, 1991).

Mathematical and numerical analysis
I am interested in the interdisciplinary approach of
combining classical analysis, numerical analysis and scientific
computing to investigate nonlinear partial differential equations
(PDEs) in the applied sciences. My main areas of expertise are the rigorous
application of
the finite element method, and techniques from mathematical
analysis to prove the global wellposedness of nonlinear PDEs.
The methodology I use for the numerical analysis of PDEs
is to mimic the properties of the continuous system in the discrete
case. Thus mathematical analysis is used as a guide to practical
computation. On the applied side I am interested in mathematical
biology and nonlinear population dynamics in particular.
Mathematical biology
I am mainly concerned with PDEs that model applied problems,
particularly in the biological sciences. An
area that interests me is the application of mathematics to
complex spatial pattern formation phenomena, for example: the
development of mammalian coat patterns, epidemiology,
planktondynamics, and the development of scrollwaves in cardiac
tissue. Many pattern formation phenomena in biology can be modeled by
nonlinear systems of reactiondiffusion equations, and thus numerical
methods have an important part to play in finding practical
solutions.
Doctoral work
During my doctoral studies I analyzed a
system of coupled nonlinear reactiondiffusion equations of lambdaomega type
that were originally devised in order to understand pattern formation in the
BelousovZhabotinskii reaction (see for example
J.D. Murray's classic book 'Mathematical Biology'). My doctoral work consisted of three parts: establishing
the wellposedness of
the PDEs using the FaedoGalerkin method of J.L Lions (see also
James Robinson's book 'InfiniteDimensional Dynamical Systems'); numerical analysis using the 'lumped mass'
finite element
method with continuous piecewise linear basis functions; and scientific computing and
simulation of spiral waves and 'target patterns' with Matlab and Fortran.
My adviser was
Dr. James Blowey. Used together these techniques provide a powerful and unified framework for the analysis of any
nonlinear elliptic or parabolic system of PDEs in the applied sciences.
For further details see my papers.
Spatiallyextended predatorprey models
My recent work involved the mathematical and numerical study of computational algorithms (finite
element/difference) for nonlinear reactiondiffusion systems modeling
predatorprey interactions. I have been focusing on a welldocumented class of predatorprey models
with the Holling Type II functional response, which in the absense of predators has
logistic growth of the prey.
We have proved for the first time the rigorous wellposedness
of these nonlinear systems and rigorously proved stability and convergence results
for two fullypractical finite element approximations of these models.
For further details see my papers and for some free Matlab software
that I have written see my webpage below:
Our results cover an important example recently reviewed
in the context of marine plankton dynamics (SIAM Review, Vol.44, No.3, p.311).
The model reactiondiffusion system has the form
where 'u(x,t)' represents the phytoplankton density, 'v(x,t)' is the zooplankton
density, and the parameters 'alpha', 'beta', and 'gamma' are positive. In our
simulations we use the zeroflux boundary conditions, which reflects our
assumption that the plankton cannot leave the boundary of interest.
For various initial conditions, our computer simulations in 2D reveal that
the evolution of this system can lead to the formation of spiral patterns (see Figure 1),
followed by irregular 'patchy' structures in the whole domain, namely
spatiotemporal 'chaos' (see Figure 2). If you have an mpeg movie player you can see an
animation (23.4 MB file) of these dynamics for T=0
to T=1000. The movie is quite fast, so you may wish to use the options on your player to slow
the animation down. For information and sample source codes for making independently playable
movies with Matlab see
The results of the simulations have important implications for
understanding the role of interspecific interactions in
the observed patchy distribution of plankton in marine
environments. In terrestrial environments spatiotemporal patterns
resembling periodic traveling waves have been observed recently in several natural
populations, for example, field voles and red grouse. The reactiondiffusion systems
are useful models for investigating possible mechanisms for this behavior (see Figure 3).
A regularized spatiallyextended predatorprey system
Another classical predatorprey system involving space and time has the form
where 'p(x,t)' and 'h(x,t)' are population densities of predators and prey, and 'd_{1}' and
'd_{2}' are the diffusion coefficients. For appropriate choices of the positive parameters 'A',
'B', and 'C' the kinetics possess a limit cycle. From a rigorous mathematical point of view
these equations are challenging because of the singular term in the 2nd equation when 'h' is zero.
To overcome this difficulty we are undertaking the analysis and numerical analysis of the
associated regularised problem, where we modify the singular term via
for some arbitrary, but small epsilon.
Unstructured grid generation
Understanding the relationship between spatial patterns in population densities and environmental heterogeneity
is crucial to the understanding of population dynamics and for the management of species in communities.
My previous simulations of predatorprey interactions
were done on a square domain (habitat). This project
extends previous work as it involves generating unstructured grids on
arbitrary shaped habitat geometries, for
example, in the simulation of plankton dynamics in a lake (see Figure 4). This work is
in collaboration with a Research Associate, John Burkardt.
There are four main issues in this project. Firstly, there is the question of
what form the input data should take for the lake geometry, i.e.,
whether we use explicit boundary data or an implicit description (e.g., level
set). One possible solution we have used is to create a polygonal domain using the Linux drawing
package Inkscape. Secondly, there is the task of generating a 'good' mesh on the
computational domain. An example of a good mesh in 2D is one where all
triangles are approximately equilateral, and there are many element quality indicators for
measuring this. There are a large number of software packages available for unstructured grid generation in 2D.
One particularly simple set of routines in Matlab that we have used is called
Distmesh.
Once we have generated our mesh for the computational domain, we need to apply the Finite Element Method discussed
above, assemble the resulting linear system, and solve the linear system in an efficient manner. See below for some
webpages that I have written for the purpose of tackling these first three tasks:
Finally, there are ecological questions concerning how the
geometry of the habitat, including possible landscape features (e.g., islands),
influences the spatiotemporal dynamics of the solutions. Some preliminary results of our work
are shown in Figures 5 and 6 for a ficticious lake with an island. If you have an mpeg movie player you can see an
animation (21.8 MB file) of these dynamics for T=0
to T=1000.
Fishnutrientplankton dynamics
I am collaborating with a Postdoctoral Research
Associate,
Catalin Trenchea (Florida State University), on the optimal control of a system that is
related to the predatorprey models discussed above (a nonlinear reactiondiffusion system).
The predatorprey model is modified to account
for the predation rate 'f' (the control variable) of fish on
zooplankton 'v'. The kinetics are given by
where 'u' represents the phytoplankton and the parameters 'r','a','b', 'm',
and 'g' are positive (see Figure 7). The positive control variable 'f' is a function
of space and time and represents the biomass of zooplankton consumed at each
small volume of water per unit time (average fish predation rate x density of fish).
The model was originally
formulated as an ordinary differential system by Scheffer (OIKOS Vol. 62, pp.271282), and has
since been generalised to include diffusion. The basic aim here is to enhance zooplankton, by reducing
planktivorous fish, thereby, reducing algal biomass. In practice, planktivorous fish
can be reduced by fish removal, or by piscivore stocking.
This manipulation of the foodweb is called biomanipulation, and is an
important approach for improving water quality in eutrophic lakes.
From a mathematical point of view it is
important to choose a cost functional that can be analysed, but is
still relevant to the field situation. We consider the problem of manipulating the
densities of plankton by minimizing the quadratic functional
where the desired densities of phytoplankton and zooplankton are denoted with a 'bar' over u and
v, and Q represents the space  time domain of interest. In Figure 8 we show some numerical results
for the optimal control problem that illustrates the ability of the control to drive the system
from a chaotic regime to an ordered one (a rotating spiral wave).
Figure 8: Snapshot at T = 100
Approximate phytoplankton densities for the controlled and uncontrolled problem evolving
from an initial chaotic state, with the target and control 'f' (bounding boxes are [0,100]^{2}).
A uniform 101by101 grid was used with a constant timestep of 1/12.
The standard ('lumped mass') Galerkin FEM was employed to solve the state equations
and the adjoint equations, with zero flux boundary conditions. A 'steepest descent' algorithms
was employed to update the control.
An animation corresponding to Figure 8 from t = 0 to t = 100 is given
here (17.7 MB file) in the MP4 format.
