A numerical framework for computing steady states of size-structured population models and their stability

Structured population models are a class of general evolution equations which are widely used in the study of biological systems. Many theoretical methods are available for establishing existence and stability of steady states of general evolution equations. However, except for very special cases, finding an analytical form of stationary solutions for evolution equations is a challenging task. In the present paper, we develop a numerical framework for computing approximations to stationary solutions of general evolution equations, which can also be used to produce existence and stability regions for steady states. In particular, we use the Trotter-Kato Theorem to approximate the infinitesimal generator of an evolution equation on a finite dimensional space, which in turn reduces the evolution equation into a system of ordinary differential equations. Consequently, we approximate and study the asymptotic behavior of stationary solutions. We illustrate the convergence of our numerical framework by applying it to a linear Sinko-Streifer structured population model for which the exact form of the steady state is known. To further illustrate the utility of our approach, we apply our framework to nonlinear population balance equation, which is an extension of well-known Smoluchowksi coagulation-fragmentation model to biological populations. We also demonstrate that our numerical framework can be used to gain insight about the theoretical stability of the stationary solutions of the evolution equations. Furthermore, the open source Python program that we have developed for our numerical simulations is freely available from our Github repository (github.com/MathBioCU).


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 combined with numerical integration of the ODE.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 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 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.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.2.In Section 2.2, 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-structure population model (also known as population balance equations in the engineering literature) described in [5,9].Moreover, in Section 4, we show that 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.

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.2, we illustrate the convergence of our approach by applying it to linear Sinko-Streifer equations, for which exact stationary solutions are known.

Infinitesimal generator approximation
Given a Banach space X , consider an abstract evolution equation, where F : D(F) ⊆ X → X is some operator defined on its domain D(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 D(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 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 ≥ 0 2. T (0) = I, the identity operator on X 3.For each fixed u 0 ∈ X , lim Moreover, showing that the operator 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 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 F n on the approximate spaces X n such that C 0 -semigroups T n (•) generated by F n approximate the C 0 -semigroup T (t) generated by F. Although there are multiple ways to approximate the infinitesimal generator F, for our purposes we use the approximation scheme based on Galerkin-type projection of the state space X [6,32,2].For the convenience of readers, we will summarize the approximation scheme here.Let X n , n = 1, 2, . . .be a sequence of subspaces of X with dim(X n ) = n and define projections π n : X → X n and canonical injections ι n : X n → X .For each subspace X n we choose basis elements {β n i } n i=1 such that each element v of subspace X n can be written in the form v = n i=1 α i β n i .Moreover, for each subspace X n we define the bijective mappings and the norm on R n by v R n = p −1 n v X .Consequently, we define bounded linear operators P n : X → R n and E n : R n → X by and respectively.Finally, we define approximate operators F n on R n by for all z ∈ R n .Accordingly, the Trotter-Kato Theorem states that the semigroup generated by the discrete operator F n converges to the semigroup generated by the infinitesimal generator F. For the convenience of the reader, we state the theorem here and refer readers to [32] for a proof.
Theorem.(Trotter-Kato) Let (T (t)) t≥0 and (T n (t)) t≥0 , n ∈ N, be strongly continuous semigroups on X and R n with generators F and F n , respectively.Furthermore, assume that they satisfy the estimate for some constants M ≥ 1, w ∈ R. Then the following are equivalent: 1.There exists a λ 0 ∈ ρ(F) as n → ∞, 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, Consequently, the solution of the IVP is mapped onto the infinite dimensional Banach space X and one has the following convergence 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 → ∞.

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 with a McKendrick-von Foerster type renewal boundary condition at x = 0 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 x < ∞, respectively.The function g(x) represents the average growth rate of the size class x and the coefficient µ(•) represents a size-dependent removal rate due to death or predation.The renewal function q(•) represents the number of new individuals entering the population due to birth.
Setting the right side of the equation ( 7) to zero and integrating over the size on (0, x) yields the exact stationary solution Multiplying both sides of (8) by q(x) integrating over the size on (0, x), we obtain a necessary condition for existence of a stationary solution, We also note that if u * is a stationary solution satisfying (8), then any multiple of u * is also stationary solution of (7).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 X n of the state space X = L 1 (0, x) are defined as otherwise for positive integer n with {x n i } n i=0 a uniform partition of [0, x], and ∆x = x n j − x n j−1 for all j.The functions β n form an orthogonal basis for the approximate solution space and accordingly, we define the orthogonal projections π n : X → X n , where α j = 1 ∆x Moreover, since the evolution equation defined in ( 7) is a linear partial differential equation, the approximate operator G n on R n is given by the following n × n matrix At this point, one can use numerical techniques to calculate zeros of the linear system For the purpose of illustration, we set the model rates to For the purpose of illustration, we arbitrarily choose a = 1/ ln 2 and b = c = 1 as these values satisfy the necessary condition (9) for the existence of the steady states of the Sinko-Streifer model.Since the approximate operator G n is an n × n matrix, we can compute the nullspace of G n using standard tools.The results of the numerical simulations are depicted in Figure 1. Figure 1a illustrates that as the dimension of the approximate space X n increases the absolute error between the exact stationary solution and the approximate stationary solution decreases.Moreover, 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].Furthermore, 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.14).
To further 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) ∈ [0, 1] × [0, 1] × [0, 1] is discretized with ∆a = ∆b = ∆c = 0.01.Consequently, we checked for the existence of a positive steady state at each of these discrete points, the resulting existence region is depicted in Figure 1c forming a nontrivial three-dimensional surface.Moreover, the existence and stability regions of the Sinko-Streifer model coincide for the chosen model rates in (12).

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 where A denotes aggregation and B denotes breakage The boundary condition is traditionally defined at the smallest size 0 and the initial condition is defined at t = 0 where the renewal rate q(x) represents the number of new flocs entering the population.
The population balance equation, presented in (13), 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 ≡ 0, 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 ( 13) 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 ( 13) 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].

Numerical implementation and results
For the numerical implementation we adopt the scheme developed in Section 2.1.Therefore, the approximate formulation of (13) becomes the following system of n nonlinear ODEs for where the matrix G n is defined as in Section 2.2, The convergence of the approximate scheme ( 18)-( 19) has been established in [1].Therefore, the stationary solutions of the microbial flocculation model ( 13) can be systematically approximated by the stationary solutions of the system of nonlinear ODEs given in (18).We used Powell's hybrid root finding method [48] as implemented in Python 2.7.10 1 to find zeros of the steady state equation [42].For faster convergence rate, we provided the solver with the exact Jacobian of F n (u n ) ( see Section 4, Eqn (24) for the formulation of the Jacobian.For the purpose of illustration, for a post-fragmentation density function we chose the well-known Beta distribution 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) Other model rates were chosen arbitrarily as where a, b and c are some positive real numbers.
The main advantage of this approximation scheme ( 18)-( 19) 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 (13).To illustrate the existence regions of the steady states of the PBE, we discretized the intervals a ∈ [0, 15] , b ∈ [0, 1] and c ∈ [0, 5] with ∆a = ∆b = ∆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., u 0 (x) = 2 i | i = 0, 1, . . ., 9 , 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, the 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 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 ( 18)-( 19) for t ∈ [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 Figures 4c and 4d, convergence is also reflected in the evolution of the total number of flocs (zeroth moment), 1 scipy.optimize.fsolve 2 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.
Moreover, to confirm that the steady state solution is not changing with increasing dimension of approximate spaces X n , we simulated our numerical scheme for different values of n. Figure 5a illustrates that stationary solutions converge to exact stationary solution of (13) as n → ∞.Furthermore, one can observe, in Figure (5)b that difference between approximate steady states for different values of n is considerably small.

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 has 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 can also be used for deriving stability conditions for stationary solutions of the general evolution equations.Towards this end, we prove the following stability result for general evolution equations.
Corollary 4.1.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 a locally asymptotically stable stationary solution of the evolution equation (1).In particular, for every closed finite time interval [0, t f ] and > 0, there exists δ > 0 such that a unique solution of (1), u(t, x), with initial condition Proof.Since the infinitesimal generator approximation scheme, presented in Section 2.1, is convergent, for every given > 0 and finite time interval [0, t f ] there exist n ∈ N such that for n ≥ n , for all t ∈ [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 ≥ n .That is, for given > 0 there is δ > 0 such that for all t > 0 and for all u 0 such that P M u 0 − P M u * ) R M = u 0 − u * X < δ (see for instance [3, §23]).Consequently, combining (21) and (22) yields for all t ∈ [0, t f ] and for all u 0 such that u 0 − u * < δ.
Having the required corollary in hand, in subsequent sections we apply it to two different examples from size-structured population modeling.In Section 4.1, we derive conditions for local stability of the stationary solutions of linear Sinko-Streifer equation (7).In Section 4.2, we derive conditions for local stability of the stationary solution of the nonlinear population balance equation defined in (13).

Sinko-Streifer model
To prove the first statement of Corollary 4.1, 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 ith diagonal element and whose radius is Consequently, by applying the Gershgorin theorem to columns of the matrix G n , the approximate ODE system for Sinko-Streifer derived in Section 2.2, yields that its eigenvalues are located in the following Gershgorin disks Therefore, a sufficient condition for eigenvalues of the matrix G n to be strictly negative is for each i = 1, • • • , n.Hence, the condition q(x) − µ(x) < 0 for all x ∈ (0, x) ensures the local stability of the non-trivial stationary solution of Sinko-Streifer equation, which is also in an agreement with the stability conditions derived in [44, Condition 2].We summarize the results of this subsection in the following proposition.
Proposition 4.2.Let u * be a stationary solution of the Sinko-Streifer model defined in (7).Moreover, assume that q(x) − µ(x) < 0 for all x ∈ [0, x], then the stationary solution u * is locally stable in the sense of Corollary 4.1.

Population balance equation
Since the approximate system for the microbial flocculation model is nonlinear, we linearize the system around its stationary solutions.Let u * ∈ L 1 (0, x) be a stationary solution of ( 13) and denote the projection of the stationary solution u * onto R n by α = P n u * = [α 1 , • • • , α n ] T , then the Jacobian of the approximate operator F n defined in ( 18) can be written as where G n is defined in (10), To bound the eigenvalues of J F (α) again we use Gershgorin theorem.Consequently, the centers and the radii of Gershgorin disks are given by respectively.Consequently, if we can show that then each of the Gershgorin disks lie strictly on the left side of the complex plane.To this end, inequality (25) can be simplified as for each i ∈ {1, • • • , n}.Accordingly, taking the limit of (26) as n → ∞ yields for all x ∈ [0, x] and together with the number conservation requirement (17) implies x−x 0 k a (x, y)u * (y) dy < 0 for all x ∈ [0, x].Conversely, note that the integral approximations in (26) are right Reimann sums.Therefore, if the functions Γ(y, x) and k a (x, y)u * (y) are decreasing in y then integral approximations in (26) are under-approximations of the integrals in (27).Thus, the inequality stated in (27) ensures that the eigenvalues of the Jacobian J F (α) 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.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 Γ is not monotonically decreasing, and thus it does not satisfy the conditions of Proposition 4.3.However, Figure 6a-c illustrates that the model rates satisfy the conditions of Corollary 4.1.In fact, Figure 6d depicts that for the steady state illustrated in Figure 3a the eigenvalues of J F (P n u * ) have strictly negative real part for n ≥ 5. Therefore, as it has already been established in Figure 4, this steady state solution is locally asymptotically stable in the sense of Corollary 4.1.

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 the existence and stability regions of the steady states of the 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 local stability 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 (20) holds for all finite time intervals, behavior of the solutions as t → ∞ is unclear.Hence, we note that the local stability of the stationary solutions specified in Corollary 4.1 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 ≥ 0. 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.

Figure 1 :
Figure 1: Results of the numerical simulations.a) Absolute error between exact stationary solution and approximate stationary solution decreases as the dimension of approximate subspaces X n increase.b) Comparison of exact stationary solution with approximate stationary solution for n = 100.c) Existence and stability surface for the steady states of the Sinko-Streifer model.
A floc is assumed to have a maximum size x < ∞.The function g(x) represents the average growth rate of the flocs of size x due to proliferation, and the coefficient µ(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 Γ(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, y 0 Γ(x; y) dx = 1 for all y ∈ (0, x].

Figure 2 :
Figure 2: Existence and stability regions for the steady states of the PBE a) Existence region for the steady states of the PBE forms a wedge like shape.b) Stability region for b = 0.1, a ∈ [0, 15] and c ∈ [0, 5].c) Stability region for b = 0.5, a ∈ [0, 15] and c ∈ [0, 5].d) Stability region for b = 1.0, a ∈ [0, 15] and c ∈ [0, 5].Color bar represents the real part of rightmost eigenvalue of the Jacobian matrix evaluated at each steady state.Yellow regions represents the region for which a positive steady state does not exists.

Figure 3 :
Figure 3: a) An example steady-state solution of the PBE for b = 0.5, a = c = 1.b) Steady states for increasing renewal rate and b = c = 1

Figure 4 :
Figure 4: Time evolution of the flocculation model with arbitrary initial conditions.a) Four different initial conditions are chosen close to the steady state.b) Solution of the PBE for those initial conditions at t = 10.c) Evolution of the total number M 0 (t) of the flocs for t ∈ [0, 10].d) Evolution of the total mass M 1 (t) of the flocs for t ∈ [0, 10].

Figure 5 :
Figure 5: Change in zeroth and first moments with increasing dimension of the approximate space X n .a) Change in the total number and the total mass of the flocs with respect to increasing dimension n.Dashed red lines and dotted green lines corresponds to the total number and the total mass of the flocs of the steady state for n = 1000, respectively.b) Steady state solution for n = 100 and n = 500.

Figure 6 :
Figure 6: Eigenvalues of the Jacobian J F (α) multiplied by ∆x for the steady state illustrated in Figure 3a.a) Eigenvalues of the Jacobian plotted in the complex plane for n = 20.b) Eigenvalues of the Jacobian plotted in the complex plane for n = 50.c) Eigenvalues of the Jacobian plotted in the complex plane for n = 200.d) Change in the rightmost eigenvalue for increasing n.Dashed black line corresponds to the rightmost eigenvalue of the Jacobian for n = 1000.