POSITIVE SYMPLECTIC INTEGRATORS FOR PREDATOR-PREY DYNAMICS

. We propose novel positive numerical integrators for approximating predator-prey models. The schemes are based on suitable symplectic proce- dures applied to the dynamical system written in terms of the log transformation of the original variables. Even if this approach is not new when dealing with Hamiltonian systems, it is of particular interest in population dynamics since the positivity of the approximation is ensured without any restriction on the temporal step size. When applied to separable M -systems, the resulting schemes are proved to be explicit, positive, Poisson maps. The approach is generalized to predator-prey dynamics which do not exhibit an M -system structure and successively to reaction-diﬀusion equations describing spatially extended dynamics. A classical polynomial Krylov approximation for the diﬀusive term joint with the proposed schemes for the reaction, allows us to propose numerical schemes which are explicit when applied to well established ecological models for predator-prey dynamics. Numerical simulations show that the considered approach provides results which outperform the numerical approximations found in recent literature.


1.
Introduction. Non-linear models for studying the predator-prey populations are sufficiently rich in their spatiotemporal dynamics to allow for the description of phenomena like species invasion, species persistence, and biological chaos [18]. As pointed out in [10,11,12], the numerical approximation of such models is a challenging aspect of the problem that deserves great attention, since whenever the theoretical solution is not available, biological conclusions rely on the obtained numerical results. Generally speaking, numerical methods that are able to mimic specific properties of the continuous models they are approximating, should be preferred to general purpose algorithms since they show, in the long run, a gain in accuracy [13].
In the recent literature, the attention on preserving positivity throughout the numerical approximation of dynamical system solutions is increased [14,1,19]. In computational ecology, where the variables may represent populations, densities, concentrations, positivity is an important feature to be preserved. However, the preservation of positivity is not at all a trivial task for standard numerical methods. Usually, it implies a severe step size restriction [1]. An interesting solution of the problem is given by splitting methods which are positive integrators independently of the time step size, provided the defining semiflows are positive. Building on that, a second-order positivity-preserving scheme for semilinear parabolic problems has been proposed in [14]. Positivity-preserving nonstandard finite difference schemes for advection-diffusion reaction equations have been proposed in [19].
In the context of population dynamics, there exists a class of systems called Msystems [17] (with the Lotka-Volterra (LV) model providing a classical example), for which population dynamics can be described by a phase space theory analogous to the Hamiltonian formulation of classical mechanics. For an effective long term approximation of the solution of this class of systems, discrete Poisson maps have to be used (see [13]).
In this paper we exploit the analogy between M -systems and Poisson systems in order to consider numerical schemes that are both Poisson and positive maps regardless the chosen temporal step size. The rationale behind the generation of these schemes is to transform the original Poisson system into the equivalent Hamiltonian system obtained via the log transform. Then, the symplectic Euler scheme is applied to the transformed system and finally, the approximated values in terms of the original variables are obtained by inverse transforming the solution via the exp transform. The approach, already reported in the Hamiltonian dynamic literature (see, for example, [2,13]), is analyzed here in the context of population dynamics: in the case of separable M -systems (e.g. LV model) we will show that, unlike methods based on composition of symplectic Euler schemes, which are implicit, we get a Poisson structure for the numerical map as well as the positivity of the approximated solution and its explicit functional form. Moreover we will show that, although the linear stability analysis for the first-order scheme, as applied to the LV system, imposes a restriction on the value of the temporal step size, this constraint is indeed less demanding than requiring the positivity of the symplectic Euler scheme.
We expand the approach to approximate general predator-prey dynamics and we analyze cases when explicit integrators can be obtained. We test the method on models of Z-controlled LV dynamics described in [16]. We show that, differently from the classical Euler explicit scheme, our approach retains the positivity even in the transient phase.
Furthermore, we consider spatially explicit population dynamics described by a class of reaction-diffusion systems. We compose positive schemes for the approximation of the reaction term with a classical polynomial Krylov approach, which consists in approximating the action of the exponential of the discretized Laplace operator on a vector. To test the approach on models of ecological interest, we perform numerical simulation on the phytoplankton-zoo-plankton dynamics provided in [18]. This example is classical in the literature in that, for sufficiently large perturbations of homogeneous equilibria, it exhibits a chaotic behavior. The numerical results obtained with our approach have been compared to the ones obtained by the authors with IMEX (implicit-explicit) [15,10,12] and IMSP (implicit-symplectic) methods [6,7]. We show that, even for large time step sizes, the novel proposed method, unlike the ones mentioned above, is not affected at all by numerical artifacts.
The paper is organized as follows. We recall the definition and the Poisson structure of M -systems in Section 2 providing some ecological examples. In Section 3 we analyze and compare some Poisson integrators for M -systems in terms of preservation of positivity and their explicit/implicit form; in the case of LV systems, the comparison is performed also by means of a linear stability analysis. In Section 4 we extend our approach to general predator-prey dynamics: the case of dynamics not depending on space is analyzed in Subsection 4.1 while the case of reactiondiffusion equations is considered in Subsection 4.2. Finally, conclusions are drawn in Section 5.
2. M -systems. A class of ecological models is described by M -systems, which are defined in [17].
1. An M -system is a dynamical system defined by a given resource function M (u, v; t) and equations of motion given bẏ Let us illustrate three simple examples of M -systems in population dynamics. The first is the Lotka-Volterra model which represents the first attempt at simulating the observed oscillating behavior of species abundances through predator-prey interactions. For dimensionless variables it can be written aṡ It is easy to check that Eq.(2) is an example of M -systems with resource function The simplest example of an M -system is a simplified description of a direct disease transmission in a total population of M = u + v individuals where v is the number of infected individuals and u is the number of individuals who are susceptible to infection. It may be obtained from (2) by setting a = b = 0 : Notice that both (2) and (3) are examples of autonomous separable M -systems as they do not explicitly depend on time and they are separable in two parts: one depending only on u and the other depending only on v. An example of still autonomous but not separable M -systems, as suggested in [17], is described by the equationsu The second equation defines the logistic evolution, where k is the growth rate for the population v; µ v represents the density dependent mortality rate. In the first equation u s is a saturation constant. The corresponding M -function is defined by ∂b ij (y) ∂y l b lk (y) + ∂b jk (y) ∂y l b li (y) + ∂b ki (y) ∂y l b lj (y) = 0 for all i, j, k.
By setting y = (u, v), it results that the resource function M (y; t), which defines an M -system in Def. 2.1, may be seen as the Hamiltonian function of the Poisson system corresponding to the skew-symmetric matrix From the previous result many other known properties may be deduced.   In case of separable M -systems, the symplectic Euler (SE) scheme and its explicit variant (EVSE) provide two examples of first-order Poisson integrators [1].
The preservation of the Poisson structure is not the solely requirement for a numerical integrator of a M -system which describes population dynamics. Indeed, in order to have physically meaningful approximated solutions, the phase space R 2 + should be preserved too. However, in general, Poisson integrators do not ensure this property. For example, the approximated solutions provided by the symplectic Euler method applied to the LV model starting from initial conditions u 0 , v 0 > 0, are not guaranteed to be positive at each step. If we choose ∆t carefully i.e., ∆t < min 1 a , 1 b the numerical solutions u n and v n stay positive for all n ∈ N (see [1]). The question we are going to answer is the following: how to construct Poisson integrators of arbitrary order which preserve both the Poisson structure of the (possibly non-separable) is an Hamiltonian system. Thus, the above change of variables allows us to construct an associated Hamiltonian system starting from an M -system evolving in the phase space R 2 + = (u, v ∈ R 2 , u, v ≥ 0). Conversely, the relations comprise the Hamilton phase space R 2 into the first octant R 2 + . From the previous relation we get a hint for the generation of numerical schemes i.e., to transform first the original Poisson system into the equivalent Hamiltonian system obtained via log function. Then, apply symplectic schemes to the transformed system; finally obtain the approximated values in terms of the original variables by inverse transforming the solution via the exp transform. The simplest first-order method is constructed by applying symplectic Euler on the transformed Hamiltonian system. It advances in time according to: The adjoint discrete flow Φ * ∆t ( [13]) defines also a first-order procedure that advances in time according to the following steps Notice that scheme (7), when applied to the LV model, is not new in the Hamiltonian dynamics literature (it may be seen, for example, as a first-order splitting procedure [2]); however, for more general population dynamics it is an original result which deserves a great attention because of its nice qualitative properties. The following theorem holds.
Theorem 3.2. The first-order method (7) (resp. (8)) is a Poisson integrator for M- . It provides positive solutions for all ∆t > 0 whenever v 0 , u 0 are positive. In case of separable M-systems it is also explicit. 1 Proof. The positivity of the solution follows by construction. It is also trivial to check that, in case of separable M-systems i.e., whenever M (u, v) = E(u) + F (v), the resulting scheme (7) (and (8), respectively) has an explicit functional form. It remains to prove that, by setting We differentiate (7) with respect to (u n , v n ) and write the results as a matrix equation: , v n ) and use the same notation for the second order partial derivatives. Under the hypothesis that 1 + ∆t u n+1 v n M uv = 0 we can invert the first matrix in (11) and evaluate the entries of the Jacobian matrix A as , Once the matrix A has been obtained, the proof follows by the straightforward calculation of the matrix product in (10). A similar reasoning can be followed to   (7) and (8)  3.2. Higher order schemes. By symmetric composition of the discrete flows we can increase the accuracy of the proposed schemes [13]. Since the composition of a Poisson scheme is itself a Poisson integrator, by composing the steps of the firstorder Poisson integrator Φ ∆t (7) with its adjoint Φ * ∆t we may construct positive Poisson integrators of arbitrary order of accuracy: As an example, in Figure 1, on the right, we show the numerical accuracy of different methods applied to (2) with a = 0.5, and b = 0.1, starting from u 0 = v 0 = 0.2 in the temporal horizon [0, 50]. As can be appreciated, the methods based on Strang composition (s = 1, β 1 = α 1 = 1/2) and Yoshida composition (s = 3, have numerical orders which reflect the theoretical second and fourth-order ones, respectively. We notice that, among the first-order integrators, the method (7) and its adjoint have the best performance when compared to the symplectic Euler method and its explicit variant.
3.3. Lotka-Volterra dynamics. A linear stability analysis. Let us compare first-order procedures for solving an example of the well-known Lotka-Volterra predator-prey dynamics. In Figure 1, on the left, we show the results of the methods (7) and (8) at T = 8.3 applied to (2) with ∆t = 1.1. The value of ∆t is greater then the value which ensures positive solutions for the symplectic Euler method (6). Notice that, unlike the symplectic Euler (as well as for its explicit variant), the two proposed first-order schemes (7) and (8)  As shown, no restriction on the temporal step size is required to enforce the positivity of the solution of the proposed methods. However, restrictions may arise from the stability requirements. Let us focus on the LV model (2) and perform a linear stability analysis of the first-order integrator (7). Firstly, we recall the stability features of the continuous system (2). The model has two steady states, the trivial steady state (u * 1 , v * 1 ) = (0, 0) and the internal steady state The Jacobian matrix evaluated at the trivial state, has two real eigenvalues a and −b, one positive and one negative, with this revealing that the trivial steady state is stable against perturbations in the predator density v, and unstable against perturbations in the prey density u. It is therefore a saddle point for all the values of the parameters. The Jacobian matrix (14) evaluated at (13) has a pair of pure imaginary eigenvalues ±i √ a b. Since in the (u, v)-phase plane the solutions describe elliptic curves through the initial state y 0 = (u 0 , v 0 ) which obey the invariance property of the Hamiltonian i.e., H(y(t)) = H(y 0 ) for all t > 0, the internal steady state is a neutral center and all the solution orbits are neutrally stable.
We study now the linear stability of the map (7). Taking into account that M uv = 0 the Jacobian matrix (9) can be rewritten in terms of solely (u n , v n ): At the trivial steady state A is diagonal with entries e ∆t a and e −∆t b . It is easy to see that 0 < e −∆t b < 1 < e ∆t a : the trivial steady state is a saddle point, attractive along v and repulsive along u, as it is in the analytical model, without restriction on the step size (as it is the case of the symplectic Euler method, see [1]).
The evaluation of the Jacobian A at the steady state (13) leads to the following In order to obtain a neutral center for the linearized scheme, it is sufficient to assume that the two eigenvalues given by with β = trace(A) = 2 − ∆t 2 a b, have a non-zero imaginary part. Indeed, in this case, for i = 1, 2 , |λ i | = λ iλi = √ λ 1 λ 2 = det(A) = 1. In order to ensure that λ 1,2 have a non-zero imaginary part, it is sufficient to require β > −2 from which the condition ∆t < 2 √ a b follows. We conclude that Theorem 3.3. Provided that ∆t < 2 √ a b , the numerical scheme (7), linearized around the equilibria, has the same stability properties of the Lotka-Volterra model i.e., it has a saddle point at the origin and a neutral center at the internal equilibrium (13).
Notice that the requirement ∆t < 2 √ a b is less demanding than ∆t < min 1 a , 1 b i.e., the condition to be satisfied in order to have positive and (linearly stable) solutions by means of the symplectic Euler method.

Space independent dynamics. Let us consider general predator-prey dynamics with equations of motion given bẏ
with u, v scalar functions of time. The existence of a continuous M -function for the previous system yields to the following necessary condition (see [17]) In this case f and g are called the components of a M-vector field.
In general, the condition (16) is not verified by the main ecological models and consequently they do not exhibit the Poisson structure defined by the matrix B given in (5). However, it is of a great importance to take into account other possible physical constraints. A crucial property that needs to be transferred to the numerical solution is the preservation of the phase space R 2 + . Numerical methods such as Runge-Kutta or linear multistep schemes do not explicitly incorporate the requirement for the solution to be positive. Even though these schemes guarantee the convergence of the discrete solution to the exact one, it occurs that qualitative properties, as positivity, are retained with an accuracy which depends on the chosen temporal step size. In order to build integrators which provide positive solutions independent on the temporal grid, one way is to follow the approach in Section 3.
Consider the system in terms of variablesũ = ln(u) andṽ = ln(v): apply a symplectic Euler scheme as follows: and perform an inverse log transform u = exp(ũ) and v = exp(ṽ) to obtain a positive approximation of (15): φ ∆t :

FASMA DIELE AND CARMELA MARANGI
Of course, the adjoint of (17) advances in time according to the steps: Needless to say, the method (17) (resp.(18)) reduces to (7) (resp. (8)) in case we are dealing with M -systems. We will refer to schemes (17) and (18) as positive symplectic Euler (PSE) schemes. In general, PSE schemes are positive implicit integrators. However, in some cases (as for the LV model), we get the explicit form.
Theorem 4.1. Consider the system (15) and suppose that and then the PSE schemes correspond to first-order splitting procedures.
Proof. Indeed, they can be obtained as composition φ ∆t = ϕ It is easy to see that we may get an explicit scheme with weaker requirements on the vector fields: Theorem 4.2. In (15) suppose that g(u, v) = v G(u). Then the PSE scheme (18) is an explicit positive integrator of (15). Similarly, if the condition (19) is solely satisfied, then the PSE scheme in (17) assumes an explicit form.

4.1.1.
Example: Z-control of Lotka-Volterra dynamics. Let us consider the Zcontrolled Lotka-Volterra system in [16] The theoretical temporal evolution is given by the functions where We set the parameters as follows: α = δ = 0.6, β = γ = 0.01, u d = 100; with this choice, the system has the equilibrium (u * , p * ) = (100, 60) as the unique attractor of the dynamics. Moreover, in [16] it has been proved that for 0.1333 ≤ λ ≤ 1.4049 the solution in (22) stays strictly positive.
We compare the performance of the positive symplectic Euler (17)  40. Notice that, since condition (19) is satisfied, (17) has an explicit functional form. The numerical solutions, obtained with step size ∆t = 0.1, are compared with the theoretical one (22) and plotted in the phase space in Figure 2 (left). In the same Figure, on the right, we plot the approximated values of the predator function v(t) against the time. Notice that, differently from the positive symplectic Euler, the explicit Euler method loses positivity in the transient phase.

4.2.
Spatially extended dynamics. PSE schemes (17) and (18) for spatially implicit dynamics may be applied in order to approximate predator-prey spatially extended phenomena described by reaction-diffusion systems (RD). The general form of a RD dynamics is given by where the real valued functions u = u(x, t) and v = v(x, t) represent the prey and predator densities at time t ≤ T in a bounded domain Ω ∈ R M (M ≤ 3) and L is the Laplace operator. Positive diffusion coefficients D u and D v and suitable initial and boundary conditions are supposed to be provided. To simplify the notation let equation (23) represent the spatial discretization of a reaction-diffusion problem. In this case L is a stiff matrix with large norm and the application of f and g is intended componentwise.
One of the most popular approaches for integrating (23) is to apply an implicitexplicit (IMEX) Runge-Kutta method i.e., a partitioned Runge-Kutta scheme consisting of a diagonally implicit method joint with an explicit one (see [15,10,12]). This approach is motivated by the fact that the stiff term L has to be suitably approximated by an implicit scheme, while an explicit method can be used to discretize the non-stiff (or mildly stiff) reaction term. Recently, in [6] and [7] the authors proposed a different approach based on a partitioned Runge-Kutta scheme which approximates the diffusive dynamics by a diagonally implicit method, while the reaction term is approximated by a partitioned symplectic Runge-Kutta method (IMSP schemes).
The separate treatment of the linear diffusion and the non linear local reaction is also at the basis of splitting methods. Splitting techniques are widely used in literature for a wide class of evolutionary problems [2] [4]. More recently, they seem to be a promising approach also for the numerical treatment of reaction-diffusion problems [14]. The convergence behaviour and the geometric properties of first and second order splitting procedures applied to semilinear evolution equations in an abstract Banach setting have been analyzed in [14].
However, it is possible to construct integrators by composing exact flows with numerical approximations (see [13]). Here we follow this approach and, by means of the PSE schemes (17) and (18), novel positive methods for spatially explicit models (23) are built. Let us denote by ϕ ∆t the exact solution of the pure diffusive problem and consider the PSE scheme (18) for solving the reaction dynamics. Then, the Lie composition advances in time according to the steps: Scheme (24) and its adjoint Φ * ∆t = ϕ ∆t • φ ∆t will be referred to as positive Lie splitting (PLS) methods.
We remark that positive Lie Splitting in (24) results to be a first-order integrator characterized by an explicit functional form in case of dynamics governed by reaction term g satisfying (20). This condition is fulfilled by some widely used models in ecological applications, evolving in spatially extended domains. Indeed, dynamics defined as in (15) with known as Rosenzweig-MacArthur (RM) model [22], or dynamics governed by both satisfy condition (20). Similarly, the positive Lie Splitting Φ * ∆t = ϕ ∆t • φ ∆t results to be a an explicit first-order integrator for dynamics governed by a reaction term f which satisfies (19).
Finally, let us remark that the positive Lie Splittings above considered, result in effective positive numerical integrators if the space discretization also preserves positivity i.e., if e ∆t Du L and e ∆t Dv L are positive operators. In this regard, we refer to [23], where the author studies the positivity of some piecewise linear finite element discretization of the initial boundary value problem for the homogeneous heat equation. In our numerical experiments we used both a suitable implicit in time procedure for the correct treatment of the stiff matrix L and a polynomial Krylov procedure to approximate the action of the discretized Laplace operator on a vector. Moreover, other approaches can be investigated to increase the efficiency, as the use of rational Krylov procedures [5], or fast Fourier methods. In any case, it is crucial the requirement that the space discretization also preserves positivity. 4.2.1. One dimensional example: a simple isothermal chemical reaction system. A simple isothermal chemical reaction system is described, in dimensionless form, by the following nonlinear system: where k is a positive parameter. It has been shown in [24] that both explicit and implicit finite difference schemes can introduce chaotic behavior in the numerical solution of (26), although such a behavior is not inherent in the original PDEs. In particular, it has been shown in [20] that (26) has traveling wave solutions for values of 0 < k < 1. We set k = 0.1 and the initial and boundary conditions considered both in [21] and [20] i.e., u( 200,200] and u(±200, t) = 1, v(±200, t) = 0.
The approximated values u n m ≈ u(x m , t n ) and v n m ≈ v(x m , t n ) given by the positive Lie splitting φ ∆t • ϕ ∆t are defined as follows: (27) Notice that the pure diffusive flows have been approximated by using the finite difference formula which are implicit in time, with this assuring positive entries of the vectors u n1 , v n1 .
For a comparison, we apply the nonstandard positivity preserving method described in [21], which evolves in space and in time according to the following scheme:  of species, the system stays homogeneous forever and no spatial pattern emerges. For a very weakly perturbed initial distribution, a smooth pattern arises that is not persistent and gradually evolves to the homogeneous distribution. For somehow stronger initial perturbations, the system evolves to the formation of a jagged spatial pattern which is persistent in time.
We solve numerically the equations when the initial conditions describe a prey (phytoplankton) patch in a domain with a constant-gradient predator (zoo-plankton) distribution: For the spatial discretization we used the five-point central difference approximation of the Laplacian in two dimensions on a grid with step size h = 1. For the time stepping procedure, the positive (explicit) Lie Splitting (24), written as Φ (RM ) : where α = 0.4, β = 2, γ = 0.6 has been implemented.
For the chosen values of the initial condition the spirals appear with their centers located in the vicinity of the equilibrium u * = 6/35 and v * = 116/245 of the spatially implicit dynamics. The disruption of the spirals, which also begins near the critical points, leads to the formation of two growing embryos of the patchy spatial pattern and finally to the appearance of irregularities in the whole domain.
In order to detect the beginning of the formation of the spirals we followed the dynamics evolution in a smaller domain [0, 300] × [0, 300] and we considered the prey dynamics at t = 120. For IMEX and IMSP schemes we start with a temporal step size ∆t = 1/3 (in order to have convergence) and then we reduce the step size. In Figure 4 (first column) we observe that as ∆t is reduced, the initial spiral patterns, generated by the IMEX scheme with large temporal step size, disappear, thus confirming that the spiral pattern of Fig. 11(b) in Medvinsky et al. (2002) is a numerical artifact. We notice that this is not the case of positive Lie splitting method, which exhibit no artifacts even for large step size i.e., ∆t = 1 ; by reducing the step size the solution does not change significantly from the one obtained with ∆t = 1 (Figure 4, third column).

4.3.
Higher order schemes. To increase the accuracy of positive Lie splitting methods, we may take as a basis the positive Lie splitting method Φ ∆t in (24) and perform successive compositions with its adjoint as in (12). As an example, the positive Strang splitting procedure (s = 1, β 1 = α 1 = 1/2) is defined as Φ ∆t/2 • Φ * ∆t/2 = φ * ∆t/2 • ϕ ∆t • φ ∆t/2 and advances according to the steps: Notice that the above procedure, in general, results to be implicit even when condition (19) is satisfied and the basis of the composition, the positive Lie splitting method Φ ∆t , turns out to be explicit. Indeed, the adjoint of an explicit integrator is always implicit. When conditions (19) and (20) are both satisfied the positive Lie splitting methods reduce to classical splitting schemes obtained by the composition of exact flows. In this case the resulting algorithm is an explicit integrator.
In the general case, besides the implicit functional form, other drawbacks arise in increasing the accuracy order of the positive Lie splitting procedures for solving spatio-temporal models described by RD systems. Without entering into details, we mention the order reduction due to the conditions imposed on the boundary of the spatial domain [8,9] and the necessity for using complex coefficients (with positive real part) for incrementing the second order accuracy [3].
However, let us remark that in spatially extended ecological models it is a common believe that numerical errors are always overwhelmed by uncertainties in the data and in the governing equations. Hence, the gain in accuracy reached by an higher order method does not justify the higher computational costs. For this reason, for many ecological applications, the positive Lie splitting method, by retaining the positivity of the solution, may represent an effective 'first-order' alternative to the widely used IMEX first-order procedures.

5.
Conclusion. Symplectic procedures applied to dynamical systems written in terms of the log transformation of the original variables allowed us to propose novel positive numerical integrators for approximating predator-prey models. The Poisson structure and the positivity of the approximation is ensured without any restriction on the temporal step size for the class of M -systems in population dynamics. When applied to separable M-systems, the resulting schemes are proved to be explicit, positive, Poisson maps. The approach is firstly extended to general predator-prey dynamics non depending on the spatial dimension and successively to reaction-diffusion equations describing spatially extended temporal phenomena. The resulting schemes turn out to be explicit integrators for some widely used predator-prey dynamics in ecological application as the RM model which is characterized by logistic growth of the prey and Holling II type functional response. Our numerical simulations provide results which outperform the solutions found in the recent literature.