Coagulation Fragmentation and Sectional Population Balance Numerics

For finite signed Borel measures $\mu ,\nu$ on $E=[a,\mathrm{\infty })$, define their bounded-Lipschitz discrepancy by

$d_{BL}(\mu ,\nu )={\sup }_{\mathrm{\parallel }\phi {\mathrm{\parallel }}_{\mathrm{\infty }}\le 1,\mathrm{Lip}(\phi )\le 1}\mid {\int }_E\phi d(\mu -\nu )\mid ,$

where the supremum is over all real functions on $E$ with supremum norm at most $1$ and Lipschitz constant at most $1$. For a finite measure $\eta$, the notation $\mathrm{\parallel }\eta {\mathrm{\parallel }}_{BL}$ means $d_{BL}(\eta ,0)$. This discrepancy only tests bounded Lipschitz observables; it does not control an unbounded moment without a separate moment estimate.

Let $x_0<x_1<\cdots <x_M=R$ be a finite grid of positive pivot masses with integer $M\ge 1$, and let $\mathrm{\Delta }={\max }_{0\le i<M}(x_{i+1}-x_i)$. For $z\in [x_i,x_{i+1}]$, define the fixed-pivot probability measure

$P_{\mathrm{\Delta }}(z)=\frac{x_{i+1}-z}{x_{i+1}-x_i}{\delta }_{x_i}+\frac{z-x_i}{x_{i+1}-x_i}{\delta }_{x_{i+1}}\mathrm{.}$

At a grid node this rule returns the single Dirac mass at that node, including at both endpoints. For $z>R$, the rule returns an overflow record of mass $z$ and no last-bin particle. For $z<x_0$ it returns an invalid-grid status. The domain convention is therefore part of the method: every retained aggregation product and every daughter size must be at least $x_0$, and a product above $R$ is handled by the overflow rule.

For error analysis, the same barycentric formula is also used on a semi-infinite grid $x_0<x_1<\cdots$ with $x_i\to \mathrm{\infty }$ and ${\sup }_i(x_{i+1}-x_i)\le \mathrm{\Delta }$. On that grid it is defined for every $z\ge x_0$ and never overflows; this is an auxiliary unbounded-pivot dynamics, distinct from the finite-cutoff method above.

Fix a basic mass $a>0$ and an integer $N\ge 2$. Species $i\in \{1,\dots ,N\}$ has mass $m_i=ia$ and number concentration $n_i\ge 0$. A symmetric nonnegative array $K_{ij}=K_{ji}$ specifies binary aggregation. In the ordered-pair convention, the event rate for $(i,j)$ is $K_{ij}{n}_i{n}_j\mathrm{/}2$; summing over ordered pairs counts an unequal pair twice and an equal pair once with its required factor $1\mathrm{/}2$.

If $i+j\le N$, an aggregation event has stoichiometric increment $e_{i+j}-e_i-e_j$, where $e_i$ is the $i$th coordinate vector. If $i+j>N$, it removes the two parents from the retained bins and adds mass $m_i+m_j$ to a scalar overflow reservoir $L$; the reservoir has no return dynamics. Thus this branch is a truncation model and is not a physical state law for particles that later fragment.

For fragmentation, each species $i\ge 2$ has event rate ${\beta }_i{n}_i$, where ${\beta }_i\ge 0$. The monomer does not fragment: set ${\beta }_1=0$ and impose no daughter normalization for $i=1$. Its expected daughter counts are $b_i(k)\ge 0$ for $1\le k<i$, with ${\sum }_{k<i}k{b}_i(k)=i$. One event has increment $-e_i+{\sum }_{k<i}{b}_i(k)e_k$. The mean equations are the finite sums of these event rates times their increments, together with the stated overflow-reservoir equation. When $b_i$ is specified as an actual random daughter law, these equations describe its expected event increment.

Fix $a>0$, put $E=[a,\mathrm{\infty })$, and define the weight $w(x)=1+x^2$. Let ${\mathcal{M}}_w(E)$ be the real vector space of finite signed Borel measures $\sigma$ on $E$ for which

$\mathrm{\parallel }\sigma {\mathrm{\parallel }}_w:={\int }_{E}w(x)\mathrm{\mid }\sigma \mathrm{\mid }(dx)<\mathrm{\infty },$

where $\mathrm{\mid }\sigma \mathrm{\mid }$ is the total-variation measure. The norm is $\mathrm{\parallel }\cdot {\mathrm{\parallel }}_w$. A positive population measure is a nonnegative member $\mu \in {\mathcal{M}}_w(E)$. Its number, mass, and second moment are $M_0(\mu )={\int }_{E}1d\mu$, $M_1(\mu )={\int }_{E}xd\mu$, and $M_2(\mu )={\int }_E{x}^{2}d\mu$. The space and its completeness are established in The weighted total-variation measure space is complete.

Let $r_0>0$ and let $(a_n{)}_{n\ge 0}$ be real coefficients such that for every $0<R<r_0$ both ${\sum }_{n\ge 0}\mathrm{\mid }{a}_n\mathrm{\mid }{R}^n$ and ${\sum }_{n\ge 1}n\mathrm{\mid }{a}_n\mathrm{\mid }{R}^{n-1}$ converge. Then $F(x)={\sum }_{n\ge 0}{a}_n{x}^n$ is differentiable for $\mathrm{\mid }x\mathrm{\mid }<r_0$, with $F^{\mathrm{\prime }}(x)={\sum }_{n\ge 1}n{a}_n{x}^{n-1}\mathrm{.}$

Input a grid $a=x_0<\cdots <x_M=R$, a nonnegative grid population $n_i$, a nonnegative reservoir mass $L$, symmetric nonnegative aggregation rates $K_{ij}$, nonnegative parent fragmentation rates ${\beta }_i$, and daughter-count measures $B_{x_i}$ supported in $[a,x_i]$ with $\int z{B}_{x_i}(dz)=x_i$. For each ordered pair $(i,j)$ use event rate $K_{ij}{n}_i{n}_j\mathrm{/}2$. If $x_i+x_j\le R$, remove one particle from each parent slot and add the fixed-pivot measure $P_{\mathrm{\Delta }}(x_i+x_j)$; if $x_i+x_j>R$, remove both parents and add mass $x_i+x_j$ to $L$, with no retained last-bin particle. For each parent $i$, use fragmentation rate ${\beta }_i{n}_i$, remove its parent particle, and add $\int P_{\mathrm{\Delta }}(z)B_{x_i}(dz)$ to the grid population.

The resulting finite-dimensional autonomous ODE is the sum of these event rates times their stated increments, with $\dot{L}$ equal to the total overflow-mass rate. Invalid grids, negative input populations, or a daughter measure with support outside $[a,x_i]$ return invalid-input. The reservoir is diagnostic truncation mass and has no later fragmentation or coagulation law.

Fix $a>0$ and $E=[a,\mathrm{\infty })$. A continuous-size coagulation-fragmentation model consists of a symmetric nonnegative aggregation kernel $K:E^2\to [0,K_{\ast }]$ and finite constants $L_K,L_{\beta }\ge 0$ such that $\mathrm{\mid }K(x,y)-K(x^{\mathrm{\prime }},y^{\mathrm{\prime }})\mathrm{\mid }\le L_K(\mathrm{\mid }x-x^{\mathrm{\prime }}\mathrm{\mid }+\mathrm{\mid }y-y^{\mathrm{\prime }}\mathrm{\mid }),\mathrm{\mid }\beta (x)-\beta (y)\mathrm{\mid }\le L_{\beta }\mathrm{\mid }x-y\mathrm{\mid },$ where $\beta :E\to [0,{\beta }_{\ast }]$ is the fragmentation rate. It also has a Borel kernel $x\mapsto B_x$ of finite nonnegative daughter-count measures on $E$. Fix a finite constant $D\ge 1$ and require $B_x(E)\le D$, $\mathrm{supp}(B_x)\subseteq [a,x]$, and ${\int }_{E}z{B}_x(dz)=x$. For $x<2a$, take $\beta (x)=0$ and $B_x={\delta }_x$. Assume $d_{BL}(B_x,B_y)\le L_B\mathrm{\mid }x-y\mathrm{\mid }$, where $d_{BL}$ is defined in Bounded-Lipschitz discrepancy for finite measures.

For signed $\mu ,\nu \in {\mathcal{M}}_w(E)$ define the symmetric bilinear coagulation operator and linear fragmentation operator by

$\mathcal{C}(\mu ,\nu )=\frac{1}{2}{\iint }_{E^2}K(x,y)({\delta }_{x+y}-{\delta }_x-{\delta }_y)\mu (dx)\nu (dy),\mathcal{F}\mu ={\int }_E\beta (x)(B_x-{\delta }_x)\mu (dx)\mathrm{.}$

On a time interval $[0,T]$, a weak solution is a path $\mu \in C([0,T];{\mathcal{M}}_w(E))$ of positive measures satisfying the following identity in integral form for each continuous $\phi$ with $\mathrm{\mid }\phi (x)\mathrm{\mid }\le C(1+x^2)$:

$\frac{d}{dt}{\int }_E\phi d{\mu }_t=\frac{1}{2}{\iint }_{E^2}K(x,y)[\phi (x+y)-\phi (x)-\phi (y)]{\mu }_t(dx){\mu }_t(dy)+{\int }_E\beta (x)[{\int }_E\phi (z)B_x(dz)-\phi (x)]{\mu }_t(dx)$

Equivalently, $\dot{\mu }=\mathcal{C}(\mu ,\mu )+\mathcal{F}\mu$ in ${\mathcal{M}}_w(E)$. The first term counts ordered coagulation pairs with factor $1\mathrm{/}2$. The second term removes one parent at each fragmentation event and adds its expected daughter measure; its first-moment normalization is the displayed identity. A binary split $x\mapsto (u,x-u)$ has $B_x={\delta }_u+{\delta }_{x-u}$.

For $\mu ,\nu \in {\mathcal{M}}_w(E)$,

$d_{BL}(\mu ,\nu )\le \mathrm{\parallel }\mu -\nu {\mathrm{\parallel }}_{TV}\le \mathrm{\parallel }\mu -\nu {\mathrm{\parallel }}_w\mathrm{.}$

In the finite discrete model of Finite discrete coagulation and fragmentation with a mass reservoir, an in-range coagulation event changes the retained number by $-1$ and retained mass by $0$. An overflow event changes retained number by $-2$, retained mass by $-(m_i+m_j)$, and reservoir mass by $m_i+m_j$. A fragmentation event of species $i$ changes number by ${\sum }_{k<i}{b}_i(k)-1$ and mass by $0$.

For every $z\in [x_i,x_{i+1}]$, the two weights in Conservative fixed-pivot redistribution on a finite positive grid are nonnegative, sum to $1$, and have weighted pivot mass $z$. Their expected absolute displacement is at most $\mathrm{\Delta }\mathrm{/}2$, and their weighted second moment satisfies

$\int u^2{P}_{\mathrm{\Delta }}(z)(du)-z^2=(z-x_i)(x_{i+1}-z)\le \frac{{\mathrm{\Delta }}^2}{4}\mathrm{.}$

At a grid endpoint the same conclusions hold for the single-pivot endpoint rule, with zero second-moment excess.

Let $T(z)={\sum }_{k\ge 1}{k}^{k-1}{z}^k\mathrm{/}k!$ and $U(z)=T(z)-T(z{)}^2\mathrm{/}2$. These series converge absolutely for $\mathrm{\mid }z\mathrm{\mid }<e^{-1}$. On that disk they satisfy $T=z{e}^T$, $z{T}^{\mathrm{\prime }}=T\mathrm{/}(1-T)$, and $z{U}^{\mathrm{\prime }}=T$. Consequently $U(z)={\sum }_{k\ge 1}{k}^{k-2}{z}^k\mathrm{/}k!$. For real $0\le z<e^{-1}$, $T(z)$ is the unique $u\in [0,1)$ satisfying $u{e}^{-u}=z$.

The normed vector space ${\mathcal{M}}_w(E)$ of Weighted total-variation space for population measures is complete.

Input a valid grid state $(n,L)$ with $n_i\ge 0$, a step $h>0$, and the finite sectional event rates of Finite sectional coagulation-fragmentation balance with reported overflow. The total per-particle loss rate is

${\mathrm{\Lambda }}_i(n)={\sum }_j{K}_{ij}{n}_j+{\beta }_i\mathrm{.}$

First validate the grid, state, rates, and daughter measures against the input conditions of Finite sectional coagulation-fragmentation balance with reported overflow, and require $h>0$; if any condition fails, return invalid-input without updating. For admissible inputs use exact arithmetic to evaluate the finite event sums defining $F$ and ${\mathrm{\Lambda }}_i$. If $h{\mathrm{\Lambda }}_i(n)>1$ for any coordinate, return cfl-failure without updating. Otherwise return success together with $(n^{+},L^{+})=(n,L)+hF(n,L)$. Exact-arithmetic positivity and conservation are stated separately in The sectional forward-Euler step is positive and conserves augmented mass ↗.

Adopt all kernel and daughter-measure hypotheses of Continuous-size coagulation-fragmentation weak model: in particular, $K$ is symmetric, nonnegative, bounded and Lipschitz, $\beta$ is nonnegative, bounded and Lipschitz, and $B_x$ is a nonnegative Borel kernel with $B_x(E)\le D$, support in $[a,x]$, ${\int }_{E}z{B}_x(dz)=x$, and the stated bounded-Lipschitz dependence and convention below $2a$. Let $G=\{x_i:i\ge 0\}$ be a strictly increasing grid with $x_0=a$, $x_i\to \mathrm{\infty }$, and $x_{i+1}-x_i\le \mathrm{\Delta }\le 1$. Use the no-cutoff pivot rule $P_{\mathrm{\Delta }}$ from Conservative fixed-pivot redistribution on a finite positive grid on every cell of this semi-infinite grid. For every signed measure $\mu \in {\mathcal{M}}_w(E)$, write $P_{\mathrm{\Delta }}\mu =\int P_{\mathrm{\Delta }}(x)\mu (dx)$ for the linear projection of a measure on $E=[a,\mathrm{\infty })$ onto $G$.

For a signed grid measure $\nu \in {\mathcal{M}}_w(E)$, define $F_{\mathrm{\Delta }}(\nu )=\frac{1}{2}\iint K(x,y)(P_{\mathrm{\Delta }}(x+y)-{\delta }_x-{\delta }_y)\nu (dx)\nu (dy)+\int \beta (x)(\int P_{\mathrm{\Delta }}(z)B_x(dz)-{\delta }_x)\nu (dx)\mathrm{.}$ All products are projected, with no overflow; positive states are the population dynamics. The finite-cutoff algorithm instead sends products above its last pivot to the reported reservoir. By the barycentric identities, each projection has total number one, reproduces its first moment, and adds at most ${\mathrm{\Delta }}^2\mathrm{/}4$ to the second moment of each projected size.

For the coagulation bilinear operator $\mathcal{C}$ of Continuous-size coagulation-fragmentation weak model, every $\mu ,\nu \in {\mathcal{M}}_w(E)$ satisfies

$\mathrm{\parallel }\mathcal{C}(\mu ,\nu ){\mathrm{\parallel }}_w\le \frac{3}{2}{K}_{\ast }\mathrm{\parallel }\mu {\mathrm{\parallel }}_w\mathrm{\parallel }\nu {\mathrm{\parallel }}_w\mathrm{.}$

Assume the daughter-count bound $B_x(E)\le D$ of Continuous-size coagulation-fragmentation weak model with $D\ge 1$. The fragmentation operator there satisfies

$\mathrm{\parallel }\mathcal{F}\mu {\mathrm{\parallel }}_w\le {\beta }_{\ast }(D+1)\mathrm{\parallel }\mu {\mathrm{\parallel }}_w(\mu \in {\mathcal{M}}_w(E))\mathrm{.}$

For fixed grid, aggregation rates, fragmentation rates, and daughter measures, the augmented right-hand side $F$ in Finite sectional coagulation-fragmentation balance with reported overflow is locally Lipschitz on the full finite-dimensional state space.

For the ODE of Finite sectional coagulation-fragmentation balance with reported overflow, every solution starting from $n_i(0)\ge 0$ and $L(0)\ge 0$ remains nonnegative for as long as it exists. The augmented mass

$m_{\mathrm{t}\mathrm{o}\mathrm{t}}(t)={\sum }_{i=0}^M{x}_i{n}_i(t)+L(t)$

is constant. The retained-bin number and retained mass need not be constant when overflow occurs.

Define concentrations $c_k(t)\ge 0$ by the discrete product-kernel equations ${\dot{c}}_k=\frac{1}{2}{\sum }_{i+j=k}ij{c}_i{c}_j-k{c}_k{\sum }_{j\ge 1}j{c}_j,k\ge 1,$ with monomer initial data $c_1(0)=1$ and $c_k(0)=0$ for $k\ge 2$. For $r=1,2$, define $M_r(t)={\sum }_{k\ge 1}{k}^r{c}_k(t)$. The pre-gelation solution is $c_k(t)=\frac{k^{k-2}}{k!}{t}^{k-1}{e}^{-kt},0\le t<1,$ where the value at $t=0$ is understood by continuity. Its first and second moments satisfy $M_1(t)=1$ and $M_2(t)=1\mathrm{/}(1-t)$. Thus the second moment becomes unbounded at the finite time $t=1$ although the pre-gelation first moment is constant.