Numerical Methods for a Two-Species Competition-Diffusion Model with Free Boundaries

: The systems of reaction-diffusion equations coupled with moving boundaries deﬁned by Stefan condition have been widely used to describe the dynamics of spreading population and with competition of two species. To solve these systems numerically, new numerical challenges arise from the competition of two species due to the interaction of their free boundaries. On the one hand, extremely small time steps are usually needed due to the stiffness of the system. On the other hand, it is always difﬁcult to efﬁciently and accurately handle the moving boundaries especially with competition of two species. To overcome these numerical difﬁculties, we introduce a front tracking method coupled with an implicit solver for the 1D model. For the general 2D model, we use a level set approach to handle the moving boundaries to efﬁciently treat complicated topological changes. Several numerical examples are examined to illustrate the efﬁciency, accuracy and consistency for different approaches.


Introduction
The reaction-diffusion equations over a changing domain to describe the dynamics of a two-species competition-diffusion model usually take the following form, The nonlinear functions f 1 (U, V) and f 2 (U, V) are assumed to be C 1 functions satisfying f i (0) = 0 , i = 1,2, and in the literature they are often taken to be two-species Lotka-Volterra type competition functions. In the rest of this paper, we will take two-species Lotka-Volterra type competition functions as an example to demonstrate the numerical methods.
The evolution of the moving domains Ω i (t) ⊂ R N , i = 1, 2, or rather their boundaries ∂Ω i (t), i = 1,2 is determined by the one phase Stefan condition which, in the case ∂Ω i (t), i = 1, 2 are C 1 manifolds in R N . For example, the evolution of species U can be described as follows: For any point x ∈ ∂Ω 1 (t), it moves with velocity µ 1 |∇ x U(t, x)|n(x), where n(x) is the unit outward normal of Ω 1 (t) at x, and µ is a given positive constant.
The moving boundaries ∂Ω 1 (t), i = 1, 2 is generally called the "free boundary", and it is well known that, in general, their smoothness is not guaranteed, even if the initial function U(0, x), V(0, x) and initial domain Ω i (0), i = 1, 2 are both smooth.
Mathematically, the two-species competition-diffusion model with two free boundaries has been intensively studied for 1D [1][2][3][4] and for high dimensions with radial symmetry [5,6] through profound mathematical analysis. For instance, the existence of the positive traveling wave solutions connecting different constant equilibria has been addressed in [7][8][9]. The asymptotic spreading speed associated with the Cauchy problem has been studied in [10][11][12]. Many other theoretical results for general models have been achieved in [13][14][15][16][17] and references cited therein.
In contrast, very few numerical methods have been developed to solve such free boundary problems. Most recently, some efficient numerical methods have been introduced to solve the single species model [18]. In general, extremely small time steps are required due to the stiffness of the system. On the other hand, it is always difficult to efficiently and accurately handle the moving boundaries especially for two species. To efficiently handle the moving boundaries, level set methods [19][20][21][22][23][24] and front tracking methods [25][26][27][28] are two popular numerical approaches. One distinct feature of front tracking [29][30][31][32][33][34] is using a pure Lagrangian approach to explicitly track locations of interfaces, but it is difficult to handle topological bifurcations in high dimensions with interaction of two free boundaries, while the level set method can efficiently overcome such difficulties. In this paper, we will introduce a front tracking framework to solve the system for the one-dimensional case, and a level set approach is employed for two dimensions.
The rest of paper is organized as follows. In Section 2, we introduce two approaches for one-dimensional two-species competition system , i.e., a front tracking approach and a front fixing approach. In Section 3, a level set method is discussed for a more general two-dimensional case. In Section 4, numerical examples are performed to show the efficiency, accuracy and consistency for these different approaches. Finally, a brief conclusion is drawn in Section 5.

Numerical Methods for 1D Two-Species Competition-Diffusion Model
A two-species competition-diffusion model with two free boundaries for the density of population of the competing species U and V depending on time t and spatial variable x states as follows: S 1 (t) = −µ 1 U x (t, S 1 (t)), S 2 (t) = −µ 2 V x (t, S 2 (t)), t ≥ 0, U(x, t) = 0 f or x ≥ S 1 (t), V(x, t) = 0 f or x ≥ S 2 (t), t ≥ 0, U(x, 0) = U 0 (x), V(x, 0) = V 0 (x), f or x ∈ [0, ∞), S 1 (0) = S 0 1 > 0, S 2 (0) = S 0 2 > 0. (9) where U(x, t) and V(x, t) represent the population densities of the two species at the position x and time t. D 1 and D 2 are the diffusion rates of species U and V. γ 1 and γ 2 are net birth rates of species U and V. K 1 and K 2 are the competition coefficients of species U and V. The parameters µ 1 and µ 2 measure the intention to spread into the new territories of u and v. Here, the two free boundaries S 1 (t) and S 2 (t) describe the spreading fronts of two competing species U(t, x) and V(t, x). We envision that the two species initially occupy the interval [0, S 1 (0)] and [0, S 2 (t)], respectively, at time t.

Method 1: Front-Tracking Method for 1D Two-Species Competition-Diffusion Model
The problem lies in solving the nonlinear parabolic partial differential Equations (3)-(9) in the unbounded fixed domain (0, ∞) × (0, L) for the variables (t, x). Let us consider the step size discretization k = t, h = x = L/M, and the mesh points (t n , x j ), with t n = kn, n ≥ 0, x j = jh, 0 ≤ j ≤ M and M positive integer. Let us denote the approximate value of U(t n , x j ) and the approximate value of V(t n , x j ) at the mesh point (t n , x j ) , Step1: Track the position of the moving front S 1 (t).
According to the Stefan condition we consider using the central approximation of the spatial derivatives to approximate ∂U ∂x (t, S 1 (t)), which can be divided into the following four cases. 1 When Let us first consider the symmetric point of x i−1 respect to the position S n 1 , which is denoted bỹ We use the Lagrange interpolation to construct polynomial P L from the value of d, h, u n i−2 , u n i−1 , u n i and S n 1 , thus atx i−1 , we use the value of P L at Mathematics 2018, xx, 1 3 of 26 t n = kn, n ≥ 0, x j = jh, 0 ≤ j ≤ M and M positive integer. Let us denote the approximate value of U(t n , x j ) and the approximate value of V(t n , x j ) at the mesh point (t n , x j ) , Step1: Track the position of the moving front S 1 (t).
According to the Stefan condition we consider using the central approximation of the spatial derivatives to approximate ∂U ∂x (t, S 1 (t)), which can be divided into the following four cases. 1 When Let us first consider the symmetric point of x i−1 respect to the position S n 1 , which is denoted byx i−1 . In particular, when S n 1 = x i , x i−1 = x i+1 . We use the Lagrange interpolation to construct polynomial P L from the value of d, h, u n i−2 , u n i−1 , u n i and S n 1 , thus atx i−1 , we use the value of When 0 = x 0 < S n 1 ≤ x 1 , the central scheme approximation of the spatial derivatives to approximate ∂U ∂x (t, S 1 (t)) involves the fictitious value u n −1 at the point (t n , −h). The value u n −1 can be estimated from the second-order discretization of the boundary condition (10), which implies that u n −1 = u n 1 = 0. It is obvious that all the values of u n i on the grid points are equal to 0 except u n 0 . Numerically, we take S 1 (t) = 0, and it can be explained that the species is only located inside one grid mesh. The simulation should stop here indicating that a more refined mesh is needed. 3 When h . Let us first consider the symmetric point of x 0 respect to the position S n 1 , which is denoted byx 0 . Then we consider the value of u n −1 = u n 1 , and use the Lagrange interpolation to construct polynomial P L from the value of h, d, u n −1 , u n 0 , u n 1 and S n 1 . Then atx 0 , we use the value of P L atx 0 instead of u(x 0 ).
When 0 = x 0 < S n 1 ≤ x 1 , the central scheme approximation of the spatial derivatives to approximate ∂U ∂x (t, S 1 (t)) involves the fictitious value u n −1 at the point (t n , −h). The value u n −1 can be estimated from the second-order discretization of the boundary condition (10), which implies that u n −1 = u n 1 = 0. It is obvious that all the values of u n i on the grid points are equal to 0 except u n 0 . Numerically, we take S 1 (t) = 0, and this can be explained by the fact the species is only located inside one grid mesh. The simulation should stop here indicating that a more refined mesh is needed. 3 When x 1 < S n 1 < x 2 as shown in Figure 2, denoting d = h . Let us first consider the symmetric point of x 0 with respect to the position S n 1 , which is denoted byx 0 . Then we consider the value of u n −1 = u n 1 , and use the Lagrange interpolation to construct polynomial P L from the value of h, d, u n −1 , u n 0 , u n 1 and S n 1 . Then atx 0 , we use the value of P L atx 0 instead of u(x 0 ). When h . Let us first consider the symmetric point of x 0 respect to the position S n 1 , which is denoted byx 0 . Then we consider the value of u n −1 = u n 1 , and use the Lagrange interpolation to construct polynomial P L from the value of h, d, u n −1 , u n 0 , u n 1 and S n 1 . Then atx 0 , we use the value of P L atx 0 instead of u(x 0 ). Step2: Track the position of the moving front of S 2 (t).
Step3: Update the value of U(t n+1 , x i ) and V(t n+1 , x i ). 1 When x i = S n+1 1 and x j = S n+1 2 , then we know that U(t n+1 , We consider the central approximation of the spatial derivatives U xx at x l , for l = 0, 1, 2, ..., i − 1, and the central approximation of the spatial derivatives V xx at x m , for m = 0, 1, 2, ..., j − 1, where U and V are updated by the backward Euler Then use the Picard Iteration (or Newton Iteration) to solve the nonlinear system (12).
2 When x i < S n+1 , we use the Lagrange interpolation to construct polynomial P L 1 from the value of h, R 1 , u n+1 i−2 , u n+1 i−1 , u n+1 i and S n+1 1 and polynomial P L 2 from the value of h, R 2 , v n+1 j−2 , v n+1 j−1 , v n+1 j and S n+1 2 . Then at x i+1 and x j+1 , we use the value of P L 1 at x i+1 instead of u n+1 i+1 and the value of P L 2 at x j+1 instead of v n+1 j+1 . For the solution u at x l , for l = 0, 1, 2, ..., i − 1, a standard central approximation in space with backward Euler in time will be employed. u n+1 l = 0, for l = i + 1, ...M. For the solution v at x l , for l = 0, 1, 2, ..., j − 1, a standard central approximation in space with backward Euler in time will be employed. v n+1 l = 0, for l = j + 1, ...M. U and V is updated by the backward Euler in time Picard Iteration (or Newton Iteration) will be applied to solve the nonlinear system (13). 3 When x i = S n+1 1 and x j < S n+1 2 < x j+1 , then we know that U(t n+1 , x i ) = 0. Let u n+1 i = 0, u n+1 l = 0, for l = i + 1, ...M. We consider the central approximation of the spatial derivatives U xx at x l , for l = 0, 1, 2, ..., i − 1. Denoting R 2 = S n+1 2 −x i h , we use the Lagrange interpolation to construct polynomial P L 2 from the value of h, R 2 , v n+1 j−2 , v n+1 j−1 , v n+1 j and S n+1 2 . Then at x j+1 , we use the value of P L 2 at x j+1 instead of v n+1 j+1 , where U and V is updated by the backward Euler in time Then use the Picard Iteration (or Newton Iteration) to solve the nonlinear system (14). 4 , we use the Lagrange interpolation to construct polynomial P L 1 from the value of h, Then at x i+1 , we use the value of P L 1 at x i+1 instead of u n+1 i+1 . For the solution u at x l , for l = 0, 1, 2, ..., i − 1, a standard central approximation in space with backward Euler in time will be employed. u n+1 l = 0, for l = i + 1, ...M. For the solution v at x l , for l = 0, 1, 2, ..., j − 1, a standard central approximation in space with backward Euler in time will be employed. v n+1 l = 0, for l = j + 1, ...M, where U and V is updated by the backward Euler in time Picard Iteration (or Newton Iteration) will be applied to solve the nonlinear system (15).

Method 2: Front-Fixing Method for 1D Two-Species Competition-Diffusion Model
Here we consider transforming Equation (3) and Equation (4) into problems with a fixed domain Step1. Update the front of S 1 (t) and the value of U by front fixing method.
Let us transform Equation (3) into a problem with a fixed domain [0, 1] using the Landau transformation [35,36] Then Equation (3) turns into the form: where: Boundary conditions (5) and Stefan condition (6) take the form: and respectively, while the initial conditions (9) become: Conditions (8) for the initial function U 0 (x) are translated to W 0 (z) as follows: After the transformation, the new problem has been changed to solve the nonlinear parabolic partial differential equations (17) in the unbounded fixed domain (0, ∞) × (0, 1) for the variables (t, y). Let us consider the step size discretization k = t, h = y = 1/M, and the mesh points (t n , y j ), with t n = kn, n ≥ 0, y j = jh, 0 ≤ j ≤ M and M positive integer. Let us denote the approximate value of M(t n , y j ) at the mesh point (t n , y j ), and let H n be the approximation of H(t n ). Let us consider the forward approximation of the time derivatives, m n+1 and the central approximation of the spatial derivatives, From (24) and (25), Equation (17) is approximated by: As usual, we assume that the Equation (26) can be also approximated at j = 0. Equation (26) written for j = 0 involves the fictitious value m n −1 at the point (t n , −h). The value m n −1 is eliminated from the discretization of the boundary and initial condition (21) and (22), Transformed Stefan condition (20) is discretized using first order forward approximation for H (t) and three points backward spatial approximation of ∂M ∂y (t, 1): to preserve accuracy of order O(k) where the coefficients are given by: Step2. Update the front S 2 (t) and the value of V by front fixing method.
Let us transform Equation (4) into a problem with a fixed domain [0, 1] using the Landau transformation [35,36] Then Equation (4) turns into the form: where: Boundary conditions (5) and Stefan condition (6) take the form: and respectively, while the initial conditions (9) become: Conditions (8) for the initial function U 0 (x) are translated to W 0 (z) as follows: After the transformation, the new problem lies in solving the nonlinear parabolic partial differential equations (31) in the unbounded fixed domain (0, ∞) × (0, 1) for the variables (t, z). Let us consider the step size discretization k = t, h = z = 1/M, and the mesh points (t n , z j ), with t n = kn, n ≥ 0, z j = jh, 0 ≤ j ≤ M and M positive integer. Let us denote the approximate value of W(t n , z j ) at the mesh point (t n , z j ), and let G n be the approximation of G(t n ). Let us consider the forward approximation of the time derivatives, w n+1 and the central approximation of the spatial derivatives, From (38) and (39), Equation (31) is approximated by: As usual, we again assume that Equation (40) can be also approximated at j = 0. Equation (40) written for j = 0 involves the fictitious value w n −1 at the point (t n , −h). The value w n −1 is eliminated from the discretization of the boundary and initial condition (35) and (36), Transformed Stefan condition (34) is discretized using first order forward approximation for G (t) and three points backward spatial approximation of ∂W ∂z (t, 1): where the coefficients are given by: Step3. Update the value of WtoM(t n , y i ) with the front information G n+1 and H n+1 . H n+1 , then we know that WtoM(t n+1 , y i ) = 0. Let wtoM n+1 We consider the central approximation of the spatial derivatives ∂ 2 WtoM ∂y 2 at y j , for j = 0, 1, 2, ..., i − 1, where WtoM is updated by the backward Euler When H n+1 , We consider to use the value of P L at y i+1 instead of wtoM n+1 j+1 .
H n+1 ≥ 1 as shown in Figure 3, we consider the central approximation of the spatial derivatives ∂ 2 WtoM ∂y 2 at y j , for , we use the Lagrange interpolation to construct polynomial P from the value of z i , z i−1 , z i+1 , W n i−1 , W n i , and W n i+1 , then Mathematics 2018, xx, 1 9 of 26 Step3. Update the value of WtoM(t n , y i ) with the front information G n+1 and H n+1 .
for l = i + 1, ...M. We consider the central approximation of the spatial derivatives ∂ 2 WtoM ∂y 2 at y j , for j = 0, 1, 2, ..., i − 1, where WtoM is updated by the backward Euler wtoM n+1 j = a n j wtoM n j−1 + b n j wtoM n j + c n j wtoM n j+1 , 2 When H n+1 , We consider to use the value of P L at y i+1 instead of wtoM n+1 j+1 .
Step4. Update the value of MtoW(t n , z i ) with the front information G n+1 and H n+1 .

Level Set Method for 2D Two-Species Competition-Diffusion Model
A general 2D two-species competition-diffusive model for the densities of population of the species U(t, x, y) and V(t, x, y) depending on time t and spatial variable (x, y) states as follows: Figure 4. when H n+1 G n+1 ≥ 1, MtoW(t n+1 , Z M ) = M(t n+1 , G n+1 H n+1 ).

Level Set Method for 2D Two-Species Competition-Diffusion Model
A general 2D two-species competition-diffusive model for the densities of population of the species U(t, x, y) and V(t, x, y) depending on time t and spatial variable (x, y) states as follows: together with the boundary conditions the Stefan conditions v 1 (t, x, y) = µ 1 |∇U(t, x, y)| n 1 (t, x, y) = −µ∇U(t, x, y), t > 0, (x, y) ∈ ∂Ω 1 (t), v 2 (t, x, y) = µ 2 |∇V(t, x, y)| n 2 (t, x, y) = −µ∇V(t, x, y), t > 0, (x, y) ∈ ∂Ω 2 (t), where v 1 (t, x, y) and n 1 (t, x, y) are, respectively, the velocity vector of the boundary point (x, y) ∈ ∂Ω 1 (t), and the unit outward normal of Ω 1 (t) at (x, y) ∈ ∂Ω 1 (t), v 2 (t, x, y) and n 2 (t, x, y) are, respectively, the velocity vector of the boundary point (x, y) ∈ ∂Ω 2 (t), and the unit outward normal of Ω 2 (t) at (x, y) ∈ ∂Ω 2 (t), and the initial conditions The initial functions U 0 (x, y) and V 0 (x, y) satisfies the following properties: Here τ 1 (t) and τ 2 (t) are the unknown moving boundaries of two species U(t, x, y) and V(t, x, y) such that the population are distributed in the domain Ω 1 (t) and Ω 2 (t) separately. D 1 > 0 and D 2 > 0 are the dispersal rates. The parameters µ 1 > 0 and µ 2 > 0 involved in the Stefan conditions (51) and (52) are the proportionality constant between the population gradient at the front and the speed of the moving boundary of two species U(t, x, y) and V(t, x, y) respectively. Following the ideas of [19], here we use a level set approach to effectively capture the front at each new time step and a finite difference scheme to solve the heat equation everywhere away from the front. The idea behind the level set method is to construct a level set function φ, then move φ with the correct speed v at the front and followed by updating u(t, x, y). The new position of the front is stored implicitly in φ. With this approach, we avoid the difficulties that arise from explicitly tracking the front and thus increase the efficiency to deal with complex interfacial geometries.
We construct a level set function φ 1 , such that at any time t, the front τ 1 (t) is equal to the zero level set of φ 1 , i.e., τ 1 (t) = {(x, y) ∈ Ω 1 : φ 1 (t, x, y) = 0} Initially, U(0, x, y) = U 0 (x, y) and φ 1 is set equal to the signed distance function from the front of species U such that φ 1 is negative in Ω 1 (0) and positive in Ω c 1 (0), where d is the distance from the front τ 1 (t).
Given the normal speed, v 1 , at which the front τ 1 (t) moves, we would construct a speed function, F 1 (t, x, y), which is a continuous extension of |v 1 (t, x, y)| from the front τ 1 (t) over the whole computational domain. The governing equation of φ 1 is then given by This equation will move φ 1 with the correct speed at the front by assuring that τ 1 (t) will always coincide with the zero level set of φ 1 at time t.
We also use φ 1 to define the outward normal vector n 1 corresponding to τ 1 by From Equations (51) and (59), we can rewrite the expression for τ 1 (t) as Since F 1 is equal to |v 1 (t)| along the interface, we can combine Equations (58) and (60) to get the following equation, which of course is only valid on the zero level set of φ 1 : Next, we need to extend the velocity function V 1 to a neighborhood of τ 1 (t). Therefore, we get the velocity function over the computational domain which is of course only valid on the zero level set of φ.
Step2. Compute the velocity filed F 1 (x, y, t) of U .
According to (58), (59) and (62), the level set equation turns into One issue in computing ∇U arises from the fact that its approximation is usually in the order O(1) at points close to or on the front.
The approximation to ∇U at τ 1 (t) is based upon approximations to the derivatives of U in four coordinate directions to cut down on grid orientation effects (please see Figure 5 for illustration). Each approximation to a derivative of U can be continuously extended away from the front by the advection equations where u 1 = ∂U/∂x, u 2 = ∂U/∂y, u 3 = ∂U/∂η and u 4 = ∂U/∂ζ on τ 1 (t). Here S is equal to the sign function. Equation (64) through Equation (67) have the effect of continuously extending u 1 , u 2 , u 3 , u 4 away from the front by advecting these fields in the proper upwind direction. Note that these equations will not degrade the value of V 1 on the front because φ 1 is zero on τ 1 (t), hence, so are S(φ 1 where u 1 = ∂U/∂x, u 2 = ∂U/∂y, u 3 = ∂U/∂η and u 4 = ∂U/∂ζ on τ 1 (t). Here S is equal to the sign function. Equation (64) through Equation (67) have the effect of continuously extending u 1 , u 2 , u 3 , u 4 away from the front by advecting these fields in the proper upwind direction. Note that these equations will not degrade the value of V 1 on the front because φ 1 is zero on τ 1 (t), hence, so are S(φ 1 From (63), we end up solving for the right hand side of the equation The spatial first derivatives of φ 1 are approximated by a second-order ENO scheme. We update φ 1 by solving (68) with a third-order Runge-Kutta scheme.
Step3: Update φ 1 to be a signed distance function for one time step.
From Equation (58) and (59), it is clear that the computation of the normal velocity, and normal vector at the front are all dependent upon the level set function φ 1 . However, by Equation (58), the level set function will cease to be an exact distance function even after one time step. In order to keep the accuracy of n 1 , and F 1 , we need to avoid having steep or flat gradients developed in φ 1 . One way to avoid these numerical difficulties is to reinitialize φ 1 to be an exact distance function from the evolving front τ 1 (t) at each time step. In order to reinitialize the level set function, we use the reinitialization scheme of Sussman [37] where φ 1 (0, x, y) = φ 0 1 and S again denotes the sign function. As in [37], the sign function S is smoothed by the equation. From (63), we end up solving for the right hand side of the equation The spatial first derivatives of φ 1 are approximated by a second-order ENO scheme. We update φ 1 by solving (68) with a third-order Runge-Kutta scheme.
Step3: Update φ 1 to be a signed distance function for one time step.
From Equation (58) and (59), it is clear that the computation of the normal velocity, and normal vector at the front are all dependent upon the level set function φ 1 . However, by Equation (58), the level set function will cease to be an exact distance function even after one time step. In order to keep the accuracy of n 1 , and F 1 , we need to avoid having steep or flat gradients developed in φ 1 . One way to avoid these numerical difficulties is to reinitialize φ 1 to be an exact distance function from the evolving front τ 1 (t) at each time step. In order to reinitialize the level set function, we use the reinitialization scheme of Sussman [37] where φ 1 (0, x, y) = φ 0 1 and S again denotes the sign function. As in [37], the sign function S is smoothed by the equation.
The basic idea behind this method is that given a function φ 0 that is not a distance function, one can evolve it into a function φ that is an exact signed distance function from the zero level set of φ 0 . This can be accomplished by iterating (69) to a steady state. As in [37], the sign function S is smoothed by the equation to avoid numerical difficulties while implemented. By using this approach, we avoid having to explicitly find the contour φ 0 1 = 0 and then resetting values of the front φ 0 1 at grid points. From Equation (69), it is clear that the original position of the front will not change, but at points away from τ 1 (t), φ 1 will be evolved into a distance function.
After moving φ 1 by the correct velocity at the front and then reinitializing φ 1 to be an exact signed distance function from τ 1 (t) in Step 3, next we update U(t, x, y). Updating U(t, x, y) essentially boils down to solving the nonlinear parabolic partial difference equation (48) over the whole computational domain in the following three cases: • At points away from the front, which means the nearby four grid points are all inside the domain Ω 1 (t), we solve the nonlinear parabolic partial difference equation (53) by combining the forward Euler method and the five-point stencil scheme.
For example, we use the scheme (71) to update U(t, x, y) at the grid point (i + 1, j) in Figure 6.
to avoid numerical difficulties while implemented. By using this approach, we avoid having to explicitly find the contour φ 0 1 = 0 and then resetting values of the front φ 0 1 at grid points. From Equation (69), it is clear that the original position of the front will not change, but at points away from τ 1 (t), φ 1 will be evolved into a distance function.
After moving φ 1 by the correct velocity at the front and then reinitializing φ 1 to be an exact signed distance function from τ 1 (t) in Step 3, next we update U(t, x, y). Updating U(t, x, y) essentially boils down to solving the nonlinear parabolic partial difference equation (48) over the whole computational domain in the following three cases:.

•
At points away from the front, which means the nearby four grid points are all inside the domain Ω 1 (t), we solve the nonlinear parabolic partial difference equation (53) by combining forward Euler method and five-point stencil scheme.
For example, we use the scheme (71) to update U(t, x, y) at the grid point (i + 1, j) in Figure 6. • For points near the front τ 1 (t), some special care should be taken. We effectively capture the front using the level set function φ 1 . We can use one-sided different sign of φ 1 to incorporate the distances between a point on the front and grid points neighboring it in either the vertical or • For points near the front τ 1 (t), some special care should be taken. We effectively capture the front using the level set function φ 1 . We can use the one-sided different sign of φ 1 to incorporate the distances between a point on the front and grid points neighboring it in either the vertical or horizontal direction. For example, Y f = (− L 2 + (i − 1)h, y f ) ∈ τ 1 (t), we consider two grid points (i, j + 1) and (i, j) which border Y f . In y-direction, we have y j ≤ y f ≤ y j+1 . We introduce and use u n i,j , u n i,j−1 , u n i,j−2 , r and U(n t, − L 2 + (i − 1)h, y f ) = 0 to construct interpolating polynomial P. When updating u n+1 i,j , we use a standard five-point stencil combing forward Euler method by employing P(− L 2 + jh) instead of u n i,j+1 , i.e., For the case when front interacts with x-axis, we use the same process in x-direction. In the special case where we cannot find enough grid points inside the domain to construct interpolating polynomial P, we employ the nearby grid points and intersect points of the front and x and y-axis to construct quadratic polynomial or straight line as the interpolating polynomial P to update U.
For the extreme configuration, where there are only intersect points of the front and x and y-axis near the grid point, we update U = 0 at the grid point.
• If a grid point lies on the front, we set the value U = 0 at that point (in view of (53)). For example, we set U n+1 i−1,j =0 for the grid point (i − 1, j).
Step6. Repeat Step 2 through Step 6 to update φ 1 , φ 2 , U and V for the next time step.

Convergence test of front-fixing method
In the 1D two-species competition-diffusion model (3)-(9) with parameters values . Here we test the order of convergence in space with very refined temporal step size. In Tables 1 and 2 the error (both L 2 and L ∞ ) and the convergence to the solution of front-fixing method is examined, with final time t end = 0.01. The error is computed by the difference of the numerical solution with the exact solution. For all the examples below when the exact solution is not given, the solution with a very fine resolution will be considered as reference or "exact" solution. As expected, a second-order convergence in space for both u and v can be observed.

Convergence test of front-tracking method
In the 1D two-species competition-diffusion model (3)-(9) with parameters values . Here we test the order of convergence in space with very refined temporal step size. In Tables 3 and 4 the error (both L 2 and L ∞ ) and the convergence to the solution of front-tracking method is examined, with final time t end = 0.01. The error is computed by the difference of the numerical solution with the exact solution. For all the examples below, when the exact solution is not given, the solution with a fine resolution will be considered as reference or "exact" solution. As expected, a second-order convergence in space for both u and v can be observed.   The Comparison of Front-fixing with Front-tracking for 1D model In Figures 7 and 8, we use the front-fixing method and front-tracking method to simulate the 1D two-species competition-diffusion model (3)

Example 1.
In the 2D two-species competition-diffusion model (48)-(56) with parameters values(D 1 , µ 1 , γ 1 , K 1 ) = (4, 10, 1, 0.6), (D 2 , µ 2 , γ 2 , K 2 ) = (0.4, 5, 3, 0.5), the initial boundary of species U is set to be an equilateral triangle which centers at the origin point (0, 0) with side length 1, while the initial boundary of species V is set to be a circle which centers at the origin point (0, 0) with radius = 1.5. The initial values U 0 (x, y), V 0 (x,y) and the initial level set functions φ 0 1 (x, y), φ 0 2 (x,y) are set as follows (76) Figure 9 shows the simulation of the evolvement of two species and their moving boundaries along time with an equilateral triangle as the initial boundary of U and a circle as the initial boundary of V. In the figure of boundary line, the red curves represent the initial boundaries, and the blue curves simulate the evolvement of free boundaries.
From Figure 9, we can see that the triangle evolves into a circle during the simulation.

Example 2.
In the 2D two-species competition-diffusion model (48)-(56) with parameters values(D 1 , µ 1 , γ 1 , K 1 ) = (4, 20, 1, 0.6), (D 2 , µ 2 , γ 2 , K 2 ) = (1, 5, 2, 0.5), the initial boundary of species U is set to be a square with side length = 1, centered at (0,0), while the initial boundary of species V is set to be a circle which centers at the origin point (0, 0) with radius = 2. The initial values U 0 (x, y), V 0 (x,y) and the initial level set functions φ 0 1 (x, y), φ 0 2 (x,y)are set as following V 0 (x, y) = 50cos( √ x 2 +y 2 π 4 ), (x, y) ∈ Ω 2 (0), 0 (x, y) ∈ Ω c 2 (0). Figure 10 shows the simulation of the evolvement of two species and their moving boundaries along time with a square as the initial boundary of U and a circle as the initial boundary of V. In the figure of boundary line, the red curves represent the initial boundaries, and the blue curves simulate the evolvement of free boundaries.
From Figure 10, we can see that the square evolves into a circle during the simulation.   Figure 10. The simulated dynamics where initial boundary of U is a square and initial boundary of V is a circle. The snapshots are taken at the times t = 0, 0.05, 0.1, respectively.
(84) Figure 11 shows the simulation of the evolvement of two species and their moving boundaries along time with circles as the initial boundary of U and V. In the figure of boundary line, the red curves represent the initial boundaries, and the blue curves simulate the evolvement of free boundaries.
From Figure 11, we can see that the circles propagate as circles during the simulation.

Conclusions
The system of reaction-diffusion equations with moving boundaries has been intensively studied analytically in recent years, however, very little numerical work has been done in this field due to numerical challenges in tracking free boundaries. In this paper, we first introduce a front tracking framework for 1D model, and compare it with a front-fixing method. Numerical experiments demonstrate that these two methods are consistent with each other. For 2D models, to overcome the difficulty of handling complicated topologically changes, we apply a level set approach to handle the moving boundaries. Numerical examples with different initial configurations demonstrate that the level set approach is able to robustly and efficiently capture different complicated geometries.
Although the level set method is very robust in handling topological changes, sometimes it is very hard to achieve high order accuracy, especially near the fronts. Currently we are extending the front tracking method to more accurately deal with topological changes for general 2D models, and to the systems of two competing species in which each species has its own moving boundary. The front will become more complicated and more challenging once two moving fronts are tangled together, and we would apply the reconstruction strategy to overcome these difficulties.
Author Contributions: X.L. and S.L. proposed and designed the numerical methods and S.L. performed the numerical experiments. X.L. and S.L. together wrote the paper.