1. Introduction

Many natural phenomena can be formulated as the differential law of the development (evolution) in time of a physical system. The resulting differential equation combined with boundary conditions affecting the system are called evolution equations. Evolution equations are a popular framework for studying the dynamics of biological populations. For example, they have proven useful in understanding the dynamics of biological invasions [51], bacterial flocculation in activated sludge tanks [7], and the spread of parasites and diseases [30]. Since many biological populations converge to a time-independent state, many researchers have used theoretical tools to investigate long-term behavior of these models. Analytical and fixed point methods have been used to show the existence of stationary solutions to size-structured population models [44,29] and semigroup theoretic methods have been used to investigate the stability of these stationary solutions [28,43,4]. For many models in the literature, the principle of linearized stability [55,34] can be used to show that the spectral properties of the infinitesimal generator (IG) of the linearized semigroup determines the stability or instability of a stationary solution. Moreover, using compactness arguments, spectral properties of the infinitesimal generator can be determined from the point spectrum of the IG, which in turn can be written as the roots of a characteristic equation.

Despite this theoretical development, the derived existence and stability conditions are oftentimes rather complex, and accordingly the biological interpretation of these conditions can be challenging. To overcome this difficulty several numerical methods for stability analysis of structured population models have been developed [12,27,37,21]. For instance, Diekmann et al. [37,22,21] present a numerical method for equilibrium and stability analysis of physiologically structured population models (PSPM) or life history models, where individuals are characterized by a (finite) set of physiological characteristics (traits such as age, size, sex, energy reserves). In this method a PSPM is first written as a system of integral equations coupled with each other via interaction (or feedback) variables. Consequently, equilibria and stability boundary of the resulting integral equations are numerically approximated using curve tracing methods. Later, de Roos [19] presented the modification of the curve tracing approach to compute the demographic characteristics (such as population growth rate, reproductive value, etc) of a linear PSPM. For a fast and accurate software for theoretical analysis of PSPMs we refer interested reader to a software package by de Roos [20]. An alternative method for stability analysis of physiologically structured population models, developed by Breda and coworkers [11,13,10], uses a pseudospectral approach to compute eigenvalues of a discretized infinitesimal generator. This method (also known as infinitesimal generator (IG) approach) has been employed to produce bifurcation diagrams and stability regions of many different linear evolution equations arising in structured population modeling [13,14,15]. Unfortunately, not all structured population models fit into the framework of PSPMs and thus there is a need for a more general numerical framework for asymptotic analysis of structured population models.

In this paper we develop a numerical framework to guide theoretical analysis of structured population models. We demonstrate that our methodology can be used for numerical computation and stability analysis of positive stationary solutions of both linear and nonlinear size-structured population models. Moreover, we illustrate the utility of our framework to produce approximate existence and stability regions for steady states of size-structured population models. We also provide an open source Python program used for the numerical simulations in our GitHub repository [42]. Although, the examples presented in this paper are size-structured population models, in Section 2, we show that the framework is applicable to more general evolution equations.

The main idea behind the numerical framework is first to write a structured population model in the form of an evolution equation and then use the well-known Trotter-Kato Theorem [53,35] to approximate the infinitesimal generator of the evolution equation on a finite dimensional space. This in turn allows one to approximate solutions (or spectrum) of the evolution equation with the solution (or spectrum) of system of differential equations. Consequently, we approximate the stationary solutions of an actual model with stationary solutions of the approximate infinitesimal generator on a finite dimensional space. Approximate local stability of the approximate steady states are then computed from the spectrum of the Jacobian of ODE system evaluated at this steady states. Our method is similar to the IG approach in [13,14,15], in a sense that we also approximate infinitesimal generator and analyze the spectrum of the resulting operator to produce existence and stability regions. However, in contrast to IG approach, our framework also computes actual steady states and is easily applied to nonlinear evolution equations arising in structured population dynamics.

The rest of the paper is structured as follows. We describe the theoretical development of our framework for general evolution equations in Section 2.1. Note that readers with more biological background can skip Section 2.1 and directly jump into the application of the framework in Section 2.3. In Section 2.3, we illustrate the convergence of the approximation method by applying it to linear Sinko-Streifer model, for which the exact form of the stationary solutions is known. To further illustrate the utility of our approach, in Section 3, we apply our framework to a nonlinear size-structured population model (also known as population balance equations in the engineering literature) described in [5,9]. Moreover, in Section 2.2, we show that approximate local stability conditions for a stationary solution can be derived from the spectral properties of the approximate infinitesimal generator. Finally, we conclude with some remarks and a summary of our work in Section 5.

---

2. Numerical framework

In this section, we demonstrate our numerical methodology for general evolution equations. We first present the numerical scheme used to approximate the infinitesimal generator of a semigroup. Subsequently, in Section 2.3, we illustrate the convergence of our approach by applying it to linear Sinko-Streifer equations, for which exact stationary solutions are known.

---

2.1. Infinitesimal generator approximation

Given a Banach space $\mathcal{X}$, consider an abstract evolution equation,

$u_{t}=\mathcal{F}(u), u(0,\bullet)=u_{0}\in\mathcal{X}\,,\label{eq:abstract evolution equation}$

(1)

where $\mathcal{F}\,:\,\mathcal{D}(\mathcal{F})\subseteq\mathcal{X}\to\mathcal{X}$ is some operator defined on its domain $\mathcal{D}(\mathcal{F})$ and $u_{0}$ is an initial condition at time $t=0$. Note that any boundary condition belonging to a given partial differential equation can be included in the domain $\mathcal{D}(\mathcal{F})$. The solution to (1) can be related to the initial state $u_{0}$ by some transition operator $T(t)$ so that

$u(t,\,x)=T(t)u_{0}(x)\,.$

The transition operator $T(t)$ is said to be a strongly continuous semigroup (or simply $C_{0}$ -semigroup) if satisfies the following three conditions:

1. $T(s)T(t)=T(s+t)$ for all $s,\,t\ge0$

2. $T(0)=I$, the identity operator on $\mathcal{X}$

3. For each fixed $u_{0}\in\mathcal{X}$,

$\lim\limits_{t\to0^{+}}\left\Vert T(t)u_{0}-u_{0}\right\Vert =0\,.$

Moreover, showing that the operator $\mathcal{F}$ generates a $C_{0}$-semigroup is equivalent to establishing well-posedness of the abstract evolution equation given in (1). In general, finding the explicit form of the transition operator is challenging. Hence, approximation methods, e.g. Yosida approximants, are used to study a more complicated evolution equation and the semigroups they generate. One of the famous approximation theorems is due to Trotter [53] and Kato [35] (see [36] for the classical literature on the approximation of generators of semigroups). The goal is to construct approximating generators $\mathcal{F}_{n}$ on the approximate spaces $\mathcal{X}_{n}$ such that $C_{0}$-semigroups $T_{n}(\cdot)$ generated by $\mathcal{F}_{n}$ approximate the $C_{0}$-semigroup $T(t)$ generated by $\mathcal{F}$. Although there are multiple ways to approximate the infinitesimal generator $\mathcal{F}$, for our purposes we use the approximation scheme based on Galerkin-type projection of the state space $\mathcal{X}$ [6,32,2]. For the convenience of readers, we will summarize the approximation scheme here.

Let $\mathcal{X}_{n}$, $n=1,2,\dots$ be a sequence of subspaces of $\mathcal{X}$ with $\dim(\mathcal{X}_{n})=n$ and define projections $\pi_{n}\,:\,\mathcal{X}\to\mathcal{X}_{n}$ and canonical injections $\iota_{n}\,:\,\mathcal{X}_{n}\to\mathcal{X}$. Assume that the projections $\pi_{n}$ are bounded, i.e., there exists $\tilde{M}>0$ such that

$\left\Vert \pi_{n}\right\Vert \le\tilde{M}$

(A1)

for all $n=1,\dots,n$. Moreover, assume that

$\mathop {\lim }\limits_{n \to \infty } {\pi _n}v = v$

(A2)

for all $v\in\mathcal{X}$. Consequently, for each subspace $\mathcal{X}_{n}$ we choose basis elements $\left\{ \beta_{i}^{n}\right\} _{i=1}^{n}$ such that each element $v$ of subspace $\mathcal{X}_{n}$ can be written in the form $v=\sum_{i=1}^{n}\alpha_{i}\beta_{i}^{n}$. Moreover, for each subspace $\mathcal{X}_{n}$ we define the bijective mappings $p_{n}\,:\,\mathcal{X}_{n}\to\mathbb{R}^{n}$ by

$p_{n}v=(\alpha_{1},\cdots,\alpha_{n})^{T}$

and the norm on $\mathbb{R}^{n}$ by

$\left\Vert v\right\Vert _{\mathbb{R}^{n}}=\left\Vert p_{n}^{-1}v\right\Vert _{\mathcal{X}}\,.$

Consequently, we define bounded linear operators $P_{n}\,:\,\mathcal{X}\to\mathbb{R}^{n}$ and $E_{n}\,:\,\mathbb{R}^{n}\to\mathcal{X}$ by

$P_{n}v=p_{n}\pi_{n}v, v\in\mathcal{X}\label{eq: space to reals}$

(2)

and

$E_{n}z=\iota_{n}p_{n}^{-1}z, z\in\mathbb{R}^{n}\,,\label{eq:reals to space}$

(3)

respectively. Finally, we define approximate operators $\mathcal{F}_{n}$ on $\mathbb{R}^{n}$ by

$\mathcal{F}_{n}(z)=P_{n}\mathcal{F}\left(E_{n}z\right)\label{eq:approximate LHS}$

(4)

for all $z\in\mathbb{R}^{n}$.

Accordingly, the Trotter-Kato Theorem states that the semigroup generated by the discrete operator $\mathcal{F}_{n}$ converges to the semigroup generated by the infinitesimal generator $\mathcal{F}$. For the convenience of the reader, we state the theorem here and refer readers to [32] for a proof.

Theorem. (Trotter-Kato) Assume that ($\mathbf{A}1$) and ($\mathbf{A}2$) are satisfied. Let $\left(T(t)\right)_{t\ge0}$ and $\left(T_{n}(t)\right)_{t\ge0}$, $n\in\mathbb{N}$, be strongly continuous semigroups on $\mathcal{X}$ and $\mathbb{R}^{n}$ with generators $\mathcal{F}$ and $\mathcal{F}_{n}$, respectively. Furthermore, assume that they satisfy the estimate

$\left\Vert T(t)\right\Vert _{\mathcal{X}},\,\left\Vert T_{n}(t)\right\Vert _{\mathbb{R}^{n}}\le Me^{wt}\text{for all }t\ge0,\,n\in\mathbb{N}\,,$

for some constants $M\ge1,\,w\in\mathbb{R}$. Then the following are equivalent:

1. There exists a $\lambda_{0}\in\rho(\mathcal{F})\cap\bigcap\limits _{i=1}^{n}\rho(\mathcal{F}_{i})$ such that for all $v\in\mathcal{X}$

$\left\Vert E_{n}\left(\lambda_{0}I_{n}-\mathcal{F}_{n}\right)^{-1}P_{n}v-\left(\lambda_{0}I-\mathcal{F}\right)^{-1}v\right\Vert _{\mathcal{X}}\to\;0\;\text{as }n\to\infty\,.$

2. For all $v\in\mathcal{X}$ and $t\ge0$,

$\left\Vert E_{n}T_{n}(t)P_{n}v-T(t)v\right\Vert _{\mathcal{X}}\to0$

as $n\to\infty$, uniformly on compact $t$ intervals.

In general, one establishes the first statement for a Trotter-Kato approximation and then uses the second statement to approximate an abstract evolution equation on a finite dimensional space. In their paper, Ito and Kappel [32] present the standard ways to establish the first statement of the theorem (see also [6,2,1]). Therefore, here we assume that for a particular problem the first statement in the theorem has already been established and thus the evolution equation in (1) can be approximated by the following system of ODEs,

$u_{n}'(t)=\mathcal{F}_{n}\left(u_{n}(t)\right), u_{n}(0)=P_{n}u(0,\,\bullet)\,.\label{eq:approximate IVP}$

(5)

Consequently, the solution of the IVP is mapped onto the infinite dimensional Banach space $\mathcal{X}$ and one has the following convergence

$\mathop {\lim }\limits_{n \to \infty } {\left\| {{E_n}{u_n}(t) - u} \right\|_{\cal X}} = 0$

(6)

for $t$ in compact intervals.

In general, finding explicit stationary solutions of abstract evolution equations is a challenging task. Conversely, many efficient root finding methods have been developed for finding steady states of a system of ODEs. For large-scale nonlinear systems, many efficient methods have been developed as well. Hence, we propose a numerical framework that utilizes those efficient root finding methods to approximate steady state solutions of general evolution equations. The idea is to use an efficient and accurate root finding method to approximate a stationary solution of the evolution equation (1) with the stationary solutions of the IVP in (5). Thus, as a consequence of the Trotter-Kato Theorem, the steady states of (5) converge to the steady states of (1) as $n\to\infty$.

---

2.2. Stability of stationary solutions

Studying the asymptotic behavior of solutions is a fundamental tool for exploring the evolution equations which arise in the mathematical modeling of real world phenomena. To this end, many mathematical methods have been developed to describe long-term behavior of evolution equations. For instance, for long-time behavior of linear evolution equations, linear semigroup theoretic methods can be used to formulate physically interpretable conditions. Furthermore, for nonlinear evolution equations, the principle of linearized stability can be used to relate the spectrum of the linearized infinitesimal generator to the local stability or instability of the stationary solution. Nevertheless, investigating the spectrum of the linearized infinitesimal generator is cumbersome and requires advanced functional analysis techniques. In contrast to general evolution equations, the asymptotic behavior of ordinary differential equations are determined by the eigenvalues of the Jacobian and well-studied. Hence, in this section we demonstrate that the approximation scheme, presented in Section 2.1, can also be used to give some insights about the stability of stationary solutions of the general evolution equations.

Stability results discussed in this section are not in a traditional Lyapunov sense. In particular, since stationary solutions discussed in this paper are only approximations to actual stationary solutions, the stability results only hold for finite time intervals. Therefore, we refer to this kind of stability as approximate local stability of stationary solutions as this stability is deduced from numerical approximation of an evolution equation. In mathematical terms, the local numerical stability of a stationary solution is defined as follows.

Definition 2.1. Stationary solution $u_{*}$ of an abstract evolution (1) is called approximately locally stable, if for every closed finite time interval $[0,\,t_{f}]$ and $\epsilon>0$, there exists $\delta>0$ such that a unique solution of (1), $u(t,\,x)$, with initial condition $u_{0}$ fulfilling $\left\Vert u_{0}-u_{*}\right\Vert _{\mathcal{X}}<\delta$ satisfies

$\left\Vert u(t,\,\cdot)-u_{*}\right\Vert _{\mathcal{X}}<\epsilon\label{eq:stability inequality}$

(7)

for all $t\in[0,\,t_{f}]$.

Having the required definition in hand, we now prove the following stability result for general evolution equations.

Corollary 2.2. Let $u_{*}$ denote a stationary solution of the abstract evolution (1) and $J_{A}(u_{n})$ denote the Jacobian of the approximate system of ODEs defined in (5). If for all sufficiently large $n$ the eigenvalues of $J_{A}(P_{n}u_{*})$ are strictly negative, then $u^{*}$ is approximately locally stable in the sense of Definition 2.1.

Proof. Since the infinitesimal generator approximation scheme, presented in Section 2.1, is convergent, for every given $\epsilon>0$ and finite time interval $[0,\,t_{f}]$ there exist $n_{\epsilon}\in\mathbb{N}$ such that for $n\ge n_{\epsilon}$,

$\left\Vert u(t,\,\cdot)-E_{n}u_{n}(t)\right\Vert _{\mathcal{X}}<\epsilon/2\label{eq:smth one}$

(8)

for all $t\in[0,\,t_{f}]$ (where the bounded linear function $E_{n}$ is defined as in (3)). Moreover, the eigenvalues of $J_{A}(P_{n}u_{*})$ are strictly negative for all sufficiently large $n$. This in turn implies that $P_{M}u_{*}$ is a locally asymptotically stable solution of (5) for some $M\ge n_{\epsilon}$. That is, for given $\epsilon>0$ there is $\delta>0$ such that

$\left\Vert u_{M}(t,\,\cdot)-P_{M}u_{*}\right\Vert _{\mathbb{R}^{M}}=\left\Vert E_{M}u_{M}(t,\,\cdot)-u_{*}\right\Vert _{\mathcal{X}}<\epsilon/2\label{eq:smth two}$

(9)

for all $t>0$ and for all $u_{0}$ such that $\left\Vert P_{M}u_{0}-P_{M}u_{*})\right\Vert _{\mathbb{R}^{M}}=\left\Vert u_{0}-u_{*}\right\Vert _{\mathcal{X}}<\delta$ (see for instance [3,ξ 23]). Consequently, combining (8) and (9) yields

$\left\Vert u(t,\,\cdot)-u_{*}\right\Vert _{\mathcal{X}}\le\left\Vert u(t,\,\cdot)-E_{M}u_{M}(t,\,\cdot)\right\Vert _{\mathcal{X}}+\left\Vert E_{M}u_{M}(t,\,\cdot)-u_{*}\right\Vert _{\mathcal{X}}<\epsilon$

for all $t\in[0,\,t_{f}]$ and for all $u_{0}$ such that $\left\Vert u_{0}-u_{*}\right\Vert _{\mathcal{X}}<\delta$.

We note that although the stability result of Corollary 2.2 holds for arbitrarily large finite time intervals, the Corollary does not guarantee Lyapunov stability of stationary solutions.

---

2.3. Numerical convergence results

To verify convergence of the proposed approximation scheme, we apply the framework to the linear Sinko-Streifer model [52] for which an exact form of the stationary solution is available. The model describes the dynamics of single species populations and takes into account the physiological characteristics of animals of different sizes (and/or ages). The mathematical model reads as

$u_{t}=\mathcal{G}(u)=-(gu)_{x}-\mu u, t\ge0,0\le x\le\overline{x}<\infty\label{eq:sinko-streifer}$

(10)

with a McKendrick-von Foerster type renewal boundary condition at $x=0$

$g(0)u(t,\,0)=\int_{0}^{\overline{x}}q(y)u(t,\,y)\,dy$

and initial condition

$u(0,\,x)=u_{0}(x)\,.$

The variable $u(t,\,x)$ denotes the population density at time $t$ with size class $x$. The population is assumed to have a minimum and a maximum size $0$ and $\overline{x}<\infty$, respectively. The function $g(x)$ represents the average growth rate of the size class $x$ and the coefficient $\mu(\bullet)$ represents a size-dependent removal rate due to death or predation. The renewal function $q(\bullet)$ represents the number of new individuals entering the population due to birth.

Setting the right side of the equation (10) to zero and integrating over the size on $(0,\,x)$ yields the exact stationary solution

$u_{*}(x)=\frac{1}{g(x)}\exp\left(-\int_{0}^{x}\frac{\mu(s)}{g(s)}\,ds\right)\int_{0}^{\overline{x}}q(y)u_{*}(y)\,dy\,.\label{eq:exact stationary solution}$

(11)

Multiplying both sides of (11) by $q(x)$ integrating over the size on $(0,\,\overline{x})$, we obtain a necessary condition for existence of a stationary solution,

$1=\int_{0}^{\overline{x}}\frac{q(x)}{g(x)}\exp\left(-\int_{0}^{x}\frac{\mu(s)}{g(s)}\,ds\right)\,dx\,.\label{eq:necessary condition}$

(12)

The convergence of the approximation scheme presented in Section 2.1 for Sinko-Streifer models has already been established in [6]. Using the basis functions for $n$-dimensional subspace $\mathcal{X}_{n}$ of the state space $\mathcal{X}=L^{1}(0,\,\overline{x})$ are defined as

$\beta_{i}^{n}(x)=\left\{ \begin{array}{cc} 1;\, & x_{i-1}^{n}<x\le x_{i}^{n}\,;\,i=1,\ldots,n\\ 0;\, & \mbox{otherwise} \end{array}\right.$

for positive integer $n$ with $\{x_{i}^{n}\}_{i=0}^{n}$ a uniform partition of $[0,\overline{x}]$, and $\Delta x=x_{j}^{n}-x_{j-1}^{n}$ for all $j$. The functions $\beta^{n}$ form an orthogonal basis for the approximate solution space

${{\cal X}_n} = \left\{ {h \in {\cal X}\;|\;h = \sum\limits_{i = 1}^n {{\alpha _i}} \beta _i^n,\;{\alpha _i} \in }\mathbb{R} \right\},$

and accordingly, we define the orthogonal projections $\pi_{n}:\,\mathcal{X}\to\mathcal{X}_{n}$

${\pi _n}h(x) = \sum\limits_{j = 1}^n {\mkern 1mu} {\alpha _j}\beta _j^n(x),{\rm{where}}\:{\alpha _j} = \frac{1}{{\Delta x}}{\mkern 1mu} \int_{x_{j - 1}^n}^{x_j^n} {\mkern 1mu} h(x){\mkern 1mu} dx.$

Moreover, since the evolution equation defined in (10) is a linear partial differential equation, the approximate operator $\mathcal{G}_{n}$ on $\mathbb{R}^{n}$ is given by the following $n\times n$ matrix

$[\mathcal{G}_{n}]_{ij}= \begin{cases} -\frac{1}{\Delta x}g(x_{i}^{n})-\mu(x_{i}^{n})+q(x_{j}^{n}) & i=j=1\\ q(x_{j}^{n}) & i=1, j\ge2\\ -\frac{1}{\Delta x}g(x_{i}^{n})-\mu(x_{i}^{n}) & i=j\ge2\\ \frac{1}{\Delta x}g(x_{i-1}^{n}) & i=j+1, j\ge1\\ 0 & \text{otherwise} \end{cases}\,.\label{eq:approximate G}$

(13)

At this point, one can use numerical techniques to calculate zeros of the linear system

$\mathcal{G}_{n}u_{n}=0\,.\label{eq:steady state equation}$

(14)

For the purpose of illustration, we set the model rates to

$q(x)=a(x+1), g(x)=b(x+1),\mu(x)=c\,.\label{eq:sinko region}$

(15)

Plugging this rates into the necessary condition (12) yields

$a=\begin{cases} \frac{\ln2}{b} & b=c\\ (b-c)/(2^{1-c/b}-1) & b\ne c \end{cases}\,.$

For $b=c$ this is a curve in 3D coordinate system and when $b\ne c$ the surface is illustrated in Figure 1a. To illustrate the utility of our approach, we used the numerical scheme described in this section to generate existence and stability regions for Sinko-Streifer model. Particularly, the interval $(a,b,\,c)\in[0,\,1]\times[0,\,1]\times[0,\,1]$ is discretized with $\Delta a=\Delta b=\Delta c=0.01$. Consequently, we checked for the existence of a positive steady state at each of these discrete points. Since the approximate operator $\mathcal{G}_{n}$ is an $n\times n$ matrix, we can check if $\mathcal{G}_{n}$ is a singular matrix using standard tools. The resulting existence region is depicted in Figure 1a forming a nontrivial three-dimensional curve (red) on the surface defined by the necessary condition for existence of the steady states. Moreover, the existence and stability regions of the Sinko-Streifer model coincide for the chosen model rates in (15).

For the purpose of illustration, we arbitrarily choose $b=0.5$ and $a=c=1$ (marked with a red star in Figure 1a) as for these values a positive steady state of the Sinko-Streifer model exists. Figure 1b indicates that even for $n=100$ the fit between approximate and actual stationary solution is satisfactory (the infinity norm of the error is $0.004$). Moreover, Figure 1c illustrates that as the dimension of the approximate space $\mathcal{X}_{n}$ increases the absolute error between the exact stationary solution and the approximate stationary solution decreases. Furthermore, the numerical algorithm has a linear convergence rate. This is due to the fact that we chose zeroth order functions as basis functions for approximate subspaces. In general, if one desires a higher order convergence for Galerkin-type approximations, choosing higher order basis functions gives higher convergence rate [33].,

---

3. Application to nonlinear population balance equation

In aerosol physics and environmental sciences, studying the flocculation of particles is widespread. The process of flocculation involves disperse particles in suspension combining into aggregates (i.e., a floc) and separating. The mathematical model used to study flocculation process is the well-known population balance equation (PBE) which describes the time-evolution of the particle size number density. The equations for the flocculation model track the time-evolution of the particle size number density $u(t,\,x)$ and can be written as

${u_t} = {\cal F}(u)$

(16)

where

$\mathcal{F}(u):=\mathcal{G}(u)+\mathcal{A}(u)+\mathcal{B}(u),$

$\mathcal{G}$ denotes growth

$\mathcal{G}(u):=-\partial_{x}(gu)-\mu(x)u(t,\,x)\,,\label{eq:Growth}$

(17)

$\mathcal{A}$ denotes aggregation

$\begin{array}{*{20}{l}} {{\cal A}(u)} & {: = \frac{1}{2}\int_0^x {{k_a}} (x - y,{\mkern 1mu} y)u(t,{\mkern 1mu} x - y)u(t,{\mkern 1mu} y){\mkern 1mu} dy}\\ {} & { - u(t,{\mkern 1mu} x)\int_0^{\bar x - x} {{k_a}} (x,{\mkern 1mu} y)u(t,{\mkern 1mu} y){\mkern 1mu} dy{\mkern 1mu} ,} \end{array}$

(18)

and $\mathcal{B}$ denotes breakage

$\mathcal{B}(u):=\int_{x}^{\overline{x}}\Gamma(x;\,y)k_{f}(y)u(t,\,y)\,dy-\frac{1}{2}k_{f}(x)u(t,\,x)\,.\label{eq:Breakage}$

(19)

The boundary condition is traditionally defined at the smallest size $0$ and the initial condition is defined at $t=0$

$g(0)u(t,\,0)=\int_{0}^{\overline{x}}q(x)u(t,\,x)dx, u(0,\,x)=u_{0}(x)\in L^{1}(0,\,\overline{x})\,,$

where the renewal rate $q(x)$ represents the number of new flocs entering the population. A floc is assumed to have a maximum size $\overline{x}<\infty$. The function $g(x)$ represents the average growth rate of the flocs of size $x$ due to proliferation, and the coefficient $\mu(x)$ represents a size-dependent removal rate due to gravitational sedimentation and death. The function $k_{a}(x,\,y)$ is the aggregation kernel, which describes the rate with which the flocs of size $x$ and $y$ agglomerate to form a floc of size $x+y$. The fragmentation kernel $k_{f}(x)$ calculates the rate with which a floc of size $x$ fragments. The integrable function $\Gamma(x;y)$ represents the post-fragmentation probability density of daughter flocs for the fragmentation of the parent flocs of size $y$. In other words, all the fractions of daughter flocs formed upon the fragmentation of a parent floc sum to unity,

$\int_{0}^{y}\Gamma(x;\,y)\,dx=1\text{ for all }y\in(0,\,\overline{x}].\label{eq:number conservation}$

(20)

The population balance equation, presented in (16), is a generalization of many mathematical models appearing in the size-structured population modeling literature and has been widely used, e.g., to model the formation of clouds and smog in meteorology [49], the kinetics of polymerization in biochemistry [56], the clustering of planets, stars and galaxies in astrophysics [39], and even schooling of fish in marine sciences [47]. For example, when the fragmentation kernel is omitted, $k_{f}\equiv0$, the flocculation model reduces to algal aggregation model used to describe evolution of a phytoplankton community [2]. When the removal and renewal rates are set to zero, the flocculation model simplifies to a model used to describe the proliferation of Klebsiella pneumoniae in a bloodstream [9]. Furthermore, the flocculation model, with only growth and fragmentation terms, was used to investigate the elongation of prion polymers in infected cells [17,26,18].

The equation (16) has also been the focus of considerable mathematical analysis. Well-posedness of the general flocculation model was first established by Ackleh and Fitzpatrick [1,2] in an $L^{2}$-space setting and later by Banasiak and Lamb [5] in an $L^{1}$-space setting. Moreover, asymptotic behavior of the equation (16) has been a challenging task because of the nonlinearity introduced by the aggregation terms. Nevertheless, under suitable conditions on the kernels, the existence of a positive steady state has been established for the pure aggregation and fragmentation case [38]. For a review of further mathematical results, we refer readers to review articles by Menon and Pego [41], and Wattis [54] and the book by Ramkrishna [50]. Lastly, although the population balance equation has received substantial theoretical work, the derivation of analytical solutions for many realistic aggregation kernels has proven elusive. Towards this end, many discretization schemes for numerical simulations of the PBEs have been proposed. For instance, to approximate steady state solutions of PBEs, numerical schemes based on the least squares spectral method [23,24,25] and the finite element method [45,46,31] have been developed. For the further review of approximation methods we refer interested readers to the review by Bortz [8].

---

3.1. Numerical implementation and results

For the numerical implementation we adopt the scheme developed in Section 2.1. Therefore, the approximate formulation of (16) becomes the following system of $n$ nonlinear ODEs for $u_{n}=(\alpha_{1},\cdots,\alpha_{n})^{T}\in\mathbb{R}^{n}$ :

${\dot u_n} = {{\cal F}_n}({u_n}) = {{\cal G}_n}{u_n} + {P_n}{\cal A}({E_n}{u_n}) + {P_n}{\cal B}({E_n}{u_n}),$

(21)

${u_n}(0,x) = {P_n}{u_0}(x){\mkern 1mu} ,$

(22)

where the matrix $\mathcal{G}_{n}$ is defined as in Section 2.3,

$P_{n}\mathcal{A}(E_{n}u_{n})=\begin{pmatrix}-\alpha_{1}\sum\limits_{j=1}^{n-1}k_{a}(x_{1}^{n},\,x_{j}^{n})\alpha_{j}\Delta x\\ \frac{1}{2}k_{a}(x_{1}^{n},\,x_{1}^{n})\alpha_{1}\alpha_{1}\Delta x-\alpha_{2}\sum\limits_{j=1}^{n-2}k_{a}(x_{2}^{n},\,x_{j}^{n})\alpha_{j}\Delta x\\ \vdots\\ \frac{1}{2}\sum\limits_{j=1}^{n-2}k_{a}(x_{j}^{n},\,x_{n-1-j}^{n})\alpha_{j}\alpha_{n-1-j}\Delta x-\alpha_{n-1}k_{a}(x_{n-1}^{n},\,x_{1}^{n})\alpha_{1}\Delta x\\ \frac{1}{2}\sum\limits_{j=1}^{n-1}k_{a}(x_{j}^{n},\,x_{n-j}^{n})\alpha_{j}\alpha_{n-j}\Delta x \end{pmatrix}$

and

$P_{n}\mathcal{B}(E_{n}u_{n})=\begin{pmatrix}\sum\limits_{j=2}^{n}\Gamma(x_{1}^{n};\,x_{j}^{n})k_{f}(x_{j}^{n})\alpha_{j}\Delta x-\frac{1}{2}k_{f}(x_{1}^{n})\alpha_{1}\\ \sum\limits_{j=3}^{n}\Gamma(x_{2}^{n};\,x_{j}^{n})k_{f}(x_{j}^{n})\alpha_{j}\Delta x-\frac{1}{2}k_{f}(x_{2}^{n})\alpha_{2}\\ \vdots\\ \Gamma(x_{n-1}^{n};\,x_{n}^{n})k_{f}(x_{n}^{n})\alpha_{n}\Delta x-\frac{1}{2}k_{f}(x_{n-1}^{n})\alpha_{n-1}\\ -\frac{1}{2}k_{f}(x_{n}^{n})\alpha_{n} \end{pmatrix}\,.$

The convergence of the approximate scheme (21)-(22) has been established in [1]. Therefore, the stationary solutions of the microbial flocculation model (16) can be systematically approximated by the stationary solutions of the system of nonlinear ODEs given in (21). We used Powell's hybrid root finding method [48] as implemented in Python 2.7.10¹ to find zeros of the steady state equation (see available code at [42]). For faster convergence rate, we provided the solver with the exact Jacobian of $\mathcal{F}_{n}(u_{n})$ (see Section 2.2, Eqn 23) for the formulation of the Jacobian). For the purpose of illustration, for a post-fragmentation density function we chose the well-known Beta distribution ² with $\alpha=\beta=2$,

$\Gamma(x,\,y)={1}_{[0,\,y]}(x)\frac{6x(y-x)}{y^{3}}\,,$

¹ scipy.optimize.fsolve

² Although normal and log-normal distributions are mostly used in the literature, Byrne et al. [16] have provided evidence that the Beta density function describes the fragmentation of small bacterial flocs.

where ${1}_{I}$ is the indicator function on the interval $I$. The aggregation kernel was chosen to describe flow within laminar shear field (i.e., orthokinetic aggregation)

$k_{a}(x,\,y)=\left(x^{1/3}+y^{1/3}\right)^{3}$

Other model rates were chosen arbitrarily as

$q(x)=a(x+1), g(x)=b(x+1)\mu(x)=cx k_{f}(x)=x\,,$

where $a,\,b$ and $c$ are some positive real numbers.

The main advantage of this approximation scheme (21)-(22) is that it can be initialized very fast using Toeplitz matrices [40]. Fast initialization of the discretization scheme allows one to check the existence of the steady states at many discrete points efficiently. This in turn allows for the generation of the existence and stability regions of the steady states of the PBE in (16). To illustrate the existence regions of the steady states of the PBE, we discretized the intervals $a\in[0,\,15]\,,b\in[0,\,1]$ and $c\in[0,\,5]$ with $\Delta a=\Delta b=\Delta c=0.1$. We note that for faster convergence the root finding method needs an initial seed close to the steady state solution. Since we have no information about the existing steady state, we seed the root finding method with $10$ different uniform initial guesses i.e.,

$\left\{ u_{0}(x)=2^{i}\,|\,i=0,1,\dots,9\right\} \,,$

before we conclude a positive steady state does not exist for a given point $(a,\,b,\,c)$. Consequently, we checked for the existence of a positive steady state at each of these discrete points. As depicted in Figure 2a, approximate existence region of positive steady states of the PBE forms a three dimensional wedge like region. Moreover, in Figures 2b-2d, to deduce stability of each steady state solution, we checked the spectrum of the Jacobian matrix evaluated at each steady state. Particularly, if the real part of rightmost eigenvalue of the Jacobian matrix is negative, the steady state is identified as locally stable (blue region). Conversely, if the real part of the rightmost eigenvalue of the Jacobian matrix is positive the steady state is identified as unstable (red region). One can observe that growth ($b$) and removal ($c$) rates can balance the smaller renewal rates ($a$), and thus locally stable steady states exist. However, as the renewal rate gets larger steady states first become unstable and then cease to exist (yellow region). This is also illustrated in Figure 3b, where steady states start diverging for the larger renewal rates ($a$).

Figure 3a illustrates an example stationary solution for $b=0.5,\,a=c=1$. To confirm that the function depicted in Figure 4 is indeed a locally stable steady state, we simulated the system of ODEs in (21)-(22) for $t\in[0,\,10]$ with a collection of arbitrary initial conditions (Figure 4a) close to the steady state solution. One can observe in Figure 4 that the stationary solution is indeed locally stable and thus initial conditions, Figure 4, converge to the steady state depicted in Figure 4b. As depicted in Figure 4c and Figure 4d, convergence is also reflected in the evolution of the total number of flocs (zeroth moment),

${M_0}(t) = \int_0^{\bar x} u (t,{\mkern 1mu} x){\mkern 1mu} dx \approx \sum\limits_{i = 1}^n {\int_{x_{i - 1}^n}^{x_i^n} {{\alpha _i}} } \beta _i^n(x){\mkern 1mu} dx = \Delta x\sum\limits_{i = 1}^n {{\alpha _i}} {\mkern 1mu} ,$

and total mass of the flocs (first moment),

${M_1}(t) = \int_0^{\bar x} x u(t,{\mkern 1mu} x){\mkern 1mu} dx \approx \sum\limits_{i = 1}^n {\int_{x_{i - 1}^n}^{x_i^n} {{\alpha _i}} } x\beta _i^n(x){\mkern 1mu} dx = \frac{{\Delta x}}{2}\sum\limits_{i = 1}^n {{\alpha _i}} \left( {x_i^n + x_{i - 1}^n} \right){\mkern 1mu} .$

Moreover, to confirm that the steady state solution is not changing with increasing dimension of approximate spaces $\mathcal{X}_{n}$, we simulated our numerical scheme for different values of $n$. Figure 5a illustrates that stationary solutions converge to exact stationary solution of (16) as $n\to\infty$. Furthermore, one can observe, in Figure (5) b that difference between approximate steady states for different values of $n$ is considerably small.

---

3.2. Conditions for numerical stability of positive steady states

In this section, we derive conditions for approximate local stability of the stationary solution of the nonlinear population balance equation defined in (16). In particular, we impose conditions on the model rates of the population balance equation for which the first statement of Corollary 2.2 holds. Towards this end, we use the well-known Gershgorin theorem for locating eigenvalues of a matrix. The Gershgorin theorem states that each eigenvalue of $A$ lies in one of the the circular areas, called Gershgorin disks, in the complex plane that is centered at the $i$th diagonal element and whose radius is $R_{i}$,

$\left\{ z\in\mathbb{C}\,:\,\left|z-a_{ii}\right|\le R_{i}\right\} \,,$

where

${R_i} = \sum\limits_{j = 1,{\kern 1pt} j \ne i}^n | {a_{ji}}|{\mkern 1mu} .$

Since the approximate system for the microbial flocculation model is nonlinear, we linearize the system around its stationary solutions. Let $u_{*}\in L^{1}(0,\,\overline{x})$ be a stationary solution of (16) and denote the projection of the stationary solution $u_{*}$ onto $\mathbb{R}^{n}$ by $\alpha=P_{n}u_{*}=[\alpha_{1},\cdots,\alpha_{n}]^{T}$, then the Jacobian of the approximate operator $\mathcal{F}_{n}$ defined in (21) can be written as

$J_{\mathcal{F}}(\alpha)=\mathcal{G}_{n}+J_{\mathcal{A}}(\alpha)+J_{\mathcal{B}}(\alpha)\,,\label{eq:jacobian of F}$

(23)

where $\mathcal{G}_{n}$ is defined in (13),

$\begin{array}{l} {J_A}(\alpha ) = \left( {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} { - {\alpha _1}{k_a}(x_1^n,{\kern 1pt} x_1^n)\Delta x} & { - {\alpha _1}{k_a}(x_1^n,{\kern 1pt} x_2^n)\Delta x} & \cdots & { - {\alpha _1}{k_a}(x_1^n,{\kern 1pt} x_{n - 1}^n)\Delta x} & 0\\ { - {\alpha _2}{k_a}(x_2^n,{\kern 1pt} x_1^n)\Delta x} & \cdots & { - {\alpha _2}{k_a}(x_2^n,{\kern 1pt} x_{n - 2}^n)\Delta x} & 0 & 0\\ \vdots & \vdots & \ddots & \vdots & \vdots \\ { - {\alpha _{N - 1}}{k_a}(x_{n - 1}^n,{\kern 1pt} x_1^n)\Delta x} & 0 & \cdots & 0 & 0\\ 0 & 0 & \cdots & 0 & 0 \end{array}} \end{array}} \right)\\ \quad \quad \quad + \left( {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} { - \sum\limits_{j = 1}^{n - 1} {{k_a}} (x_1^n,{\kern 1pt} x_j^n){\alpha _j}\Delta x} & 0 & \cdots & 0 & 0\\ {{\alpha _1}{k_a}(x_1^n,{\kern 1pt} x_1^n)\Delta x} & { - \sum\limits_{j = 1}^{n - 2} {{k_a}} (x_2^n,{\kern 1pt} x_j^n){\alpha _j}\Delta x} & 0 & \cdots & 0\\ {{\alpha _2}{k_a}(x_1^n,{\kern 1pt} x_2^n)\Delta x} & {{\alpha _1}{k_a}(x_2^n,{\kern 1pt} x_1^n)\Delta x} & \ddots & 0 & \vdots \\ \vdots & \vdots & \ddots & { - \sum\limits_{j = 1}^1 {{k_a}} (x_{n - 1}^n,{\kern 1pt} x_j^n){\alpha _j}\Delta x} & 0\\ {{\alpha _{n - 1}}{k_a}(x_1^n,{\kern 1pt} x_{n - 1}^n)\Delta x} & {{\alpha _{n - 2}}{k_a}(x_2^n,{\kern 1pt} x_{n - 2}^n)\Delta x} & \cdots & {{\alpha _1}{k_a}(x_{n - 1}^n,{\kern 1pt} x_1^n)\Delta x} & 0 \end{array}} \end{array}} \right), \end{array}$

and

$J_{\mathcal{B}}(\alpha)=\begin{pmatrix}\begin{smallmatrix}-\frac{1}{2}k_{f}(x_{1}^{n}) & \Gamma(x_{1}^{n};x_{2}^{n})k_{f}(x_{2}^{n})\Delta x & \Gamma(x_{1}^{n};x_{3}^{n})k_{f}(x_{3}^{n})\Delta x & \cdots & \Gamma(x_{1}^{n};x_{n}^{n})k_{f}(x_{n}^{n})\Delta x\\ 0 & -\frac{1}{2}k_{f}(x_{2}^{n}) & \Gamma(x_{2}^{n};x_{3}^{n})k_{f}(x_{3}^{n})\Delta x & \cdots & \Gamma(x_{2}^{n};x_{n}^{n})k_{f}(x_{n}^{n})\Delta x\\ \vdots & 0 & \ddots & \ddots & \vdots\\ 0 & \cdots & 0 & -\frac{1}{2}k_{f}(x_{n-1}^{n}) & \Gamma(x_{n-1}^{n};x_{n}^{n})k_{f}(x_{n}^{n})\Delta x\\ 0 & 0 & \cdots & 0 & -\frac{1}{2}k_{f}(x_{n}^{n}) \end{smallmatrix}\end{pmatrix}$

To bound the eigenvalues of $J_{\mathcal{F}}(\alpha)$ we use Gershgorin theorem. Consequently, the centers and the radii of Gershgorin disks are given by

${a_{ii}} = - \frac{1}{{\Delta x}}g(x_i^n) - \mu (x_i^n) - \frac{1}{2}{k_f}(x_i^n) - {\alpha _i}{k_a}(x_i^n,{\mkern 1mu} x_i^n)\Delta x - \sum\limits_{j = 1}^{n - i} {{k_a}} (x_i^n,{\mkern 1mu} x_j^n){\alpha _j}\Delta x$

and

${R_i} \le \frac{1}{{\Delta x}}g(x_i^n) + q(x_i^n) + \sum\limits_{j = 1}^{i - 1} \Gamma (x_j^n;x_i^n){k_f}(x_i^n)\Delta x + \sum\limits_{j = 1}^{n - i} {{\alpha _j}} {k_a}(x_j^n,{\mkern 1mu} x_i^n)\Delta x$

$+ \sum\limits_{j = 1,{\kern 1pt} j \ne i}^{n - i} {{\alpha _j}} {k_a}(x_i^n,{\mkern 1mu} x_j^n)\Delta x{\mkern 1mu} ,$

respectively. Consequently, if we can show that

$\left|a_{ii}\right|>R_{i}\text{for each }i\in\left\{ 1,\dots,n\right\} \,,\label{eq:center larger than radius}$

(24)

then each of the Gershgorin disks lie strictly on the left side of the complex plane. To this end, inequality (24) can be simplified as

$\mu (x_i^n) + \frac{1}{2}{k_f}(x_i^n) > q(x_i^n) + \sum\limits_{j = 1}^{i - 1} \Gamma (x_j^n;x_i^n){k_f}(x_i^n)\Delta x + \sum\limits_{j = 1}^{n - i} {{\alpha _j}} {k_a}(x_j^n,{\mkern 1mu} x_i^n)\Delta x$

(25)

for each $i\in\left\{ 1,\cdots,n\right\}$. Accordingly, taking the limit of (25) as $n\to\infty$ yields

$\mu(x)+\frac{1}{2}k_{f}(x)>q(x)+\int_{0}^{x}\Gamma(y,\,x)k_{f}(x)\,dy+\int_{0}^{\overline{x}-x}k_{a}(x,\,y)u_{*}(y)\,dy\label{eq:precise inequality}$

(26)

for all $x\in[0,\,\overline{x}]$ and together with the number conservation requirement (20) implies

$q(x)+\frac{1}{2}k_{f}(x)-\mu(x)+\int_{0}^{\overline{x}-x}k_{a}(x,\,y)u_{*}(y)\,dy<0$

for all $x\in[0,\,\overline{x}]$. Conversely, note that the integral approximations in (25) are right Riemann sums. Therefore, if the functions $\Gamma(y,\,x)$ and $k_{a}(x,\,y)u_{*}(y)$ are decreasing in $y$ then integral approximations in (25) are under-approximations of the integrals in (26). Thus, the inequality stated in (26) ensures that the eigenvalues of the Jacobian $J_{\mathcal{F}}(\alpha)$ are strictly negative for all sufficiently large $n$. Now, we are in a position to summarize the results of this section in the following proposition.

Proposition 3.1. Let $u_{*}$ be a stationary solution of the microbial flocculation model (16). Moreover, if

$q(x)+\frac{1}{2}k_{f}(x)-\mu(x)+\int_{0}^{\overline{x}-x}k_{a}(x,\,y)u_{*}(y)\,dy<0\label{eq:first condition}$

(27)

for all $x\in[0,\,\overline{x}]$ and

$\partial_{y}\left(k_{a}(x,\,y)u_{*}(y)\right)\le0\text{ and }\partial_{y}\Gamma(y,\,x)\le0\label{eq:second condition}$

(28)

for all $x\in[0,\,\overline{x}]$ and $y\in[0,\,\overline{x}]$, then stationary solution $u_{*}$ is approximately locally stable in the sense of Definition 2.1.

To illustrate the utility of the framework developed in this section we applied our approach to the model rates given in Section 3.1 for generation of Figure 3a. The Beta distribution used for the post-fragmentation function $\Gamma$ is not monotonically decreasing, and thus it does not satisfy the conditions of Proposition 3.1. However, Figure 6a-c illustrates that the model rates satisfy the conditions of Corollary 2.2. In fact, Figure 6d depicts that for the steady state illustrated in Figure 3a the eigenvalues of $J_{\mathcal{F}}(P_{n}u_{*})$ have strictly negative real part for $n\ge5$. Therefore, as it has already been established in Figure 4, this steady state solution is locally asymptotically stable in the sense of Corollary 2.2.

---

4. Concluding remarks

Our primary motivation in this paper was to show that available numerical tools in the literature can facilitate theoretical analysis of evolution equations. Towards this end we developed a numerical framework for theoretical analysis of evolution equations arising in population dynamical models. The main idea behind this framework is to approximate generators of semigroups of evolution equations and use numerical tools to study stability of steady states of evolution equations. We demonstrated the utility of our approach by applying the numerical framework to both linear and nonlinear size-structured population models. In particular, we generated approximate existence and stability regions of the steady states of both models (which can be difficult to obtain by using conventional analytical tools). We showed that our numerical framework can also be used to gain insight about the approximate local stability (see Definition 2.1) of stationary solutions. Furthermore, code used for the numerical simulations can be obtained from our GitHub repository under the open source MIT License (MIT) [42].

Although the stability inequality in (7) holds for arbitrarily large finite time intervals, behavior of the solutions as $t\to\infty$ is unclear. Hence, we note that the local stability of the stationary solutions specified in Corollary 2.2 is not in a usual Lyapunov sense. In order to prove Lyapunov stability of stationary solutions using the approximation scheme presented in Section 2.1, one has to prove uniform convergence of the approximation scheme for all $t\ge0$. Hence, as a future plan we wish to utilize the numerical framework presented here to establish Lyapunov stability of stationary solutions of general evolution equations.

---

Acknowledgments

Funding for this research was supported in part by grants NSF-DMS 1225878 and NIH-NIGMS 2R01GM069438-06A2. This work utilized the Janus supercomputer, which is supported by the National Science Foundation (award number CNS-0821794) and the University of Colorado Boulder. The Janus supercomputer is a joint effort of the University of Colorado Boulder, the University of Colorado Denver and the National Center for Atmospheric Research.

---