Coagulation-fragmentation model for animal group-size statistics

We study coagulation-fragmentation equations inspired by a simple model proposed in fisheries science to explain data for the size distribution of schools of pelagic fish. Although the equations lack detailed balance and admit no $H$-theorem, we are able to develop a rather complete description of equilibrium profiles and large-time behavior, based on recent developments in complex function theory for Bernstein and Pick functions. In the large-population continuum limit, a scaling-invariant regime is reached in which all equilibria are determined by a single scaling profile. This universal profile exhibits power-law behavior crossing over from exponent $-\frac23$ for small size to $-\frac32$ for large size, with an exponential cut-off.


1.
Introduction. Most animals in nature aggregate in groups of different sizes. These sizes vary in their frequency and obviously depend on the species. So the question arises whether and how typical distributions of group sizes emerge. Related questions are: Can we find adequate models for these distributions? How do the distributions evolve over time? Is there an (or several) equilibrium distribution(s)? Can one say something about the trend towards these equilibria?
Various models of describing the coagulation and fragmentation of groups of animals have been suggested and analysed in the past (cf. e.g. [1,2,8,9,19]). The model this work rests upon was introduced by Hiro-Sato Niwa in 2003 [17] related to studies in [15,16,18] and has turned out to hold for data from pelagic fish and mammalian herbivores in the wild. The model can be formalized into coagulationfragmentation integral equations where the coagulation rate is a constant independent from the group sizes and the fragmentation rate is also a constant independent from the fragment. By analogy with an Itô Stochastic Differential Equation Niwa shows that the equilibrium must be given by (1.1) W (N ) is the stationary probability density function of group sizes and N P is the average of the population distribution among group sizes, i.e. the expected size of the groups which an arbitrary individual is part of. For a continuum of cluster sizes, this is defined as In the discrete setting the integrals are replaced by sums.
In [17] Niwa shows that the proposed equilibrium distribution (1.1) matches empirical data of several species of pelagic fish very well. Ma et al. [12] provide a critical discussion of Niwa's result and point out some obscurities in the analysis. Due to the appealing simplicity of Niwa's model and the good empirical match to the data, mathematical clarification is important. Degond et al. [5] have pursued with Niwa's model and given a rigorous description of the equilibria for continuous (model C) and discrete (model D) cluster sizes, which differ from (1.1). The lack of a detailed balanced condition has made the analysis difficult. However, by introducing the so called Bernstein transformation, they have shown that there exists a unique equilibrium, under a suitable normalization condition, for both the discrete and the continuous cluster size case.
The task of the present work is a numerical investigation of both models and their equilibria. The continuous equilibrium is approximated numerically using three different methods whose accuracy will be examined. One of them is a recursive algorithm derived from model D in [12] which enables a transition from the discrete to the continuous equilibrium. The other two, a Newton and a time-dependent method, operate within a discretized truncated model, denoted by D', of the continuous model C. There is an abundant amount of literature about discretizations of coagulation (and fragmentation) integral equations using finite volume methods (e.g. [3,6,7,10,11,21]) or finite element methods (e.g. [13,14,20]). In our case, the discretization scheme is already predetermined by model D.
It is investigated how well the numerically generated equilibria match the analytically predicted decay rate and the small-size asymptotic behaviour of the model C equilibrium. We find all three methods to be very accurate apart from small deviations of the large-size behaviour in the case of the Newton and the time-dependent method due to truncation. The Newton method turns out to be extremely fast, providing a very close approximation of the equilibrium after five iterations. The recursive algorithm is the best numerical approach to this particular model with respect to a couple of aspects: it is numerically cheap, doesn't require truncation, is completely accurate for the discrete model D and approximates the continuous case properly without any aberrations. However, the other two methods are far more flexible regarding changes of the models since, in principal, they don't require constant coagulation-fragmentation parameters p and q as opposed to the recursive algorithm. The Newton scheme as an approach to prove the existence and uniqueness of the equilibrium, as introduced in this work for model C', has the advantage of not depending on fixed parameters as contrasted with the Bernstein method (see [5]) which needs p and q to be equal to one.
Hence, the truncated model and the associated numerical methods provide the tools to work in more sophisticated models with the coagulation and fragmentation depending on the group sizes and/or time. In this context, the model under investigation serves as a toy model to show the accuracy of the suggested schemes.
In addition to that, the Euler scheme helps to examine the convergence of timedependent solutions to the stationary one, something that hasn't been understood comprehensively in the analysis in [5]. The numerical approach indicates uniform convergence on finite intervals with super-exponential convergence rates independent from the group sizes.
We introduce model C and model D in Section 2. In Section 3 we summarise the analytical results concerning equilibria in models C and D. We introduce our own truncated model C' and the constructive approximation (Newton) method to its equilibrium in Section 4. Section 5 provides a description of the different numerical algorithms whose validations and insights are shown in Section 6.
2. The governing equations: The continuous and the discrete version. In weak form it reads, for all test functions ϕ in a class to be specified later: The coagulation rate a(x, y) and fragmentation rate b(x, y) are both nonnegative and symmetric. The coagulation and fragmentation reactions can be written schematically By a change of variables, (2.1) can be transformed into Note that by taking ϕ(x) = x, one obtains the conservation of mass The intuition behind (2.1) becomes clearer when we consider the strong form. In the following, Q C shall denote the coagulation operator and Q F the fragmentation operator. They both have a gain and a loss component and build up the strong form of the equation as ∂f ∂t The case where the cluster sizes form a discrete set can be described analogously. So consider a system of clusters with discrete sizes i ∈ N. Merging and splitting with the coagulation rate a i,j and fragmentation rate b i,j are ruled by the following coagulation-fragmentation reactions The system is described by the number density f i (t) of clusters of size i at time t evolving according to the discrete coagulation-fragmentation equation. Written in weak form the equation reads for any test sequence ϕ i in a class to be specified later The equation can also be written similarly to (2.2) as (2.8) If one takes ϕ k = k, it can be seen immediately that mass is conserved: Let Q Ci and Q F i denote the coagulation resp. fragmentation operator for cluster size i. Then the strong form can be written as The equations based on Niwa's model. According to Niwa's model, we assume s different zones of space on which Φ individuals move. The number of individuals is conserved through time. At each time step every group moves towards a randomly selected site with equal probability. When i-and j-sized groups meet at the same site, they aggregate to a group of size i + j. So the coagulation rate is independent from the group sizes and can be written as a i,j = 2q for any i, j > 0 where q > 0 is the fixed coagulation parameter. The fragmentation rate b i,j expresses the fact that at each time step each group with size k ≥ 2 splits with probability p independent of k, and that if it does split, it breaks into one of the pairs with sizes (1, k − 1), (2, k − 2), . . . , (k − 1, 1) with equal probability. As the actually distinct pairs are counted twice in such an enumeration, one gets for all Summarizing, we can express Niwa's model in the discrete system of equations introduced above by choosing As already indicated, Ma et al. [12] have studied the coagulation-fragmentation system with these rates. Gueron and Levin [9] had proposed coagulation and fragmentation rates that satisfied a detailed balance condition. That means that their choice of a and b was such that there exists an equilibrium distribution f fulfilling The detailed balance condition is not satisfied in Niwa's model (cf. [5,Section 7]). Degond et al. have chosen the same fragmentation and coagulation rates as Niwa in the continuous case but slightly different ones in the discrete case. The results of these steps are the discrete model D and the continuous model C, as described below: • Model D (Discrete): : (2.14) The fragmentation of a group of size k in Model D can now be understood as breaking into the pairs (0, i), . . . , (i, 0) with equal probability 1/k + 1. This means that we also consider the cases in which actually nothing changes. This results in a significantly simpler analysis. To summarize, we will consider the following models: Model D. The weak form for Model D (derived from (2.8)) reads, for all bounded test sequences ϕ i , The test space is chosen to be l ∞ because of the finite mass assumption on the solutions which is specified in Section 3.2. The strong form becomes Model C. The continuous model C can be written in weak form, for any test function ϕ ∈ C([0, ∞)), as x + y dxdy.
(2. 19) or as (2.20) The test space is chosen to be C([0, ∞)) because of the finite mass assumption on the solutions which is specified in Section 3.1. The strong form can be written as By introducing the method of Bernstein transformations, the existence and uniqueness of an equilibrium can be shown. The following section summarizes the important findings of [5], and prepares us for the numerical investigation.
3. Preliminary findings in the analysis of the coagulation-fragmentation model. This section contains a summary of important results from [5].
3.1. Equilibrium in the continuous case. Let k ∈ N and f : x ∈ R + → f (x) ∈ R + . The kth moment m k (f ) is given by For initial condition There is a scaling invariance for model C ((2.19) -(2.23)): with m 1 (f 0 ) =: m 1 < ∞ and let f p,q (x, t) be the solution of (2.21) with parameters p and q. Then, one has Due to this proposition, we can assume p = 1, q = 1 and m 1 = 1. The problem in strong form becomes In weak form it reads as The main theorem can be stated as follows: There is a unique equilibrium distribution function f ∞ of (3.1) and (3.3) satisfying (3.2). It can be written as where γ is a completely monotone function and has the following asymptotic behaviour: γ(x) ∼ 1 3Γ(4/3) x −2/3 as x → 0, (3.4) γ(x) ∼ 9 16Γ(3/2) x −3/2 as x → ∞. 3.2. Equilibrium in the discrete case. One can show that there is a scaling invariance for model D ((2.15)-(2.18)) as well (cf. [5,Section 2.3]). Hence, we will work with p = q = 1 in the following. In the discrete setting the kth moment of a sequence f = (f i ) i∈N is given by Let us further introduce the sets Let F t (x) denote a time-dependent solution of the continuous model C ((2.19)-(2.23)). For a transition from the continuous to the discrete model, we introduce a grid size h > 0 and the approximation

PIERRE DEGOND AND MAXIMILIAN ENGEL
for the number of clusters with sizes in the interval I h i . For a smooth test function ϕ(x), we can write ϕ i = ϕ(ih) and require that Letting h → 0, leads to an approximation of the continuous model by the discrete one. Since the physically attainable scenario in the case of animal group size statistics is given by h = 1, one might wonder why smaller grid sizes are of interest. Clearly this is the case not just for mathematical reasons but also due to the fact that the coagulation-fragmentation equation can be used to model a variety of phenomena where the cluster sizes are not necessarily integer numbers. Understanding the relation between the continuous and discrete model for small h is therefore an important part of this paper.
We define the zeroth and first moment of an equilibrium distribution by The following theorem tells us that such an equilibrium actually exists and gives details about the asymptotic behaviour (cf.
where γ is a completely monotone sequence with the asymptotic behaviour Further, the following mass-number relation holds: Complete monotonicity in the discrete context means that where the difference operator ∆ is given by Let Let again F t (x) be a solution of the continuous model C ((2.19)-(2.23)). One can show that for a certain correspondence of initial data, for each t > 0, we have F h t → F t narrowly as h → 0 (for a proof see [5,Theorem 16.1]). In the following, we want to approximate these equilibria numerically. We are going to apply three different methods. The one based on model D ((2.15)-(2.18)) will rest upon a recursive algorithm introduced in Section 5.1. The other two, a Newton and a time-dependent method, require for a truncation in model C ((2.19)-(2.23)) onto a compact interval of R. This new model C' will be treated in the next section.
4. Model C': A truncated version of model C.

4.1.
The time-dependent problem. We introduce a truncation of the weak formulation of model C to the interval [0, L]. Let ϕ be a test function. The truncation is chosen as follows: for all t > 0 and test functions ϕ.
Note that, indeed, by chosing ϕ(x) = x mass conservation is still obtained: Proposition 4.2. Let Q C T denote the coagulation operator and Q F T the fragmentation operator. Then the strong form of model C' can be written down as: Proof. Obvious calculation.
Further, we can state the following local existence and uniqueness result: Then there is an α > 0 such that the initial value problem corresponding with (4.2) Proof. This is an immediate application of the Cauchy-Lipschitz Theorem for initial value problems in Banach spaces as Q C T is a continuous quadratic and Q F T is a continuous linear operator from L 1 ([0, L]) to itself (cf. Lemma 4.4).

4.2.
The equilibrium: A constructive approximation method. We present a constructive approach to find the equilibrium in model C' ((4.1)-(4.4)). It relies on a Newton method. The stationary version of (4.2) is This equation can also be written as T is a linear operator whereas q is a bilinear form with Starting with an appropriate f 0 , we want to find a recursive scheme giving a convergent sequence (f n ) n∈N with limit f ∞ , the equilibrium. Observe the following: If f n+1 was an equilibrium, we'd have Hence, the following Newton scheme rests upon neglecting this quadratic term and defines a sequence (f n ) n∈N by iteratively solving the following linear problem: (4.8) Introducing δf = f n+1 − f n , by adding −T f n and 2Q C T (f n ) on both sides of equation (4.8), we get We introduce the notation where W f n is a linear operator and G n is a function. W φ can be written as where In the following, R (W φ ) denotes the range of W φ and N (W φ ) its null space.
. Then W φ , as given in equations (4.10), (4.11), is a bounded, linear operator from Proof. The Lemma follows immediately from the definitions.
In addition, we can find out the following about the range of W f n (We choose the index f n instead of φ in order to build on equation (4.9)): Proof. For any test function ϕ we have By adding and subtracting T f n (x) and Q C T (f n )(x), one can see that Since this is true for any δf and the first two summands can be cancelled due to (4.12), for any λ > 0 it holds that Extracting the λ and dividing by λ leaves the factor λ on the right hand side of the equation. Due to arbitrariness of λ, it can be chosen arbitrarily small which shows that the left hand side is zero. Now, we conjecture the following based on Fredholm theory (cf. [4]): Proving this conjecture allows to single out the solution of W φ f = g by imposing xf dx = 1. This is the subject of current work. A natural question is how large L has to be such that the equilibrium of model C' ((4.1) -(4.4)) is reasonably close to the equilibrium profile f ∞ of model C ((2.19) -(2.23)). Recall from Theorem 3.3 that where γ is a completely monotone function with The analysis in [5] suggests that for L ∼ 10 the function γ is already approximated very well by these asymptotics. This means that for L > 10 If we choose L = 100, as will be done in the following numerical experiments, the neglected mass of the equilibrium is of order e −15 . This indicates that the truncation error is indeed very small.
5. Numerical methods. This section contains three numerical methods to approach an equilibrium distribution. The first one concerns a recursive computation of the equilibrium sequence for model D ((2.15)-(2.18)) already proposed in [12] and [5]. The other approaches rely on model D', a discretised version of truncated model C' ((4.1)-(4.4)). The first one simulates the evolution of the size distribution in time via an explicit Euler scheme and shall reach the steady state after a certain time span. The other one follows the Newton method theoretically outlined in Section 4.2. Note that the second method provides also an approximation of the time-dependent problem while the first and third methods only allow for the computation of the equilibrium. For a test function ϕ with ϕ i = ϕ(ih), the equilibrium profile satisfies Define now Taking ϕ j ≡ 1 yields Further, with taking ϕ k = 1 if k = i and 0 otherwise, we get Note that this is the stationary version of the strong form of model D given by (2.16) -(2.18) which is solved by the equilibrium profile. The specific choice of the test sequence ϕ j ≡ 1 in the weak form has lead to the relation (5.1) which is used to obtain the following recursive algorithm: Choose m h 0 ∈ (0, 1) and set Then for i = 1, 2, 3, . . . : for all test sequences φ i .
Observe that mass is preserved over time according to this equation.

5.2.2.
Time discretization of the time-evolution scheme. The explicit Euler scheme in time is applied with time step size ∆t. Let t k = k∆t. The sequence {f k } k∈N denotes an approximation of {f (t k )} k∈N and is defined by the following recursive scheme: where for any point ih with 1 ≤ i ≤ N , is given by A linear stability analysis of the scheme would mean to study the eigenvalues of the following operator L f , which is the Jacobian of the right hand side in (5.6): Since a spectral analysis of this operator is analytically impossible and numerically very costly, we restrain ourselves to adjusting the time-step recursively. Starting with ∆t = 10 −2 , the time step size is increased by ten per cent as long as the distribution stays non-negative and monotone. If one of these criteria is violated, the step size is reduced by ten per cent. The maximal time step size given by that scheme is ∆t = 1.1.

Equilibrium in model
Analogously to Eq. (4.5) involving the operators T and q, the discretized problem can be written as S is a linear operator and p is a bilinear form. Write P (f ) = p(f, f ). Hence, the task is to find f such that its image under the linear operator S equals its image under the quadratic form P derived from the bilinear form p. Following our considerations in Section 4.2, we apply the Newton method expressed by Eq. (4.8). Starting with an appropriate f 0 the following recursive scheme is applied: Sf n+1 − 2p(f n+1 , f n ) + P (f n ) = 0. The limit of this sequence, if it exists, satisfies the stationary equation (5.9).
Analogously to (4.9), the recursive scheme can be written as where we introduce the notation This equation can be written explicitly as We transfer our considerations concerning the invertibility of W f n in Section 4.2 to the discretised version V f n . Let x = (1, . . . , N ). The range of the operator is restricted to span{x} ⊥ , i.e. to N − 1 dimensions, and, hence, consider the above equation just for 1 ≤ i ≤ N − 1. Thereby we win a degree of freedom to implement the mass conservation in form of This scheme provides us with an algorithm to approximate numerically the solution of the stationary problem (4.5). As always, the performance of Newton's method crucially depends on the choice of the initialization. Here, we choose with m 1 > 0 denoting the mass to be chosen which will lead to convergence.
6. Numerical investigations. The numerical methods introduced in Section 5 shall now be applied. In the first subsection we check if the computed equilibrium distributions actually show the behaviour analytically predicted in [5]. Hence, we have to account for non-negativity and the predicted asymptotics for small and large sizes. We supplement the validation of the schemes by a comparison of the large-size asymptotics in model D ((2.15)-(2.18)) and model C ((2.19)-(2.23)). Further, we exploit the codes to gain new insights into the small-size behaviour in model D and the convergence rates to equilibrium in time. In the following, it will be appropriate to display the distributions mainly in a log scale using the decadic logarithm if not declared otherwise.
6.1. Validation of the numerical schemes.
6.1.1. The Newton method. First, we want to check the accuracy of the Newton method presented in the previous section. In particular, we will compare the predicted asymptotic behaviour with the asymptotic behaviour displayed by the computed equilibrium distribution. Recall from Theorem 3.3 that according to equation The following plots show that the approximation of the equilibrium generated by the Newton method matches the predicted asymptotic behaviour very well. First, we are interested in the asymptotic behaviour for large sizes. We choose m 1 = 1, truncation size L = 100 and h = 0.01. We perform five iterations. In Fig. 1a, the solid blue line shows the logarithmic distribution as a function of the group sizes whereas the dashed red line shows the predicted asymptotic behaviour for x → ∞ in (6.2). The distribution is chosen in a log scale while the group size is shown in a linear scale in order to illustrate the leading behaviour for the logarithmic distribution, − 4 27 x log 10 e, in a linear shape. Second, we focus on the small-size behaviour x → 0. We truncate at L = 5 and take h = 0.0005. Since for the case of m 1 = 1 and a calculation up to L = 100, the mass concentrated in [0, 5] equals 0.5676, we take this as our starting value for the mass. Again, we perform five iterations. In Fig. 1b the blue solid graph shows the log of the distribution as a function of the log of the group sizes whereas the red dashed graph shows the predicted asymptotic behaviour close to 0. These graphs show a linear behaviour consistent with the leading order term being given by −(2/3) log 10 x (see Eq. (6.1)).
Note that the distribution as shown in Fig. 1a tends to zero very quickly (already f ∞ (x) < 10 −2 at x = 10) but never becomes negative as intended. Observe the perfect convergence of both graphs for the group sizes becoming higher and higher. This means that the large-size asymptotic behaviour of the equilibrium generated by the Newton scheme is utterly accurate. There is a very small kink at the cut-off   Fig. 1a, we take mass m1 = 1, grid size h = 0.01 and the cut-off at L = 100. The plot shows the generated distribution (blue solid line) in a log scale against the group sizes in a linear scale and presents the theoretically found large-size asymptotic behaviour (red dashed line) in a log scale for the sake of comparison. The group sizes are taken in a linear scale in order to illustrate the leading behaviour for large group sizes as a straight line. For Fig. 1b, the equilibrium distribution is approximated by the Newton scheme taking mass m1 = 1, grid size h = 0.0005 and the cut-off at L = 5. The plot shows the generated distribution (blue solid line) in a log scale and presents the theoretically found asymptotic small-size behaviour (red dashed line) in a log scale for the sake of comparison. The group sizes are taken in a log scale as well in order to illustrate the leading behaviour close to zero as a straight line.
at L = 100. This is a consequence of the truncation. In model C' ((4.1)-(4.4)) the groups of size L = 100 cannot be part of coagulation into a group of bigger size, as opposed to model C ((2.19)-(2.23)) which is defined on [0, ∞). Also groups with sizes slightly smaller than 100 are concerned as they are involved in significantly less coagulation than in the case without truncation. Summarizing, the cut-off leads to a small overestimate of the probability of occurrence for group sizes in a small neighbourhood of 100 compared to model C. Varying h in the range (0, 0.1) doesn't make a visible difference regarding the kink. For 0.1 ≤ h ≤ 1 the kink becomes much smaller. This indicates that the missing coagulation concerns mainly a neighbourhood of L with radius 0.1. Group sizes outside that range are not visibly affected by not being able to merge into groups of size bigger than 100.
Note the approach of both graphs for x → 0 in Fig. 1b. We can see a high similarity to the predicted small-size behaviour but no real convergence. This divergence close to 0 can be explained by the fact that model C ((2.19)-(2.23)) is continuous and has a singularity at 0 whereas the numerical equilibrium is discrete. Further recall from equations (2.13) and (2.14) that we have chosen the discrete fragmentation rate to be b i,j = 2 i+j+1 whereas the continuous rate is given by b x,y = 2 x+y . Hence, the fragmentation probability is smaller in the discrete setting than in model C. This explains that the generated distribution lies beneath the asymptotic behaviour of model C.
If we choose the computed equilibrium distributions shown in Fig. 1 as initial distributions for the time-dependent scheme described in Section 5.2.2, they actually stay the same over an arbitrary long period of time (taking time step size ∆t ≤ 1.1). This confirms that the computed equilibrium is indeed a proper approximation of the stationary solution of (4.2)-(4.4).
6.1.2. The Euler scheme. Let us now turn to the convergence to the equilibrium in the time evolution scheme. In the following we start with a uniform distribution. We take the time step size ∆t = 1 (which is accurate due to the remark in Section 5.2.2) and work with m 1 = 1. We observe in Fig. 2 that there is actually convergence to the equilibrium. Again, start with the large sizes and take the truncation size L = 100 and the grid size h = 0.01. The stationary distribution reached after time length T = 30 has exactly the same shape as Fig. 1a. As we can see in Fig. 2a, the predicted large-size asymptotics are reached. As in the case of the Newton algorithm, one can also observe the kink at the cut-off due to the reason explained above. For the investigation of the small-size behaviour, we truncate at L = 5 and take h = 0.0005. As in the case of the Newton algorithm for generating the equilibrium, we choose 0.5676 as starting value for the mass to simulate the process for an overall mass of m 1 = 1. For generating the small-size behaviour accurately enough, we have to choose ∆t = 0.5. After T = 6 we get the small-size behaviour displayed in the following Fig. 2b. It seems to equal the predicted asymptotics up to a point very close to 0 where it diverges slightly from the theoretical prediction. This is exactly the same observation as in the Newton scheme. The possible reasons are obviously the same.   Fig. 2a, we take mass m1 = 1, grid size h = 0.01 and the cut-off at L = 100. The plot shows the generated distribution (blue solid line) in a log scale as a function of the group sizes in a linear scale and presents the theoretically found large-size asymptotic behaviour (red dashed line) in a log scale for the sake of comparison. The group sizes are taken in a linear scale in order to illustrate the leading behaviour for large group sizes as a straight line. For Fig. 2b, the equilibrium distribution is approximated by the Euler scheme taking mass m1 = 1, grid size h = 0.0005 and the cut-off at L = 5. The plot shows the generated distribution (blue solid line) in a log scale and presents the theoretically found asymptotic small-size behaviour (red dashed line) in a log scale for the sake of comparison. The group sizes are taken in a log scale as well in order to illustrate the leading behaviour close to zero as a straight line.
we do not have to care about truncation. For the sake of comparison with the continuous model, we will look at f h i as hf (ih) in accordance with equation (3.6). Again, we want to compare the predicted asymptotic behaviour with the asymptotic behaviour displayed by the computed equilibrium distribution: recall from Theorem (3.5) that the equilibrium f h n for mass m h 1 satisfies the large-size asymptotic behaviour given by log 10 f h n ∼ log 10 C − n log 10 z − (3/2) log 10 n as n → ∞, There is no theoretical prediction for the small-size behaviour since the recursive scheme was derived from the discrete model which obviously doesn't have an equilibrium with singularity at zero as opposed to the continuous case. However, we will discuss the possibility of a small-size analysis in Section 6.2. The plots in Fig. 3 indicate that the distribution generated by the algorithm matches very well the predicted asymptotic behaviour for the equilibrium for any h > 0. Again for the sake of comparison with the continuous setting, we choose m h 1 = 1 and compute the terms of the sequence until L = 100. In Fig. 3a, we choose the grid size h = 1 which gives the actual realistic distribution with integer group sizes. The plot compares the predicted asymptotic behaviour given by Eq.(6.3) with the one given by our computed equilibrium. In Fig. 3b, we do the same for h = 0.01. Observe that in both cases the equilibrium is non-negative. Note that  . The equilibrium distribution is approximated by the recursive scheme (Section 5.1). In Fig. 3a, the mass is m h 1 = 1, the grid size h = 1 and the equilibrium sequence is computed till L = 100. It shows the generated distribution (blue solid line) in a log scale and presents the theoretically found asymptotic behaviour (red dashed line) in a log scale (Eq. 6.3) for the sake of comparison. One can observe perfect agreement for large sizes. In Fig. 3b, exactly the same is done for grid size h = 0.01. Again, one can observe that the generated distribution shows the predicted asymptotics. the asymptotics are perfectly matched for both choices of h. As opposed to the truncated discretisation of the continuous model, one cannot observe any kink at the right-hand side of the graph. Obviously, this is the case since we don't need any truncation for the recursive algorithm. Additionally, one can observe that the large-size asymptotics differ for h = 1 and h = 0.01. We are going to investigate this phenomenon more precisely in the next section where we compare the large-size asymptotics of model D ((2.15)-(2.18)) and model C ((2.19)-(2.23)).   equilibrium f h n . The discretization error between f h n and f ∞ (nh), where f ∞ is the continuous equilibrium profile, accumulates with increasing n. Hence, the solutions of the two models separate for larger sizes if h is fixed. We illustrate that in Fig. 6 where we compare the predicted asymptotic behaviour for model D ((2.15)-(2.18)) and model C ((2.19)-(2.23)) at large group sizes x. The plots show the asymptotic behaviour close to x = 200, x = 1000, x = 2000 for fixed h = 0.01. One can see how the difference increases which means that for fixed h the continuous and discrete equilibrium diverge as x → ∞.
6.2. Small-size behaviour for model D. We turn towards the asymptotics of the equilibrium sequence in the case h → 0. Although this seems not significant for the application of the model to animal group sizes, we investigate the question for the sake of a more thorough comparison of model C ((2.19) -(2.23)) and model D ((2.15) -(2.18)). In applications with non-integer cluster sizes this will definitely be relevant.
First, we need to investigate m h 0 for h → 0. As pointed out in [5, Section 15] we can immediately see from Eq. (3.8) that the leading behaviour for h → 0 is given by Obviously, f h 1 is a good indicator of the sought behaviour since it is the first term of the sequence. With the above and using (5.2)-(5.4), one gets (for taking m h 1 = 1 in the end) Observe that for x becoming greater, the difference between both graphs increases significantly.
due to Eq. (3.6). One can see that -except for the factor 1 Γ(4/3) ≈ 1.12 -the discrete case actually has the same leading behaviour as the continuous one: h −2/3 as h → 0 (see (3.5)). (6.6) In Fig. 7a we compare 1 h f h 1 for h ∈ [5 * 10 −5 , 1] with the small-size behaviour of model C. We observe an approximation for decreasing h due to the converging leading behaviour but the preservation of a small gap between the two graphs due to the different constants as seen in (6.6). In Fig. 7b we look at the equilibrium sequence given by the recursive algorithm for h = 5 * 10 −5 , just in the interval [h, 1], and compare it to the behaviour predicted for the continuous case. We note that the two curves show a very close approximation for decreasing group sizes with the very first members of the sequence exhibiting the gap explained above. So the slope close to 0 becomes the same but diverges slightly for the first few members of the sequence. Again, this can be explained by model D ((2.15)-(2.18)) providing a smaller fragmentation rate than model C ((2.19)-(2.23)), in connection with the fact that whereas the continuous equilibrium is defined on (0, ∞) and has a singularity at 0, the discrete equilibrium is a sequence. 6.3. Determination of convergence rates. Degond et al. have proven in [5] that model C ((2.19)-(2.23)) exhibits weak convergence to equilibrium as time goes to ∞. However, there is no finding about convergence almost everywhere. We want to show that the time-dependent solution f (x, t) of Eq. (3.1) converges uniformly to the equilibrium f ∞ if we start with a uniform distribution or also an exponential distribution. We also investigate the convergence rates for different group sizes. For simulating the convergence process, we work with the Euler method in the discretized version D' ((5.5)-(5.6)) of the truncated model C'((4.1)-(4.4)). Denote the discrete approximation of the time-dependent solution by f i (t) ( ∼ f (ih, t)) and the discrete approximation of the equilibrium by f ∞ i . Let's again choose the cut-off at L = 100, grid size h = 0.01, mass m 1 = 1 and time step size ∆t = 1. As initial distribution we first take the uniform distribution (Table 1) as described in Section 5.2.2 and then the exponential distribution (Table 2) as for the Newton method, given by Eq. (5.11). The discretized equilibrium distribution f ∞ i is approximated by conducting the Euler scheme until t = 30. Further, we calculate f i (25), f i (20), f i (15), f i (10) and f i (5) representing f (x, 25), . . . , f (x, 5). We evaluate the distributions at i = 500, 3500, 6500, 9500 (representing x = 5, 35, 65, 95) and consider the relative distance to the equilibrium 10,15,20,25 and i = 500, 3500, 6500, 9500. Table 1 gives an overview of the results for starting with a uniform distribution and Table 2 for starting with an exponential distribution. The tables indicate that the convergence is uniform on a    Table 2. Starting with an exponential distribution the time-dependent solution of model C ((2.19)-(2.23)), f (x, t), is approximated via the Euler scheme for model D' ((5.5)-(5.6)), taking L = 100, h = 0.01, ∆t = 0.5 (smaller than in the previous case due to stabilisation problems for small sizes) and m1 = 1. This approximation, fi(t), is evaluated at t = 5, 10, 15, 20 and the equilibrium distribution is approximated via following the Euler scheme until t = 30. The table shows the relative distances to the equilibrium, for t = 5, 10, 15, 20, 25 and i = 500, 3500, 6500, 9500.
bounded interval since the distance to equilibrium decreases in time monotonically for any i (resp. x). Taking a uniform initial distribution effects in the relative distances being on a much smaller scale for small i than for large i. The impact of the initial distribution vanishes on the long run. The convergence rates seem to approach each other for different group sizes (with the exception of x = 35 in the second table) which can be seen by comparing the quotients of subsequent entries in the different columns.
We are investigating the speed of convergence depending on the sizes and time more thoroughly. Consider the following approach for determining the exponential convergence rate δ x,t where x stands for the group size and t for time: one can express f (x, t) as f (x, t) = f ∞ (x)(1 − e −tδx,t ) + f 0 (x)e −tδx,t .
Thus, if the convergence rate is the same for t 2 and t 1 , it can be expressed as µ(x, t 2 ). (6.8) We have estimated δ x,t2 numerically for x = 5 and x = 95 by calculating the relative distances µ(x, t 1 ), µ(x, t 2 ) as for Table 1 and Table 2. The points of time t 1 , t 2 were taken to be t 1 = 20, . . . , 28 and t 2 = t 1 + 1. We have started with a uniform distribution (Fig. 8a) and with an exponential distribution (Fig. 8b) and observed -as expected -the same limit behaviour for the convergence rates. Note that in (b) Convergence rates for exponential initial Figure 8. Starting with a uniform distribution (Fig. 8a) and with an exponential distribution (Fig. 8b), the time-dependent solution of model C ((2.19)-(2.23)), f (x, t), is approximated via the Euler scheme for model D' ((5.5)-(5.6)), taking L = 100, h = 0.01, m1 = 1 and ∆t = 1 for uniform initial and ∆t = 0.5 for exponential initial (due to stability issues for small sizes). The approximation, fi(t), is evaluated at t = 20, . . . , 29 and the equilibrium distribution is approximated via following the Euler scheme until t = 30. Calculating the relative distances to the equilibrium, µi(t) = |f ∞ i − fi(t)| /f ∞ i , for i = 500 and i = 9500 (representing x = 5 and x = 95), we estimate the exponential convergence rate δx,t 2 (∼ δx,t 1 ) for t1 = 20, . . . , 28 and t2 = t1 + 1 according to Eq. (6.8).
both cases the estimated convergence rates become the same for the small and the large size. The increase in time indicates super-exponential convergence rates. 7. Conclusion. In this work, we have investigated numerically the coagulationfragmentation model for animal group size distributions theoretically discussed by Degond et al. in [5]. The central point of this work was to approximate the equilibria numerically and investigate convergence to equilibrium. We have worked with three different numerical methods: a recursive algorithm -first introduced by Ma et al. in [12] -and a Newton and a time-dependent method -developed in this paper. We have validated our numerical methods by checking the accordance with the predicted asymptotic behaviour and used the time-dependent scheme to show that there is super-exponential convergence to equilibrium in time on finite intervals.
We have seen that the Newton method provides a very fast approximation of the equilibrium after just five iterations. We suggest that the algorithm could be used in more complicated models with coagulation and fragmentation rates depending on the group sizes and/or time. Further, the Newton scheme could be deployed to prove the existence and uniqueness of the equilibrium in such models where the Bernstein method -used in [5] -fails as it solely works for fixed coagulation and fragmentation parameters. Another topic of possible future work is to analyse the indicated super-exponential convergence more precisely and determine the convergence rates analytically.