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

Abstract

An algorithm is proposed for finding numerical solutions of a kinetic equation that describes an infinite system of point particles 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 numerical algorithm which we use to solve it is based on a space-time discretization, boundary conditions, composite Simpson and trapezoidal rules, Runge-Kutta methods, and adjustable system-size schemes. We show that, for special choices of the model parameters, the solutions manifest unusual time behavior. A numerical error analysis of the obtained results is also carried out.

1 Introduction

The stochastic motion of infinite populations dwelling in a continuous habitat may include 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 $\mathbb {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 in $\mathbb {R}$ up to their collision, then they coalesce and continue moving as single particles. Recently, an alternative model of this kind was proposed [7, 8]. Here, analogously as in the Kawasaki model [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]). In particular, they are used to study phytoplankton dynamics [13,14,15]. Phytoplankton cells dispersed in water constitute the basis for the vast majority of oceanic and freshwater food chains. Another example can be modeling melanoma cells migrating and coalescing to form tumors [16, 17].

Quite recently [18], first numerical results on the jump-coalescence kinetic equation were reported, where the role of repulsion in the appearance of a spatial heterogeneity was elucidated. There it was shown also how the kinetic equation can rigorously be obtained from the corresponding microscopic evolution equations obtained and studied in [8]. However, in [18], the numerical consideration was restricted to a few particular examples of the model parameters. Moreover, very little was said about the algorithm used to obtain the numerical solutions.

In the present paper, we derive an algorithm enabling to solve the repulsion-jump-coalescence kinetic equation. It combines different techniques and methods allowing for studying infinite systems on the basis of their 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.

2 A kinetic equation of a population model

In a broader sense, a kinetic equation is an integro-differential equation describing the evolution in time of the density (or similar aggregated quantities) of a large population. This scale of description is often called mesoscopic as a bridge between microscopic and macroscopic ones. Unlike the microscopic description, at the mesoscale, the individual particles and their interactions are not considered directly. However, it still allows to efficiently describe the local structure of the population. The idea of such description, probably best known from the famous Boltzmann equation, is employed in different areas, such as vehicular traffic, human crowds (see, e.g., [19]), or ecological systems (see, e.g., [20]). Usually, kinetic equations are introduced by phenomenological/heuristic arguments. Recently, considerable attention is directed to multiscale modeling based upon a consistent description at different scales. The kinetic equation studied in this work is derived by a certain procedure from the corresponding chain of evolution equations obtained at a microscopic level (see [8, 18] for more detail).

We consider an infinite population dwelling in $\mathbb {R}^{d}$, d ≥ 1, that undergoes the following random events: 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]:

$$ \begin{array}{@{}rcl@{}} \frac{\partial n(x,t)}{\partial t} &= &- \int a(x,y) \exp \left( - \int {\varPhi}(y-u) n(u,t) du \right) n(x,t) dy \\ &&+ \int a(x,y) \exp \left( - \int {\varPhi}(x-u) n(u,t) du \right) n(y,t) dy \\ &&- \int \int \left[ b(x,y,z) + b(y,x,z) \right] n(x,t) n(y,t) dy dz \\ &&+ \int \int b(y,z,x) n(y,t) n(z,t) dy dz , \end{array} $$

(1)

where a,Φ,b ≥ 0. The first term in the rhs of (1) represents jumping from x to y with intensity a(x,y) modulated by an exponential factor. This results in a decrease of n(x,t). The exponential factor corresponds to the repulsion of the particle jumping to y caused by the rest of the system, where Φ(y − u) denotes the corresponding repulsion 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 describes free coalescence where 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 (1) are performed in $\mathbb {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 \leftrightarrow 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 ${\int \limits } a(x, y)dx = {\int \limits } a(x, y)dy = \mu _{a}$, ${\int \limits } {\varPhi }(x) = \mu _{{\varPhi }}$, and ${\int \limits } b(x_{1}, x_{2}, x_{3} )dx_{i} dx_{j} = \mu _{b}$, where i,j = 1,2,3 with i≠j. 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|).

The coalescence intensity will be taken in the form:

$$ b(x,y,z) = b(x-y) \delta((x + y) / 2 - z), $$

(2)

where b(x − y) = b(|x − y|) and δ is the Dirac function. That is, the resulting point is located in the middle z = (x + y)/2. This is quite natural for identical particles. Instead of the δ-function, one can choose smoother functions, e.g., Gaussian-like ones.

The δ-function is used here to simplify the calculations as follows. By integrating over z one transforms kinetic equation (1) to:

$$ \begin{array}{@{}rcl@{}} \frac{\partial n(x,t)}{\partial t} &= &- \int a(x-y) \exp \left( - \int {\varPhi}(y-u) n(u,t) du \right) n(x,t) dy \\ &&+ \int a(x-y) \exp \left( - \int {\varPhi}(x-u) n(u,t) du \right) n(y,t) dy \\ &&- \int \left( 2 b(x-y) n(x,t) -2^{d} b(2(x-y)) n(2x - y, t) \right) n(y,t) dy,\\ \end{array} $$

(3)

where the symmetric properties of the kernels have been taken into account. Imposing an initial condition n(x,0), (3) leads to a complicated partial integro-differential equation with respect to the unknown function n(x,t) at t > 0. Unfortunately, we cannot apply the convolution theorem to avoid the spatial integration in the last coalescence term of (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.

3 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 $\mathbb {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 (3) is:

$$ \begin{array}{@{}rcl@{}} \frac{dn_{i}}{dt} &= &h \sum\limits_{j} \left( a_{i-j} \exp \left[ -h \sum\limits_{k} {\varPhi}_{i-k}n_{k} \right] n_{j} - a_{i-j} \exp \left[ -h \sum\limits_{k} {\varPhi}_{j-k} n_{k} \right] n_{i} \right) \\ &&- 2h \sum\limits_{j} b_{i-j} n_{i} n_{j} + 2h \sum\limits_{j} b_{2i - 2j} n_{j} n_{2i - j}, \end{array} $$

(4)

where n_i ≡ n_i(t) = n(x_i,t) with a_i−j = a(x_i − x_j) = a((i − j)h), b_i−j = b(x_i − x_j) = b((i − j)h), Φ_i−k = Φ(x_i − x_k) = Φ((i − k)h), and the infinite sums over j and k represent the spatial integrals. It is obvious that in the limit $h \rightarrow 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:

$$ \lambda_{i} = \exp \left[ -h \sum\limits_{k} {\varPhi}_{i-k} n_{k} \right], $$

(5)

one obtains that (4) can be cast more compactly in the following equivalent form:

$$ \frac{dn_{i}}{dt} = h \sum\limits_{j} \left( a_{j} (\lambda_{i} n_{i-j} - \lambda_{i-j} n_{i}) - 2b_{j} n_{i} n_{i-j} + 2b_{2j} n_{i-j} n_{i+j} \right) $$

(6)

where a_j = a(jh), b_j = b(jh), and b_2j = b(2jh).

In computer simulations, we cannot operate with infinite-size samples leading to the infinite summation over j in (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 \in ] - \infty , \infty [$ by finite area [−L/2,L/2]. In view of the aforesaid, (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} = \pm (L/2 - {\max \limits } \sigma )$, where ${\max \limits } \sigma$ is the maximal radius of the kernels (see Section 4). In numerical calculations, this can be expressed by the condition $n(x_{B}, t) < \varepsilon \ \max \limits _{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) > \varepsilon \ \max \limits _{x} n(x,t)$, we should enlarge L (e.g., gradually doubling it) until the required first condition is satisfied, use DBC again, and continue the simulation for t > T. A special case, which is discussed in Section 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 (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^{\prime }$ 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^{\prime }} = x_{j} \pm KL$, where $K = 1, 2, \ldots , \infty$ and $n_{j^{\prime }} = 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 interparticle distance |x_j|, the summation over j in (6) can be actually restricted to a finite number of terms for which |x_j|≤ R_a,b < L/2. The truncation radii 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^{\prime }} = n(t)$ for all $x_{j^{\prime }} < -L/2$ and $n_{j^{\prime }} = n(t)$ for all $x_{j^{\prime }} > 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^{\prime }} = n(-L/2,t)$ for all $x_{j^{\prime }} < -L/2$ while $n_{j^{\prime }} = n(L/2,t)$ for all $x_{j^{\prime }} > 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, (6) transforms to:

$$ \begin{array}{@{}rcl@{}} \frac{dn_{i}}{dt} &= &h \sum\limits_{j=1}^{j_{a}} \xi_{j}^{(a)} a_{j} \left( \lambda_{i} (n_{i-j} + n_{i+j}) - n_{i} (\lambda_{i-j} + \lambda_{i+j}) \right) \\ &&- 2h \xi_{0}^{(b)} b_{0} {n_{i}^{2}} -2h n_{i} \sum\limits_{j=1}^{j_{b}} \xi_{j}^{(b)} b_{j} (n_{i-j} + n_{i+j}) \\ &&+ 2h \xi_{0}^{(2b)} b_{0} {n_{i}^{2}} +4h \sum\limits_{j=1}^{j_{b}/2} \xi_{j}^{(2b)} b_{2j} n_{i-j} n_{i+j} \end{array} $$

(7)

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 (5) that:

$$ \lambda_{i} = \exp \left[ -h \sum\limits_{j}^{BC} \xi_{j}^{({\varPhi})} {\varPhi}_{j} n_{i-j} \right] = \exp \left[ -h \xi_{0}^{({\varPhi})} {\varPhi}_{0} n_{i} - h \sum\limits_{j=1}^{j_{{\varPhi}}} \xi_{j}^{({\varPhi})} {\varPhi}_{j} (n_{i-j} + n_{i+j}) \right], $$

(8)

where i = 1,2,…,N and j_Φ = R_Φ/h with truncation radius RΦ < L/2. It is evident that in the limits $L, N, R \rightarrow \infty$ provided $h \rightarrow 0$, the discretized equation (7) with (8) coincide with its original, continuous counterpart (3).

Weights $\xi _{j}^{(a,b,2b,{\varPhi })}$ at kernel values a_j,b_j,b_2j, and Φ_j were introduced in (7) and (8) to improve precision of the numerical integration over coordinate space. They are determined according to the chosen method [21]. These weights satisfy the normalization condition $\xi _{0} + 2 {\sum }_{j=1}^{j^{*}} \xi _{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 $\xi _{0} = \xi _{1} = \xi _{2} = {\ldots } = \xi _{j^{*}-1} = 1$, while $\xi _{j^{*}} = 1/2$. The composite Simpson method of the fourth order yields $\xi _{j^{*}} = 1/3, \xi _{j^{*} - 1} = 4/3, \xi _{j^{*} - 2} = 2/3, \xi _{j^{*} - 3} = 4/3, \xi _{j^{*} - 4} = 2/3, \ldots , \xi _{0} = 4/3$ if j^∗ is odd and ξ₀ = 2/3 if j^∗ is even. The latter is more accurate than the former and involve numerical uncertainties of order of ${\mathscr{O}}(h^{4})$ versus ${\mathscr{O}}(h^{2})$.

We mention that n_i are explicitly defined at i = 1,2,…,N, so that the corresponding BC should be applied in (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 x_i = −L/2 + (i − 1/2)h, where i = 1,2,…,N with even N. 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:

$$ n_{i-j} = \left\{ \begin{array}{ll} n_{i-j+N}, & i - j < 1\\ n_{i-j}, & 1 \leq i - j \leq N \end{array} \right. ,\quad n_{i+j} = \left\{ \begin{array}{ll} n_{i+j-N} & i + j > N\\ n_{i+j} & 1 \leq i - j \leq N \end{array} \right. . $$

(9)

Note that i = 1,2,…,N and j = 1,2,…,j^∗ < N/2. The applications of DBC and ABC result in:

$$ n_{i-j} = \left\{ \begin{array}{ll} 0, & i - j < 1\\ n_{i-j}, & 1 \leq i - j \leq N \end{array} \right. ,\quad n_{i+j} = \left\{ \begin{array}{ll} 0 & i + j > N\\ n_{i+j} & 1 \leq i - j \leq N \end{array} \right. . $$

(10)

and

$$ n_{i-j} = \left\{ \begin{array}{ll} n_{1}, & i - j < 1\\ n_{i-j}, & 1 \leq i - j \leq N \end{array} \right. ,\quad n_{i+j} = \left\{ \begin{array}{ll} n_{N} & i + j > N\\ n_{i+j} & 1 \leq i - j \leq N \end{array} \right. . $$

(11)

respectively. In the cases d = 2 and d = 3, the boundary transformations can be implemented similarly to those given by (9)–(11).

Time integration

In the most general form, our coupled system of N autonomous ordinary differential equations can be given as:

$$ \frac{d n_{i}}{d t} = \dot{n}_{i} = s_{i}(n_{1}, n_{2}, {\ldots} n_{N}), $$

(12)

where i = 1,2,…,N and s_i(n₁,n₂,…,n_N) represents the rhs of (7) with taking into account (8). Introducing the set Γ(t) = {n_i(t)} of N dynamical variables and their time derivatives Φ(t) = {s_i(t)}, (12) can be compactly cast as $\dot {{\varGamma }} = {\varPhi }({\varGamma })$. A common practice to solve differential equations of such a type is to use classical Runge-Kutta schemes [21]. The scheme of the fourth order (RK4) reads:

$$ {\varGamma}(t + {\varDelta} t) = {\varGamma}(t) + \frac{{\varDelta} t}{6} \left( {\varPhi}_{1} + 2 {\varPhi}_{2} + 2 {\varPhi}_{3} + {\varPhi}_{4}\right) + \mathscr{O}({\varDelta} t^{5}), $$

(13)

where Δt is the time step and Φ₁ = Φ(Γ(t)), Φ₂ = Φ(Γ(t) + Φ₁Δt/2), Φ₃ = Φ(Γ(t) + Φ₂Δt/2), Φ₄ = Φ(Γ(t) + Φ₃Δt) are the derivatives in intermediate stages. 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 ${\varGamma } (t + {\varDelta } t) = {\varGamma } (t) + [{\varPhi }({\varGamma } (t)) + {\varPhi }({\varGamma } (t) + {\varPhi } ({\varGamma }(t)) {\varDelta } t)] {\varDelta } t/2 + {\mathscr{O}}({\varDelta } t^{3})$. The second version (RK2’) relates to a middle point scheme, ${\varGamma } (t + {\varDelta } t) = {\varGamma } (t) + {\varPhi }({\varGamma } (t) + {\varPhi }({\varGamma } (t)) {\varDelta } t/2) {\varDelta } t + {\mathscr{O}}({\varDelta } t^{3})$. 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 ${\mathscr{O}}({\varDelta } t^{5})$ versus ${\mathscr{O}}({\varDelta } 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 (13). The ${\mathscr{O}}({\varDelta } 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². 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.

4 Performed simulations and analysis of the results

The results presented in this section show how the discussed model describes a variety of interesting phenomena when its parameters are changed. Probably the most spectacular one is the appearance of the persisting heterogeneity (see Fig. 6). One of its possible applications is the description of the stochastic dynamics of a pelagic marine population [11] or phytoplankton [13,14,15]. Important feature of the discussed model is that the population constituents are indistinguishable apart from the described treat $x \in \mathbb {R}^{d}$, which can be e.g. their spatial location or in the case of d = 1 also their mass. It is important to choose appropriate kernels to fit the desired application. In general, they do not need to be symmetric, as we assumed within this work. We present some simple cases that shed light on how the dynamics depend on the parameter functions, as well as on the initial conditions.

Numerical details

The numerical simulations were carried out for one-dimensional 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 $\mathbb {R}$ was done with a mesh of h = 0.0125 – 0.2 dependent 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 $\varepsilon \sim 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:

$$ G_{\mu, \sigma} (x) = \frac{\mu}{(2 \pi \sigma^{2})^{\frac{1}{2}}} \exp \left( - \frac{x^{2}}{2 \sigma^{2}} \right) $$

(14)

or rectangle

$$ C_{\mu, \sigma} (x) = \left\{ \begin{array}{lll} \frac{\mu}{2\sigma}, & & |x| \leq \sigma\\ 0, & & |x| > \sigma \end{array} \right. $$

(15)

functions, where μ = μ_a, μ_b, or μ_Φ and σ = σ_a, σ_b, or σ_Φ are the intensities and ranges of the corresponding interactions, respectively. The kernels are normalized, ${\int \limits } \{a, b, {\varPhi }\}(x)dx = \mu _{a}, \mu _{b}, \mu _{{\varPhi }}$. A symmetrical pair of shifted Gaussian or rectangle functions was involved as well:

$$ F_{\mu, \sigma, s}(x) = \frac{1}{2} \left( F_{\mu, \sigma}(x-s) + F_{\mu, \sigma}(x+s) \right), $$

(16)

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_{\mu ,\sigma } (Q \sigma )/G_{\mu ,\sigma } (0) \sim 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 the 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₀(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.

Rectangle initial density profiles

The first example relates to the initial condition in the form of periodic rectangle function n(x,0) = C_1,1,L(x). The infinite system is reproduced by repeating the single rectangle segment C_1,1(x) at x ∈Ω = [−L/2,L/2] on the interval $]-\infty , \infty [$ with a period of L and applying PBC. The jump $a(x) = G_{\mu _{a}, \sigma _{a}}(x)$, repulsion ${\varPhi }(x) = G_{\mu _{{\varPhi }},\sigma _{{\varPhi }}} (x)$, and coalescence $b(x) = G_{\mu _{b}, \sigma _{b}} (x)$ kernels are modelled by the Gaussians with intensities μ_a = 1, μ_Φ = 20, and μ_b = 1. Three sets of kernel ranges were considered, namely, σ_a = σ_Φ = σ_b = 1, σ_a < σ_Φ < σ_b with σ_a = 0.5, σ_Φ = 1, σ_b = 2 and vice versa, σ_a > σ_Φ > σ_b with σ_a = 2, σ_Φ = 1, σ_b = 0.5. The corresponding time evolution of spatial structure n(x,t) is presented in Fig. 1 for the cases of pure free jumps, repulsive jumps, pure coalescence, and repulsive jumps with coalescence with σ_a = σ_Φ = σ_b = 1, parts (a), (b), (c), and (d), respectively, as well, with σ_a < σ_Φ < σ_b and σ_a > σ_Φ > σ_b for coalescing repulsive jumps, parts (e) and (f).

Fig. 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)

From Fig. 1a, 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 \gtrsim 10$. For longer times, $t \gtrsim 40$, the solution n(x,t) tends to more and more homogeneous densities. In the limit $t \rightarrow \infty$, we expect absolutely flat function $\lim _{t \rightarrow \infty } n(x,t) = n$ independent on x, where $n = \frac {1}{L} {\int \limits }_{{\varOmega }} n(x, 0) dx = 1/L = 0.05$. Moreover, for any time t, the total number ${\int \limits }_{{\varOmega }} 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 \sim 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 \lesssim 40$ in the ranges x ≈± 2 due to the repulsion between particles. We believe that all the maxima disappear at $t \rightarrow \infty$ with the same asymptotic behaviour $\lim _{t \rightarrow \infty } 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 \rightarrow \infty$. 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 \gtrsim 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 become more visible with respect to those for equal-range 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 Fig. 1a–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 ${\mathscr{N}}_{0}$ shifted single rectangle functions $C_{\nu _{k}, \sigma _{k}} (x + s_{k})$ with intensities ν_k and ranges σ_k, namely:

$$ n(x,0) = C_{\nu_{k},\sigma_{k}, s, \mathscr{N}_{0}} (x) = \sum\limits_{k=1}^{\mathscr{N}_{0}} C_{\nu_{k}, \sigma_{k}} (x + s_{k}) $$

(17)

where $s_{k} = (-{\mathscr{N}}_{0} + 2k - 1)s$ are shifting parameters. We use $s = L/(2 {\mathscr{N}}_{0})$ and therefore $s_{k} = -L/2 + (k - 1/2)L/{\mathscr{N}}_{0}$. Then, repeating (17) with period L, we should apply PBC to deal with the infinite system. We used a particular case of (17) with ${\mathscr{N}}_{0} = 3$ and L = 20 as well as three different amplitudes ν_1,2,3 generated at random in the interval ]0,1[. The corresponding result on this is shown in Fig. 2.

Fig. 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

Looking at Fig. 2 and comparing it with Fig. 1, we see that behavior 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_{\nu _{k}, \sigma _{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 \rightarrow \infty$ flattens in x and seems tend to the nonzero constant $n = \frac {1}{L} {\int \limits }_{{\varOmega }} n(x, 0) dx = (\nu _{1} + \nu _{2} + \nu _{3})/L$ for any initial distributions. We mention that at μ_b = 0, the total amount ${\int \limits }_{{\varOmega }} n(x,t)dx$ of particles in Ω is unchanged and equal to its initial value ${\int \limits }_{{\varOmega }} n(x, 0)dx$. For μ_b > 0, the density function n(x,t) seems to take its zeroth asymptotics $\lim _{t \rightarrow \infty } n(x,t) = 0$ at $t \rightarrow \infty$ [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, \ldots , \infty$, 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 ${\int \limits }_{{\varOmega }} n(x,t)dx \rightarrow 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 [(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_{\nu _{k},\sigma _{k},s',{\mathscr{N}}_{0}} (x)$ of two (${\mathscr{N}}_{0} = 2$) single rectangle functions [ (17)] with the same amplitude ν_k = μ_b = 1 and ranges σ_k = σ_b ≡ σ = 1 as those of the coalescence kernel b(x) but with different shifting parameter $s^{\prime } \neq s$ satisfying the constraint $s^{\prime } = 2s - 2 \sigma$. Choosing s = 5 one finds at σ = 1 that $s^{\prime } = 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. 3a. 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, 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 \gtrsim 2 \cdot 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^{\prime } - \sigma = 7 < |x| < 9 = s^{\prime } + \sigma$, so that the absence of coalescence at a given 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. 3b 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 \lesssim 4$ achieving a maximum at $t \sim 4$, and then decreases at $t \gtrsim 4$.

Fig. 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

Trigonometric initial density profiles

Another interesting case is the initial density in the form of a trigonometric function:

$$ n(x,0) = T_{n_{0}, \mu_{0}, k} (x) = n_{0} \left( 1 + \mu_{0} \cos (2 \pi k x / L) \right), $$

(18)

where 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₀ = 1, μ₀ = 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.

Fig. 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

From Fig. 4a, 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₀(1 + μ₀)/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, $\lim _{t \rightarrow \infty } n(x,t) = 1$. Again, the total number ${\int \limits }_{{\varOmega }} n(x,t)dx$ of particles in Ω ∈ [−L/2,L/2] does not change in time. The repulsion makes the simple trigonometric (cosine) form much more complicated with the existence of additional maxima and minima inside Ω; see Fig. 4b. Moreover, the homogeneity is being achieved here much slower than in the case of free jumps (compare density at $t \gtrsim 5000$ versus for $t \gtrsim 10$ in Fig 4a). This can be explained by the strong intensity (μ_Φ = 8) of repulsion potential Φ(x) = G_8,1(x).

For pure coalescence in Fig. 4c, 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 \rightarrow \infty$. The inclusion of repulsive jumps changes the behavior of n(x,t) near its minima in a characteristic way. Indeed, in Fig. 4c, these minima remain to be equal to zero, while in Fig. 4d 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 \lesssim 10$ despite their strength. Moreover, the shape of the density profile at long $t \gtrsim 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:

$$ n(x,0) = H_{n_{0}}(x) = \left\{ \begin{array}{lll} n_{0}, & & x \leq 0\\ 0, & & x > 0 \end{array} \right. . $$

(19)

Here, the system is initially (t = 0) considered on the finite interval [−L/2,L/2] with no PBC. Further (t > 0) its size is gradually increased to the infinity on the unbounded interval $x \in ]-\infty , \infty [$ at $t \rightarrow \infty$ 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 \rightarrow \infty } n(x,t) = 0$ for all t. From the left, where $\lim _{x \rightarrow - \infty } n(x,t) = n(t) \neq 0$ with n(0) = n₀, we must employ ABC. This means that we should exploit the DBC rule for n_i+j given by the second equality of (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 (11), i.e., replace n_i−j by n₁ 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) > \varepsilon \max \limits _{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₀ 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 $\varepsilon \max \limits _{x} n(x,t)$ we should enlarge L according to the automatically adjustable system-size approach. The asymptotic value n₁(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 _{{\varDelta } t \rightarrow 0, L \rightarrow \infty } n_{1} (t) = n^{h} (t)$, we come to the limiting ABC-DBC scheme: $\lim _{x \rightarrow - \infty } n(x,t) = n^{h} (t)$ and $\lim _{x \rightarrow \infty } 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 \rightarrow - \infty$, 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₁(x) is presented in Fig. 5 using Gaussian jump $G_{\mu _{a},\sigma _{a}}$, repulsion $G_{\mu _{{\varPhi }}, \sigma _{{\varPhi }}}$, and coalescence $G_{\mu _{a}, \sigma _{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. 5a for pure free jumps, the discontinuous step function n(x,0) = H₁(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 ${\int \limits }_{-\infty }^{0}[1 - n(x,t)]dx = {\int \limits }_{0}^{\infty } 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 \lesssim 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 \sim 1$. The latter disappears at $t \gtrsim 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.

Fig. 5

Time evolution of density profile for initial condition H₁(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 \in ]-\infty , \infty [$ at $t \rightarrow \infty$ according to the automatically adjusted approach. Other notations are the same as for Fig. 1

For pure coalescence we can observe in Fig. 5c that the initial step function n(x,0) = H₁(x) changes to a step-like dependence n(x,t) at t > 0 containing one small peak in $x \sim -1$. The density in $]-\infty , 0[$ decreases from 1 at t = 0 to zero at larger t > 0. No particles appear at all in $]0, \infty [$ for any t. Moreover, the initial discontinuity does not vanish even for relatively long times. The inclusion of repulsive jumps, see Fig. 5d, prevents the appearance of peaks at $x \sim -1$, while at x > 0 the profile n(x,t) is more flat when compared to that in Fig. 5b. Here n(x,t) decreases to zero at sufficiently large times $(t \gtrsim 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. 5e), the peaks in n(x,t) appear again at $x \sim -\sigma _{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 \sim 1$ shift their positions to $x \sim 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. 5f). 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₁(−x), the results n(x,t) presented in Fig. 5 for n(x,0) = H₁(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₁(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^{\prime } = 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^{\prime } = 8$. For n(x,0) = H₁(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₁(−x) the pattern is inverse). The inclusion of coalescence even with a slight intensity of μ_b = 0.05 drastically changes the situation (see Fig. 6b). Here, the oscillations are not so strong, especially at x > 0, and they quickly disappear with increasing time.

Fig. 6

Dynamics of the system starting from unit step function H₁(x) in the cases of pure repulsive jumps with different kernel shifts: as = 2, $s^{\prime } = 4$; c s = 1, $s^{\prime } = 2$; and d s = 4, $s^{\prime } = 8$; as well as of b repulsive jumps and coalescence with s = 2 and $s^{\prime } = 4$. Initially (t = 0) the system is considered on the finite interval [− 20,20] with no periodic conditions and further (t > 0) its size gradually increases to the infinity on $x \in ]-\infty , \infty [$ at $t \rightarrow \infty$ 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^{\prime }}$ and coalescence G_0.05,1,2 kernels (see the legends inside)

Furthermore, we decreased and increased the shifting parameters two times to values s = 1 and $s^{\prime } = 2$ as well as to s = 4 and $s^{\prime } = 8$. The results for these two cases are shown in parts (c) and (d) of Fig. 6, respectively, at the absence of coalescence. As can be seen from Fig. 6c, the decrease of s and $s^{\prime }$ 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^{\prime } = 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^{\prime } = 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” data $\breve {n}_{i}$ at time t using the relation:

$$ \varTheta(h, {\varDelta} t) = \frac{{\sum}_{i=1}^{\breve{N}} |n_{i} - \breve{n}_{i}|}{{\sum}_{i=1}^{\breve{N}} \breve{n}_{i}} \times 100 [<percent>]. $$

(20)

The total number $\breve {N} \leq N$ of grid points involved in summation (20) depends on the spatial region considered. The “exact” (or rather reference) values $\breve {n}_{i}$ were obtained at high enough space and time resolutions with a tiny mesh of h = 0.005 and a tiny time step of Δt = 0.01 in order to be entitled to ignore the numerical uncertainties. The spatial and time integrations were performed with the help of the composite Simpson rule and RK4 algorithm, respectively. The numerical error analysis was done by carrying out a series of simulations at different h = 0.01,0.02,0.04,0.1,0.2 and Δt = 0.1,0.2,0.4,1.0,2.0,2.5. The numerical deviations were then estimated by (20) to plot Θ(h,Δt) in a wide range of varying h and Δt. For the purpose of comparison, the simulation runs with the composite trapezoidal rule and RK2 algorithm were performed as well. Error results presented below are related to one of the situations considered in Figs. 1, 2, 3, 4, 5, and 6 when choosing initial condition n(x,0), namely, in the case of Fig. 5d at t = 50. There, the homogeneous and inhomogeneous intervals were x ∈ [− 80,− 20] and x ∈ [− 20,20], respectively. Similar results were observed for all the rest choices of n(x,0).

The numerical uncertainties Θ cased by spatial discretization and spatial integration are shown in Fig. 7a 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 behavior 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⁴ or h². Indeed, the log-log plots are lines with slopes 2 or 4, i.e., $\log \varTheta = c_{2} + 2 \log h$ or $\log \varTheta = c_{4} + 4 \log h$, where c₂ and c₄ denote some constants. This confirms the fact [21] that the composite Simpson and trapezoidal rules are accurate to the fourth ${\mathscr{O}}(h^{4})$ and second ${\mathscr{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.

Fig. 7

a Uncertainties Θ of the simulations related to single composite Simpson (SMPS) or trapezoidal (TRPZ) spatial integrations and the overall spatial discretization obtained in homogeneous (SMPSh) and inhomogeneous (SMPSi) regions at different mesh h; b Dependencies of Θ on the size of the time step Δt for the RK2, RK2’, and RK4 time integrations

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 ${\mathscr{O}}(h^{3})$. Taking into account that the total number of numerical operations is proportional to N², the global errors appear as an accumulation of local ones leading to the overall uncertainties Θ of order of $N^{2} {\mathscr{O}}(h^{3}) \sim {\mathscr{O}}(h)$, since h = L/N. In other words, the SMPSh and SMPSi deviations behave as $\log \varTheta = c_{1} + \log h$. Despite the linear dependence ${\mathscr{O}}(h)$, even the SMPSi errors are relatively small and do not exceed about 0.01–0.03% at h = 0.01–0.02 (see 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. 7b. It can be seen clearly that for sufficiently small values of Δt, the relative deviations Θ are proportional to Δt² or Δt⁴ 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 \varTheta ({\varDelta } t) = c_{2} + 2 \log {\varDelta } t$ or $\log \varTheta ({\varDelta } t) = c_{4} + 4 \log {\varDelta } t$ (similarly as for dependence of Θ on h for the single Simpson and trapezoidal spatial integrations, see Fig. 7a). We mention that the local (single-step) errors of the RK4 algorithm are of order of ${\mathscr{O}}({\varDelta } t^{5})$ (see (13)). So that the numerical integration over long times t ≫Δt leads to the global errors of order of $t / {\varDelta } t {\mathscr{O}} ({\varDelta } t^{5}) = {\mathscr{O}} ({\varDelta } t^{4})$ as this requires repeating of the single-step integration by t/Δt times. Analogously, the RK2 and RK2’ local uncertainties ${\mathscr{O}}({\varDelta } t^{3})$ modify the global ${\mathscr{O}}({\varDelta } t^{2})$-errors. Since Δt⁴ decreases with decreasing Δt much faster than Δt², 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 $\varTheta \sim 10^{-9}\%$ at ${\varDelta } t \sim 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. 7b). 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. 7a). For this reason, the RK4 algorithm can be applied with much larger time steps of ${\varDelta } t \sim 1$ – 2, where the time-discretization uncertainties did not exceed $\sim 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 $\varTheta \sim 10^{-5} \%$) with respect to the RK2 and RK2’ schemes at Δt = 0.1 (where a worse precision of $\varTheta \sim 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%.

5 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 $\mathbb {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 us 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 from persisting. 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 exploited 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 [22] or Fisher-Kopeliovich [23] 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:

$$ \frac{\partial n(x, t)}{\partial t} = - n(x, t) \int a(x - y) g(y, t) dy + g(x, t) \int a(x - y) n(y, t) dy, $$

(A1)

where

$$ g(x, t) = \exp \left( - \int {\varPhi} (x - u) n(u, t) du \right) . $$

(A2)

To simplify notation, consider the case d = 1 and introduce the discrete direct and inverse one-dimensional spatial Fourier transforms [24]:

$$ \tilde{f}_{k} = \sum\limits_{i=0}^{N-1} f_{i} \exp (-2 \pi \imath k i / N) \equiv \{f_{i}\}_{k}, \quad f_{i} = \frac{1}{N} \sum\limits_{k=0}^{N-1} \tilde{f}_{k} \exp (2 \pi \imath k i / N) \equiv \{f_{k}\}_{i}^{-1} $$

(A3)

where f_i = f(x_i) and x_i = −L/2 + (i − 1/2)h ∈] − L/2,L/2[ with even N and h = L/N. Now, applying the convolution theorem, the discrete counterpart of (A1) in view of (A2) can be presented as:

$$ \dot{n}_{i} = \frac{dn_{i}}{dt} = -h n_{i} \{\tilde{a}_{k} \tilde{g}_{k}\}_{i}^{-1} + h g_{i} \{\tilde{a}_{k} \tilde{n}_{k}\}_{i}^{-1}, \quad g_{i} = \exp (-h \{ \tilde{{\varPhi}}_{k} \tilde{n}_{k} \}_{i}^{-1}) $$

(A4)

where n_i = n(x_i,t) with $\tilde {n}_{k}$, $\tilde {a}_{k}$, and $\tilde {{\varPhi }}_{k}$ being the discrete Fourier transforms of n_i,a_i = a(x_i), and Φ_i = Φ(x_i) according to (A3) at f_i ≡ n_i,a_i,g_i and $\tilde {f}_{k} \equiv \tilde {n}_{k}, \tilde {a}_{k}, \tilde {{\varPhi }}_{k}$.

Because the kernel function a(x) and repulsive potential Φ(x) do not depend on time, their Fourier components $\tilde {a}_{k}$ and $\tilde {{\varPhi }}_{k}$ can be calculated once at the very beginning. Furthermore, having the current values n_i for i = 0, 1,…,N − 1, we compute their Fourier duplicates $\tilde {n}_{k}$ for k = 0, 1,…,N − 1. On the basis of the known products $\tilde {a}_{k} n_{k}$ and $\tilde {{\varPhi }}_{k} \tilde {n}_{k}$, we evaluate their coordinate counterparts $\{\tilde {a}_{k} \tilde {n}_{k} \}_{i}^{-1}$ and $\{{\varPhi }_{k} N_{k}\}_{i}^{-1}$ by exploiting the inverse Fourier transform. Constructing $g_{i} = \exp (-\{\tilde {{\varPhi }}_{k} \tilde {n}_{k} \}_{i}^{-1})$, we calculate $\tilde {g}_{k}$ by means of the direct Fourier transform and use again the inverse transform to obtain $\{\tilde {a}_{k} \tilde {g}_{k} \}_{i}^{-1}$. In such a way, we have all the necessary components to form the time derivatives $\dot {n}_{i}$ according to (A4). Despite this requiring several direct and inverse transformations, the order of total number of operations can be reduced from N² (when the spatial integration is performed directly, see Section 3) to $N \ln N$ when exploiting the fast discrete direct and inverse Fourier transforms. This is a very important feature because N ≫ 1 is large (in practice $N \sim 10^{2}$ – 10³ at least). However, such an approach can work not so well for discontinuous functions. Moreover, it assumes that all involved functions are periodic on the interval [−L/2,L/2].

The generalization to any higher dimensionality d > 1 can also be done by replacing i and k on the d-dimensional vectors i = (i₁,i₂,…,i_d) and k = (k₁,k₂,…,k_d) as well as their multiplication ki on the scalar product k ⋅i with summation in (A3) from 0 to N − 1 for each dimension, where $f_{\mathbf {i}} = f(x_{i_{1}}, x_{i_{2}}, \ldots , x_{i_{d}})$. As a result, the total number of operations increases from N² to N^2d or from $N \ln N$ to $N^{d} \ln N^{d}$ when applying the product of d one-dimensional usual or fast Fourier transform in Cartesian coordinates. If the functions possess radial symmetry, we can consider the Fourier transform in polar or spherical coordinates at d = 2 or 3, for instance. This allows reducing d-dimensional integrations to discrete Hankel [25] or Fourier-Bessel [26] transforms in one (radial) dimension by integrating analytically out all angle dependencies. The chief advantage of this approach is the fact that then the total number of operations becomes independent of d and remains the same (of order of $N \ln N$) for any dimensionality [27], just as at d = 1.

In the spatially homogeneous case, i.e., when n(x,t) = n(t) does not change on x, we can carry out the spatial integrations in (4) and (5) over y explicitly. This significantly simplifies the kinetic equations to the form dn/dt = −μ_bn² which gives the analytical solution n(t) = n(0)/[1 + μ_bn(0)t]. It tends to zero with $t \rightarrow \infty$ for any initial density n(0) and μ_b≠ 0. In this case, the jump contribution completely disappears. At μ_b = 0, this leads to the time-independent solution n(t) = n(0). It disappears also in the spatially inhomogeneous case for systems with finite number ${\mathscr{N}} = \int \limits n(x)dx$ of particles when calculating $d {\mathscr{N}} /dt$ (then during integration in the rhs of (3) over x, the first two terms are mutually cancelled).