Dr. Marcus R Garvie

Research:



    Figure 1: Snapshot at T = 150

    Numerical solution of spiral waves for phytoplankton on [0,400]2 in a phytoplankton-zooplankton reaction-diffusion model. The standard Galerkin FEM using continuous piecewise linear basis functions was employed with zero flux boundary conditions. Solved on a uniform 401-by-401 grid with a fixed time-step of 1/384.



    Figure 2: Snapshot at T = 1000

    Numerical solution of spatiotemporal 'chaos' for phytoplankton on [0,400]2 in a phytoplankton-zooplankton reaction-diffusion model. The standard Galerkin FEM using continuous piecewise linear basis functions was employed with zero flux boundary conditions. Solved on a uniform 401-by-401 grid with a fixed time-step of 1/384.



    Figure 3: Snapshot at T = 110

    Numerical solution of an expanding wavefront for prey on [0,400]2 in a predator-prey reaction-diffusion 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 401-by-401 grid with a fixed time-step 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.329-345, 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 phytoplankton-zooplankton reaction-diffusion 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 phytoplankton-zooplankton reaction-diffusion 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 fish-nutrients-plankton model

    Arrows indicate positive effects, circles indicate negative effects (redrawn and modified from OIKOS Vol. 62, pp.271-282, 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 well-posedness 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, plankton-dynamics, and the development of scroll-waves in cardiac tissue. Many pattern formation phenomena in biology can be modeled by nonlinear systems of reaction-diffusion 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 reaction-diffusion equations of lambda-omega type that were originally devised in order to understand pattern formation in the Belousov-Zhabotinskii reaction (see for example J.D. Murray's classic book 'Mathematical Biology'). My doctoral work consisted of three parts: establishing the well-posedness of the PDEs using the Faedo-Galerkin method of J.L Lions (see also James Robinson's book 'Infinite-Dimensional 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.

Spatially-extended predator-prey models

My recent work involved the mathematical and numerical study of computational algorithms (finite element/difference) for nonlinear reaction-diffusion systems modeling predator-prey interactions. I have been focusing on a well-documented class of predator-prey 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 well-posedness of these nonlinear systems and rigorously proved stability and convergence results for two fully-practical finite element approximations of these models. For further details see my papers and for some free Matlab software that I have written see my web-page 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 reaction-diffusion 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 zero-flux boundary conditions, which reflects our assumption that the plankton cannot leave the boundary of interest. For various initial conditions, our computer simulations in 2-D 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 reaction-diffusion systems are useful models for investigating possible mechanisms for this behavior (see Figure 3).

A regularized spatially-extended predator-prey system

Another classical predator-prey system involving space and time has the form



where 'p(x,t)' and 'h(x,t)' are population densities of predators and prey, and 'd1' and 'd2' 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 predator-prey 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 web-pages 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.

Fish-nutrient-plankton 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 predator-prey models discussed above (a nonlinear reaction-diffusion system). The predator-prey 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.271-282), 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 food-web 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 101-by-101 grid was used with a constant time-step 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.

   
   
  Home    (Best viewed with a monitor resolution of 1280x960)

Last updated: August 15, 2005 by Marcus R Garvie