Numerical approximation of a coagulation-fragmentation model for animal group size statistics

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 coagulation-fragmentation 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

$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 time-dependent 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

2.1.   General form of the equations

The continuous version of a coagulation-fragmentation equation, called also Smoluchowski equation, describes the evolution of the number density $f(x,t)$ of continuous sizes $x \geq 0$ at time $t$. In weak form it reads, for all test functions $\varphi$ 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 $\varphi(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

The case where the cluster sizes form a discrete set can be described analogously. So consider a system of clusters with discrete sizes $i \in \mathbb{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 $\varphi_i$ in a class to be specified later

The equation can also be written similarly to (2.2) as

If one takes $\varphi_k = k$, it can be seen immediately that mass is conserved:

Let $Q_{Ci}$ and $Q_{Fi}$ denote the coagulation resp. fragmentation operator for cluster size $i$. Then the strong form can be written as

2.2.   The equations based on Niwa's model

According to Niwa's model, we assume $s$ different zones of space on which $\Phi$ 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} = 2 q$ 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 \geq 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), \dots, (k-1,1)$ with equal probability. As the actually distinct pairs are counted twice in such an enumeration, one gets for all $1 \leq i,j < k$ with $i + j =k$: $b_{i,j} = \frac{p}{(i+j-1)/2} = \frac{2p}{i+j-1}$. 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 $\overline{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):

● Model C (Continuous):

The fragmentation of a group of size $k$ in Model D can now be understood as breaking into the pairs $(0,i), \dots, (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 $\varphi_i$,

The test space is chosen to be $l^{\infty}$ 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 $\varphi \in C([0, \infty))$, as

or as

The test space is chosen to be $C([0, \infty))$ 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 \in \mathbb{N}$ and $f: x \in \mathbb{R}_+ \mapsto f(x) \in \mathbb{R}_{+}$. The $k$th moment $m_k(f)$ is given by

For initial condition $f_0$ with $m_1(f_0) < \infty$, we know from (2.3) that $m_1(f(t)) = m_1(f_0) := m_1$. There is a scaling invariance for model C ((2.19) -(2.23)):

Proposition 3.1. Let $f_0: x \in \mathbb{R}_{+} \mapsto f_0(x) \in \mathbb{R}_{+}$ be an initial condition for (2.21) with $m_1(f_0) =: m_1 < \infty$ and let $f_{\overline{p},\overline{q}}(x,t)$ be the solution of (2.21) with parameters $\overline{p}$ and $\overline{q}$. Then, one has

where $f_{1,1}$ is associated with the initial condition $\tilde{f}_0$ such that

Due to this proposition, we can assume $\overline{p}=1$, $\overline{q}=1$ and $m_1=1$. The problem in strong form becomes

In weak form it reads as

Definition 3.2. A function $f: x \in (0, \infty) \to \mathbb{R}$ is said to be completely monotone if it is $C^{\infty}$ and such that

The main theorem can be stated as follows:

Theorem 3.3. There is a unique equilibrium distribution function $f_{\infty}$ of (3.1) and (3.3) satisfying (3.2). It can be written as

where $\gamma$ is a completely monotone function and has the following asymptotic behaviour:

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 $k$th moment of a sequence $f = (f_i)_{i \in \mathbb{N}}$ is given by

Let us further introduce the sets

One can establish the well-posedness of the initial value problem (for a proof see [5,Theorem 12.1]):

Theorem 3.4. Let $k \geq 0$ and $f_{in} =(f_{in,i})_{i \in \mathbb{N}}$ be given in $\ell_{1,k}$. Then there exists a unique global-in-time strong solution $f \in C^1([0,\infty), l_{1,k})$ for system (2.16)-(2.18) with $f(0) = f_{in}$. If $k \geq 1$, then $m_1(f(t)) = m_1(f_{in})$ for all $t\geq 0$.

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

for the number of clusters with sizes in the interval $I_i^h$.

For a smooth test function $\varphi(x)$, we can write $\varphi_i = \varphi(ih)$ and require that $f^h(t) =(f_i^h(t))_{i \in \mathbb{N}}$ solves Model D ((2.15)-(2.18)) as a discretization of Model C ((2.19)-(2.23)):

Letting $h \to 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. [5,Sections 11 and 15]):

Theorem 3.5. For any $m_1^h \in [0, \infty)$, there is a unique equilibrium solution $f^h =(f_i^h)_{i \in \mathbb{N}}$ of model D ((2.15)-(2.18)). The solution has the form

where $\gamma$ 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 $\Delta$ is given by $(\Delta f^h)_n = f_{n+1}^h - f_n^h$ and $(\Delta^k f^h)_n = \left( \Delta\left( \Delta^{k-1} f^h \right) \right)$.

Let $F_t^h$ denote the discrete measure on the grid $\{ih: \ i= 1, \dots \}$ formed from the solution $f^h(t)$ of model D ((2.15)-(2.18)):

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_t^h \to F_t$ narrowly as $h \to 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 $\mathbb{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 $\varphi$ be a test function. The truncation is chosen as follows:

Definition 4.1. A time-dependent size distribution $f(t,x)$ in Model C' is characterised as a solution of the weak problem

for all $t > 0$ and test functions $\varphi$.

Note that, indeed, by chosing $\varphi(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:

Proposition 4.3. Let $f_0 \in L^1([0,L])$ and $R > 0$. Then there is an $\alpha >0$ such that the initial value problem corresponding with (4.2)

has a unique solution on $[0, \alpha]$ with values in $\overline{B}(f_0, R) \subset L^1([0,L])$.

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

with

$T$ is a linear operator whereas $q$ is a bilinear form with $Q_{C_T}(f) = q(f,f)$.

Starting with an appropriate $f_0$, we want to find a recursive scheme giving a convergent sequence $(f^n)_{n \in \mathbb{N}}$ with limit $f_{\infty}$, the equilibrium. Observe the following: If $f^{n+1}$ was an equilibrium, we'd have

with

when $\left|f^{n+1} - f^n\right|$ is small. Hence, the following Newton scheme rests upon neglecting this quadratic term and defines a sequence $(f^n)_{n \in \mathbb{N}}$ by iteratively solving the following linear problem:

Introducing $\delta f = f^{n+1} - f^n$, by adding $- T f^n$ and $2 Q_{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_{\phi}$ can be written as

where

In the following, ${\rm R}{(W_{\phi})}$ denotes the range of $W_{\phi}$ and ${\rm N}{(W_{\phi})}$ its null space.

For $f \in L^1([0,L])$, $g \in L^{\infty}([0,L])$ we define

and for $V \subset L^{\infty}([0,L])$ we define

Lemma 4.4. Let $\phi \in L^1([0,L])$. Then $W_{\phi}$, as given in equations (4.10), (4.11), is a bounded, linear operator from $L^1([0,L])$ to $L^1([0,L])$.

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 $\phi$ in order to build on equation (4.9)):

Lemma 4.5. For any $f^n \in L^1 ([0,L])$ it holds that ${\rm R}{(W_{f^n})} \subset {\rm {span}}{\{x\}}^{\perp}$.

Proof. For any test function $\varphi$ we have

So if we set $\varphi(x) = x$, we get

By adding and subtracting $Tf^{n}(x)$ and $Q_{C_T}(f^{n})(x)$, one can see that

Since this is true for any $\delta f$ and the first two summands can be cancelled due to (4.12), for any $\lambda > 0$ it holds that

Extracting the $\lambda$ and dividing by $\lambda$ leaves the factor $\lambda$ on the right hand side of the equation. Due to arbitrariness of $\lambda$, 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]):

Conjecture 4.6. ${\rm R}{(W_{\phi})} = {\rm {span}}{\{x\}}^{\perp}$, $\dim{ {\rm N}{(W_{\phi})}} = 1$ and ${\rm N}{(W_{\phi})} \cap {\rm {span}}{\{x\}}^{\perp} = \{0\}$.

Proving this conjecture allows to single out the solution of $W_{\phi}f = g$ by imposing $\int x f 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_{\infty}$ of model C ((2.19) -(2.23)). Recall from Theorem 3.3 that

where $\gamma$ is a completely monotone function with

The analysis in [5] suggests that for $L \sim 10$ the function $\gamma$ 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.

5.1.   A recursive algorithm for model D equilibria

The equilibrium sequence in model D ((2.15) -(2.18)), $(f_i^h)_{ i \in \mathbb{N}}$, can be computed recursively for any $h > 0$ (see [5,Section 4.2.3] and [12,Eq. (13)-(15)]).

For a test function $\varphi$ with $\varphi_i = \varphi(ih)$, the equilibrium profile satisfies

Define now

Taking $\varphi_j \equiv 1$ yields

Further, with taking $\varphi_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 $\varphi_j \equiv 1$ in the weak form has lead to the relation (5.1) which is used to obtain the following recursive algorithm:

Choose $m_0^h \in (0,1)$ and set

Then for $i = 1,2,3, \dots$:

5.2.   Model D': The discretized form of model C'

5.2.1.   Setting of the model

We consider solutions $f(x,t)$ of the truncated model C' and write $f_i(t) = f(ih,t)$ for the discretised function. Let $L >0$ be the truncation size, $h$ the grid size and $N= L/h$. Write $\phi_i = \phi(ih)$ for a test function $\phi$.

Definition 5.1. The weak form of model D', the discretisation of model C', is given by the following evolution equation for the discrete size distribution $f_i(t)$:

for all test sequences $\phi_i$.

Observe that mass is preserved over time according to this equation.

Remark 5.2. Note that the link between model C' ((4.1)-(4.4)) and model D' resembles the link between model C ((2.19)-(2.23)) and D ((2.15)-(2.18)) as discussed in Section 3.2. However, note that equation (3.6) defines $f_i^h(t)$ to be interpreted as $hF_t(ih)$, if $F_t(x)$ is a solution of model C.

Proposition 5.3. The strong form of model D' is given by

for $1\leq i \leq N$.

Proof. Obvious calculation.

5.2.2.   Time discretization of the time-evolution scheme

The explicit Euler scheme in time is applied with time step size $\Delta t$. Let $t_k = k \Delta t$. The sequence $\{f^k\}_{k \in \mathbb{N}}$ denotes an approximation of $\{f(t_k)\}_{k \in \mathbb{N}}$ and is defined by the following recursive scheme:

where for any point $ih$ with $1 \leq i \leq N$, $\left\{ \left(\frac{df_i}{d\tau}\right) ^k \right\}_{i=1, \dots, N}$ is given by

A linear stability analysis of the scheme would mean to study the eigenvalues of the following operator $\mathcal{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 $\Delta 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 $\Delta t = 1.1$.

5.2.3.   Equilibrium in model D': The Newton method

The stationary equation in the discretized setup of model C' ((4.1)-(4.4)) reads, for $1\leq i \leq N$,

Analogously to Eq. (4.5) involving the operators $T$ and $q$, the discretized problem can be written as

where for $1 \leq i \leq N$

$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:

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, \dots, N)$. The range of the operator is restricted to ${\rm {span}} \{x\}^{\perp}$, i.e. to $N-1$ dimensions, and, hence, consider the above equation just for $1 \leq i \leq 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 (3.4) the unique equilibrium $f_{\infty}$ for mass $m_1 = 1$ satisfies

Due to equation (3.5), the large-size asymptotic behaviour of $f_{\infty}$ is given by

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 \to \infty$ 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, $-\frac{4}{27}x \log_{10}{e}$, in a linear shape. Second, we focus on the small-size behaviour $x \to 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_{\infty}(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 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, \infty)$. 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 \leq h \leq 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 \to 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.3) and (2.4) that we have chosen the discrete fragmentation rate to be $b_{i,j} = \frac{2}{i+j+1}$ whereas the continuous rate is given by $b_{x,y} = \frac{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 $\Delta t \leq 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 $\Delta 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 $\Delta 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.

6.1.3.   The recursive computation of the equilibrium sequence

Now we turn to checking the accuracy of the recursive scheme introduced in Section 5.1. In the following we will choose $m_1^h$ and then $m_0^h$ such that equation (3.8) is satisfied. Using the recursive algorithm determined by equations (5.2)-(5.4), one can compute the equilibrium $(f_i^h)_{i \in \mathbb{N}}$ up to an arbitrarily large integer. As opposed to model C' ((4.1)-(4.4)), we do not have to care about truncation. For the sake of comparison with the continuous model, we will look at $f_i^h$ as $h f(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_n^h$ for mass $m_1^h$ satisfies the large-size asymptotic behaviour given by

where

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_1^h =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 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)).

6.1.4.   Link between discrete and continuous model

i) Convergence for fixed interval length L. For $m_1 = m_1^h =1$ the continuous and discrete models can be compared as follows: according to Eq. (6.2) the leading term in the asymptotics of the continuous equilibrium $f_{\infty}$ is given by $e^{-(4/27)x}$ as $x \to \infty$. Set $x=hn$. Then, due to Eq. (6.3), the leading term in the asymptotics of the discrete equilibrium $f_n^h$ is given by $[(1+\frac{4h}{27})^{-1/h}]^x$ as $n (=x/h) \to \infty$. Since we have

the leading term of the discrete equilibrium converges to the leading term of the continuous equilibrium as $h \to 0$. Deploying the Newton method and the recursive scheme, we verify numerically if the same holds true for the truncated models uniformly on a fixed interval $[0,L]$. Indeed, we can observe that for $h$ small enough and a fixed truncation size $L$, the discretized equilibrium for model C' ((4.1)-(4.4)), i.e. for model D' ((5.5)-(5.6)), as approximated by the Newton method and the equilibrium for model D ((2.15)-(2.18)) generated by the recursive algorithm are very close. We have chosen $L=100$, $h=0.01$ and $m_1 = m_1^h =1$. The equilibrium computed by the Newton scheme -the solid blue line in Fig. 4 -and the equilibrium computed by the recursive scheme -the dotted red line in Fig. 4 -are the same up to a maximal absolute error of magnitude $10^{-6}$. This can be seen as an additional validation of the Newton method.

We have verified numerically uniform convergence of model C' ((4.1)-(4.4)) and model D' ((5.5)-(5.6)) in their large-size behaviour on finite intervals as $h \to 0$. This reflects the uniform convergence of model C ((2.19)-(2.23)) and model D ((2.15)-(2.18)) on finite intervals as indicated by Eq. (6.4). We illustrate this by fixing $L=200$ and comparing the asymptotics of model C and model D for $h$ becoming smaller. Fig. 5 shows the asymptotic large-size behaviour of the discrete equilibria $(f_i^h)_{ i \in \mathbb{N}}$ generated by the recursive algorithm in Section 5.1 and the analytically predicted continuous one (Eq. (6.2)). We consider the grid sizes $h=1$, $h=0.1$ and $h=0.01$ and observe the expected convergence of both models. The equilibrium in the genuine discrete case of model D ((2.15)-(2.18)), i.e. $h=1$, differs from the stationary solution of model C ((2.19)-(2.23)) in its large-size behaviour. This difference becomes smaller for $h=0.1$ and even much smaller, invisible in the shown scale, for $h=0.01$.

ii) Divergence on increasing intervals. If we fix $h$ and increase the investigated intervals of group sizes, the large-size behaviour of the discrete and continuous model diverge. This is due to the dependence on $h$ in the asymptotics of the discrete equilibrium $f_n^h$. The discretization error between $f_n^h$ and $f_{\infty}(nh)$, where $f_{\infty}$ 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 \to \infty$.

6.2.   Small-size behaviour for model D

We turn towards the asymptotics of the equilibrium sequence in the case $h \to 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_0^h$ for $h \to 0$. As pointed out in [5,Section 15] we can immediately see from Eq. (3.8) that the leading behaviour for $h \to 0$ is given by

Obviously, $f_1^h$ 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_1^h =1$ in the end)

Let's compare this behaviour with the small-size asymptotics of the stationary solution of model C ((2.19)-(2.23)), denoted by $f$. We need to collate $f(h)$ with $\frac{1}{h} f_1^h$ due to Eq. (3.6). One can see that -except for the factor $\frac{1}{\Gamma(4/3)} \approx 1.12$ -the discrete case actually has the same leading behaviour as the continuous one:

In Fig. 7a we compare $\frac{1}{h} f_1^h$ for $h \in [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,\infty)$ 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 $\infty$. 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_{\infty}$ 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)$ ($\sim f(ih,t)$) and the discrete approximation of the equilibrium by $f_i^{\infty}$.

Let's again choose the cut-off at $L=100$, grid size $h=0.01$, mass $m_1 = 1$ and time step size $\Delta 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^{\infty}$ 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), \dots, 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 $\frac{\left|f_i^{\infty} - f_i(t)\right|}{f_i^{\infty} }$ for $t = 5,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 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 $\delta_{x,t}$ where $x$ stands for the group size and $t$ for time: one can express $f(x,t)$ as

Substracting and dividing both sides by $f_{\infty}(x)$ and taking absolute values gives

Hence, for two different points of time $t_1$ and $t_2$, one gets

Thus, if the convergence rate is the same for $t_2$ and $t_1$, it can be expressed as

We have estimated $\delta_{x, t_2}$ numerically for $x = 5$ and $x = 95$ by calculating the relative distances $\mu(x,t_1), \mu(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, \dots, 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 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 coagulation-fragmentation 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.

Acknowledgments

The authors would like to thank J-G. Liu and R. Pego for enlightening discussions. This work has been supported by the Engineering and Physical Sciences Research Council (EPSRC) under grant ref. EP/M006883/1, and by the National Science Foundation (NSF) under grant RNMS11-07444 (KI-Net). P. D. is on leave from CNRS, Institut de Mathématiques, Toulouse, France. He acknowledges support from the Royal Society and the Wolfson foundation through a Royal Society Wolfson Research Merit Award. M.E. has been supported by the German National Academic Foundation during the first part of this work and is now supported by a Roth Scholarship from the Department of Mathematics at Imperial College London.