Source: https://file.scirp.org/Html/5-1100604_83079.htm
Timestamp: 2019-04-21 16:56:07+00:00

Document:
Reaction-diffusion equations modeling Predator-Prey interaction are of current interest. Standard approaches such as first-order (in time) finite difference schemes for approximating the solution are widely spread. Though, this paper shows that recent advance methods can be more favored. In this work, we have incorporated, throughout numerical comparison experiments, spectral methods, for the space discretization, in conjunction with second and fourth-order time integrating methods for approximating the solution of the reaction-diffusion differential equations. The results have revealed that these methods have advantages over the conventional methods, some of which to mention are: the ease of implementation, accuracy and CPU time.
Numerical methods are important tools in investigating the solution’s behavior of non-linear realistic models in biology  where no closed form solutions exist. A class of these models is reaction diffusion (RD) problems that, for instant, reproduce some of the complex pattern observed on the skin of certain animals. This includes the development of coat patterns on mammals and the patterning of butterfly wings  .
Reaction diffusion models have been studied extensively since the RD theory first proposed by Turing  to describe the range of spatial patterns observed in the developing embryo.
In recent years, several theoretical models, regarding spatial pattern, have been introduced and led to a better understanding of how these patterns arise from certain mechanism. A general theoretical framework for studying pattern formation in biological systems within a growing domain was developed in  . In addition, reaction diffusion equations modeling predator-prey interactions have been paid attention and studied extensively    . In these models, the most important element is the “ functional response”, the function that describes the number of prey consumed by predator per unit time.
For the numerical solution of nonlinear reaction diffusion problems (most of which are stiff1), fully explicit temporal methods are avoided due to the sever restriction on the time step size imposed by the stiff diffusion term. Fully implicit temporal schemes are not recommended since they have the task of solving large implicit systems at each time step. Alternatively, several methods were introduced. The vast majority of scientific community uses, for the spatial discretization, lower-order centred finite differences schemes due to the ease of implementation and incorporation of the boundary conditions in a straightforward way. The time integration is carried out, usually, by using some types of low-order time step methods such as, second-order fully implicit Crank-Nicholson methods.
The Exponential Time Differencing (ETD) methods  in conjunction with spectral methods    have emerged as viable alternatives to classical schemes for a wide variety of problems. Their applications to solve reaction-diffusion problems      have been steadily growing. Spectral methods are popular for generating spectral accurate spatial derivatives, through built-in codes in Matlab  . ETD schemes recover the exact solution to the linear part (which is generally the most stiff part of the system), and integrate exactly an approximation of the non-linear terms. We may approximate the non-linear parts by some polynomial in time that may be calculated using previous steps of the integration process (producing multi-step ETD methods) or by Runge-Kutta-like stages (resulting in ETD schemes of RK type), see  for a comprehensive review. The numerical computation of the ETD coefficients, which are functions related to the matrix exponential, arises some difficulties   . Fortunately, several numerical algorithms   -  were introduced to enhance the accuracy in computing these coefficients.
The work presented in this article is motivated by the theoretical and computational study by Gavie  for the dynamical properties of the 2-component reaction-diffusion system modeling predator-prey interactions with the Holling type II functional response. The system displays a wide spectrum of ecologically relevant behavior, including chaos. The model and an illustration of the computational study conducted by Garvie  , using first-order finite difference schemes, are briefly described in §2. In §3, we numerically approximate the solution of the predator-prey model in two dimensions employing spectral methods, for the space discretization, in conjunction with high-order time integrating methods. The results of the numerical experiment are compared to Garvie’s findings. Section §4 recaps our results and points out some conclusions.
• h ( ⋅ ) is strictly increasing on [ 0, ∞ ) .
h ( η ) = η 1 + η , ( η = a u ) , with a = 1 / α , b = β , c = γ .
We consider the natural biological meaningful region u ≥ 0, v ≥ 0 with bounded positive initial data, where the solution remains positive at all time. From an environmental point of view, the choice of parameters for the numerical simulation of the full reaction-diffusion system insure that both the spatial and temporal local dynamics (the densities) of the predator and preys are oscillatory. It is important to note that the above model takes into account the invasion of the prey species by predators but does not include stochastic effects or any influences from the environment.
where Δ t = T / N and h = ( B − A ) / J represent the time and space step respectively. In addition, suppose that U i , j n and V i , j n denote the two dimensional approximation to the solutions u and v respectively at the point ( x i , y j , t n ) , where the time levels t n = n Δ t .
and zero-flux boundary conditions, for studying the dynamics of spatially extended predator-prey interaction with logistic growth of the prey (1) in two space dimensions. The “analytical solution” of the predator-prey system (1) is approximated using schemes (3) and (4) on a fine mesh and small time step. The linear system of algebraic equations, resulting from the solution of the problem using either schemes, is sparse, banded and has a simple structure (tridiagonal and block tridiagonal system in one and two dimensions respectively). In addition, the coefficient matrices of the linear system are strictly diagonally dominant, which guarantee the convergence of the GMRES algorithm (an iterative solver) used to solve the system (in two dimensions), see  .
The convergence of Schemes (3) and (4) was decided on upon reducing the time step until the difference between the approximations from either schemes is negligible. On the other hand, the space step is kept sufficiently small to see clearly the qualitative feature of the solutions.
• For various initial conditions, the evolution of the system led to the formulation of spiral patterns, see Figure 1, followed by irregular patches covering the whole domain (spatio-temporal chaos).
• Schemes (3) and (4) are simple to program, stable and convergent provided that the time step is below a (non-restrictive) critical value.
• The disadvantage of the code is that the run time can be prohibitive when using a combination of large domain size and final time T coupled with small space and time step.
From an ecological point of view, Garvie stated that in the absence of external influences, certain initial conditions can lead to spatial and temporal variants in the densities of predator and prey that persist indefinitely.
The section following gives an overview of higher-order competing methods: Exponential Time Differencing (ETD) methods §3.1, Integrating Factor (IF) methods §3.2 and Implicit-Explicit methods (IMEX) §3.3 utilized for simulating the numerical solution of the prey-predator system (1) with kinetics (2) in two dimensions.
Figure 1. The approximated prey densities U i , j n for kinetics (2) in two dimensions at T = 150 . The initial conditions are given by U i , j 0 = 6 / 35 − 2 × 10 − 7 ( x i − 0.1 y j − 225 ) ( x i − 0.1 y j − 675 ) , V i , j 0 = 116 / 245 − 3 × 10 − 5 ( x i − 450 ) − 1.2 × 10 − 4 ( y j − 150 ) with α = 0.4 , β = 2.0 , γ = 0.6 , δ = 1 , h = 1 and Δ t = 1 / 384 .
§3.1, Integrating Factor (IF) methods §3.2 and Implicit-Explicit methods (IMEX) §3.3, for approximating the solution of the prey-predator system (1) with kinetics (2) in two dimensions. The aim is to observe the effectiveness of the competing methods, taking into account the accuracy and CPU time consumed by the methods.
where Δ t is the time step and u n and F n denote the numerical approximation to u ( t n ) and F ( u ( t n ) , t n ) respectively.
The terms a n and b n approximate the values of u at t n + Δ t / 2 and the term c n approximates the value of u at t n + Δ t .
u n + 1 = ( u n + Δ t F n ) e c Δ t .
For the time discretization, this method treats the non-linear reaction term explicitly utilizing the second-order Adams-Moulton methods, while the diffusion term is treated implicitly utilizing the second-order Adams-Bashforth schemes.
IMEX methods are restricted from having an order higher than two if A-stability is required. Therefore, despite their simplicity and frequent usage, they are not extendable to higher order.
The overall efficiency of the methods in our numerical experiments is measured by three dominant testing parameters: the accuracy, the start-up overhead cost and the CPU time consumed by the methods. All the calculations presented in this paper are performed using Matlab codes.
V i , j 0 = 116 / 245 − 3 × 10 − 5 ( x i − 450 ) − 1.2 × 10 − 4 ( y j − 150 ) .
Considering the implementation of the above time discretization methods, we evaluate three coefficients2 for the ETD2RK1(6), the ETD2RK2 (7) methods and four coefficients for the ETD4RK (8) method, once at the beginning of the integration for each value of the time-step sizes, by means of contour integration3 in the complex plane approximated by the Trapezium rule. We choose circular contours, each is centred at one of the elements that are on the diagonal matrix of the linear part of the semi-discretized model (12), with radius R = 1 and sampled at 32 equally spaced points. Moreover, in the main loop of integration, the ETD2RK1, the ETD2RK2, the IFRK2 (9), the ETD4RK and the IFRK4 (10) methods perform two (for the second-order methods) and four (for the fourth-order methods) function transforms per time step. The IF schemes require, in addition, the evaluation of one or more matrix exponentials, for which acceptable algorithms are well known   .
In our tests, we measure the accuracy of the methods in terms of the relative error evaluated in the integrated error norm. The numerical “exact” solution is approximated on a fine mesh with a very small time-step size and compared with the corresponding approximated solution with a sequence of large time steps.
Figure 2. Relative errors versus (a) time step (b) CPU time for the prey-predator system (1) with kinetics (2) in two dimensions. The initial conditions are given by U i , j 0 = 6 / 35 − 2 × 10 − 7 ( x i − 0.1 y j − 225 ) ( x i − 0.1 y j − 675 ) , V i , j 0 = 116 / 245 − 3 × 10 − 5 ( x i − 450 ) − 1.2 × 10 − 4 ( y j − 150 ) with α = 0.4 , β = 2.0 , γ = 0.6 and δ = 1 .
and confirms the expected order for all of them. In the figure, we plot the numerical relative error of the integrated error norm as a function of the time step. The time-step values are selected to ensure that all methods achieve stable accurate results. The plot indicates, firstly, that the fourth-order methods take the largest time-step size, i.e. the fewest number of steps, to converge to a solution within a fixed given relative error in the figure. On the other hand, the first-order schemes (3) and (4) break down (errors are of O ( 1 ) ) for time-step size larger than Δ t ≈ 2 − 4 because of the numerical stability constraints. Secondly, a fixed reduction in the time-step size will produce a greater reduction in the error for the second and fourth order methods than for the first-order schemes (3) and (4).
For an additional preference between the methods compared, we valued the accuracy with respect to the CPU time, illustrated in Figure 2(b). For the first-order schemes (3) and (4), the computation cost is almost identical, per time step.
Second-order convergence is confirmed in Figure 2(a) for the ETD2RK1(6), the ETD2RK2 (7), the AB2AM2 (11) and the IFRK2 (9) methods. All second-order methods successfully integrate the system for time-step sizes Δ t ≈ 2 − 2 and less. In addition, the corresponding errors for all methods (except for that of the AB2AM2 method which is considered to be the least accurate) lie approximately on the same line for all values of the time-step. Regarding the CPU time, Figure 2(b) indicates that the variation in time consumption for all methods, for a given level of accuracy, is insignificant.
For the fourth-order methods, the performance of the IFRK4 (10) method resembles that of the ETD4RK (8) method, and the errors are identical and all scale with the expected ( Δ t ) 4 order, see Figure 2(a). However, according to Figure 2(b), the IFRK4 method uses a slightly longer computation time than the ETD4RK method for a given error tolerance.
Clearly the fourth-order methods have a superior performance, as the accuracy and speed are improved significantly compared to lower order methods referring to Figure 2(a) and Figure 2(b) respectively. The fourth order methods have the advantages of gaining order of magnitude of accuracy, even for larger time steps Δ t ≈ 2 − 1 , greater than lower order methods.
In general, we find that out of all comparable methods, the fourth-order methods are the most accurate, for a given time-step size, and the least time consuming for a given level of accuracy, see Figure 2(a) and Figure 2(b) respectively. However, the ETD4RK method was found to be the best scheme for providing accurate stable solution in a fast and efficient way. On the other hand, the most expensive with high computational cost and unpleasant performance are schemes (3) and (4). Thus, the greater accuracy of the ETD4RK and the IFRK4 methods rewards the extra programming efforts.
In most of the cases in using numerical computations of reaction-diffusion problems, a simple first-order semi-implicit scheme is the first choice  due to the difficulties introduced by the combination of non-linearity and stiffness of a PDE, the complexity both of analysis and implementation of higher order methods, and the increased computer storage required by these methods. However, in tow and higher space dimensions problems, simulations based upon the more conventional ideas become more time consuming. Higher order time integration methods can result in significant benefits in run time reduction while maintaining high accuracy.
In this paper, spectral methods, coupled with second and fourth order exponential integrator methods, are used for solving a reaction-diffusion problem modeling predator-prey interaction. The spectral methods offer several advantages over finite-difference methods such as, accurate discretization of spatial derivatives, ease of implementation and applicability to a wide varies of equations    .
Exponential time stepping solvers are used to remove the the stiffness associated with the diffusive terms, and hence, the severe restriction on the time-step for reasons of stability, allowing much larger time-steps to be selected and the selection is only limited by accuracy.
In our comparison experiment, we have shown that the fourth-order ETD4RK (8) and IFRK4 (10) methods have aided the simulation of model (1) in tow dimensions with improvements of the computational efficiency and speed over first and second-order methods.
Although we focus on system (1) in this paper, we believe that the numerical methods described here and our results can also be extended with modification to numerically solve other non-liner wave problems, models which include advection and reaction-diffusion-convection systems.
1. Chauhan, A., Kaith, B.S., Singha, A.S. and Pathania, D. (2010) Induction of the Morphological Changes in Hibiscus sabdariffa on Graft Copolymerization with Acrylonitrile and Co-Vinyl Monomer in Binary Mixtures. Malaysian Polymer Journal, 52, 140-150.
3. Kaith, B.S., Singha, A.S., Dwivedi, D.K., Kumar, S., Kumar, D. and Dhemeniya, D. (2003) Preparation of Polystyrene Matrix Based Composites Using Flax-Gcopolymers as Reinforcing Agent and Evaluation of the Mechanical Behavior. International Journal of Plastic Technology, 7, 119-125.
4. Kaith, B.S., Singha, A.S. and Sharma, S.K. (2004) Synthesis of Graft Copolymers of Binary Vinyl Monomer Mixtures and Flax Fiber Using FAS-KPS Redox System. International Journal of Chemical Science, 2, 37-43.
5. Kaith, B.S., Singha, A.S. and Sharma, S.K. (2003) Graft Copolymerization of Flax Fibers with Binary Vinyl Monomer Mixtures and Evaluation of Swelling, Moisture Absorbance and Thermal Behavior of the Grafted Fibers. Journal of Polymeric Materials, 20, 195.
6. Singha, A.S., Kaith, B.S. and Kumar, S. (2004) Evaluation of Physical and Chemical Properties of FAS-KPS Induced Graft Co-Polymerization of Binary Vinyl Monomer Mixtures onto Mercerized Flax. International Journal of Chemical Sciences, 2, 472.
10. Kaith, B.S., singha, A.S. and Kalia, S. (2007) Grafting MM onto Flax under the Influence of Microwave Radiation and the Use of Flax-g-poly(MMA) in Preparing PF Composites. Autex Research Journal, 7, 119-129.
17. Kaith, B.S. and Kalia, S. (2007) Synthesis and Characterization of Graft Co-Polymers of Flax Fiber with Binary Vinyl Monomers. International Journal of Polymer Analysis and Characterization, 12, 401-412.
25. Murray, J.D. (1993) Mathematical Biology. 2nd Edition, Springer-Verlag, Berlin Heidelberg.
37. Boyd, J.P. (2001) Chebyshev and Fourier Spectral Methods. 2nd Edition, Dover, New York.
40. Craster, R.V. and Sassi, R. (2006) Spectral Algorithm for Reaction-Diffusion Equations. Note del Polo 99, Universita degli Studi di Milano, Polp Didattico e di Ricerca di Crema.
42. Kassam, A.K. (2003) Solving Reaction-Diffusion Equations 10 Times Faster. Oxford University, Oxford, Numerical Analysis Group Research Report No. 16.
45. Minchev, B.V. and Wright, W.M. (2005) A Review of Exponential Integrators for First Order Semi-Linear Problems. Tech. Rep. NTNU.
47. Ashi, H.A. (2008) Numerical Methods for Stiff Systems. PhD Thesis, University of Nottingham, Nottingham.
49. Kassam, A.K. (2004) High Order Time Stepping for Stiff Semi-Linear Partial Differential Equations. PhD Thesis, Oxford University, Oxford.
50. Schmelzer, T. and Trefethen, L.N. (2007) Evaluating Matrix Functions for Exponential Integrators via Carathéodory-Fejér Approximation and Contour Integrals. Electronic Transactions on Numerical Analysis, 29, 1-18.
54. Berland, H., Skeflestad, B. and Wright, W.M. (2007) EXPINT—A Matlab Package for Exponential Integrators. ACM Transactions on Mathematical Software, 33, Article No. 4.
57. Ascher, U.M., Ruuth, S.J. and Wetton, B.T. (1995) Implicit-Explicit Methods for Time-Dependent Partial Differential Equations. SIAM Journal on Numerical Analysis, 32, 797-823.
1Stiff differential equations are categorized as those whose solutions (or different components of a single solution) evolve on very different time scales occurring simultaneously, i.e. the rates of change of the various components of the solutions differ markedly. For a comprehensive review of this phenomena see      .
2The ETD-RK methods require an accurate algorithm for evaluating the coefficients of F ( u ( t n ) , t n ) to avoid numerical difficulties, see   .
3The “Cauchy integral” approach was proposed by Kassam and Trefethen, see   for further details.
where α = 0.4 , β = 2 , γ = 0.6 , δ = 1 and k x , k y represent the wave-numbers in two dimensions. The linear (diffusion) operator in the resulting uncoupled system of ordinary differential equations (ODEs) (12) is diagonal, with elements − ( k x 2 + k y 2 ) , of which some has large negative real eigenvalues that represent decay on a time scale much shorter than that typical of the non-linear term (strong dissipation dynamics), causing system (12) to be stiff. The non-linear term is transformed to physical space and evaluated at the uniform grid points and then transformed back to spectral space.

References: §2
 §3
 §4
 §3
 §3
 §3

§3
 §3
 §3