Numerical coalescence of chaotic trajectories

Pairs of numerically computed trajectories of a chaotic system may coalesce because of finite arithmetic precision. We analyse an example of this phenomenon, showing that it occurs surprisingly frequently. We argue that our model belongs to a universality class of chaotic systems where this numerical coincidence effect can be described by mapping it to a first-passage process. Our results are applicable to aggregation of small particles in random flows, as well as to numerical investigation of chaotic systems.


Introduction
When we numerically compute two distinct trajectories of a deterministic system, after a certain number of iterations their coordinates may happen to be exactly equal, to machine precision.At this point all subsequent computations of the two trajectories yield exactly the same values.We might say that the trajectories have undergone numerical coalescence.This coalescence effect is readily observed in systems with stable dynamics, where nearby trajectories approach each other, typically with an exponentially decreasing separation.It might, however, be expected that it would be unusual to see this numerical coalescence in chaotic systems, where nearby trajectories have an exponentially growing separation (the concept of chaos in the theory of dynamical systems is discussed in [1,2]).In this paper we analyse numerical coalescence in a simple dynamical system which shows a transition to chaos, exploring its dependence upon the numerical precision of the calculation, which is denoted by δ.In our numerical work we used variable precision arithmetic implemented in the mpmath package [3] for the Python programming language, and comparable results were obtained using the Maple mathematical software system [4].These packages allow the number of decimal digits, M, to be set to an arbitrary integer value.If the typical magnitude of the numbers representing the trajectory is x, we may write δ ∼ x 10 −M . ( Numerical experiments showed that coalescence is surprisingly frequent.We were able to explain this observation, showing that our numerical results are in agreement with a calculation which maps the numerical coalescence to a first-passage problem [5] for a simple stochastic process.We argue that our results are representative of a physically significant universality class of chaotic systems. As well as being relevant to the numerical exploration of complex dynamical systems, the effect can also serve as a model for aggregation of small particles in complex flows (reviewed in [6,7]).If the particles aggregate when their separation falls below a threshold, δ, which is analogous to computing trajectories with a finite numerical precision.
Our explanation of the numerical coalescence phenomenon is related to earlier studies [8,9], where it was shown that trajectories of a chaotic system can display a clustering effect, which is characterised by a power-law distribution of their separations.We review the relationship between the results in this paper and those earlier works in our conclusions, section 5.

An example
We start with a simple example illustrating the effect that we wish to explain and analyse.Consider trajectories of a one-dimensional map in the case where F is a function which is continuous and differentiable almost everywhere, and which has unit period, F (x + 1) = F (x).The φ n are random numbers, chosen independently at each time step, with a distribution which is uniform on [0, 1].We also found it convenient to require that the mean value of F (x) is equal to zero.Note that, while a single trajectory of (2) executes a random walk, this equation describes a differentiable dynamical system, so that it is possible to calculate the Lyapunov exponent [1].A system with similar properties to (2) is discussed in some detail in [10].
In our study the function F (x) was piecewise linear, having a sawtooth form illustrated by figure 1, with gradients ±g on intervals of length 1/2, With this choice of F (x), an individual trajectory x n executes a random walk for which the variance of the step sizes is g 2 /48, implying that where • denotes expectation value.
Because the map is periodic, with period unity, if two trajectories x n and x ′ n differ by an integer, then x n − x ′ n remains equal to the same integer for all subsequent iterations.Note that the map has multiple pre-images when g > 1.We shall assume that g > 1 throughout our discussion of the map defined by (2) and (3).
x In section 3 we demonstrate that the map is chaotic (because the Lyapunov exponent, λ, is greater than zero [1]) for g > √ 2. When λ > 0, nearby trajectories should separate exponentially [1]: if ∆x n is the separation of two trajectories after n iterations, then approaches the Lyapunov exponent λ as n → ∞, provided ∆x 0 is sufficiently small that |∆x n | ≪ 1.When λ < 0, trajectories converge to the same point (if their initial  1.Two trajectories of a realisation of (2) after n = 2.5 × 10 5 iterations, computed with M = 10 digit arithmetic.One of the initial conditions was random in [0, 1], and the initial separation was ∆x 0 = 0.25 in all cases.
separation is small) or else to integer separations due to the periodic nature of the onedimensional map.When λ > 0, it might be expected that a very small initial separation of trajectories eventually grows to be of order unity, and the similarity of (2) to the equation of a random walk leads one to expect that subsequent growth of separation is diffusive, ∆x ∼ √ n.
However, numerical investigations of (2) show a different behaviour.Table 1 lists the separations for pairs of trajectories with x 0 random and an initial separation of ∆x 0 = 0.25 after n = 2.5 × 10 5 iterations of the map, for a range of values of g.Here we used an arithmetic precision which was defined by specifying M = 10 decimal places.The numerically computed separations of the two trajectories, x n and x ′ n , are either exactly zero or exactly integers, for two values of g which are in excess of the threshold for chaos, at g = √ 2 = 1.41421 . ... It is easy to understand why this numerical coalescence effect is possible for values of g beyond the threshold for chaos.When λ > 0, it is still possible for trajectories to coalesce if two trajectories happen to reach exactly the same point (this is possible because the calculation uses finite precision arithmetic, and because points may have multiple pre-images).However, we might expect that this occurrence is rare.Above the threshold for chaos, nearby points are separating exponentially.If these separations were to fill the unit interval with constant density, there would be a probability P ∼ δ to coalesce at each iteration of the map.If the probability to coalesce were δ upon each iteration, then after n iterations the probability P c for the paths to have undergone coalescence would be In our numerical investigation illustrated by table 1, n = 2.5 × 10 5 , implying (from equation (4), using g = 2 as a representative value) that x = √ 2Dn ≈ 100.This would give an estimate for the probability for coalescence of P c ∼ 2.5 × 10 5 × 100 × 10 −10 = 2.5×10 −3 .According to this estimate, the results displayed in table 1 appear to be highly unlikely.The remainder of this paper will explain and quantify the effect illustrated in table 1, (section 3), and discuss the extent to which it is a manifestation of a universal phenomenon (section 4).

Numerical investigation
When two trajectories, started from randomly chosen initial conditions, reach exactly the same coordinate (or else two coordinates with an exactly integer separation), their values remain locked together for all subsequent iterations.If the arithmetic of a numerical iteration of the map has a finite precision, say M decimal places, then this phenomenon of numerical coalescence provides an explanation for the effect illustrated in table 1.The principal difficulty is to explain why the effect happens so much more frequently than the simple estimate, equation (6).
We can take two trajectories and determine the 'time' N (that is, number of iterations) for them to coalesce.This will be different for different initial conditions.We should therefore look at the probability P (N), or else at statistics such as the moments N j .The simplest statistic is just the mean time for coalescence, N , and the remainder of this section is concerned with estimating this quantity.
We investigated the mean number of iterations for coalescence of trajectories N as a function of g and of the number of decimal digits, M. The results are presented in table 2. The initial separation was ∆x 0 = 0.1 in all cases, and for all of the data points we averaged over 1000 realisations, with initial conditions distributed randomly with uniform density on [0, 1].
Note that most of these values of N are sufficiently small that x ∼ 2D N is not a very large number: for example at g = 1.4 and M = 20, we found N ≈ 2500, implying that x ≈ 10, which indicates that only one of the M = 20 digits is required to store the integer part of x n leaving 19 digits after the decimal point.For this reason we our analysis of these data will make the simplifying assumption that the precision of the calculation, δ, is given by δ = 10 −M , where M is the number of digits.

Theory: relation to mean first-passage time
Next we consider how to analyse the results in table 2. Because the map (2) has unit period, two trajectories which have an exactly integer separation maintain the same separation for all subsequent iterations.In this sense, integer separations of trajectories are equivalent to a zero separation, and accordingly we define ∆x n as the magnitude of the separation of two trajectories modulo unity.In order to facilitate the calculation of N , the separation of two trajectories, ∆x, is transformed into a logarithmic variable, If trajectories were computed with arbitrary precision, the variable Z would occupy the interval Z ∈ [0, ∞).When Z is large, the separation of two trajectories is very small,  2. Tabulation of the mean number of iterations before coalescence, N , determined numerically for different values of the coefficient g and of the number of digits, M .The Lyapunov exponent λ, given by the first line of equations (9), is positive when g > √ 2 = 1.41 . ... and (at points where the derivative of the map exists) the iteration of Z is described by linearisation of the map, so that For the map (2), the dynamics of ( 8) is Markovian, with displacements ∆Z ± = − ln(1±g) occurring with random choices of the sign, having probabilities P ± = 1/2.Because ( 8) is the equation of a random walk, over many iterations the motion of Z can be modelled by an advection-diffusion equation, with a drift velocity v and a diffusion coefficient D (see [11] for a discussion of the relationships between random walks or Langevin processes and diffusion or Fokker-Planck equations).Note that, from the definition of the Lyapunov exponent λ [1], we have λ = −v.The values of v and D are determined from the statistics of ∆Z: if time is measured by the number of iterations, the drift velocity v and diffusion coefficient D are obtained from The first of these relations shows that λ = ln g 2 − 1, which is positive when g > √ 2. When ∆x is small, implying that Z is large, the behaviour of Z is determined by the linearised equation of motion, leading to (8).In the vicinity of Z = 0, however, Z has complex dynamics which could, in principle, be determined from the equation of motion of x, (that is, from equation ( 2)).However we can observe that the representative point Z does not pass Z = 0, and eventually will leave the region close to Z = 0.This implies that the diffusion process can be modelled as having a reflecting barrier at Z = 0.If we take account of the coalescence of trajectories due to finite-precision representation of arithmetic, this implies that the separation of trajectories becomes zero when ∆x = δ = 10 −M where δ is the floating point precision.At this point, the diffusive model ceases to be appropriate.We can represent this by introducing an absorbing barrier at a position The mean number of iterations before coalescence, N , is the same as the mean time for first contact with the absorbing barrier in the diffusion process.This quantity is known as the mean first-passage time.

Calculation of mean time to coalescence
There is an extensive literature on the mean first-passage time for diffusive processes [5].We shall use a standard formula for the mean first-passage time for the Langevin equation where V (Z) is a potential and η(t) is white noise.We assume that there is a reflecting barrier at Z 1 , absorbing barrier at Z 0 , and particles initially released from Z i .The mean first-passage time is (see [12,13]) As Z approaches ∞, the mean velocity v is becomes independent of Z, and is equal to −λ, where λ is the Lyapunov exponent.The dynamics of Z is also subject to fluctuations about this mean motion (which are quantified by the diffusion coefficient D).If the potential is V (x) = −vx, then the drift velocity, −V ′ (x), is a constant, v, as required.The true boundary condition at Z = 0 could be modelled more faithfully by introducing a delay time, but we shall assume that using a reflecting boundary at Z = 0 is sufficient.We could, in principle, average over the initial conditions by averaging Z i over the steady-state distribution of Z, but again we adopt a simplifying assumption, assuming all of the representative points are injected at Z i = 0 (we expect that this approximation will cause our calculation to overestimate the time taken to reach the absorbing barrier).Defining and setting V (x) = −vx in equation (12), we obtain x 0 dy exp(αy) .
Evaluating the integrals gives Noting that the expected number of iterations for path coalesce is N = T , this can be written in the form where and Because this approximation assumes |Z 0 | ≫ 1, the predicted values of N are unobservably large or small unless g is close enough to g 0 = √ 2 that we can make the following two approximations: where g 0 = √ 2 is the parameter value for the threshold of chaos, and D 0 is the value of the diffusion coefficient at g 0 .Using these approximations we find The function f (X) is positive for all real X.It implies the following limiting behaviours: The first of these is what would be expected from the definition of the Lyapunov exponent.The value of N should not exceed δ −1 = 10 M , so that exponential growth in the limit as λ → ∞ which is predicted by (21) will only be correct when g − g 0 is sufficiently small.Equations (20) are expected to give an asymptotic approximation to N , which is valid when M ≫ 1 and g− √ 2 ≪ 1.We compared this theory against the data tabulated in table 2 by testing whether it would collapse onto a single curve, representing the function f (X).Because of the exponential growth of f (X) for positive X, it is more convenient to graph ln f (X).Accordingly, in order to test the prediction contained in equations ( 16) to (18), we made a plot of as a function of X, as defined by equation (20).In figure 2 we display our plot of Y versus X for the data in table 2. We use different point styles (and colours, online) to distinguish the values of M, and included some additional data for more densely sampled values of the coefficient g.The points with different colours all collapse onto the same curve, which is well approximated by the function Y (X) = ln f (X), plotted as a solid curve.This validates the theory described by ( 16), ( 17) and (18) as a description of the finite-precision path-coalescence effect for map (2) and (3).because a stable periodic orbit either disappears, or else undergoes a bifurcation, the ∆Z n form a periodic sequence when g is on the stable side of the transition.In this case the diffusion coefficient D is equal to zero in the stable phase, and we must have D 0 = 0, so that the theory is not applicable.

Conclusions
We have exhibited an example of a chaotic system where trajectories coalesce at a surprisingly high rate in the vicinity of the threshold of chaos due to arithmetic truncation errors.The effect was explained by a transformation of the separation of two trajectories to a logarithmic variable, which leads to analogy with a first-passage problem for a diffusive process.This effect is a consequence of transient convergence of chaotic trajectories, which was previously analysed in [8] and [9].Our analysis of the numerical coalescence phenomenon is based upon a principle that was used in those earlier works, namely that the linearised equation of motion is mapped to an advectiondiffusion equation if we make a logarithmic transformation, equation (8) in this present work.Here we have shown that the numerical coalescence effect may be understood using a surprisingly simple treatment of the associated first-passage problem.
Our analysis is also relevant to understanding the aggregation of particles in complex flows where the particles aggregate when their separation is equal to δ.In particular, our theory for N describes the mean time for collision of particles which are advected by a flow which has a Lyapunov exponent which is close to zero.
The theory that we have presented is valid when the stability factors ∆Z n have random or pseudo-random fluctuations.This includes almost all real-world applications of chaotic dynamics, where noise will be present.In section 4 we pointed out that the effect will may be absent in systems where the transition to chaos arises from perturbation of a periodic orbit.This is consistent with a viewpoint that lowdimensional autonomous dynamical systems may not be a good models for physical applications of chaos.