Algorithm for numerical solutions to the kinetic equation of a spatial population dynamics model with coalescence and repulsive jumps
aa r X i v : . [ m a t h . D S ] J u l ALGORITHM FOR NUMERICAL SOLUTIONS TO THE KINETICEQUATION OF A SPATIAL POPULATION DYNAMICS MODEL WITHCOALESCENCE AND REPULSIVE JUMPS
IGOR OMELYAN, YURI KOZITSKY, AND KRZYSZTOF PILORZA
BSTRACT . An algorithm is proposed for finding numerical solutions of a ki-netic equation that describes an infinite system of point articles placed in R d ( d ≥ ) . The particles perform random jumps with pair wise repulsion, in the courseof which they can also merge. The kinetic equation is an essentially nonlin-ear and nonlocal integro-differential equation, which can hardly be solved an-alytically. The derivation of the algorithm is based on the use of space-timediscretization, boundary conditions, composite Simpson and trapezoidal rules,Runge-Kutta methods, adjustable system-size schemes, etc. The algorithm isthen applied to the one-dimensional version of the equation with various initialconditions. It is shown that for special choices of the model parameters, the so-lutions may have unexpectable time behaviour. A numerical error analysis of theobtained results is also carried out.
1. I
NTRODUCTION
Stochastic dynamics of infinite populations placed in a continuous habitat mayinclude such aspects as their motion accompanied by merging. A pioneering workin this direction was published by Arratia in 1979, in which an infinite system ofcoalescing Brownian particles in R was proposed and studied [1]. It was thenextended to self-repelling motion of merging particles [2]. Over the next twodecades, different mathematical aspects of the Arratia model were intensively stud-ied, see [3–6] and the references therein. In the Arratia flow, an infinite number ofBrownian particles move independently up to their collision, then they coalesceand continue moving as single particles.Recently, an alternative model of this kind has been proposed [7,8]. Here, analo-gously to the Kawasaki approach [9, 10], particles make random jumps with repul-sion acting on the target point. Additionally, two particles merge into one particleplaced elsewhere with probability per time dependent on all the three locationsand independent of the remaining particles. Thereafter, this new particle partici-pates in the motion. The intensities and magnitude of jumps and coalescence aredetermined 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 Mathematics Subject Classification.
Key words and phrases.
Integro-differential equations, Population dynamics, Coalescence, Hop-ping, Spatial inhomogeneity, Numerical methods . study the phytoplankton dynamics [13–15]. Phytoplankton cells are dispersed inthe water, leading to a patchy distribution of entities. The latter is the basis for thevast majority of oceanic and freshwater food chains. Another example concerns themodeling of cancer mechanisms according to one of which melanoma cells migrateand 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 equa-tion have been reported. As a consequence, the role of repulsion interactions inthe possible appearance of a spatial heterogeneity in the system has been eluci-dated. There it was shown also that the kinetic equation can rigorously be obtainedfrom the corresponding microscopic evolution equations using previous theoreticalapproaches developed earlier [19–24] for other spatial population models. How-ever, the numerical consideration was restricted to a few particular examples whenchoosing 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 therepulsion-jump-coalescence kinetic equation. The derivation combines differenttechniques and develops methods allowing to reproduce properties of infinite sys-tems on the basis of finite samples. The algorithm is applied to population dynam-ics simulations of one-dimensional systems with various initial spatially inhomo-geneous densities and forms of the jump, coalescence, and repulsion kernels. Acomprehensive analysis of the obtained numerical results is also provided.2. K
INETIC 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 coa-lescence and repulsive jumps. Time evolution of the population will be describedin terms of the particle density n ( x , t ) . Performing the passage from the micro-scopic individual-based dynamics to the mesoscopic description by means of aVlasov-type scaling, one derives the following kinetic equation for the density inthe Poisson approximation with the jump-coalescence model [7, 8, 18]: ∂ n ( x , t ) ∂ t = − Z a ( x , y ) exp (cid:18) − Z ϕ ( y − u ) n ( u , t ) du (cid:19) n ( x , t ) dy + Z a ( x , y ) exp (cid:18) − Z ϕ ( x − u ) n ( u , t ) du (cid:19) n ( y , t ) dy − Z Z h b ( x , y , z ) + b ( y , x , z ) i n ( x , t ) n ( y , t ) dydz + Z Z b ( y , z , x ) n ( y , t ) n ( z , t ) dydz . (1)The first term in the rhs of Eq. (1) represents the process of jumping of agentsfrom 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 exponen-tial factor determines the strength of repulsion between the particle jumping to y LGORITHM FOR NUMERICAL SOLUTIONS TO THE KINETIC EQUATION OF COALESCING JUMPS3 and the rest of the system, where ϕ ( y − u ) denotes the corresponding interparticlepotential. The intensity contribution and repulsion strength depend on the configu-ration of all particles, so that the spatial integration over y and u is carried out. Thesecond term relates to the reverse process when particles located anywhere in thecoordinate space can appear at x due to their random jumps to the latter point. Thisincreases the density n ( x , t ) . The third term defines the process of free coalescencewhen two particles, located at x and y , merge into a single particle located at z withintensities b ( x , y , z ) and b ( y , x , z ) , resulting in a decrease of n ( x , t ) . Finally, mergedparticles can be randomly created in x owing to the coalescence of arbitrary twoother particles located in y and z with intensity b ( y , z , x ) , leading to an increase ofthe 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 coalescenceson the change of n ( x , t ) in x at time t . These integrations are invariant with re-spect to the mutual replacement y ↔ z in b ( x , y , z ) , b ( y , x , z ) , and b ( y , z , x ) . Thejump and coalescence kernels are nonnegative functions symmetric w.r.t. first twoarguments, a ( x , y ) = a ( y , x ) ≥ b ( x , y , z ) = b ( y , x , z ) ≥
0. The kernels sat-isfy the integrability conditions R a ( x , y ) dx = R a ( x , y ) dy = µ a , R ϕ ( x ) = µ ϕ , and R b ( x , x , x ) dx i dx j = µ b , where i , j = , , i = j . The quantities µ a , µ b , and µ ϕ can be treated as parameters of the model. Moreover, according to the transla-tion 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 choicefor the jump kernel is a ( x , y ) = a ( x − y ) = a ( | x − y | ) .The coalescence intensity will be modeled by the form b ( x , y , z ) = b ( x − y ) δ (( x + y ) / − z ) , (2)where b ( x − y ) = b ( | x − y | ) and the Dirac δ -function is applied. It implies thatthe resulting point of the coalescence of two particles in x and y is at the middle z = ( x + y ) / δ -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 lowerdimensionality of the integration. Indeed, in view of Eq. (2) and integrating over z to eliminate δ -function, the kinetic equation (1) transforms to ∂ n ( x , t ) ∂ t = − Z a ( x − y ) exp (cid:18) − Z ϕ ( y − u ) n ( u , t ) du (cid:19) n ( x , t ) dy + Z a ( x − y ) exp (cid:18) − Z ϕ ( x − u ) n ( u , t ) du (cid:19) n ( y , t ) dy − Z (cid:16) b ( x − y ) n ( x , t ) − d b ( ( x − y )) n ( x − y , t ) (cid:17) n ( y , t ) dy , (3)where the symmetrical properties of the kernels have been taken into account. Im-posing an initial condition n ( x , ) , Eq. (3) leads to a complicated partial integro-differential equation with respect to the unknown function n ( x , t ) at t >
0. Because
IGOR OMELYAN, YURI KOZITSKY, AND KRZYSZTOF PILORZ of the presence of spatial integrals and nonlinearity, we doubt it can be solved an-alytically in general. Moreover, we cannot apply the convolution theorem to avoidthe spatial integration in the last coalescence term of Eq. (3) due to the existence ofthree functions under the integrand which all depend on y . That is why we developa numerical approach for this equation. The convolution method in the absence ofcoalescence is presented in Appendix. Analytical solutions in the spatially homo-geneous case are also given there.3. N UMERICAL ALGORITHM
Spatial discretization. — In order to solve numerically the kinetic equation it isnecessary, first of all, to perform its discretization in coordinate space. Let x i be thegrid points uniformly distributed over R d with mesh h along all dimensions. Forsimplification of our further presentation we will restrict ourselves in this study toa particular case of d = d > dn i dt = h ∑ j (cid:16) a i − j exp h − h ∑ k ϕ i − k n k i n j − a i − j exp h − h ∑ k ϕ j − k n k i n i (cid:17) − h ∑ j b i − j n i n j + h ∑ j b i − j n j n i − j , (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 →
0, the discretizedkinetic 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 infinitenumber of terms, and introducing the auxiliary quantities λ i = exp h − h ∑ k ϕ i − k n k i , (5)one obtains that Eq. (4) can be cast more compactly in the following equivalentform dn i dt = h ∑ j (cid:16) a j ( λ i n i − j − λ i − j n i ) − b j n i n i − j + b j n i − j n i + j (cid:17) (6)where a j = a ( jh ) , b j = b ( jh ) , and b j = b ( jh ) .In computer simulations we cannot operate with infinite-size samples leading tothe infinite summation over j in Eqs. (4), (5), and (6). Because of this we considera finite number N of grid points x i uniformly distributed over the area [ − L / , L / ] with spacing h = L / N , where i = , , . . . , N . This area will represent an interval,a square, or a cube in cases d =
1, 2, or 3, respectively. The finite length L shouldbe 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 discretizationnoise. Then h will be sufficiently small to provide a high accuracy of the spatialintegration. The finite-size effects can be reduced by applying the correspondingboundary conditions (BC) when mapping infinite range x ∈ ] − ∞ , ∞ [ by finite area LGORITHM FOR NUMERICAL SOLUTIONS TO THE KINETIC EQUATION OF COALESCING JUMPS5 [ − L / , L / ] . In view of the aforesaid, Eq. (6) presents a coupled system of N autonomous ordinary differential equations, where i = , , . . . , N and summationover 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 choicedepends on initial function n ( x , ) and properties of solution n ( x , t ) . For example,if n ( x , ) takes nonzero values only within a narrow interval [ − l / , l / ] with l ≪ L , we can apply the DBC by letting n j = | x j | > L /
2. This means thatduring the finite simulation time 0 ≤ t ≤ T , the nonzero values of n ( x , t ) do notapproach the boundaries x B = ± ( L / − max σ ) , where max σ is the maximal radiusof the kernels (see Sect. 4). In numerical calculations this can be expressed by thecondition n ( x B , t ) < ε max x n ( x , t ) , where 0 < ε ≪ n ( x B , t ) > ε max x n ( x , t ) , we shouldenlarge 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 isdiscussed in Sect. 4, where members of infinite configuration are initially absentin one half-space, requires a modified BC that combines DBC and ABC with anaddition of adjustable system-size approach.If n ( x , ) and, thus, n ( x , t ) are periodic functions, it is necessary to apply PBCwith fixed finite size L . According to PBC, the summation in Eq. (6) for eachcurrent i = , , . . . , N is performed not only over all j = , , . . . , N but also overall infinite number of images j ′ of j . The images are obtained by repeating thebasic interval [ − L / , L / ] by the infinite number of times to the left and to theright of it, so that x j ′ = x j ± KL , where K = , , . . . , ∞ and n j ′ = n j . This repro-duces the periodicity n ( x ± KL , ) = n ( x , ) , where x ∈ [ − L / , L / ] . The solution n ( x ± KL , t ) = n ( x , t ) will also be periodic for any time t > 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 interpar-ticle distance | x j | , the summation over j in Eq. (6) can be actually restricted toa 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 | ) ≈ b ( | x | ) ≈ | 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 / 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 spatiallyinhomogeneous solutions n ( x , t ) which are flat for x < − L / x > L / t where they take nonzero constant values. Then n j ′ = n ( − L / , t ) for all x j ′ < − L / n j ′ = n ( L / , t ) for all x j ′ > L /
2. If in the course of time theflatness is violated at a current L , the basic length should be enlarged using theautomatically adjustable system-size approach mentioned above.Taking into account the properties a − j = a j , b − j = b j and the fact that the in-fluence of a ( x ) and b ( x ) can be neglected at large distances x > R a , b , i.e. assuming IGOR OMELYAN, YURI KOZITSKY, AND KRZYSZTOF PILORZ that a ( x ) = | x | > R a and b ( x ) = | x | > R b , Eq. (6) transforms to dn i dt = h j a ∑ j = ξ ( a ) j a j (cid:16) λ i ( n i − j + n i + j ) − n i ( λ i − j + λ i + j ) (cid:17) − h ξ ( b ) b n i − hn i j b ∑ j = ξ ( b ) j b j ( n i − j + n i + j )+ h ξ ( b ) b n i + h j b / ∑ j = ξ ( b ) j b j n i − j n i + j (7)where i = , , . . . , N and the summations over j are performed already with finitepositive nonzero integers up to j a = R a / h or j b = R b / h . Quite similarly, using theproperties ϕ − j = ϕ j and ϕ ( x ) = | x | > R ϕ , one finds from Eq. (5) that λ i = exp h − h BC ∑ j ξ ( ϕ ) j ϕ j n i − j i = exp h − h ξ ( ϕ ) ϕ n i − h j ϕ ∑ j = ξ ( ϕ ) j ϕ j ( n i − j + n i + j ) i , (8)where i = , , . . . , N and j ϕ = R ϕ / h with truncation radius R ϕ < L /
2. It is evidentthat in the limits L , N , R → ∞ provided h →
0, the discretized equation (7) with Eq.(8) coincide with its original, continuous counterpart (3).Weights ξ ( a , b , b , ϕ ) j at kernel values a j , b j , b 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 satisfythe normalization condition ξ + ∑ j ∗ j = ξ j = j ∗ , where j ∗ = j a , b , ϕ or j ∗ = j b / ξ = ξ = ξ = . . . = ξ j ∗ − =
1, while ξ j ∗ = /
2. The composite Simpson methodof the fourth order yields ξ j ∗ = / , ξ j ∗ − = / , ξ j ∗ − = / , ξ j ∗ − = / , ξ j ∗ − = / , . . . , ξ = / j ∗ is odd and ξ = / j ∗ is even. The latter is more accu-rate than the former and involve numerical uncertainties of order of O ( h ) versus O ( h ) .We mention that n i are explicitly defined at i = , , . . . , N , so that the corre-sponding BC should be applied in Eqs. (7) and (8) to n i − j and/or n i + j whenever i − j < i + j > N , because j = , , . . . , j ∗ < N /
2. The uniform knot distri-bution over [ − L / , L / ] can be chosen in the form x i = − L / + ( i − / ) h , where i = , , . . . , N with even N . This provides the symmetricity of knot positions withrespect to x =
0. Then according to PBC the calculations of n i ± j should be per-formed as n i − j = (cid:26) n i − j + N , i − j < n i − j , ≤ i − j ≤ N , n i + j = (cid:26) n i + j − N i + j > Nn i + j ≤ i − j ≤ N . (9)Note that i = , , . . . , N and j = , , . . . , j ∗ < N /
2. The applications of DBC andABC result in n i − j = (cid:26) , i − j < n i − j , ≤ i − j ≤ N , n i + j = (cid:26) i + j > Nn i + j ≤ i − j ≤ N . (10) LGORITHM FOR NUMERICAL SOLUTIONS TO THE KINETIC EQUATION OF COALESCING JUMPS7 and n i − j = (cid:26) n , i − j < n i − j , ≤ i − j ≤ N , n i + j = (cid:26) n N i + j > Nn i + j ≤ i − j ≤ N . (11)respectively. In the cases d = d =
3, the boundary transformations can beimplemented similarly to those given by Eqs. (9) – (11).
Time integration. — In the most general form, our coupled system of N au-tonomous ordinary differential equations can be given as dn i dt = ˙ n i = s i ( n , n , . . . n N ) , (12)where i = , , . . . , N and s i ( n , n , . . . , n N ) represents the rhs of Eq. (7) with takinginto account Eq. (8). Introducing the set Γ ( t ) = { n i ( t ) } of N dynamical variablesand their time derivatives Φ ( t ) = { s i ( t ) } , Eq. (12) can be compactly cast as ˙ Γ = Φ ( Γ ) . A common practice to solve differential equations of such a type is to useclassical Runge-Kutta schemes [25]. The scheme of the fourth order (RK4) reads Γ ( t + ∆ t ) = Γ ( t ) + ∆ t (cid:16) Φ + Φ + Φ + Φ ) + O ( ∆ t ) , (13)where ∆ t is the time step and Φ = Φ ( Γ ( t )) , Φ = Φ ( Γ ( t ) + Φ ∆ t / ) , Φ = Φ ( Γ ( t ) + Φ ∆ t / ) , Φ = Φ ( Γ ( t ) + Φ ∆ t ) are the derivatives in intermediate stages.The Runge-Kutta scheme of the second order (RK2) can also be applied. There aretwo versions of RK2. The first one (will be referred to as simply RK2) is known asthe Heun method and can be related to the trapezoidal or Verlet-like integration. Ithas the form Γ ( t + ∆ t ) = Γ ( t ) + [ Φ ( Γ ( t )) + Φ ( Γ ( t ) + Φ ( Γ ( t )) ∆ t )] ∆ t / + O ( ∆ t ) .The second version (RK2’) relates to a middle point scheme, Γ ( t + ∆ t ) = Γ ( t ) + Φ ( Γ ( t ) + Φ ( Γ ( t )) ∆ t / ) ∆ t + O ( ∆ t ) . The RK2 and RK20 approaches are two-stageschemes. The RK4 integration consists of four stages, so that at the same ∆ t it willrequire twice larger number of operations to cover the same time interval T of theobservation with respect to RK2. However, RK4 is more accurate than RK2 with O ( ∆ t ) versus O ( ∆ t ) local truncation errors.In such a way, the numerical solution n i ( t ) can be found for any t = p ∆ t , where p = , , . . . , P over the interval t ∈ [ , T ] with T = P ∆ t by sequentially repeating P times the transformations given by Eq. (13). The O ( ∆ t ) -uncertainties can bereduced to an arbitrary small level by decreasing the length ∆ t of the time step. Wemention that the finite-size effects are minimized by choosing a sufficiently largesize L of the system and applying the corresponding boundary conditions. Theuncertainties 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 totalnumber 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 > d , namely as N d . For convolution solutions this number is lowered to order of N d ln N d (seeAppendix). In the case of simple rectangle kernels the total number of operationscan be decreased to be proportional to N d . IGOR OMELYAN, YURI KOZITSKY, AND KRZYSZTOF PILORZ
4. A
PPLICATIONS AND ANALYSIS OF THE RESULTS
Numerical details. — The numerical simulations were carried out for one-dimensional systems ( d =
1) of free or repulsive jumping particles which can coa-lesce. The kinetic equation (3) was solved by using the algorithm described in thepreceding section. The discretization in R was done with a mesh of h = . . L =
20 within either the fixed (for PBC) or adjustable-size(for DBC and ABC) regime. Spatial integration was performed by employing thecomposite Simpson rule. Time integration was done with the help of the RK4scheme at a step of ∆ t = .
1. The relative tolerable level was chosen to be equal to ε ∼ − . Further increasing space and time resolutions did not noticeably affectthe solutions (see the last subsection).The jump a ( x ) , coalescence b ( x ) , and repulsion ϕ ( x ) kernels were modelled byGaussian G µ , σ ( x ) = µ ( πσ ) exp (cid:18) − x σ (cid:19) (14)or rectangle C µ , σ ( x ) = (cid:26) µ σ , | x | ≤ σ , | x | > σ (15)functions, where µ = µ a , µ b , or µ ϕ and σ = σ a , σ b , or σ ϕ are the intensities andranges of the corresponding interactions, respectively. The kernels are normalized, R { a , b , ϕ } ( x ) dx = µ a , µ b , µ ϕ . A symmetrical pair of shifted Gaussian or rectanglefunctions was involved as well, F µ , σ , s ( x ) = (cid:16) F µ , σ ( x − s ) + F µ , σ ( x + s ) (cid:17) , (16)where F ≡ G or C and s is the shifting interval. The truncation radius for Gaussiankernels was R = Q σ with Q = G µ , σ ( Q σ ) / G µ , σ ( ) ∼ − , so that theircontribution at | x | > R can be neglected. In the case of discontinuous rectanglekernels we have Q = R = Q σ + s .In view of the above, we have six kernel parameters, µ a , µ b , and µ ϕ , as wellas σ a , σ b , and σ ϕ for the repulsion-jump-coalescence model. This leads to a largenumber of all possible combinations. Because of this we will deal with most char-acteristic examples when choosing values for these parameters for single Gaussianand rectangle kernels or their ± s -shifted double counterparts. In order to consis-tently analyze their influence on the dynamics of populations the following foursituations will be considered: (i) pure free jumps; (ii) repulsive jumps; (iii) purecoalescence; and (iv) repulsive jumps with coalescence. It is worth emphasiz-ing also that solution n ( x , t ) depends cardinally on the choice of initial condition n ( x ) ≡ n ( x , ) . Like kernels, the initial condition can be selected in the form ofGaussians (14) or rectangles (15), trigonometric or step functions, etc. Rectangle initial density profiles. — The first example relates to the initialcondition in the form of periodic rectangle function n ( x , ) = C , , L ( x ) . The in-finite system is reproduced by repeating the single rectangle segment C , ( x ) at LGORITHM FOR NUMERICAL SOLUTIONS TO THE KINETIC EQUATION OF COALESCING JUMPS9 x ∈ Ω = [ − L / , L / ] on the interval ] − ∞ , ∞ [ with a period of L and applyingPBC. The jump a ( x ) = G µ a , σ a ( x ) , repulsion ϕ ( x ) = G µ ϕ , σ ϕ ( x ) , and coalescence b ( x ) = G µ b , σ b ( x ) kernels are modelled by the Gaussians with intensities µ a = µ ϕ =
20, and µ b =
1. Three sets of kernel ranges were considered, namely, σ a = σ ϕ = σ b = σ a < σ ϕ < σ b with σ a = . σ ϕ = σ b = σ a > σ ϕ > σ b with σ a = σ ϕ = σ b = .
5. The corresponding time evolu-tion of spatial structure n ( x , t ) is presented in Fig. 1 for the cases of pure freejumps, repulsive jumps, pure coalescence, and repulsive jumps with coalescencewith σ 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).From Fig. 1(a) we see that for pure free jumps, the discontinuity of C , ( x ) at x = ± t &
10. For longer times, t &
40, the solution n ( x , t ) tendsto more and more homogeneous densities. In the limit t → ∞ we expect absolutelyflat function lim t → ∞ n ( x , t ) = n independent on x , where n = L R Ω n ( x , ) dx = / L = .
05. Moreover, for any time t the total number R Ω n ( x , t ) dt of particles on the in-terval Ω = [ − L / , L / ] is constant because of the absence of coalescence. Forstrong repulsive jumps [part (b)] the density function n ( x , t ) remains discontinuousat x = ± t ∼
40 and its shape at shorter times is more complicated. In par-ticular, apart from the existence of the main maximum at x = t , two symmetric secondary maxima appear additionallyat 0 < t .
40 in the ranges x ≈ ± t → ∞ with the same asymptotic behaviourlim t → ∞ n ( x , t ) = n = / L = .
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 parti-cles remain to be located exclusively within the initial interval [ − , ] and they areabsent outside of it at any t . In other words, no density propagation is indicatedbecause the particles do not move apart from coalescence.When the repulsive jumps are carried out in the presence of coalescence at equalinteraction ranges σ a = σ ϕ = σ b =
1, see part (d) of Fig. 1, the pattern is somewhatsimilar to that of part (b). However, the three-maximum structure dissipates nowmuch faster, leading to the zeroth asymptotics already at t &
40. For short-rangedjumps, where σ a = . σ ϕ = , σ b =
2, the central peaks at x = x ≈ ± σ a = , σ ϕ = , σ b = .
5, the central peaks transform into a more complicated structure withone central minimum at x = x ≈ ± .
5, see part (f). Thesecondary maxima at x ≈ ± n ( − x , t ) = n ( x , t ) , like the initial condition, n ( − x , ) = n ( x , ) .This follows from the symmetry of the kinetic equation.The second choice deals with asymmetric initial condition n ( − x , ) = n ( x , ) inthe form of N shifted single rectangle functions C v k , σ k ( x + s k ) with intensities v k F IGURE
1. Time evolution of density n ( x , t ) as dependent onspatial coordinate x at several moments of time t for rectangle ini-tial condition n ( x , ) = C , , L ( x ) in the cases: (a) pure free jumps;(b) repulsive jumps; and (c) pure coalescence; as well as repul-sive jumps and coalescence for (d) equal interaction ranges; (e)short-ranged jumps; and (f) short-ranged coalescence. The infi-nite system is periodic ( L =
20) on the interval x ∈ [ − , ] withjump, repulsion, and coalescence Gaussian kernels of different in-tensities and ranges (see the legends inside). LGORITHM FOR NUMERICAL SOLUTIONS TO THE KINETIC EQUATION OF COALESCING JUMPS11 and ranges σ k , namely, n ( x , ) = C v k , σ k , s , N ( x ) = N ∑ k = C v k , σ k ( x + s k ) (17)where s k = − L / + ( k − / ) L / N are shifting parameters. Repeating (17) withperiod L , we should apply PBC to deal with the infinite system. We used a par-ticular case of (17) with N = L =
20 as well as three different amplitudes v , , generated at random in the interval ] , [ . The corresponding result on this isshown 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 so-lutions obtained for single rectangle initial densities C v k , σ k ( x + s k ) . With increasingtime, a coherence between the separate solutions appears. Again, in the absenceof 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 = L R Ω n ( x , ) dx = ( v + v + v ) / L forany initial distributions. We mention that at µ b =
0, the total amount R Ω n ( x , t ) dx ofparticles in Ω is unchanged and equal to its initial value R Ω n ( x , ) dx . For µ b > n ( x , t ) seems to take its zeroth asymptotics lim t → ∞ n ( x , t ) = 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 , ) = n ( x , ) with x ∈ [ − L / , L / ] , where K = , , . . . , ∞ , the solution n ( x ± KL , t ) = n ( x , t ) will also be periodic for any time t > L . Then, in particular, n ( − L / , t ) = n ( L / , t ) as is confirmed in Fig. 2.Investigations show that the increase of the strength µ and range σ of the jump andcoalescence kernels accelerates this process, see for instance, parts (e) and (f) ofFigs. 1 and 2. Using the Gaussian (instead of rectangle) initial conditions and rec-tangle (instead of Gaussian) kernels leads to results (not shown) which are similarto those of Figs. 1 and 2.As visualized in Figs. 1 and 2, the presence of coalescence ( b ( x ) =
0) seems tolead to the zeroth asymptotic R Ω n ( x , t ) dx → b ( x ) = C µ , σ , s ( x ) of a pair of two shifted rectangle functions [Eq. (16)] with ap-propriate shifting parameter s to avoid the zeroth density limit. The initial inho-mogeneous density n ( x , ) over the basic interval [ − L / , L / ] with PBC at period L should also be chosen correspondingly. For instance, we can consider n ( x , ) inthe form C v k , σ k , s ′ , N ( x ) of two ( N =
2) single rectangle functions [Eq. (17)] withthe same amplitude v k = µ b = σ k = σ b ≡ σ = b ( x ) but with different shifting parameter s ′ = s satisfying the con-straint s ′ = s − σ . Choosing s = σ = s ′ =
8. The correspond-ing result for the case n ( x , ) = C , , , ( x ) with L =
20 and b ( x ) = C , , ( x ) in theabsence ( a =
0) of jumps is depicted in Fig. 3(a). We see that at the beginning, therectangles soon transform into triangle-shaped peaks centered at x = ± ( + KL ) ,while the additional triangles appear exactly in the middle of them at x = ± KL (below we will omit the terms ± KL to simplify notation). With increasing time, F IGURE
2. Time evolution of solution n ( x , t ) for asymmetric ini-tial density in the form of three rectangle functions with differentamplitudes. Other notations are the same as in Fig. 1.the widths of the triangle-shaped peaks decrease but their maxima at x = ± x =
0. At sufficiently long times t & · themodification of the density profile slows down to a level suggesting that the sys-tem approaches a non-trivial stationary state in which ∂ n ( x , t ) / ∂ t =
0. In otherwords, the shape of the density profile is transformed to such a form at which anycoalescence processes become impossible in view of the specific form of the co-alescence kernel b ( x ) = C , , ( x ) . The latter accepts nonzero values only in theinterval s ′ − σ = < | x | < = s ′ + σ , so that the absence of coalescence at a given LGORITHM FOR NUMERICAL SOLUTIONS TO THE KINETIC EQUATION OF COALESCING JUMPS13 F IGURE
3. Dynamics of spatial structure n ( x , t ) for initial den-sity in the form of two rectangle functions, n ( x , ) = C , , , ( x ) in the cases: (a) pure coalescence and (b) coalescence with freejumps. The interactions are modelled by shifted pair rectangle co-alescence C , , ( x ) and Gaussian jump G . , ( x ) kernels.spatial configuration means that there exists no pair of particles with interparti-cle separations | x | lying in the interval [ , ] . Allowing particles to jump changesthe situation radically, as is demonstrated in Fig. 3(b) for Gaussian jump kernel G . , ( x ) . Even for relatively small jump intensity µ a = . σ a =
1, thedensity quickly decreases to zero for each x , after initial period of time, when newpeaks are formed. It is interesting to remark that the monotonic decrease of mainmaxima in x ± x = ± KL ). This magnitude first increases at 0 < t . t ∼
4, and then decreases at t & Trigonometric initial density profiles. — Another interesting case is the initialdensity in the form of a trigonometric function, n ( x , ) = T n , µ , k ( x ) = n (cid:16) + µ cos ( π kx / L ) (cid:17) , (18)where 0 < µ < k ≥ L / k in coordinate. Then PBC should be used to reproduce the infinite system. Theobtained solutions n ( x , t ) for n ( x , ) = T , , ( x ) with n = µ = k =
3, and L =
20 are plotted in Fig. 4 when jump, repulsion, and coalescence kernels areGaussians a ( x ) = G , ( x ) , ϕ ( x ) = G , ( x ) , and b ( x ) = G . , ( x ) , respectively. Forconvenience, we presented n ( x , t ) only on the right-hand side of coordinate spaceat x ≥ k and periodicity L . In particular, the density continues to oscillate around the samelevel n ( + µ ) / = t , so that the particle density profile seems to becomeflat at long times, lim t → ∞ n ( x , t ) =
1. Again, the total number R Ω n ( x , t ) dx ofparticles in Ω ∈ [ − L / , L / ] does not change in time. The repulsion makes the F IGURE
4. Density profile n ( x , t ) starting from trigonometricinitial condition n ( x , ) = T , , ( x ) . The jump, repulsion, and co-alescence kernels are Gaussians G , ( x ) , G , ( x ) , and G . , ( x ) ,respectively. Other notations are the same as for Fig. 1.simple trigonometric (cosine) form much more complicated with the existence ofadditional maxima and minima inside Ω , see Fig. 4(b). Moreover, the homogeneityis being achieved here much slower than in the case of free jumps (compare densityat t & t &
10 in Fig 4(a)). This can be explained by the strongintensity ( µ ϕ =
8) of repulsion potential ϕ ( x ) = G , ( x ) .For pure coalescence in Fig. 4(c), the distribution n ( x , t ) is no longer of atrigonometric form at t > t → ∞ . Theinclusion of repulsive jumps changes the behaviour of n ( x , t ) near its minima in acharacteristic 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 timepractically does not change at t .
10 despite their strength. Moreover, the shapeof the density profile at long t &
10 is modified as well, making it more flat withsmaller oscillations. As a result, spatial homogeneity is obtained here faster.
LGORITHM FOR NUMERICAL SOLUTIONS TO THE KINETIC EQUATION OF COALESCING JUMPS15
Step initial density profiles. — A special case presents initial condition in theform of the Heaviside step function n ( x , ) = H n ( x ) = (cid:26) n , x ≤ , x > . (19)Here the system is initially ( t =
0) considered on the finite interval [ − L / , L / ] with no PC. Further ( t >
0) its size is gradually increased to the infinity on theunbounded interval x ∈ ] − ∞ , ∞ [ at t → ∞ according to the automatically adjustedsystem-size (AASS) approach. Then, a modified BC should be applied by combi-nation of DBC and ABC together with AASS when analyzing the densities on theboundaries ± L /
2. The DBC are used from the right, where lim x → ∞ n ( x , t ) = t . From the left, where lim x →− ∞ n ( x , t ) = n ( t ) = n ( ) = n , we must em-ploy ABC. This means that we should exploit the DBC rule for n i + j given by thesecond equality of Eq. (10), i.e., replace n i + j by 0 if i + j > N , and the ABC rulefor n i − j given by the first equality of Eq. (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 / , t ) > ε max x n ( x , t ) , the system size L is enlarged (in two times) and the sim-ulations are continued.By monitoring from the left the difference between the actual values of n ( x , t ) at x = − L / n hRK ( t ) [obtained by solving numeri-cally the kinetic equation for the spatially homogeneous initial condition n ( x , ) = 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 exceedsthe predefined level ε max x n ( x , t ) we should enlarge L according to the automati-cally adjustable system-size approach. The asymptotic value n ( t ) will differ fromthe exact solution n h ( t ) even at very large L because of the approximate char-acter of time integration. In the limit of sufficiently small step sizes ∆ t whenlim ∆ t → , L → ∞ n ( 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 valuesof n ( x , t ) at x = − L / n h ( t ) at x → − ∞ we can estimatethe influence of two effects in one fashion, namely, those caused by the finite-ness of L (should be significantly large) and ∆ t (should be significantly small). Insuch a way, both ABC and DBC deviations are analyzed when deciding on thesize enlargement within the ABC-DBC scheme. This completes the automaticallyadjustable system-size approach.Time evolution of n ( x , t ) for n ( x , ) = H ( x ) is presented in Fig. 5 using Gauss-ian jump G µ a , σ a , repulsion G µ ϕ , σ ϕ , and coalescence G µ a , σ b kernels with differentintensities of µ a = µ ϕ =
8, and µ b = .
1, respectively, as well as with equal[ σ a , ϕ , b =
1, parts (a)–(d)] or different [ σ a = . σ ϕ = σ b =
2, part (e), and σ a = σ ϕ = σ b = .
5, part (f)] ranges. As can be seen from Fig. 5(a) forpure free jumps, the discontinuous step function n ( x , ) = H ( x ) transforms into acontinuous S-shaped curve at finite times t >
0. All the curves intersect each otherin the same point ( , / ) , where 1 / the decrease of the amount of particles for x < x >
0, that can be written as R − ∞ [ − n ( x , t )] dx = R ∞ n ( x , t ) dx . The samestatement concerns repulsive jumps, but here the curves are intersected in a pointwhich lies below ( , / ) . 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 morecomplicated, including the existence of maximum in x ∼
1. The latter disappearsat 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 , ) = H ( x ) changes to a step-like dependence n ( x , t ) at t > x ∼ −
1. The density in ] − ∞ , [ decreases from 1 at t = t >
0. No particles appear at all in ] , ∞ [ for any t . Moreover, the initialdiscontinuity does not vanish even for relatively long times. The inclusion of re-pulsive jumps, see Fig. 5(d), prevents the appearance of peaks at x ∼ −
1, whileat x > 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 & ) as in the case of purecoalescence. As the range of jumps becomes shorter, σ a = .
5, and the range ofcoalescence increases, σ b =
2, see Fig. 5(e), the peaks in n ( x , t ) appear again at x ∼ − σ b = − σ b = x ∼ x ∼ . σ a =
2, and the coalescencerange decreases, σ b = .
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 ap-proaches zero more slowly here. Note that for the inverse initial unit step function n ( x , ) = H ( − x ) , the results n ( x , t ) presented in Fig. 5 for n ( x , ) = H ( x ) shouldalso be inversed by n ( − x , t ) to obtain solutions without resolving the kinetic equa-tion. This statement is quite general and remains in force not only for step func-tions, but for any other asymmetric initial conditions. Namely, when picking theinitial condition n ∗ ( x , ) = n ( − x , ) , we automatically obtain n ∗ ( x , t ) = n ( − x , t ) ,where n ( x , t ) is the solution for n ( x , ) . This follows from the structure of kineticequation (3) and the symmetricity of kernel functions.Consider, finally, the dynamics of the initial step distribution n ( x , ) = H ( x ) inthe presence of pairs of shifted Gaussian kernels. First, we used shifting parameter s = µ a =
1) jump kernel a ( x ) = G , , ( x ) and weak ( µ b = .
05) coa-lescence interaction b ( x ) = G . , , ( x ) , while s ′ = µ ϕ =
10) repulsionpotential ϕ ( x ) = G , , ( x ) (the ranges were σ a = σ b = σ ϕ = t > s = n ( x ) = s ′ =
8. For n ( x , ) = H ( x ) the inhomogeneous front propagates to theleft with increasing amplitude at not too large x and to the right with dumping oscil-lations (for n ( x , ) = H ( − x ) the pattern is inverse). The inclusion of coalescenceeven with a slight intensity of µ b = .
05 drastically changes the situation, see Fig.
LGORITHM FOR NUMERICAL SOLUTIONS TO THE KINETIC EQUATION OF COALESCING JUMPS17 F IGURE
5. Time evolution of density profile for initial con-dition H ( x ) . The jump, repulsion, and coalescence kernels areGaussians with different intensities and ranges (see the legendsinside). Initially ( t =
0) the system is considered on the finite in-terval [ − , ] with no periodic conditions and further ( t >
0) itssize gradually increases to the infinity on x ∈ ] − ∞ , ∞ [ at t → ∞ according to the automatically adjusted approach. Other notationsare the same as for Fig. 1.6(b). Here, the oscillations are not so strong, especially at x >
0, and they quicklydisappear with increasing time.Further, we decreased and increased the shifting parameters in two times to val-ues s = s ′ = s = s ′ =
8. The results for these two cases F IGURE
6. Dynamics of the system starting from unit step func-tion H ( x ) in the cases of pure repulsive jumps with different ker-nel shifts: (a) s = s ′ =
4; (c) s = s ′ =
2; and (d) s = s ′ = s = s ′ =
4. Initially ( t =
0) the system is considered on the finite in-terval [ − , ] with no periodic conditions and further ( t >
0) itssize gradually increases to the infinity on x ∈ ] − ∞ , ∞ [ at t → ∞ ac-cording to the automatically adjusted approach. The interactionsare described by the shifted pair Gaussian jump G , , s , repulsion G , , s ′ and coalescence G . , , 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 ofdensity propagation by lowering its speed and oscillation amplitude. Moreover, thedistance between peaks in n ( x , t ) and their width also decrease correspondingly to2 s ′ = s =
1. This is in a contrast to the opposite case when the kernel shiftsincrease. Such an increase leads to stimulation of density propagation with highamplitude of the oscillations. Then the distance between peaks and their widthincrease to 2 s ′ =
16 and s =
4, respectively. The propagation speed grows propor-tionally as well.
Precision of the simulations.
Accuracy of our numerical simulations was mea-sured in terms of relative absolute deviations of density values n i = n ( x i , t ) in grid LGORITHM FOR NUMERICAL SOLUTIONS TO THE KINETIC EQUATION OF COALESCING JUMPS19 points x i from “exact” data ˘ n i at time t using the relation Θ ( h , ∆ t ) = ∑ ˘ Ni = | n i − ˘ n i | ∑ ˘ Ni = ˘ n i × [ % ] . (20)The total number ˘ N ≤ N of grid points involved into summation (20) depends onthe spatial region considered. The “exact” (or rather reference) values ˘ n i were ob-tained at high enough space and time resolutions with a tiny mesh of h = . ∆ t = .
01 in order to be entitled to ignore the numeri-cal uncertainties. The spatial and time integrations were performed with the helpof the composite Simpson rule and RK4 algorithm, respectively. The numeri-cal error analysis was done by carrying out a series of simulations at different h = . , . , . , . , . ∆ t = . , . , . , . , . , .
5. The numerical de-viations were then estimated by Eq. (20) to plot Θ ( h , ∆ t ) in a wide range of varying h and ∆ t . For the purpose of comparison, the simulation runs with the compositetrapezoidal rule and RK2 algorithm were performed as well. Error results presentedbelow are related to one of the situations considered in Figs. 1–6 when choosinginitial condition n ( x , ) , namely, to the case of Fig. 5(d) at t =
50. There, the ho-mogeneous and inhomogeneous intervals were x ∈ [ − , − ] and x ∈ [ − , ] ,respectively. Similar results were observed for all the rest choices of n ( x , ) .The numerical uncertainties Θ cased by spatial discretization and spatial inte-gration are shown in Fig. 7(a) as functions of the length h of space mesh (at a fixedRK4 time step of ∆ t = . Θ ( h , ∆ t ) at small h and ∆ t in more detail. Here we distin-guish three kinds of errors. The first is related to uncertainties due to single spatialintegrations. The corresponding results for them obtained by the composite Simp-son 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 Θ = c + h or log Θ = c + h , where c and c denote some con-stants. This confirms the fact [25] that the composite Simpson and trapezoidal rulesare accurate to the fourth O ( h ) and second O ( h ) orders in h , respectively. Theaccuracy is lowered significantly when considering the overall uncertainties causedby spatial discretization of the kinetic equation. Such uncertainties observed in thehomogeneous and inhomogeneous interval with the Simpson integration are pre-sented by curves labeled correspondingly as SMPSh and SMPSi. The TRPZh andTRPZi data of the trapezoidal integration appear to be approximately the same (seebelow the reason) and are not shown in the figure.Maximal uncertainties are achieved in the inhomogeneous region because thedependence 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 adiscretization are of order of O ( h ) . Taking into account that the total num-ber of numerical operations is proportional to N , the global errors appear asan accumulation of local ones leading to the overall uncertainties Θ of order of F IGURE
7. (a) Uncertainties Θ of the simulations related to sin-gle composite Simpson (SMPS) or trapezoidal (TRPZ) spatial in-tegrations and the overall spatial discretization obtained in homo-geneous (SMPSh) and inhomogeneous (SMPSi) regions at differ-ent mesh h ; (b) Dependencies of Θ on the size of the time step ∆ t for the RK2, RK2’, and RK4 time integrations. N O ( h ) ∼ O ( h ) , since h = L / N . In other words, the SMPSh and SMPSi de-viations behave as log Θ = c + log h . Despite the linear dependence O ( h ) , eventhe SMPSi errors are relatively small and do not exceed about 0 .
01 – 0 .
03% at h = .
01 – 0 .
02 (see Fig. 7). They are, however, much larger than those of the sin-gle spatial integration. Therefore, the errors caused by spatial discretization of thekinetic equation significantly prevail over the spatial integration uncertainties (theboth 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 theSMPSh/SMPSi and TRPZh/TRPZi data are very close to each other.Dependencies of Θ on the size of the time step ∆ t , obtained in the simulationswith the RK2, RK2’, and RK4 time integrations (at a fixed spatial mesh of h = . x ∈ [ − , ]) 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 or ∆ t for the algorithms of the second (RK2 and RK2’)or fourth (RK4) orders, respectively. Here, the log-log plots are lines with slopes2 or 4, i.e., log Θ ( ∆ t ) = c + ∆ t or log Θ ( ∆ t ) = c + ∆ t [similarly as fordependence 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 algorithmare of order of O ( ∆ t ) , see Eq. (13). So that the numerical integration over longtimes t ≫ ∆ t leads to the global errors of order of t / ∆ t O ( ∆ t ) = O ( ∆ t ) as thisrequires repeating of the single-step integration by t / ∆ t times. Analogously, theRK2 and RK2’ local uncertainties O ( ∆ t ) modifies to the global O ( ∆ t ) -errors.Since ∆ t decreases with decreasing ∆ t much faster than ∆ t , the fourth-order RK4algorithm produces much more precise solutions than the second-order integratorsRK2 and RK2’. Thus, the former can be recommended when a very high accuracyis required. For instance, the RK4 algorithm is precise up to a negligible level of LGORITHM FOR NUMERICAL SOLUTIONS TO THE KINETIC EQUATION OF COALESCING JUMPS21 order of Θ ∼ − % at ∆ t ∼ .
1. On the other hand, the RK2 and RK2’ integra-tors can provide here only an accuracy of 10 − [RK2 is slightly better than RK2’,see Fig. 7(b)]. Remember, however, that the overall errors include also the spatialdiscretization uncertainties which are of order of 10 − % at h = .
01, see Fig. 7(a).For this reason, the RK4 algorithm can be applied with much larger time steps of ∆ t ∼ ∼ − andare comparable with the space-discretization errors. This significantly acceleratesthe simulations since the computational time is inverse proportional to ∆ t . In par-ticular, the RK4 speedup at t = Θ ∼ − %) with respect to the RK2 andRK2’ schemes at ∆ t = . Θ ∼ − % is observeddespite the smaller time step) is of order of 1 / . / = / = ∆ t , all the integratorsbecome unstable and cannot be used at ∆ t >
2, where Θ > ONCLUSION
In this paper we have derived an efficient algorithm to obtain numerical solutionsfor the time-differential kinetic equation which approximates nonlocal stochasticevolution 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 thepresence of complicated spatial integrals with nonlinear and nonlocal terms andtherefore requires numerical analysis. The proposed algorithm is based on a setof techniques including time-space discretization, periodic, Dirichlet, and asymp-totic boundary conditions, composite Simpson and trapezoidal rules, Runge-Kuttaschemes, and automatically adjustable system-size approaches. This has allowedas to carry out simulations of dynamics for one-dimensional systems with variousinitial inhomogeneous densities and different forms of the coalescence, jump andrepulsion kernels, giving a comprehensive study of the jump-coalescence dynam-ics. A numerical error analysis of the obtained results has also been performed.For some specific choices of the model parameters and initial densities, a non-trivial dynamics has been revealed. In the case of pure free jumps the systemalways tends to a nonzero homogeneous density for any initial conditions. In con-trast, the presence of strong repulsion potentials can result in the appearance ofpersistent wave-like density propagations when the repulsion and jump kernels arechosen in the form of a sum of shifted single Gaussian functions. The shifting pa-rameters define the shape and spatial period of density structures. The introductionof even relatively weak coalescence prevents them to persist. In the case of purecoalescence, the population eventually goes extinct, except for a special choice ofthe initial density profile and coalescence kernel. For instance, if the particles areinitially located exclusively on an archipelago of islands and the minimal distancebetween them is equal to the shifting parameter of the coalescence kernel consist-ing of two simple rectangle functions with the ranges which do not exceed the sizes of the islands, then a stationary state with inhomogeneous particle distribu-tions can arise. The inclusion of any jumps even with extremely small intensityradically changes the situation, leading to the collapse of the system.The proposed algorithm is implemented in a Fortran program code which canbe received by request and exploit by any researcher free of charge. The code canbe adapted to more complex multicomponent models of jumping and coalescenceparticles of different types. The numerical simulations can be extended to systemsof higher dimensions. The Poisson approximation used for the derivation of thekinetic equation can be improved by advancing to the Kirkwood [26] or Fisher-Kopeliovich [27] ansatz like for birth-death models. Mass and size of particlescan also be taken into account. These and other topics as well as possible appli-cations of the obtained results to real populations will be the subject of our furtherinvestigations. A PPENDIX
In the absence of coalescence [ b ( x ) =
0] we can apply the convolution theoremto avoid the direct spatial integration. Indeed, then the kinetic equation (4) can berewritten in the form ∂ n ( x , t ) ∂ t = − n ( x , t ) Z a ( x − y ) g ( y , t ) dy + g ( x , t ) Z a ( x − y ) n ( y , t ) dy , (A1)where g ( x , t ) = exp (cid:16) − Z ϕ ( x − u ) n ( u , t ) du (cid:17) . (A2)To simplify notation consider the case d = f k = N − ∑ i = f i exp ( − π ıki / N ) ≡ { f i } k , f i = N N − ∑ k = ˜ f k exp ( π ıki / N ) ≡ { f k } − i (A3)where f i = f ( x i ) and x i = − L / + ( i − / ) h ∈ ] − L / , L / [ with even N and h = L / N . Now, applying the convolution theorem, the discrete counterpart of Eq. (A1)in view of Eq. (A2) can be presented as˙ n i = dn i dt = − hn i { ˜ a k ˜ g k } − i + hg i { ˜ a k ˜ n k } − i , g i = exp ( − h { ˜ ϕ k ˜ n k } − i ) (A4)where n i = n ( x i , t ) with ˜ n k , ˜ a 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 and ˜ f k ≡ ˜ n k , ˜ a k , ˜ ϕ k .Because the kernel function a ( x ) and repulsive potential ϕ ( x ) do not depend ontime, their Fourier components ˜ a k and ˜ ϕ k can be calculated once at the very begin-ning. Further, having the current values n i for i = , , . . . , N − n k for k = , , . . . , N −
1. On the basis of the known products˜ a k n k and ˜ ϕ k ˜ n k we evaluate their coordinate counterparts { ˜ a k ˜ n k } − i and { ϕ k N k } − i by exploiting the inverse Fourier transform. Constructing g i = exp ( −{ ˜ ϕ k ˜ n k } − i ) we calculate ˜ g k by means of the direct Fourier transform and use again the inverse LGORITHM FOR NUMERICAL SOLUTIONS TO THE KINETIC EQUATION OF COALESCING JUMPS23 transform to obtain { ˜ a k ˜ g k } − i . In such a way we have all necessary components toform the time derivatives ˙ n i according to Eq. (A4). Despite this requires severaldirect and inverse transformations, the order of total number of operations can bereduced from N (when the spatial integration is performed directly, see Sect. 3)to N ln N when exploiting the fast discrete direct and inverse Fourier transforms.This is very important feature because N ≫ N ∼ – 10 at least). However, such an approach can work not so well for discontinuous func-tions. Moreover, it assumes that all involved functions are periodic on the interval [ − L / , L / ] .The generalization to any higher dimensionality d > 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 Eq.(A3) from 0 to N − f i = f ( x i , x i , . . . , x i d ) . As a result,the total number of operations increases from N to N d or from N ln N to N d ln N d when applying the product of d one-dimensional usual or fast Fourier transformin Cartesian coordinates. If the functions possess radial symmetry, we can con-sider the Fourier transform in polar or spherical coordinates at d = d -dimensional integrations to discrete Hankel [29]or Fourier-Bessel [30] transforms in one (radial) dimension by integrating analyt-ically out all angle dependencies. The chief advantage of this approach is the factthat then the total number of operations becomes independent of d and remains thesame (of order of N ln N ) for any dimensionality [31], just as at d = n ( x , t ) = n ( t ) does not changeon x , we can carry out the spatial integrations in Eqs. (4) and (5) over y explic-itly. This significantly simplifies the kinetic equations to the form dn / dt = − µ b n which gives the analytical solution n ( t ) = n ( ) / [ + µ b n ( ) t ] . It tends to zero with t → ∞ for any initial density n ( ) and µ b =
0. In this case the jump contribu-tion completely disappears. At µ b = n ( t ) = n ( ) . It disappears also in the spatially inhomogeneous case for systemswith finite number N = R n ( x ) dx of particles when calculating d N / dt (then dur-ing integration in the rhs of Eq. (3) over x , the first two terms are mutually can-celled). R EFERENCES [1] Arratia, R.A.: Coalescing Brownian motion on the line. PhD thesis, University of Wisconsin,Madison, ProQuest LLC, Ann Arbor, MI (1979)[2] T´oth, B., Werner W.: The true self-repelling motion. Probab. Theory Relat. Fields 111, 375–452(1998)[3] Le Jan, Y., Raimond, O.: Flows, coalescence and noise. Ann. Probab. 32, 1247–1315 (2004)[4] Konarovskii, V.V.: On an infinite system of diffusing particles with coalescing. Teor. Veroyatn.Primen. 55, 157–167 (2010)[5] Berestycki, N., Garban, Ch., Sen, A.: Coalescing Brownian flows: a new approach. Ann. Probab.43, 3177–3215 (2015)[6] Konarovskii, V.V., von Renesse, M.: Modified massive Arratia flow and Wasserstein diffusion.Comm. Pure Appl. Math. 72, 764–800 (2019) [7] Pilorz, K.: A kinetic equation for repulsive coalescing random jumps in continuum. Ann. Univ.Mariae Curie-Skłodowska Sect. A 70, 47–74 (2016)[8] Kozitsky, Yu., Pilorz, K.: Random jumps and coalescence in the continuum: evolution of statesof an infinite system. arXiv:1807.07310 (2018)[9] Bara´nska, J., Kozitsky, Yu.: The global evolution of states of a continuum Kawasaki model withrepulsion. IMA J. Appl. Math. 83, 412–435 (2018)[10] Berns, C., Kondratiev, Y., Kozitsky, Y., Kutoviy, O.: Kawasaki dynamics in continuum: micro-and mesoscopic descriptions. J. Dynam. Differential Equations 25, 1027–1056 (2013)[11] Capit´an, J.A., Delius, G.W.: Scale-invariant model of marine population dynamics. Phys. Rev.E 81, 061901 (2010)[12] Law, R., Plank, M.J., James A., Blanchard J.L.: Size-spectra dynamics from stochastic preda-tion and growth of individuals. Ecology 90, 802–811 (2009)[13] Weisse, T.: Dynamics of autotrophic picoplankton in marine and freshwater ecosystems. In:Jones, J.G.(ed.) Advances in Microbial Ecology, vol. 13, pp. 327–370. Plenum Press, New York(1993)[14] Rudnicki, R., Wieczorek, R.: Fragmentation-coagulation models of phytoplankton. Bull. Pol.Acad. Sci. Math. 54, 175–191 (2006)[15] Rudnicki, R., Wieczorek, R., Phytoplankton dynamics: from the behaviour of cells to a trans-port equation. Math. Model. Nat. Phenom. 1, 83–100 (2006)[16] Ambrose J., Livitz M., Wessels D., Kuhl S., Lusche D.F., Scherer A., Voss E., Soll D.R.:Mediated coalescence: a possible mechanism for tumor cellular heterogeneity. Am. J. Cancer Res.5, 3485–3504 (2015)[17] Wessels D., Lusche D.F., Voss E., Kuhl S., Buchele E.C., Klemme M.R, Russell K.B., AmbroseJ., Soll B.A., Bossler A., Milhem M., Goldman C., Soll D.R.: Melanoma cells undergo aggressivecoalescence in a 3D Matrigel model that is repressed by anti-CD44. PLoS One. 12, e0173400 (2017)[18] Kozitsky, Yu., Omelyan I., Pilorz K.: Jumps and coalescence in the continuum: a numericalstudy. Appl. Math. Comput. [submitted] (2019)[19] Finkelshtein, D., Kondratiev, Y., Oliveira, M.J.: Markov evolutions and hierarchical equationsin the continuum. I. One-component systems. J. Evol. Equ. 9, 197–233 (2009)[20] Finkelshtein, D., Kondratiev, Y., Kutoviy, O.: Vlasov scaling for stochastic dynamics of con-tinuous systems. J. Stat. Phys. 141, 158–178 (2010)[21] Banasiak, J., Lachowicz, M.: Methods of Small Parameter in Mathematical Biology.Birh¨auser/Springer, Cham (2014)[22] Finkelshtein, D., Kondratiev, Yu., Kozitsky, Yu., Kutoviy, O.: The statistical dynamics of aspatial logistic model and the related kinetic equation. Math. Models Methods Appl. Sci. 25, 343–370 (2015)[23] Finkelshtein, D., Kondratiev, Y., Kutoviy, O.: Statistical dynamics of continuous systems: per-turbative and approximative approaches. Arabian Journal of Mathematics 4, 255–300 (2015)[24] Kondratiev, Yu., Kozitsky, Yu.: The evolution of states in a spatial pupulation model. J. Dynam.Differ. Equ. 30, 135–173 (2018)[25] Chapra, S.C., Canale, R.P.: Numerical Methods for Engineers, 7th edn. McGraw-Hill Educa-tion, Penn Plaza, New York (2015)[26] Omelyan, I., Kozitsky, Yu.: Spatially inhomogeneous population dynamics: beyond the meanfield approximation. J. Phys. A.: Math. Theor. 52, 305601 (2019)[27] Omelyan, I.: Spatial population dynamics: beyond the Kirkwood superposition approximationby advancing to the Fisher-Kopeliovich ansatz. arXiv:1907.00223, Physica A [submitted] (2019)[28] Press, W.H., Teukolsky S.A., Vetterling W.T., Flannery B.P.: Numerical Recipes, The Art ofScientific Computing, 3rd edn., Cambridge University Press, Cambridge (2007)[29] Baddour, N.: Operational and convolution properties of two-dimensional Fourier transforms inpolar coordinates. J. Opt. Soc. Am. A 26, 1767–1777 (2009)[30] Baddour, N.: Operational and convolution properties of three-dimensional Fourier transformsin spherical polar coordinates. J. Opt. Soc. Am. A 27, 2144–2155 (2010)
LGORITHM FOR NUMERICAL SOLUTIONS TO THE KINETIC EQUATION OF COALESCING JUMPS25 [31] Nikitin, A.A., Nikolaev M.V.: Equilibrium integral equations with Kurtosian kernels in spacesof various dimensions. Mosc. Univ. Comput. Math. Cybern. 42, 105–113 (2018)I
NSTITUTE FOR C ONDENSED M ATTER P HYSICS , N
ATIONAL A CADEMY OF S CIENCES OF U KRAINE , 1 S
VIENTSITSKII S TREET , UA-79011 L
VIV , U
KRAINE
E-mail address : [email protected] I NSTYTUT M ATEMATYKI , U
NIWERSYTET M ARII C URIE -S KODOWSKIEJ , 20-031 L
UBLIN ,P OLAND
E-mail address : [email protected] I NSTYTUT M ATEMATYKI , U
NIWERSYTET M ARII C URIE -S KODOWSKIEJ , 20-031 L
UBLIN ,P OLAND
E-mail address ::