A Moving-Mesh Finite-Difference Method for Segregated Two-Phase Competition-Diffusion

: A moving-mesh ﬁnite-difference solution of a Lotka-Volterra competition-diffusion model of theoretical ecology is described in which the competition is sufﬁciently strong to spatially segregate the two populations, leading to a two-phase problem with a coupling condition at the moving interface. A moving mesh approach preserves the identities of the two species in space and time, so that the parameters always refer to the correct population. The model is implemented numerically with a variety of parameter combinations, illustrating how the populations may evolve in time.


Introduction
Ecological competition is a widely studied concept in both theoretical and experimental ecology. In particular, the interspecific competition has long been one of ecology's most prevailing pursuits as it is one of the factors that affect the evolution of species and can alter populations and community structure.
Many researchers have made efforts to develop models to investigate species competition from the viewpoint of mathematical ecology. For a theoretical understanding of spatial patterns arising in population dynamics for competitive species, several free boundary problems have been proposed [1][2][3][4][5] which provide useful theoretical results. Various numerical methods have been applied to solve such free boundary problems, for instance, in [6] the authors use a front tracking approach and a front fixing approach to a two-species competition-diffusion model with two free boundaries, whereas in [7] a moving mesh finite element method is used solve a two competitive segregated species that cannot coexist in space.
Due to the complexity of the equations numerical approximation is useful both for extracting quantitative solutions and for achieving a qualitative understanding of the behaviour of the solution. However, special attention must be paid to the moving interface, whose location usually requires higher resolution than the rest of the domain while the solution may exhibit singular behaviour there.
To avoid this difficulty adaptive methods have been used which modify the resolution of the domain as the solution evolves in response to changes in the dependent variable. Adaptive methods become preferable to fixed mesh methods when the areas of interest are a fraction of the domain, such as the moving interface. With adaptive methods greater precision can be achieved locally without having to increase the resolution everywhere in the domain, which would lead to a very computationally expensive scheme.
There are three most used adaptive mesh methods, namely h-, p-, and r-refinement. The first involves repeated subdivision of intervals of a fixed mesh around the areas of interest. The second is often used in combination with h-refinement and includes higherorder polynomials in each interval between the mesh points, so that functions are better approximated. Finally, r-refinement is the movement of existing mesh points at each time-step to track a feature of interest.
Constructions of r-refinement moving mesh methods vary considerably. In [8] moving mesh methods are classified as velocity-based or location-based approaches. In the latter case the solution is found from a moving form of the PDE which is solved simultaneously with a mesh equation to generate the node position and the solution together. This moving mesh approach (MMPDE) is described in [9,10], see also the Moving Finite Element method [11,12]. In the present paper a velocity-based r-refinement scheme is used in which the mesh is generated by one Eulerian conservation principle while the solution is determined algebraically from a Lagrangian form of conservation [13][14][15]. We first briefly describe the moving mesh method based on conservation.
The conservation method uses an integral to preserve a desired conserved quantity within each patch of elements from which the velocities are constructed. For a mass conserving problem, in which the global mass is constant this quantity is the conserved local mass. However, for a non-mass conserving problem the theory uses the concept of relative mass. By the Leibniz Integral Rule we can then construct an equation from which the velocity of each node is derived. For a unique solution of that equation, the population density must be known at a node or a velocity must be applied to a node which may be thought of as an 'anchor' point.
The approach used in this paper is similar to the finite difference method used in [15] where the method was applied to one-dimensional moving boundary problems such as the mass conserving porous medium equation, Richard equation in hydrology, and the Crank-Gupta problem that does not conserve mass. Finite element versions of the velocity-based moving mesh approach have been used by Baines, Hubbard and Jimack in [13,14] and by Baines, Hubbard, Jimack and Mahmood in [16].
In this paper we apply the moving mesh finite difference method based on conservation to a PDE system of Lotka-Volterra competition model, first proposed by [17], which describes a two phase segregated reaction-diffusion system with a high competition limit where the species are completely spatially segregated and only interact through an interface condition. It is shown in [17] that where the competition is strong enough to spatially segregate the two populations the Lotka-Volterra system can be reduced to a form similar to a Stefan problem in physics [18]. The two major differences are firstly, that there are additional logistic growth terms in the Lotka-Volterra model and secondly, there is a parameter in the Lotka-Volterra model of the interface condition (the equivalent of the latent heat coefficient of the Stefan problem) which is set equal to zero. Unlike the Stefan problem, one species does not transform into another which means that the competition system has an interface condition that specifies the interface velocity only implicitly.
In [19] the authors considered the segregation problem due to high competition with inhomogeneous Dirichlet boundary condition while similar studies in the case of Neumann boundary conditions are presented in [20,21].
The system of equations presented in this paper is suitable to describe concepts in ecology when two species with similar ecological niches cannot co-exist, known as the competitive exclusion principle [22]. One will always overcome the other, so the more competitive species will stay and the subordinate one will either adapt or be excluded by either emigration or extinction.
The layout of the paper is as follows. Section 2 gives details of the Lotka-Volterra competition model with a high competitive rate and describes the relative conservation principle approach and its finite difference implementation, together with the algorithm of the moving mesh finite difference method. In Section 3 illustrations are given for a variety of parameter combinations, observing the various behaviours that dominate as the species evolve through time. Finally, Section 4 gives a brief discussion of the results and potential research directions.

The Lotka-Volterra System
The Lotka-Volterra system is the two-component reaction-diffusion system where u 1 (x, t) and u 2 (x, t) are the population densities of two competing species in abutting regions R 1 (t) and R 2 (t), the parameters δ 1 , δ 2 are constant diffusion coefficients, and are reaction terms in which K 1 , K 2 are species-specific competition rates, k 1 , k 2 are the carrying capacities of the species, and r 1 , r 2 > 0 are reproductive rate parameters.
In [17] it is demonstrated that for two species completely segregated the reaction terms can be reduced to so that Equations (1) and (2) become The resulting system represents the limit in which the carrying capacities k 1 , k 2 values are very large, i.e., the competition rate is high enough that the two species cannot coexist in space and interact only through the interface boundary.
Initial conditions on u 1 and u 2 are selected such that one species is in growth and the other in decline. These are shown in Figure 1.
Zero Neumann boundary conditions ∂u 1 /∂x = 0 and ∂u 2 /∂x = 0 are applied at fixed external boundaries away from the interface.

The Interface Conditions
At the interface between the two species there is a condition that gives the relationship between their fluxes. In essence, the species both flow into the interface and annihilate each other with a rate determined by the ratio of the interspecific competition coefficients . This condition is given as Because the annihilation is complete we also have zero Dirichlet conditions u 1 = u 2 = 0 at the interface. Initial conditions for the competition system, with population density u 1 of species 1 and u 2 of species 2 presented by the graph in 0 ≤ x ≤ 0.5 and in 0.5 ≤ x ≤ 1 respectively. The interface node has zero population and must always satisfy the interface condition.

A Relative Conservation Principle
Define the total population of species p as . Then by Leibnitz' Integral Rule, The final term vanishes by the boundary and interface conditions, sȯ From (3) and (4), which can be integrated in time to give θ p . We now suppose that population fractions c(Ω p ) in each moving subdomain Ω p (t) are independent of time, so that θ p (t) and u p (x, t) satisfy the relative conservation principle Note that u p θ p are conserved quantities. Since the population fractions c(Ω p ) are constant in time, they are determined by the conditions at the initial time t 0 , i.e., Writing (8) as and differentiating the left hand side of (9) with respect to time using Leibnitz Integral Rule, where v p is the velocity of points of the domain. Therefore, by (9), given the population fractions c(Ω p ) of the total massθ p , the velocity v p and u p satisfy the equations where theθ p are given by (3) or (4), leading to We let the subdomains Ω 1 (t) in the region R 1 (t) consist of the interval (a, x(t)) where a is a fixed boundary and let x(t) be any point in the region R 1 (t). Similarly the subdomains Ω 2 (t) in the region R 2 (t) consist of the interval (x(t), b) where b is a fixed boundary and x(t) is any point in the region R 2 (t).
The boundary conditions at the external boundaries a and b are ∂u 1 /∂x = ∂u 2 /∂x = 0, and also v 1 = v 2 = 0 because the boundaries are fixed. Together with the condition that u 1 = u 2 = 0 at the interface boundary, Equation (10) for the velocities v 1 and v 2 and the rates of change of the total massθ 1 andθ 2 satisfy and respectively, where From (6),θ Since the population density u = 0 at the interface and the population densities either side of the interface are positive, the density function is 'V' shaped at the interface.
From [17] the interface condition is given by (5). Whilst the interface velocity is not given explicitly by (5) this equation does determine the location of the interface implicitly. Thus, if we know ∂u/∂x adjacent to the interface in each region we may use the condition that u = 0 at the interface to infer an interface position such that the values of δ p ∂u p /∂x either side of the interface are in the ratio −µ.
We now describe a finite difference numerical method for the solution of the problem.

Numerical Solution
Let the domain (a, b) be (0, 1). At time level t = t n define time-dependent mesh points where x n m is the node at the moving interface, and let u n i , (0 ≤ i ≤ N + 1), approximate u(x, t) by u n i at these points. The initial values θ 0 1 and θ 0 2 of the total mass approximations θ n 1 ≈ θ 1 (t) and θ n 2 ≈ θ 2 (t) of (13) and (14) are estimated by the composite trapezium rule and the constant-in-time relative masses c 1,i and c 2,i in the interval (x n i−1 , x n i ) by at the initial time t = t 0 .
For the initial conditions we take the x 0 i to be equally spaced and the u 0 i pointwise from an initial function chosen to resemble the one in [7] (see Figure 1). We remark that in the case of the chosen initial conditions, a uniform mesh x 0 i can be obtained only for an odd N.

Rates of Change of the Total Populations
The rates of change of the total populationsθ 1 ,θ 2 of (6) are approximated by composite trapezium rules, in region 1 from (13), and in region 2, from (14),θ

Approximating the Velocities
From (11), using the composite trapezium rule, the velocity v n i in region 1 satisfies, where we have taken the subdomain Ω n 1 to be the interval (x n , x n m ). Similarly, from (12), the velocity v n 2 in region 2 satisfies where we have taken the subdomain Ω n 2 to be the interval (x n m , x n ).

Time-Stepping
For the time integration we adopt an explicit Euler time-stepping approach. Given the u i , we update the total masses θ p from the equationθ p = dθ p /dt, (p = 1, 2) using (18) and (19) by θ n+1 p = θ n p + ∆tθ n p (20) (p = 1, 2), where ∆t is the time step, and the mesh points x n i are updated from the equation The updates are first-order accurate in time and subject to limitations on the time step to preserve node ordering.
The explicit Euler method is adequate since the time step used is small and the main purpose of the paper is the treatment of the interface. Implicit schemes for the nonlinear Equations (3) and (4) require convergence of an iteration, details of which are not central to the method. Higher order explicit schemes can be found in [23,24].
Note that in case of a zero velocity there is the following well-known sufficient condition on a time step ∆t in the explicit scheme to prevent the u n+1 i (and hence the local mass in an interval) going unstable, Here we take (22) as a guide for a safe time step in the moving mesh case.

The Population Densities
In order to determine the approximate population densities u i at the new time step t = t n+1 from the θ n+1 p and x n+1 i we approximate the relative conservation principle (8) as where from (16) and (17) the constants are dependent only on initial values. Thus, once the x n+1 i have been found, in region 1, the approximate population density u n+1 i is given by and in that region 2 by while u n+1 m = 0 from the interface condition. Note that the values of u n+1 m±1 determined by (23) depend on x n+1 m , which is not yet known at t n+1 . This value can however be found using the one-sided approximations where from (16) and (17)

Approximating the Interface
The interface condition (5) is approximated by where the subscript m denotes the interface node and the x m±1 , u m±1 are adjacent node positions and solution values. Since u m = 0, from (27) an approximation to the position of the interface node x n+1 m in terms of adjacent nodal values at m ± 1 is Thus, once the other x n+1 i , u n+1 i have been updated, x n+1 m can be found from (28).

Algorithm
In summary, the moving mesh finite difference solution of the competition-diffusion problem given by Equations (3) and (4) with the interface condition (5) on the moving mesh in 1-D generated by (8) is given by the following algorithm.
From the initial mesh and the initial condition compute the initial values θ p (0), (p = 1, 2) of the total populations of the species from (15) and the values of the relative masses c p,i and c p,i from (16), (17) and (24).
Then for each time step: 1.
Generate the nodes x i at the next time-step from the v i using the explicit Euler scheme (21), 5.
Update the population densities u i at the next time level in each region from (25) and (26), 6.
Update the new position of the interface node x m at the next time level from (28).

Results
We find that the model is stable and robust for a variety of parameter choices even when using the explicit Euler integration scheme when using a sufficiently small time step. We observe minimal oscillations affecting the smoothness of the results.

A Parameter Choice
In the body of work concerning Lotka-Volterra equations there is a vast range of parameter values in use because there are so many varied but suitable examples of the type of competition that are described here. We select a conservatively representative set of parameters, chosen to demonstrate some of the behaviour that this model is able to describe.
For the first example we choose a set of parameters that favour species 1, namely δ 1 = δ 2 = 0.01, k 1 = k 2 = 100, r 1 = r 2 = 1 and µ = 3. Even if species 2 makes territorial gains at early stages with the moving interface shifting towards left, the increase in mass of species 1 becomes its weapon to transform it from the inferior species to the superior one ( Figure 2). The moving interface changes direction, moving towards the right at an approximately constant velocity where species 1 continues to increase in density with a rate that decreases as time progresses (Figure 3). As we approach the annihilation of species 2, the interface velocity increases again. This is due to the low mass of species 2 affecting its ability to grow (Figure 4). The movement of the interface is shown in Figure 5.  Figure 4. Result of competition model at t = 8. Here we use δ 1 = δ 2 = 0.01, k 1 = k 2 = 100, r 1 = r 2 = 1 and µ = 3. We run the model with a time step of 0.00001 for 800,000 iterations and plot the results every 0.01. We observe that whilst species 2 initially grew in mass, it will now be wiped out by competition with species 1.The initial conditions are shown in red, with species 1 in blue and species 2 in green.  Figure 5. Movement of the interface position x m for the competition model with parameters δ 1 = δ 2 = 0.01, k 1 = k 2 = 100, r 1 = r 2 = 1 and µ = 3. We run the model with a time step of 0.00001 for 800,000 iterations. We see the interface velocity accelerate as we approach an annihilation event.

Carrying Capacities
We now investigate other parameter choices. We restrict the growth of species 1 by lowering its carrying capacity and observe that in this scenario neither species is dominant, even though all the competition and diffusion characteristics are unchanged. Here we use δ 1 = δ 2 = 0.01, k 1 = 50, k 2 = 150, r 1 = r 2 = 1 and µ = 3. With these differently chosen carrying capacities we find that the interface position is approximately steady and the two species are in balance. This scenario is shown in Figure 6. Density dependence can affect the ability of a species to compete. In the case of decreasing the carrying capacity of species 1 we observe that species can both exist in space, still competing for common resources with none shifting towards other regions in space or in some cases becoming extinct.

Diffusion Characteristics
Alternatively we may adjust the diffusion characteristics of the system. By allowing species 2 to diffuse at a higher rate, we observe that species 2 is able to make territorial gains due to this property alone (Figure 7). Here we use δ 1 = 0.01, δ 2 = 0.05, k 1 = k 2 = 100, r 1 = r 2 = 1 and µ = 3. Due to the growth characteristics, we can see interesting temporal effects. Here the interface velocity has actually reversed directions as the system changes from diffusion dominated to growth dominated. We observe that species 2 is able to make territory gains initially due to its high diffusion rate, even though the competition rate is unaltered. However, as time goes on, the growth and competition characteristics become increasingly important. We see species 1 becoming more dominant over time, so that the interface velocity actually reverses direction. Figure 8 shows the evolution of the system at t = 11 and Figure 9 shows the movement of the interface with the direction reversal.  Figure 7. Result of competition model at t = 1.5, considering the effect of an increased diffusion rate for species 2. Here we use δ 1 = 0.01, δ 2 = 0.05, k 1 = k 2 = 100, r 1 = r 2 = 1 and µ = 3. We run the model with a time step of 0.00001 for 150,000 iterations and plot the results every 0.01. The figure shows the rapid territorial gains of species 2 over species 1 due to its high diffusion rate. The initial conditions are shown in red, with species 1 in blue and species 2 in green. Here we use δ 1 = 0.01, δ 2 = 0.05, k 1 = k 2 = 100, r 1 = r 2 = 1 and µ = 3. We run the model with a time step of 0.00001 for 1,100,000 iterations and plot the results every 0.01. We see that the initial diffusion-driven gains by species 2 are reversed, and that the overall growth characteristics are dominating so that species 1 is gaining territory. The initial conditions are shown in red, with species 1 in blue and species 2 in green. Here we use δ 1 = 0.01, δ 2 = 0.05, k 1 = k 2 = 100, r 1 = r 2 = 1 and µ = 3. We run the model with a time step of 0.00001 for 1,100,000 iterations. Due to the growth characteristics we can see interesting temporal effects. Here the interface velocity has actually reversed direction as the system changes from diffusion-dominated to growth-dominated.

Discussion
In this paper we constructed a moving mesh finite difference method based on conservation for the Lotka-Volterra competition system with a high competition limit, such that the species are completely spatially segregated at an interface. The system of equations produced interesting and realistic behaviour in a very stable model. We were able to implement the model with a wide variety of creative parameter combinations, and observed various effects dominating in turn as the populations evolve through time.
The illustrations presented above give confidence that the model and the moving mesh finite difference approach is likely to be able to satisfy the requirements of modelling a wide variety of competition systems and is numerically stable to a large choice of set-up parameters and is able to produce complex behaviours without problems.
For a set of parameters that favour species 1 we see an increasing interface velocity in the initial stages followed by a change in direction and a long steady phase where the interface velocity is approximately constant. Although the population of species 2 initially makes territorial gains it is eventually wiped out by the competition with species 1. As the annihilation of species 2 is approached, the interface velocity increases again. This is due to the low population of species 2 affecting its ability to compete with species 1.
If the growth of species 1 is restricted by lowering its carrying capacity, interestingly, we observe that neither species is dominant, even though all the competition and diffusion characteristics are unchanged. Therefore, density dependence can affect the competitive ability of a species.
In the case of increasing the diffusion rate for species 2, this species is able to make initial territorial gains, even though the competition rate is unaltered. However, as time goes on, growth and competition characteristics become increasingly important and species 1 becomes more dominant, so the interface velocity reverses direction.
A natural extension is to two dimensions along the lines described in [18], a first attempt appearing in reference [25] which foundered only on stability issues. In further work it would be interesting to compare the behaviour of the model against an empirical data set. The model lends itself to alterations to the logistic terms and changes to parameters without the need for any further development. The aim should be to understand the requirements from both a mathematical and quantitative perspective.

Acknowledgments:
The authors wish to acknowledge the work of Watkins [7,25] using finite elements in the motivation for this work.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

MMFDM Moving Mesh Finite Difference Method MMPDE
Moving Mesh Partial Differential Equations PDE Partial Differential Equations