Algorithm for numerical solutions to the kinetic equation of a spatial population dynamics model with coalescence and repulsive jumps

An algorithm is proposed for finding numerical solutions of a kinetic equation that describes an infinite system of point articles placed in $\mathbb{R}^d (d \geq 1)$. The particles perform random jumps with pair wise repulsion, in the course of which they can also merge. The kinetic equation is an essentially nonlinear and nonlocal integro-differential equation, which can hardly be solved analytically. The derivation of the algorithm is based on the use of space-time discretization, boundary conditions, composite Simpson and trapezoidal rules, Runge-Kutta methods, adjustable system-size schemes, etc. The algorithm is then applied to the one-dimensional version of the equation with various initial conditions. It is shown that for special choices of the model parameters, the solutions may have unexpectable time behaviour. A numerical error analysis of the obtained results is also carried out.


INTRODUCTION
Stochastic dynamics of infinite populations placed in a continuous habitat may include such aspects as their motion accompanied by merging. A pioneering work in this direction was published by Arratia in 1979, in which an infinite system of coalescing Brownian particles in R was proposed and studied [1]. It was then extended to self-repelling motion of merging particles [2]. Over the next two decades, different mathematical aspects of the Arratia model were intensively studied, see [3][4][5][6] and the references therein. In the Arratia flow, an infinite number of Brownian particles move independently up to their collision, then they coalesce and continue moving as single particles.
Recently, an alternative model of this kind has been proposed [7,8]. Here, analogously to the Kawasaki approach [9,10], particles make random jumps with repulsion acting on the target point. Additionally, two particles merge into one particle placed elsewhere with probability per time dependent on all the three locations and independent of the remaining particles. Thereafter, this new particle participates in the motion. The intensities and magnitude of jumps and coalescence are determined by the spatially nonlocal kernels.
Similar models are used to describe predation in marine ecological systems, see, e.g. Refs. [11,12]. Identical processes are considered in freshwater ecosystems to study the phytoplankton dynamics [13][14][15]. Phytoplankton cells are dispersed in the water, leading to a patchy distribution of entities. The latter is the basis for the vast majority of oceanic and freshwater food chains. Another example concerns the modeling of cancer mechanisms according to one of which melanoma cells migrate and coalesce to form tumors [16,17]. Fresh tumor cells grow into clonal islands, or primary aggregates, that then coalesce to form heterogeneous formations.
Quite recently [18], first numerical results on the jump-coalescence kinetic equation have been reported. As a consequence, the role of repulsion interactions in the possible appearance of a spatial heterogeneity in the system has been elucidated. There it was shown also that the kinetic equation can rigorously be obtained from the corresponding microscopic evolution equations using previous theoretical approaches developed earlier [19][20][21][22][23][24] for other spatial population models. However, the numerical consideration was restricted to a few particular examples when choosing the parameters of the models and forms of the initial density. Moreover, very little was said about an algorithm used to obtain the numerical solutions.
In the present paper, we consistently derive the algorithm enabling to solve the repulsion-jump-coalescence kinetic equation. The derivation combines different techniques and develops methods allowing to reproduce properties of infinite systems on the basis of finite samples. The algorithm is applied to population dynamics simulations of one-dimensional systems with various initial spatially inhomogeneous densities and forms of the jump, coalescence, and repulsion kernels. A comprehensive analysis of the obtained numerical results is also provided.

KINETIC EQUATION OF A POPULATION MODEL
Consider an infinite population in the continuous space R d of dimensionality d ≥ 1. We will deal with two classes of stochastic events in the system: free coalescence and repulsive jumps. Time evolution of the population will be described in terms of the particle density n(x,t). Performing the passage from the microscopic individual-based dynamics to the mesoscopic description by means of a Vlasov-type scaling, one derives the following kinetic equation for the density in the Poisson approximation with the jump-coalescence model [7,8,18]: ∂ n(x,t) ∂t = − a(x, y) exp − ϕ(y − u)n(u,t)du n(x,t)dy + a(x, y) exp − ϕ(x − u)n(u,t)du n(y,t)dy − b(x, y, z) + b(y, x, z) n(x,t)n(y,t)dydz + b(y, z, x)n(y,t)n(z,t)dydz.
The first term in the rhs of Eq. (1) represents the process of jumping of agents from point x to location y with intensity a(x, y) modulated by an exponential factor. This results in a decrease of density n(x,t) in point x at time t. The exponential factor determines the strength of repulsion between the particle jumping to y and the rest of the system, where ϕ(y − u) denotes the corresponding interparticle potential. The intensity contribution and repulsion strength depend on the configuration of all particles, so that the spatial integration over y and u is carried out. The second term relates to the reverse process when particles located anywhere in the coordinate space can appear at x due to their random jumps to the latter point. This increases the density n(x,t). The third term defines the process of free coalescence when two particles, located at x and y, merge into a single particle located at z with intensities b(x, y, z) and b(y, x, z), resulting in a decrease of n(x,t). Finally, merged particles can be randomly created in x owing to the coalescence of arbitrary two other particles located in y and z with intensity b(y, z, x), leading to an increase of the local density.
The double integrations over y and z in the rhs of Eq. (1) are performed in R d . They take into account the influence of all possible pair-wise coalescences on the change of n(x,t) in x at time t. These integrations are invariant with respect to the mutual replacement y ↔ z in b(x, y, z), b(y, x, z), and b(y, z, x). The jump and coalescence kernels are nonnegative functions symmetric w.r.t. first two arguments, a(x, y) = a(y, x) ≥ 0 and b(x, y, z) = b(y, x, z) ≥ 0. The kernels satisfy the integrability conditions a(x, y)dx = a(x, y)dy = µ a , ϕ(x) = µ ϕ , and b( The quantities µ a , µ b , and µ ϕ can be treated as parameters of the model. Moreover, according to the translation invariance of the system, the following shifting identities should be satisfied, a(x, y) = a(x + u, y + u) and b(x, y, z) = b(x + u, y + u, z + u). The obvious choice for the jump kernel is a(x, y) = a(x − y) = a(|x − y|).
The coalescence intensity will be modeled by the form where b(x − y) = b(|x − y|) and the Dirac δ -function is applied. It implies that the resulting point of the coalescence of two particles in x and y is at the middle z = (x + y)/2 of their locations. This is quite natural for identical particles. Instead of δ -function, we can, in principle, choose smoother dependencies, for example, Gaussian-like ones. The δ -function was used to simplify the calculations as then it allows us to lower dimensionality of the integration. Indeed, in view of Eq. (2) and integrating over z to eliminate δ -function, the kinetic equation (1) transforms to where the symmetrical properties of the kernels have been taken into account. Imposing an initial condition n(x, 0), Eq. (3) leads to a complicated partial integrodifferential equation with respect to the unknown function n(x,t) at t > 0. Because of the presence of spatial integrals and nonlinearity, we doubt it can be solved analytically in general. Moreover, we cannot apply the convolution theorem to avoid the spatial integration in the last coalescence term of Eq. (3) due to the existence of three functions under the integrand which all depend on y. That is why we develop a numerical approach for this equation. The convolution method in the absence of coalescence is presented in Appendix. Analytical solutions in the spatially homogeneous case are also given there.

NUMERICAL ALGORITHM
Spatial discretization. -In order to solve numerically the kinetic equation it is necessary, first of all, to perform its discretization in coordinate space. Let x i be the grid points uniformly distributed over R d with mesh h along all dimensions. For simplification of our further presentation we will restrict ourselves in this study to a particular case of d = 1 (the generalization to any higher dimensionality d > 1 is straightforward and can be given elsewhere). Then the discrete duplicate of Eq. where , and the infinite sums over j and k represent the spatial integrals. It is obvious that in the limit h → 0, the discretized kinetic equation (4) coincides with its original continuous version (3). Replacing i − j by j, taking into account that the summation is carried out over the infinite number of terms, and introducing the auxiliary quantities one obtains that Eq. (4) can be cast more compactly in the following equivalent form where a j = a( jh), b j = b( jh), and b 2 j = b(2 jh).
In computer simulations we cannot operate with infinite-size samples leading to the infinite summation over j in Eqs. (4), (5), and (6). Because of this we consider a finite number N of grid points x i uniformly distributed over the area [−L/2, L/2] with spacing h = L/N, where i = 1, 2, . . . , N. This area will represent an interval, a square, or a cube in cases d = 1, 2, or 3, respectively. The finite length L should be sufficiently big with respect to all characteristic coordinate scales of the system. The number N of grid points must be large enough to minimize the discretization noise. Then h will be sufficiently small to provide a high accuracy of the spatial integration. The finite-size effects can be reduced by applying the corresponding boundary conditions (BC) when mapping infinite range x ∈] − ∞, ∞[ by finite area [−L/2, L/2]. In view of the aforesaid, Eq. (6) presents a coupled system of N autonomous ordinary differential equations, where i = 1, 2, . . . , N and summation over j is performed according to BC.
Boundary conditions. -Three types of BC can be employed, namely, Dirichlet (DBC), periodic (PBC), and asymptotic (ABC) boundary conditions. The choice depends on initial function n(x, 0) and properties of solution n(x,t). For example, if n(x, 0) takes nonzero values only within a narrow interval [−l/2, l/2] with l ≪ L, we can apply the DBC by letting n j = 0 for all |x j | > L/2. This means that during the finite simulation time 0 ≤ t ≤ T , the nonzero values of n(x,t) do not approach the boundaries x B = ±(L/2− max σ ), where max σ is the maximal radius of the kernels (see Sect. 4). In numerical calculations this can be expressed by the condition n(x B ,t) < ε max x n(x,t), where 0 < ε ≪ 1 is the relative tolerable level (a negligibly small quantity slightly exceeding machine zero). When the propagation front becomes too close to the boundaries, i.e., n(x B ,t) > ε max x n(x,t), we should enlarge L (e.g. gradually doubling it) until to satisfy the required first condition, use DBC again, and continue the simulation for t > T . A special case, which is discussed in Sect. 4, where members of infinite configuration are initially absent in one half-space, requires a modified BC that combines DBC and ABC with an addition of adjustable system-size approach.
If n(x, 0) and, thus, n(x,t) are periodic functions, it is necessary to apply PBC with fixed finite size L. According to PBC, the summation in Eq. (6) for each current i = 1, 2, . . . , N is performed not only over all j = 1, 2, . . . , N but also over all infinite number of images j ′ of j. The images are obtained by repeating the basic interval [−L/2, L/2] by the infinite number of times to the left and to the right of it, so that x j ′ = x j ± KL, where K = 1, 2, . . . , ∞ and n j ′ = n j . This reproduces the periodicity n(x ± KL, 0) = n(x, 0), where x ∈ [−L/2, L/2]. The solution n(x ± KL,t) = n(x,t) will also be periodic for any time t > 0 with the same (finite) period L. In such a way the infinite system can be handled by a finite-size sample. Because the kernel values a j and b j decrease to zero with increasing the interparticle distance |x j |, the summation over j in Eq. (6) can be actually restricted to a finite number of terms for which |x j | ≤ R a,b < L/2. The truncations radiuses R a,b are chosen to satisfy the conditions a(|x|) ≈ 0 and b(|x|) ≈ 0 for |x| > R a and |x| > R b , respectively.
In the spatially homogeneous case when n(x,t) = n(t), we should apply ABC, i.e. n j ′ = n(t) for all x j ′ < −L/2 and n j ′ = n(t) for all x j ′ > L/2. For this case, PBC and ABC lead to the same results. The ABC can also be used for spatially inhomogeneous solutions n(x,t) which are flat for x < −L/2 and x > L/2 at a given t where they take nonzero constant values. Then n j ′ = n(−L/2,t) for all x j ′ < −L/2 while n j ′ = n(L/2,t) for all x j ′ > L/2. If in the course of time the flatness is violated at a current L, the basic length should be enlarged using the automatically adjustable system-size approach mentioned above.
Taking into account the properties a − j = a j , b − j = b j and the fact that the influence of a(x) and b(x) can be neglected at large distances x > R a,b , i.e. assuming that a(x) = 0 for |x| > R a and b(x) = 0 for |x| > R b , Eq. (6) transforms to where i = 1, 2, . . . , N and the summations over j are performed already with finite positive nonzero integers up to j a = R a /h or j b = R b /h. Quite similarly, using the properties ϕ − j = ϕ j and ϕ(x) = 0 for |x| > R ϕ , one finds from Eq. (5) that at kernel values a j , b j , b 2 j , and ϕ j were introduced in Eqs. (7) and (8) to improve precision of the numerical integration over coordinate space. They are determined according to the chosen method [25]. These weights satisfy the normalization condition ξ 0 + 2 ∑ j * j=1 ξ j = 2 j * , where j * = j a,b,ϕ or j * = j b /2. In particular, we can use the composite Simpson or trapezoidal rules. For example, in the composite trapezoidal scheme of the second order we have that ξ 0 = ξ 1 = ξ 2 = . . . = ξ j * −1 = 1, while ξ j * = 1/2. The composite Simpson method of the fourth order yields ξ j * = 1/3, ξ j * −1 = 4/3, ξ j * −2 = 2/3, ξ j * −3 = 4/3, ξ j * −4 = 2/3, . . . , ξ 0 = 4/3 if j * is odd and ξ 0 = 2/3 if j * is even. The latter is more accurate than the former and involve numerical uncertainties of order of We mention that n i are explicitly defined at i = 1, 2, . . . , N, so that the corresponding BC should be applied in Eqs. (7) and (8) to n i− j and/or n i+ j whenever i − j < 1 and/or i + j > N, because j = 1, 2, . . . , j * < N/2. The uniform knot distribution over [−L/2, L/2] can be chosen in the form This provides the symmetricity of knot positions with respect to x = 0. Then according to PBC the calculations of n i± j should be performed as Note that i = 1, 2, . . . , N and j = 1, 2, . . . , j * < N/2. The applications of DBC and ABC result in respectively. In the cases d = 2 and d = 3, the boundary transformations can be implemented similarly to those given by Eqs. (9) - (11). Time integration. -In the most general form, our coupled system of N autonomous ordinary differential equations can be given as where i = 1, 2, . . . , N and s i (n 1 , n 2 , . . . , n N ) represents the rhs of Eq. (7) with taking into account Eq. (8). Introducing the set (12) can be compactly cast asΓ = Φ(Γ). A common practice to solve differential equations of such a type is to use classical Runge-Kutta schemes [25]. The scheme of the fourth order (RK4) reads where ∆t is the time step and The Runge-Kutta scheme of the second order (RK2) can also be applied. There are two versions of RK2. The first one (will be referred to as simply RK2) is known as the Heun method and can be related to the trapezoidal or Verlet-like integration. It has the form The second version (RK2') relates to a middle point scheme, . The RK2 and RK20 approaches are two-stage schemes. The RK4 integration consists of four stages, so that at the same ∆t it will require twice larger number of operations to cover the same time interval T of the observation with respect to RK2. However, RK4 is more accurate than RK2 with O(∆t 5 ) versus O(∆t 3 ) local truncation errors. In such a way, the numerical solution n i (t) can be found for any t = p∆t, where p = 1, 2, . . . , P over the interval t ∈ [0, T ] with T = P∆t by sequentially repeating P times the transformations given by Eq. (13). The O(∆t 5 )-uncertainties can be reduced to an arbitrary small level by decreasing the length ∆t of the time step. We mention that the finite-size effects are minimized by choosing a sufficiently large size L of the system and applying the corresponding boundary conditions. The uncertainties caused by the discretization are reduced by decreasing mesh h = L/N. The latter can be achieved at sufficiently large values of N ≫ 1. In general, the total number of operations per one time step, which are necessary to obtain solutions n i (t) for the set of N differential equations (12), is proportional to N 2 . For d > 1, the computational efforts will increase drastically with increasing d, namely as N 2d . For convolution solutions this number is lowered to order of N d ln N d (see Appendix). In the case of simple rectangle kernels the total number of operations can be decreased to be proportional to N d .

APPLICATIONS AND ANALYSIS OF THE RESULTS
Numerical details. -The numerical simulations were carried out for onedimensional systems (d = 1) of free or repulsive jumping particles which can coalesce. The kinetic equation (3) was solved by using the algorithm described in the preceding section. The discretization in R was done with a mesh of h = 0.0125 -0.2 in dependence on the choice of initial conditions and kernel parameters. The initial basic length was L = 20 within either the fixed (for PBC) or adjustable-size (for DBC and ABC) regime. Spatial integration was performed by employing the composite Simpson rule. Time integration was done with the help of the RK4 scheme at a step of ∆t = 0.1. The relative tolerable level was chosen to be equal to ε ∼ 10 −12 . Further increasing space and time resolutions did not noticeably affect the solutions (see the last subsection).
The jump a(x), coalescence b(x), and repulsion ϕ(x) kernels were modelled by Gaussian or rectangle functions, where µ = µ a , µ b , or µ ϕ and σ = σ a , σ b , or σ ϕ are the intensities and ranges of the corresponding interactions, respectively. The kernels are normalized, {a, b, ϕ}(x)dx = µ a , µ b , µ ϕ . A symmetrical pair of shifted Gaussian or rectangle functions was involved as well, where F ≡ G or C and s is the shifting interval. The truncation radius for Gaussian kernels was R = Qσ with Q = 6 for which G µ,σ (Qσ )/G µ,σ (0) ∼ 10 −8 , so that their contribution at |x| > R can be neglected. In the case of discontinuous rectangle kernels we have Q = 1 by definition. For the sums of two shifted kernels (16) the truncation distance increases to R = Qσ + s. In view of the above, we have six kernel parameters, µ a , µ b , and µ ϕ , as well as σ a , σ b , and σ ϕ for the repulsion-jump-coalescence model. This leads to a large number of all possible combinations. Because of this we will deal with most characteristic examples when choosing values for these parameters for single Gaussian and rectangle kernels or their ±s-shifted double counterparts. In order to consistently analyze their influence on the dynamics of populations the following four situations will be considered: (i) pure free jumps; (ii) repulsive jumps; (iii) pure coalescence; and (iv) repulsive jumps with coalescence. It is worth emphasizing also that solution n(x,t) depends cardinally on the choice of initial condition n 0 (x) ≡ n(x, 0). Like kernels, the initial condition can be selected in the form of Gaussians (14) or rectangles (15), trigonometric or step functions, etc.
From Fig. 1(a) we see that for pure free jumps, the discontinuity of C 1,1 (x) at x = ±1 soon disappears with increasing time, transforming the initial density into a Gaussian-like shape at t 10. For longer times, t 40, the solution n(x,t) tends to more and more homogeneous densities. In the limit t → ∞ we expect absolutely flat function lim t→∞ n(x,t) = n independent on x, where n = 1 L Ω n(x, 0)dx = 1/L = 0.05. Moreover, for any time t the total number Ω n(x,t)dt of particles on the interval Ω = [−L/2, L/2] is constant because of the absence of coalescence. For strong repulsive jumps [part (b)] the density function n(x,t) remains discontinuous at x = ±1 up to t ∼ 40 and its shape at shorter times is more complicated. In particular, apart from the existence of the main maximum at x = 0 which decreases in amplitude with increasing t, two symmetric secondary maxima appear additionally at 0 < t 40 in the ranges x ≈ ±2 due to the repulsion between particles. We believe that all the maxima disappear at t → ∞ with the same asymptotic behaviour lim t→∞ n(x,t) = n = 1/L = 0.05 as for free jumps. In contrast, for free coalescence [part (c)] we expect decay of n(x,t) to zero at t → ∞. Moreover, here the particles remain to be located exclusively within the initial interval [−1, 1] and they are absent outside of it at any t. In other words, no density propagation is indicated because the particles do not move apart from coalescence.
When the repulsive jumps are carried out in the presence of coalescence at equal interaction ranges σ a = σ ϕ = σ b = 1, see part (d) of Fig. 1, the pattern is somewhat similar to that of part (b). However, the three-maximum structure dissipates now much faster, leading to the zeroth asymptotics already at t 40. For short-ranged jumps, where σ a = 0.5, σ ϕ = 1, σ b = 2, the central peaks at x = 0 become sharper, while the secondary side maxima at x ≈ ±2 do not appear, see part (e) and compare it with subsect (d). In the case of short-ranged coalescence when σ a = 2, σ ϕ = 1, σ b = 0.5, the central peaks transform into a more complicated structure with one central minimum at x = 0 and two side maxima at x ≈ ±0.5, see part (f). The secondary maxima at x ≈ ±2 becomes more visible with respect to those for equalrange interactions [cf. part (f)]. Thus, the influence of jumps on the dynamics increases not only with increasing their intensity but range as well. The same concerns the coalescence. Note also that the density profiles in Figs. 1 (a)-(f) are symmetric, i.e., n(−x,t) = n(x,t), like the initial condition, n(−x, 0) = n(x, 0). This follows from the symmetry of the kinetic equation.
The second choice deals with asymmetric initial condition n(−x, 0) = n(x, 0) in the form of N 0 shifted single rectangle functions C v k ,σ k (x + s k ) with intensities v k FIGURE 1. Time evolution of density n(x,t) as dependent on spatial coordinate x at several moments of time t for rectangle initial condition n(x, 0) = C 1,1,L (x) in the cases: (a) pure free jumps; (b) repulsive jumps; and (c) pure coalescence; as well as repulsive jumps and coalescence for (d) equal interaction ranges; (e) short-ranged jumps; and (f) short-ranged coalescence. The infinite system is periodic (L = 20) on the interval x ∈ [−10, 10] with jump, repulsion, and coalescence Gaussian kernels of different intensities and ranges (see the legends inside). and ranges σ k , namely, where s k = −L/2 + (k − 1/2)L/N 0 are shifting parameters. Repeating (17) with period L, we should apply PBC to deal with the infinite system. We used a particular case of (17) with N 0 = 3 and L = 20 as well as three different amplitudes v 1,2,3 generated at random in the interval ]0, 1[. The corresponding result on this is shown in Fig. 2. Looking at Fig. 2 and comparing it with Fig. 1 we see that behaviour of n(x,t) at short times can be approximately presented as a sum of independent separate solutions obtained for single rectangle initial densities C v k ,σ k (x + s k ). With increasing time, a coherence between the separate solutions appears. Again, in the absence of coalescence [µ b = 0, parts (a) and (b)] the density n(x,t) at t → ∞ flattens in x and seems tend to the nonzero constant n = 1 L Ω n(x, 0)dx = (v 1 + v 2 + v 3 )/L for any initial distributions. We mention that at µ b = 0, the total amount Ω n(x,t)dx of particles in Ω is unchanged and equal to its initial value Ω n(x, 0)dx. For µ b > 0, the density function n(x,t) seems to take its zeroth asymptotics lim t→∞ n(x,t) = 0 at t → ∞ [parts (c)-(f)] with the flow of time (except special choices, see below). For populations with an infinite number of particles and periodic initial conditions n(x ± KL, 0) = n(x, 0) with x ∈ [−L/2, L/2], where K = 1, 2, . . . , ∞, the solution n(x ± KL,t) = n(x,t) will also be periodic for any time t > 0 with the same (finite) period L. Then, in particular, n(−L/2,t) = n(L/2,t) as is confirmed in Fig. 2. Investigations show that the increase of the strength µ and range σ of the jump and coalescence kernels accelerates this process, see for instance, parts (e) and (f) of Figs. 1 and 2. Using the Gaussian (instead of rectangle) initial conditions and rectangle (instead of Gaussian) kernels leads to results (not shown) which are similar to those of Figs. 1 and 2.
As visualized in Figs. 1 and 2, the presence of coalescence (b(x) = 0) seems to lead to the zeroth asymptotic Ω n(x,t)dx → 0 at long times provided the kernels are single rectangle functions with positive values around zero (the same concerns simple Gaussians). However, the coalescence kernel can be chosen in the form b(x) = C µ,σ ,s (x) of a pair of two shifted rectangle functions [Eq. (16)] with appropriate shifting parameter s to avoid the zeroth density limit. The initial inhomogeneous density n(x, 0) over the basic interval [−L/2, L/2] with PBC at period L should also be chosen correspondingly. For instance, we can consider n(x, 0) in the form C v k ,σ k ,s ′ ,N 0 (x) of two (N 0 = 2) single rectangle functions [Eq. (17)] with the same amplitude v k = µ b = 1 and ranges σ k = σ b ≡ σ = 1 as those of the coalescence kernel b(x) but with different shifting parameter s ′ = s satisfying the constraint s ′ = 2s − 2σ . Choosing s = 5 one finds at σ = 1 that s ′ = 8. The corresponding result for the case n(x, 0) = C 1,1,5,2 (x) with L = 20 and b(x) = C 1,1,8 (x) in the absence (a = 0) of jumps is depicted in Fig. 3(a). We see that at the beginning, the rectangles soon transform into triangle-shaped peaks centered at x = ±(5 + KL), while the additional triangles appear exactly in the middle of them at x = 0 ± KL (below we will omit the terms ±KL to simplify notation). With increasing time, FIGURE 2. Time evolution of solution n(x,t) for asymmetric initial density in the form of three rectangle functions with different amplitudes. Other notations are the same as in Fig. 1. the widths of the triangle-shaped peaks decrease but their maxima at x = ±5 become higher while unchanged at x = 0. At sufficiently long times t 2 · 10 4 the modification of the density profile slows down to a level suggesting that the system approaches a non-trivial stationary state in which ∂ n(x,t)/∂t = 0. In other words, the shape of the density profile is transformed to such a form at which any coalescence processes become impossible in view of the specific form of the coalescence kernel b(x) = C 1,1,8 (x). The latter accepts nonzero values only in the interval s ′ − σ = 7 < |x| < 9 = s ′ + σ , so that the absence of coalescence at a given FIGURE 3. Dynamics of spatial structure n(x,t) for initial density in the form of two rectangle functions, n(x, 0) = C 1,1,5,2 (x) in the cases: (a) pure coalescence and (b) coalescence with free jumps. The interactions are modelled by shifted pair rectangle coalescence C 1,1,8 (x) and Gaussian jump G 0.2,1 (x) kernels. spatial configuration means that there exists no pair of particles with interparticle separations |x| lying in the interval [7,9]. Allowing particles to jump changes the situation radically, as is demonstrated in Fig. 3(b) for Gaussian jump kernel G 0.2,1 (x). Even for relatively small jump intensity µ a = 0.2 and range σ a = 1, the density quickly decreases to zero for each x, after initial period of time, when new peaks are formed. It is interesting to remark that the monotonic decrease of main maxima in x ± 5 is accompanied by nonmonotonic change of the magnitude of the newly formed peaks in x = 0 (and ±KL). This magnitude first increases at 0 < t 4 achieving a maximum at t ∼ 4, and then decreases at t 4.
Trigonometric initial density profiles. -Another interesting case is the initial density in the form of a trigonometric function, n(x, 0) = T n 0 ,µ 0 ,k (x) = n 0 1 + µ 0 cos(2πkx/L) , where 0 < µ 0 < 1 is the coefficient of the modulation and k ≥ 1 defines the period L/k in coordinate. Then PBC should be used to reproduce the infinite system. The obtained solutions n(x,t) for n(x, 0) = T 1,1,3 (x) with n 0 = 1, µ 0 = 1, k = 3, and L = 20 are plotted in Fig. 4 when jump, repulsion, and coalescence kernels are Gaussians a(x) = G 1,1 (x), ϕ(x) = G 8,1 (x), and b(x) = G 0.5,1 (x), respectively. For convenience, we presented n(x,t) only on the right-hand side of coordinate space at x ≥ 0. From Fig. 4(a) we see that pure free jumps do not change the form of density profile which remains to be of the trigonometric shape with the same frequency k and periodicity L. In particular, the density continues to oscillate around the same level n 0 (1 + µ 0 )/2 = 1 for any t. However, the amplitude of these oscillations decreases with increasing t, so that the particle density profile seems to become flat at long times, limt → ∞n(x,t) = 1. Again, the total number Ω n(x,t)dx of particles in Ω ∈ [−L/2, L/2] does not change in time. The repulsion makes the FIGURE 4. Density profile n(x,t) starting from trigonometric initial condition n(x, 0) = T 1,1,3 (x). The jump, repulsion, and coalescence kernels are Gaussians G 1,1 (x), G 8,1 (x), and G 0.25,1 (x), respectively. Other notations are the same as for Fig. 1. simple trigonometric (cosine) form much more complicated with the existence of additional maxima and minima inside Ω, see Fig. 4(b). Moreover, the homogeneity is being achieved here much slower than in the case of free jumps (compare density at t 5000 versus for t 10 in Fig 4(a)). This can be explained by the strong intensity (µ ϕ = 8) of repulsion potential ϕ(x) = G 8,1 (x).
For pure coalescence in Fig. 4(c), the distribution n(x,t) is no longer of a trigonometric form at t > 0 contrary to pure free jumps and like as in the case of repulsion. In addition, here the density seems to decrease to zero at t → ∞. The inclusion of repulsive jumps changes the behaviour of n(x,t) near its minima in a characteristic way. Indeed, in Fig. 4(c) these minima remain to be equal to zero, while in Fig. 4(d) they have a tendency to increase to positive values at some t. On the other hand, the speed of decrease of maxima in n(x,t) with increasing time practically does not change at t 10 despite their strength. Moreover, the shape of the density profile at long t 10 is modified as well, making it more flat with smaller oscillations. As a result, spatial homogeneity is obtained here faster.
Step initial density profiles. -A special case presents initial condition in the form of the Heaviside step function Here the system is initially (t = 0) considered on the finite interval [−L/2, L/2] with no PC. Further (t > 0) its size is gradually increased to the infinity on the unbounded interval x ∈] − ∞, ∞[ at t → ∞ according to the automatically adjusted system-size (AASS) approach. Then, a modified BC should be applied by combination of DBC and ABC together with AASS when analyzing the densities on the boundaries ±L/2. The DBC are used from the right, where lim x→∞ n(x,t) = 0 for all t. From the left, where lim x→−∞ n(x,t) = n(t) = 0 with n(0) = n 0 , we must employ ABC. This means that we should exploit the DBC rule for n i+ j given by the second equality of Eq. (10), i.e., replace n i+ j by 0 if i + j > N, and the ABC rule for n i− j given by the first equality of Eq. (11), i.e., replace n i− j by n 1 if i − j < N.
When nonzero values of n(x,t) approach the right boundary at x = L/2, i.e., when n(L/2,t) > ε max x n(x,t), the system size L is enlarged (in two times) and the simulations are continued. By monitoring from the left the difference between the actual values of n(x,t) at x = −L/2 and their homogeneous counterpart n h RK (t) [obtained by solving numerically the kinetic equation for the spatially homogeneous initial condition n(x, 0) = n 0 in parallel to the spatially inhomogeneous case using the same RK method] we can estimate the influence of the finiteness of L. If this difference exceeds the predefined level ε max x n(x,t) we should enlarge L according to the automatically adjustable system-size approach. The asymptotic value n 1 (t) will differ from the exact solution n h (t) even at very large L because of the approximate character of time integration. In the limit of sufficiently small step sizes ∆t when lim ∆t→0,L→∞ n 1 (t) = n h (t) we come to the limiting ABC-DBC scheme: lim x→−∞ n(x,t) = n h (t) and lim x→∞ n(x,t) = 0. Calculating the difference between the actual values of n(x,t) at x = −L/2 and their exact counterpart n h (t) at x → −∞ we can estimate the influence of two effects in one fashion, namely, those caused by the finiteness of L (should be significantly large) and ∆t (should be significantly small). In such a way, both ABC and DBC deviations are analyzed when deciding on the size enlargement within the ABC-DBC scheme. This completes the automatically adjustable system-size approach.
Time evolution of n(x,t) for n(x, 0) = H 1 (x) is presented in Fig. 5 using Gaussian jump G µ a ,σ a , repulsion G µ ϕ ,σ ϕ , and coalescence G µ a ,σ b kernels with different intensities of µ a = 1, µ ϕ = 8, and µ b = 0.1, respectively, as well as with equal [σ a,ϕ,b = 1, parts (a)-(d)] or different [σ a = 0.5, σ ϕ = 1, σ b = 2, part (e), and σ a = 2, σ ϕ = 1, σ b = 0.5, part (f)] ranges. As can be seen from Fig. 5(a) for pure free jumps, the discontinuous step function n(x, 0) = H 1 (x) transforms into a continuous S-shaped curve at finite times t > 0. All the curves intersect each other in the same point (0, 1/2), where 1/2 is the arithmetic mean of two initial values (1 to the left and 0 to the right). The slope of these curves becomes smaller with increasing time, so that the system tends to the mid-value everywhere. Moreover, the decrease of the amount of particles for x < 0 is equal to the corresponding increase for x > 0, that can be written as 0 −∞ [1 − n(x,t)]dx = ∞ 0 n(x,t)dx. The same statement concerns repulsive jumps, but here the curves are intersected in a point which lies below (0, 1/2). In addition, at short times t 25 the density profile n(x,t) remains to be discontinuous in x = 0, and the shape of the curves is more complicated, including the existence of maximum in x ∼ 1. The latter disappears at t 25, and n(x,t) becomes more and more flat with the flow of time. However, the flattening process is not as fast as in the case of free jumps.
For pure coalescence we can observe in Fig. 5(c) that the initial step function n(x, 0) = H 1 (x) changes to a step-like dependence n(x,t) at t > 0 containing one small peak in x ∼ −1. The density in ] − ∞, 0[ decreases from 1 at t = 0 to zero at larger t > 0. No particles appear at all in ]0, ∞[ for any t. Moreover, the initial discontinuity does not vanish even for relatively long times. The inclusion of repulsive jumps, see Fig. 5(d), prevents the appearance of peaks at x ∼ −1, while at x > 0 the profile n(x,t) is more flat when compared to that in Fig. 5(b). Here n(x,t) decreases to zero at sufficiently large times (t 100) as in the case of pure coalescence. As the range of jumps becomes shorter, σ a = 0.5, and the range of coalescence increases, σ b = 2, see Fig. 5(e), the peaks in n(x,t) appear again at x ∼ −σ b = −2 (i.e. they shift to the left with respect to those for σ b = 1). Moreover, they become more pronounced with larger amplitudes. The right lying peaks at x ∼ 1 shifts their positions to x ∼ 0.5 and decrease their amplitudes. This is contrary to the case when the jump range increases, σ a = 2, and the coalescence range decreases, σ b = 0.5, see Fig. 5(f). Then the peaks on the left do not appear, while the maximum on the right becomes more visible. In addition, the density approaches zero more slowly here. Note that for the inverse initial unit step function n(x, 0) = H 1 (−x), the results n(x,t) presented in Fig. 5 for n(x, 0) = H 1 (x) should also be inversed by n(−x,t) to obtain solutions without resolving the kinetic equation. This statement is quite general and remains in force not only for step functions, but for any other asymmetric initial conditions. Namely, when picking the initial condition n * (x, 0) = n(−x, 0), we automatically obtain n * (x,t) = n(−x,t), where n(x,t) is the solution for n(x, 0). This follows from the structure of kinetic equation (3) and the symmetricity of kernel functions.
Consider, finally, the dynamics of the initial step distribution n(x, 0) = H 1 (x) in the presence of pairs of shifted Gaussian kernels. First, we used shifting parameter s = 2 for moderate (µ a = 1) jump kernel a(x) = G 1,1,2 (x) and weak (µ b = 0.05) coalescence interaction b(x) = G 0.05,1,2 (x), while s ′ = 4 for strong (µ ϕ = 10) repulsion potential ϕ(x) = G 10,1,4 (x) (the ranges were σ a = σ b = σ ϕ = 1 for all the kernels). The corresponding results are presented in Fig. 6. From part (a) of this figure we see that for pure repulsive jumps such a choice leads at t > 0 to the emergence of self-propagating spatial inhomogeneity in the from of persistent oscillations with thin peaks of width of order of s = 2 at level of n(x) = 1 and distance between them of order of 2s ′ = 8. For n(x, 0) = H 1 (x) the inhomogeneous front propagates to the left with increasing amplitude at not too large x and to the right with dumping oscillations (for n(x, 0) = H 1 (−x) the pattern is inverse). The inclusion of coalescence even with a slight intensity of µ b = 0.05 drastically changes the situation, see Fig.  FIGURE 5. Time evolution of density profile for initial condition H 1 (x). The jump, repulsion, and coalescence kernels are Gaussians with different intensities and ranges (see the legends inside). Initially (t = 0) the system is considered on the finite interval [−10, 10] with no periodic conditions and further (t > 0) its size gradually increases to the infinity on x ∈] − ∞, ∞[ at t → ∞ according to the automatically adjusted approach. Other notations are the same as for Fig. 1.   6(b). Here, the oscillations are not so strong, especially at x > 0, and they quickly disappear with increasing time.
Further, we decreased and increased the shifting parameters in two times to values s = 1 and s ′ = 2 as well as to s = 4 and s ′ = 8. The results for these two cases  20,20] with no periodic conditions and further (t > 0) its size gradually increases to the infinity on x ∈] − ∞, ∞[ at t → ∞ according to the automatically adjusted approach. The interactions are described by the shifted pair Gaussian jump G 1,1,s , repulsion G 10,1,s ′ and coalescence G 0.05,1,2 kernels (see the legends inside).
are shown in parts (c) and (d) of Fig. 6, respectively, at the absence of coalescence. As can be seen from Fig. 6(c), the decrease of s and s ′ suppresses the processes of density propagation by lowering its speed and oscillation amplitude. Moreover, the distance between peaks in n(x,t) and their width also decrease correspondingly to 2s ′ = 4 and s = 1. This is in a contrast to the opposite case when the kernel shifts increase. Such an increase leads to stimulation of density propagation with high amplitude of the oscillations. Then the distance between peaks and their width increase to 2s ′ = 16 and s = 4, respectively. The propagation speed grows proportionally as well.
Precision of the simulations. Accuracy of our numerical simulations was measured in terms of relative absolute deviations of density values n i = n(x i ,t) in grid points x i from "exact" datan i at time t using the relation The total numberN ≤ N of grid points involved into summation (20)  The numerical uncertainties Θ cased by spatial discretization and spatial integration are shown in Fig. 7(a) as functions of the length h of space mesh (at a fixed RK4 time step of ∆t = 0.01). Note that the log-log presentation was employed to show the behaviour of Θ(h, ∆t) at small h and ∆t in more detail. Here we distinguish three kinds of errors. The first is related to uncertainties due to single spatial integrations. The corresponding results for them obtained by the composite Simpson and composite trapezoidal rules are marked by SMPS and TRPZ, respectively. We see that for sufficiently small values of h, the relative SMPS or TRPZ errors Θ are proportional to h 4 or h 2 . Indeed, the log-log plots are lines with slopes 2 or 4, i.e., log Θ = c 2 + 2 log h or log Θ = c 4 + 4 log h, where c 2 and c 4 denote some constants. This confirms the fact [25] that the composite Simpson and trapezoidal rules are accurate to the fourth O(h 4 ) and second O(h 2 ) orders in h, respectively. The accuracy is lowered significantly when considering the overall uncertainties caused by spatial discretization of the kinetic equation. Such uncertainties observed in the homogeneous and inhomogeneous interval with the Simpson integration are presented by curves labeled correspondingly as SMPSh and SMPSi. The TRPZh and TRPZi data of the trapezoidal integration appear to be approximately the same (see below the reason) and are not shown in the figure.
Maximal uncertainties are achieved in the inhomogeneous region because the dependence of n(x,t) on x makes the functions under spatial integrals more sharp, decreasing the precision of the discretization. The local uncertainties of such a discretization are of order of O(h 3 ). Taking into account that the total number of numerical operations is proportional to N 2 , the global errors appear as an accumulation of local ones leading to the overall uncertainties Θ of order of   Fig. 7). They are, however, much larger than those of the single spatial integration. Therefore, the errors caused by spatial discretization of the kinetic equation significantly prevail over the spatial integration uncertainties (the both use the same mesh h). For this reason it is not so important what the method (Simpson or trapezoidal) is applied to spatial integration. This explains why the SMPSh/SMPSi and TRPZh/TRPZi data are very close to each other.
Dependencies of Θ on the size of the time step ∆t, obtained in the simulations with the RK2, RK2', and RK4 time integrations (at a fixed spatial mesh of h = 0.01 with the involved inhomogeneous interval x ∈ [−20, 20]) are depicted in Fig. 7(b). It can be seen clearly that for sufficiently small values of ∆t, the relative deviations Θ are proportional to ∆t 2 or ∆t 4 for the algorithms of the second (RK2 and RK2') or fourth (RK4) orders, respectively. Here, the log-log plots are lines with slopes 2 or 4, i.e., log Θ(∆t) = c 2 + 2 log ∆t or log Θ(∆t) = c 4 + 4 log ∆t [similarly as for dependence of Θ on h for the single Simpson and trapezoidal spatial integrations, see Fig. 7(a)]. We mention that the local (single-step) errors of the RK4 algorithm are of order of O(∆t 5 ), see Eq. (13). So that the numerical integration over long times t ≫ ∆t leads to the global errors of order of t/∆tO(∆t 5 ) = O(∆t 4 ) as this requires repeating of the single-step integration by t/∆t times. Analogously, the RK2 and RK2' local uncertainties O(∆t 3 ) modifies to the global O(∆t 2 )-errors. Since ∆t 4 decreases with decreasing ∆t much faster than ∆t 2 , the fourth-order RK4 algorithm produces much more precise solutions than the second-order integrators RK2 and RK2'. Thus, the former can be recommended when a very high accuracy is required. For instance, the RK4 algorithm is precise up to a negligible level of order of Θ ∼ 10 −9 % at ∆t ∼ 0.1. On the other hand, the RK2 and RK2' integrators can provide here only an accuracy of 10 −3 [RK2 is slightly better than RK2', see Fig. 7(b)]. Remember, however, that the overall errors include also the spatial discretization uncertainties which are of order of 10 −2 % at h = 0.01, see Fig. 7(a). For this reason, the RK4 algorithm can be applied with much larger time steps of ∆t ∼ 1 -2, where the time-discretization uncertainties did not exceed ∼ 10 −3 and are comparable with the space-discretization errors. This significantly accelerates the simulations since the computational time is inverse proportional to ∆t. In particular, the RK4 speedup at t = 1 (when Θ ∼ 10 −5 %) with respect to the RK2 and RK2' schemes at ∆t = 0.1 (where a worse precision of Θ ∼ 10 −3 % is observed despite the smaller time step) is of order of 1/0.1/2 = 10/2 = 5 (the factor "2" takes into account the fact that RK4 requires a twice larger number of operations per step than those of RK2 or RK2' ). With further increasing ∆t, all the integrators become unstable and cannot be used at ∆t > 2, where Θ > 1%.

CONCLUSION
In this paper we have derived an efficient algorithm to obtain numerical solutions for the time-differential kinetic equation which approximates nonlocal stochastic evolution of coalescing and repulsively jumping particles in the continuous space R d . The equation is very difficult from the analytical point of view due to the presence of complicated spatial integrals with nonlinear and nonlocal terms and therefore requires numerical analysis. The proposed algorithm is based on a set of techniques including time-space discretization, periodic, Dirichlet, and asymptotic boundary conditions, composite Simpson and trapezoidal rules, Runge-Kutta schemes, and automatically adjustable system-size approaches. This has allowed as to carry out simulations of dynamics for one-dimensional systems with various initial inhomogeneous densities and different forms of the coalescence, jump and repulsion kernels, giving a comprehensive study of the jump-coalescence dynamics. A numerical error analysis of the obtained results has also been performed.
For some specific choices of the model parameters and initial densities, a nontrivial dynamics has been revealed. In the case of pure free jumps the system always tends to a nonzero homogeneous density for any initial conditions. In contrast, the presence of strong repulsion potentials can result in the appearance of persistent wave-like density propagations when the repulsion and jump kernels are chosen in the form of a sum of shifted single Gaussian functions. The shifting parameters define the shape and spatial period of density structures. The introduction of even relatively weak coalescence prevents them to persist. In the case of pure coalescence, the population eventually goes extinct, except for a special choice of the initial density profile and coalescence kernel. For instance, if the particles are initially located exclusively on an archipelago of islands and the minimal distance between them is equal to the shifting parameter of the coalescence kernel consisting of two simple rectangle functions with the ranges which do not exceed the sizes of the islands, then a stationary state with inhomogeneous particle distributions can arise. The inclusion of any jumps even with extremely small intensity radically changes the situation, leading to the collapse of the system.
The proposed algorithm is implemented in a Fortran program code which can be received by request and exploit by any researcher free of charge. The code can be adapted to more complex multicomponent models of jumping and coalescence particles of different types. The numerical simulations can be extended to systems of higher dimensions. The Poisson approximation used for the derivation of the kinetic equation can be improved by advancing to the Kirkwood [26] or Fisher-Kopeliovich [27] ansatz like for birth-death models. Mass and size of particles can also be taken into account. These and other topics as well as possible applications of the obtained results to real populations will be the subject of our further investigations.

APPENDIX
In the absence of coalescence [b(x) = 0] we can apply the convolution theorem to avoid the direct spatial integration. Indeed, then the kinetic equation (4) can be rewritten in the form ∂ n(x,t) ∂t = −n(x,t) a(x − y)g(y,t)dy + g(x,t) a(x − y)n(y,t)dy, where g(x,t) = exp − ϕ(x − u)n(u,t)du . (A2) To simplify notation consider the case d = 1 and introduce the discrete direct and inverse one-dimensional spatial Fourier transforms [28] f k = where n i = n(x i ,t) withñ k ,ã k , andφ k being the discrete Fourier transforms of n i , a i = a(x i ), and ϕ i = ϕ(x i ) according to Eq. (A3) at f i ≡ n i , a i , g i andf k ≡ n k ,ã k ,φ k . Because the kernel function a(x) and repulsive potential ϕ(x) do not depend on time, their Fourier componentsã k andφ k can be calculated once at the very beginning. Further, having the current values n i for i = 0, 1, . . . , N − 1 we compute their Fourier duplicatesñ k for k = 0, 1, . . . , N − 1. On the basis of the known products a k n k andφ kñk we evaluate their coordinate counterparts {ã kñk } −1 i and {ϕ k N k } −1 i by exploiting the inverse Fourier transform. Constructing g i = exp(−{φ kñk } −1 i ) we calculateg k by means of the direct Fourier transform and use again the inverse