# Efficient Numerical Schemes for
Multidimensional Population Balance Models

Pavan K. Inguva  Massachusetts Institute of Technology, 77 Massachusetts
Avenue, Cambridge, MA 02139, USA  Richard D. Braatz  Massachusetts Institute
of Technology, 77 Massachusetts Avenue, Cambridge, MA 02139, USA

###### Abstract

Multidimensional population balance models (PBMs) describe chemical and
biological processes having a distribution over two or more intrinsic
properties (such as size and age, or two independent spatial variables). The
incorporation of additional intrinsic variables into a PBM improves its
descriptive capability and can be necessary to capture specific features of
interest. As most PBMs of interest cannot be solved analytically,
computationally expensive high-order finite difference or finite volume
methods are frequently used to obtain an accurate numerical solution. We
propose a finite difference scheme based on operator splitting and solving
each sub-problem at the limit of numerical stability that achieves a
discretization error that is zero for certain classes of PBMs and low enough
to be acceptable for other classes. In conjunction to employing specially
constructed meshes and variable transformations, the scheme exploits the
commutative property of the differential operators present in many classes of
PBMs. The scheme has very low computational cost – potentially as low as just
memory reallocation. Multiple case studies demonstrate the performance of the
proposed scheme.

##  1 Introduction

Multidimensional population balance models (PBM) are of significant interest
due to having the ability to describe population dynamics that vary across
multiple intrinsic variables. Examples of such populations include crystals
that vary along two independent spatial dimensions such as length and width
[1, 2, 3, 4]; granules that vary in solid, liquid, and gas content [5, 6] or
porosity, binder content, and composition [7]; and cell populations that vary
in multiple properties such as cell size, age, and intracellular
concentrations of species of interest such as enzymes [8, 9, 10]. Another
important class of problems where multidimensional PBMs are valuable are
systems where, in addition to variations in one or more intrinsic variables,
there is spatial variation such as in slug or plug flow in a tube which
introduces an additional intrinsic variable in the form of axial position or
residence time [11, 12, 13, 14, 15].

Most PBMs of interest cannot be solved analytically, and so many numerical
methods have been developed. Compared to 1D PBMs, multidimensional PBMs are
more challenging to solve numerically, as the solution process is much more
memory intensive and prone to numerical diffusion and/or dispersion. As such,
various numerical schemes have been specifically targeted at multidimensional
PBMs based on the finite difference method (FDM) [8, 16], the finite volume
method (FVM) [17, 5, 18, 19, 20], the finite element method (FEM) [21, 22],
spectral methods [23], and Lattice Boltzmann methods [24]. Often, these
methods are computationally costly and/or highly mathematically involved which
impacts their adoption and deployment.

In previous work, we demonstrated how FDM, when thoughtfully applied using
specially constructed meshes and variable transformations, can accurately and
efficiently solve 1D PBMs [15]. This article describes the extension of those
methods to multidimensional PBMs. Many classes of PBMs have differential
operators that commute which enables the use of operator splitting techniques
with no splitting error [25]. These operator splitting techniques transform a
multidimensional PBM into a series of 1D sub-problems, each of which can be
solved highly accurately and efficiently as previously shown. Although the
presented methods and case studies focus on 2D PBMs, the extension to
n-dimensional PBMs naturally follows.

##  2 Theory and Methods

This section is structured as a series of cases in which each case outlines
the development of the finite difference numerical scheme for a specific class
of PBM. Finite difference schemes are benchmarked with the SharpClaw solver
[26], a high-order weighted essentially non-oscillatory (WENO) solver from the
PyClaw package [27, 28]. Default solver settings with a maximum CFL number of
1 for all problems were used. For cases where the growth rate in the PBM is
variable (in the context of hyperbolic conservation equations, this is also
called variable velocity advection), PyClaw requires the equation to be
formulated non-conservatively which is done by employing a variable
transformation (e.g., see (16)) or by expanding out the “spatial” derivatives.

###  2.1 Case 1: PBMs with Constant Growth Rate

Consider the homogeneous PBM,

$\frac{\partial f(t,a_{1},a_{2})}{\partial t}+g_{1}\frac{\partial
f(t,a_{1},a_{2})}{\partial a_{1}}+g_{2}\frac{\partial
f(t,a_{1},a_{2})}{\partial a_{2}}=0,\quad
f(0,a_{1},a_{2})=f_{0}(a_{1},a_{2}),$ (1)

where $a_{1}$ and $a_{2}$ are the intrinsic variables and $g_{1}$ and $g_{2}$
are the constant positive growth rates. PBMs of this form arise when the
driving force for growth is constant which, for example, can occur in
crystallization where an external control system maintains a constant
supersaturation (see [15] and citations therein). The efficient and accurate
numerical solution to this class of PBMs using the techniques described in
this section is already known (e.g., see [29, 30, 25]). The solution of the
PBM (1) is considered here for pedagogical reasons, as a robust understanding
of this section simplifies the presentation of subsequent, more complicated
PBMs.

The application of the upwind finite difference scheme to (1) gives

$\frac{f_{i+1}^{j,k}-f_{i}^{j,k}}{\Delta
t}+g_{1}\frac{f_{i}^{j,k}-f_{i}^{j-1,k}}{\Delta
a_{1}}+g_{2}\frac{f_{i}^{j,k}-f_{i}^{j,k-1}}{\Delta a_{2}}=0,$ (2)

where $i$ is the time index, $j$ is the index for $a_{1}$, and $k$ is the
index for $a_{2}$. The discretized PDE can be specified either in natural
variables (such as in the Taylor series expansions) or by employing the index
previously used, i.e.,

$\displaystyle\frac{f_{i+1}^{j,k}-f_{i}^{j,k}}{\Delta
t}+g_{1}\frac{f_{i}^{j,k}-f_{i}^{j-1,k}}{\Delta
a_{1}}+g_{2}\frac{f_{i}^{j,k}-f_{i}^{j,k-1}}{\Delta a_{2}}\equiv$
$\displaystyle\,\frac{f(t+\Delta t,a_{1},a_{2})-f(t,a_{1},a_{2})}{\Delta t}$
$\displaystyle+g_{1}\frac{f(t,a_{1},a_{2})-f(t,a_{1}-\Delta
a_{1},a_{2})}{\Delta a_{1}}$
$\displaystyle+g_{2}\frac{f(t,a_{1},a_{2})-f(t,a_{1},a_{2}-\Delta
a_{2})}{\Delta a_{2}}.$ (3)

Defining $\alpha=\frac{g_{1}\Delta t}{\Delta a_{1}}$ and
$\beta=\frac{g_{2}\Delta t}{\Delta a_{2}}$ for compactness, the corresponding
upwind scheme is

$f_{i+1}^{j,k}=f_{i}^{j,k}-\alpha(f_{i}^{j,k}-f_{i}^{j-1,k})-\beta(f_{i}^{j,k}-f_{i}^{j,k-1}).$
(4)

To characterize the local truncation error of the scheme, consider the Taylor
series expansions

$\displaystyle f(t+\Delta t,a_{1},a_{2})$
$\displaystyle=\sum_{n=0}^{\infty}\frac{(\Delta
t)^{n}}{n!}\frac{\partial^{n}f}{\partial t^{n}}\bigg{|}_{t,a_{1},a_{2}},$
$\displaystyle f(t,a_{1}-\Delta a_{1},a_{2})$
$\displaystyle=\sum_{n=0}^{\infty}(-1)^{n}\frac{(\Delta
a_{1})^{n}}{n!}\frac{\partial^{n}f}{\partial
a_{1}^{n}}\bigg{|}_{t,a_{1},a_{2}},$ $\displaystyle f(t,a_{1},a_{2}-\Delta
a_{2})$ $\displaystyle=\sum_{n=0}^{\infty}(-1)^{n}\frac{(\Delta
a_{2})^{n}}{n!}\frac{\partial^{n}f}{\partial
a_{2}^{n}}\bigg{|}_{t,a_{1},a_{2}}.$ (5)

Equation 1 implies that the higher order derivatives are related by

$\displaystyle\frac{\partial^{n}f}{\partial t^{n}}$
$\displaystyle=(-1)^{n}\sum_{p=0}^{n}\frac{n!}{p!(n-p)!}g_{1}^{n-p}g_{2}^{p}\frac{\partial^{n}f}{\partial
a_{1}^{n-p}\partial a_{2}^{p}}$
$\displaystyle=(-1)^{n}\\!\left(g_{1}^{n}\frac{\partial^{n}f}{\partial
a_{1}^{n}}+g_{2}^{n}\frac{\partial^{n}f}{\partial
a_{2}^{n}}+\sum_{p=1}^{n-1}\frac{n!}{p!(n-p)!}g_{1}^{n-p}g_{2}^{p}\frac{\partial^{n}f}{\partial
a_{1}^{n-p}\partial a_{2}^{p}}\right)$ (6)

Correspondingly, the local truncation error is

Error $\displaystyle=\\!\bigg{[}\frac{\partial f}{\partial
t}+g_{1}\frac{\partial f}{\partial a_{1}}+g_{2}\frac{\partial f}{\partial
a_{2}}\bigg{|}_{t,a_{1},a_{2}}-\\!\left(\frac{f_{i+1}^{j,k}-f_{i}^{j,k}}{\Delta
t}+g_{1}\frac{f_{i}^{j,k}-f_{i}^{j-1,k}}{\Delta
a_{1}}+g_{2}\frac{f_{i}^{j,k}-f_{i}^{j,k-1}}{\Delta a_{2}}\right)$
$\displaystyle=\sum_{n=2}^{\infty}\frac{1}{n!}\bigg{[}g_{1}(-1)^{n}(\Delta
a_{1})^{n-1}\frac{\partial^{n}f}{\partial a_{1}^{n}}+g_{2}(-1)^{n}(\Delta
a_{2})^{n-1}\frac{\partial^{n}f}{\partial a_{2}^{n}}-(\Delta
t)^{n-1}\frac{\partial^{n}f}{\partial t^{n}}\bigg{|}_{t,a_{1},a_{2}}$
$\displaystyle=\sum_{n=2}^{\infty}\bigg{(}\frac{(-1)^{n}}{n!}\bigg{[}(g_{1}(\Delta
a_{1})^{n-1}-g_{1}^{n}(\Delta t)^{n-1})\frac{\partial^{n}f}{\partial
a_{1}^{n}}+(g_{2}(\Delta a_{2})^{n-1}-g_{2}^{n}(\Delta
t)^{n-1})\frac{\partial^{n}f}{\partial a_{2}^{n}}$ $\displaystyle-(\Delta
t)^{n-1}\sum_{p=1}^{n-1}\frac{n!}{p!(n-p)!}g_{1}^{n-p}g_{2}^{p}\frac{\partial^{n}f}{\partial
a_{1}^{n-p}\partial a_{2}^{p}}\bigg{|}_{t,a_{1},a_{2}}\bigg{)}$
$\displaystyle=\sum_{n=2}^{\infty}\bigg{(}\frac{(-1)^{n}}{n!}\bigg{[}\frac{1}{g_{1}^{n}}\bigg{(}\frac{(\Delta
a_{1})^{n-1}}{g_{1}^{n-1}}-(\Delta
t)^{n-1}\bigg{)}\frac{\partial^{n}f}{\partial
a_{1}^{n}}+\frac{1}{g_{2}^{n}}\bigg{(}\frac{(\Delta
a_{2})^{n-1}}{g_{2}^{n-1}}-(\Delta
t)^{n-1}\bigg{)}\frac{\partial^{n}f}{\partial a_{2}^{n}}$
$\displaystyle-(\Delta
t)^{n-1}\sum_{p=1}^{n-1}\frac{n!}{p!(n-p)!}g_{1}^{n-p}g_{2}^{p}\frac{\partial^{n}f}{\partial
a_{1}^{n-p}\partial a_{2}^{p}}\bigg{|}_{t,a_{1},a_{2}}\bigg{)}$
$\displaystyle=\sum_{n=2}^{\infty}\bigg{(}\frac{(-1)^{n}(\Delta
t)^{n-1}}{n!}\bigg{[}\frac{1}{g_{1}^{n}}\left(\frac{1}{\alpha^{n-1}}-1\right)\frac{\partial^{n}f}{\partial
a_{1}^{n}}+\frac{1}{g_{2}^{n}}\left(\frac{1}{\beta^{n-1}}-1\right)\frac{\partial^{n}f}{\partial
a_{2}^{n}}$
$\displaystyle-\sum_{p=1}^{n-1}\frac{n!}{p!(n-p)!}g_{1}^{n-p}g_{2}^{p}\frac{\partial^{n}f}{\partial
a_{1}^{n-p}\partial a_{2}^{p}}\bigg{|}_{t,a_{1},a_{2}}\bigg{)}$ (7)

It can be shown using von Neumann stability analysis that the explicit upwind
scheme is conditionally stable for

$\alpha,\beta\geq 0,\quad\alpha+\beta\leq 1.$ (8)

The above upwind scheme is not able to solve the PBM exactly irrespective of
the values of $\alpha$ and $\beta$. The selection of $\alpha$ and $\beta$ to
satisfy the second inequality in (8) results in excessive numerical diffusion
[25].

An alternative approach is to employ operator splitting (also called
fractional step methods) [25, 31] which enables the original PBM (1) to be
expressed in a manner more amenable to numerical solution. Furthermore,
operator splitting incurs no additional error penalty when the operators
commute [25]. To illustrate this point, consider the application of first-
order order splitting (also known as dimensional splitting in this instance)
to (1),

$\displaystyle\frac{\partial f^{*}}{\partial t}+g_{1}\frac{\partial
f^{*}}{\partial a_{1}}=0,\quad f^{*}(t,a_{1},a_{2})=f(t,a_{1},a_{2}),\quad
t\in[t,t+\Delta t],$ $\displaystyle\frac{\partial f^{**}}{\partial
t}+g_{2}\frac{\partial f^{**}}{\partial a_{2}}=0,\quad
f^{**}(t,a_{1},a_{2})=f^{*}(t+\Delta t,a_{1},a_{2}),\quad t\in[t,t+\Delta t],$
$\displaystyle f(t+\Delta t,a_{1},a_{2})=f^{**}(t+\Delta t,a_{1},a_{2}).$ (9)

To write the subsequent equations more compactly, $\mathcal{A}$ and
$\mathcal{B}$ are used to represent the operators
$g_{1}\frac{\partial}{\partial a_{1}}$ and $g_{2}\frac{\partial}{\partial
a_{2}}$ respectively. An expression for the splitting error for first-order
splitting can be obtained by considering the Taylor series expansion of $f$
about $t+\Delta t$ [25],

$\displaystyle f(t+\Delta t,a_{1},a_{2})$
$\displaystyle=f(t,a_{1},a_{2})+\sum_{n=1}^{\infty}\frac{(\Delta
t)^{n}}{n!}(\mathcal{A}+\mathcal{B})^{n}f(t,a_{1},a_{2})$
$\displaystyle=\sum_{n=0}^{\infty}\frac{(\Delta
t)^{n}}{n!}(\mathcal{A}+\mathcal{B})^{n}f(t,a_{1},a_{2})$
$\displaystyle=e^{\Delta t(\mathcal{A}+\mathcal{B})}f(t,a_{1},a_{2}).$ (10)

For the dimensional splitting method (9),

$\displaystyle f^{**}(t+\Delta t,a_{1},a_{2})$ $\displaystyle=e^{\Delta
t\mathcal{B}}e^{\Delta t\mathcal{A}}f(t,a_{1},a_{2})$
$\displaystyle=\left(\sum_{n=0}^{\infty}\frac{(\Delta
t)^{n}}{n!}\mathcal{B}^{n}\right)\left(\sum_{n=0}^{\infty}\frac{(\Delta
t)^{n}}{n!}\mathcal{A}^{n}\right)f(t,a_{1},a_{2})$
$\displaystyle=\sum_{n=0}^{\infty}\frac{(\Delta
t)^{n}}{n!}\left(\sum_{p=0}^{n}\frac{n!}{p!(n-p)!}\mathcal{B}^{p}\mathcal{A}^{n-p}\right)f(t,a_{1},a_{2}).$
(11)

Therefore the splitting error is

Error $\displaystyle=f(t+\Delta t,a_{1},a_{2})-f^{**}(t+\Delta t,a_{1},a_{2})$
$\displaystyle=\sum_{n=0}^{\infty}\frac{(\Delta
t)^{n}}{n!}\left[(\mathcal{A}+\mathcal{B})^{n}-\left(\sum_{p=0}^{n}\frac{n!}{p!(n-p)!}\mathcal{B}^{p}\mathcal{A}^{n-p}\right)\right]f(t,a_{1},a_{2})$
$\displaystyle=\sum_{n=2}^{\infty}\frac{(\Delta
t)^{n}}{n!}\left[(\mathcal{A}+\mathcal{B})^{n}-\left(\sum_{p=0}^{n}\frac{n!}{p!(n-p)!}\mathcal{B}^{p}\mathcal{A}^{n-p}\right)\right]f(t,a_{1},a_{2})$
$\displaystyle\approx\frac{(\Delta
t)^{2}}{2}(\mathcal{A}\mathcal{B}-\mathcal{B}\mathcal{A})f(t,a_{1},a_{2})+\mathcal{O}((\Delta
t)^{3}).$ (12)

If the operators commute, i.e.,

$\mathcal{A}\mathcal{B}-\mathcal{B}\mathcal{A}=0,$ (13)

then the second-order term is zero. The operators do commute for
$g_{1}\frac{\partial}{\partial a_{1}}$ and $g_{2}\frac{\partial}{\partial
a_{2}}$ for sufficiently smooth $f$,

$g_{1}\frac{\partial}{\partial a_{1}}g_{2}\frac{\partial}{\partial
a_{2}}=g_{1}g_{2}\frac{\partial^{2}}{\partial a_{1}\partial
a_{2}}=g_{2}\frac{\partial}{\partial a_{2}}g_{1}\frac{\partial}{\partial
a_{1}}.$ (14)

More generally, the binomial formula can be used to show that, for
$\mathcal{A}$ and $\mathcal{B}$ that commute, the splitting error (12) is
exactly zero [25]. Correspondingly, the PBM (1) can be solved exactly by
solving each 1D PBM sub-problem in (9) exactly. Each sub-problem can be
exactly solved very efficiently using the upwind finite difference scheme when
$\text{CFL}=1$ is employed as the scheme simplifies to a form that only
requires memory reallocation [15, 25].

As we will subsequently demonstrate, many PBMs have operators that commute
which enables the use of operator splitting techniques for effective solution.
In cases where the operators do not commute, more accurate splitting schemes
such as Strang splitting [32, 33] can be more accurate.

###  2.2 Case 2: PBMs with Growth Rate Dependent on Intrinsic Variables Only

####  2.2.1 PBMs with Growth Rate $G_{i}=G_{i}(a_{i})$

Consider a homogeneous PBM with a growth rate given by
$\mathbf{G}(\mathbf{a})=(G_{1}(a_{1}),G_{2}(a_{2}))$ expressed in conservative
form,

$\frac{\partial f}{\partial t}+\frac{\partial(G_{1}(a_{1})f)}{\partial
a_{1}}+\frac{\partial(G_{2}(a_{2})f)}{\partial a_{2}}=0,\quad
f(0,a_{1},a_{2})=f_{0}(a_{1},a_{2}),$ (15)

with $G_{1}(a_{1})$ and $G_{2}(a_{2})$ continuous in $a_{1}$ and $a_{2}$
respectively, bounded, and $G_{1}(a_{1}),G_{2}(a_{2})>0,\forall a_{1},a_{2}$.
PBMs of this form are relevant when the growth rate depends on the intrinsic
variable such as size-dependent growth in crystallization and precipitation or
age-dependent growth in cell population models (see [15] and citations
therein). Equation 15 can be transformed by multiplying each term in the
equation with $G_{1}(a_{1})G_{2}(a_{2})$ which gives

$\frac{\partial(G_{1}(a_{1})G_{2}(a_{2})f)}{\partial
t}+G_{1}(a_{1})\frac{\partial(G_{1}(a_{1})G_{2}(a_{2})f)}{\partial
a_{1}}+G_{2}(a_{2})\frac{\partial(G_{1}(a_{1})G_{2}(a_{2})f)}{\partial
a_{2}}=0.$ (16)

Defining a new variable $\hat{f}=G_{1}(a_{1})G_{2}(a_{2})f$ enables (16) to be
written as

$\frac{\partial\hat{f}}{\partial
t}+G_{1}(a_{1})\frac{\partial\hat{f}}{\partial
a_{1}}+G_{2}(a_{2})\frac{\partial\hat{f}}{\partial a_{2}}=0.$ (17)

Two variable transformations are introduced,

$\tilde{a}_{1}=\int_{0}^{a_{1}}\frac{1}{G_{1}(a_{1})}da_{1},\quad\tilde{a}_{2}=\int_{0}^{a_{2}}\frac{1}{G_{2}(a_{2})}da_{2},$
(18)

which can be used to transform (17),

$\frac{\partial\hat{f}}{\partial
t}+G_{1}(a_{1})\frac{\partial\hat{f}}{\partial\tilde{a}_{1}}\frac{d\tilde{a}_{1}}{da_{1}}+G_{2}(a_{2})\frac{\partial\hat{f}}{\partial\tilde{a}_{2}}\frac{d\tilde{a}_{2}}{da_{2}}=0,$
(19) $\frac{\partial\hat{f}}{\partial
t}+\frac{\partial\hat{f}}{\partial\tilde{a}_{1}}+\frac{\partial\hat{f}}{\partial\tilde{a}_{2}}=0.$
(20)

Under the assumptions, each function $\tilde{a}_{i}(a_{i})$ is invertible
[15]. Reparametrizing $\hat{f}$ in terms of $\tilde{a}_{i}$,

$\tilde{f}(t,\tilde{a}_{1},\tilde{a}_{2})=\hat{f}(t,a_{1},a_{2}),$ (21)

where $\tilde{f}$ denotes that $\hat{f}$ has been reparametrized. This
reparametrization enables (20) to be written in the form

$\frac{\partial\tilde{f}}{\partial
t}+\frac{\partial\tilde{f}}{\partial\tilde{a}_{1}}+\frac{\partial\tilde{f}}{\partial\tilde{a}_{2}}=0$
(22)

that can be solved exactly as discussed in the previous section. Then all
values for $\tilde{f}(t,\tilde{a}_{1},\tilde{a}_{2})$ are mapped to
$\hat{f}(t,a_{1},a_{2})$ by applying the inverses of $\tilde{a}_{i}(a_{i})$.

To demonstrate the effectiveness of the proposed exact method, multiple
solution strategies are considered alongside the exact method,

  1. 1.

Apply the upwind scheme to (15) on a uniform mesh (“Con-Uniform,Upwind”)

  2. 2.

Apply the WENO scheme to (16) on a uniform mesh (“Trans-Uniform,WENO”)

  3. 3.

Apply the upwind scheme to (15) on a non-uniform mesh to locally enforce CFL =
1 (“Con-Nonuniform,Upwind”)

  4. 4.

Apply the upwind scheme to (16) on a uniform mesh (“Trans-Uniform,Upwind”)

  5. 5.

Apply the upwind scheme to (16) on a non-uniform mesh to locally enforce CFL =
1 (“Trans-Nonuniform,Upwind”)

  6. 6.

Employ the exact method presented in this work (“Exact”)

When the uniform meshes are employed in strategies 1 and 4, the value of
$\Delta t$ is computed to enforce CFL $\leq 1$ which ensures numerical
stability,

$\Delta t=\frac{1.0}{\frac{\max G_{1}(a_{1})}{\Delta a_{1}}+\frac{\max
G_{2}(a_{2})}{\Delta a_{2}}}.$ (23)

To construct the nonuniform mesh for 2D problems for strategies 3 and 5, the
mesh is constructed backwards from its end using a specified $\Delta t$,

$a_{1}^{j-1}=a_{1}^{j}-\gamma G_{1}(a_{1}^{j})\Delta t,\quad
a_{2}^{k-1}=a_{2}^{k}-(1-\gamma)G_{2}(a_{2}^{k})\Delta t,\quad 0\leq\gamma\leq
1,$ (24)

where $\gamma$ is a free parameter which controls the mesh spacing in a
particular direction, with $\gamma=0.5$ used in the case studies. The
implementation of the exact method requires computing (18) and their inverses.
These steps can be done offline either analytically or numerically as
discussed in previous work [15].

####  2.2.2 PBMs with Growth Rate $G_{i}=G_{i}(\mathbf{a})$

Consider a homogeneous PBM with growth rate
$\mathbf{G}(\mathbf{a})=(G_{1}(a_{1},a_{2}),G_{2}(a_{1},a_{2}))$ expressed in
conservative form,

$\frac{\partial f}{\partial t}+\frac{\partial(G_{1}(a_{1},a_{2})f)}{\partial
a_{1}}+\frac{\partial(G_{2}(a_{1},a_{2})f)}{\partial a_{2}}=0.$ (25)

While it is not possible in general to solve this class of PBMs exactly,
employing the aforementioned strategies can improve the numerical solution.
Applying first-order splitting to (25) gives

$\displaystyle\frac{\partial f^{*}}{\partial
t}+\frac{\partial(G_{1}(a_{1},a_{2})f^{*})}{\partial a_{1}}=0,\quad
f^{*}(t,a_{1},a_{2})=f(t,a_{1},a_{2}),$ $\displaystyle\frac{\partial
f^{**}}{\partial t}+\frac{\partial(G_{2}(a_{1},a_{2})f^{**})}{\partial
a_{2}}=0,\quad f^{**}(t,a_{1},a_{2})=f^{*}(t+\Delta t,a_{1},a_{2}),$
$\displaystyle f(t+\Delta t,a_{1},a_{2})=f^{**}(t+\Delta t,a_{1},a_{2}).$ (26)

As the operators do not commute, this splitting incurs an error proportional
to $(\Delta t)^{2}$. Although it might be advantageous to use a higher order
splitting scheme, they are not considered in the present work as the intention
is to demonstrate how comparatively simple techniques can be employed to
effectively solve these PBMs. Defining the variable transformation
$\hat{f}=G_{i}(a_{1},a_{2})f$ for each sub-problem in (26) results in

$\displaystyle\frac{\partial\hat{f}^{*}}{\partial
t}+G_{1}(a_{1},a_{2})\frac{\partial\hat{f}^{*}}{\partial
a_{1}}=0,\quad\hat{f}^{*}(t,a_{1},a_{2})=G_{1}(a_{1},a_{2})f(t,a_{1},a_{2}),$
$\displaystyle\frac{\partial\hat{f}^{**}}{\partial
t}+G_{2}(a_{1},a_{2})\frac{\partial\hat{f}^{**}}{\partial
a_{2}}=0,\quad\hat{f}^{**}(t,a_{1},a_{2})=\hat{f}^{*}(t+\Delta
t,a_{1},a_{2})\frac{G_{2}(a_{1},a_{2})}{G_{1}(a_{1},a_{2})},$ $\displaystyle
f(t+\Delta t,a_{1},a_{2})=\frac{\hat{f}^{**}(t+\Delta
t,a_{1},a_{2})}{G_{2}(a_{1},a_{2})}.$ (27)

Each sub-problem in (27) can be solved exactly by defining the variable
transformation,

$\tilde{a}_{i}=\int_{0}^{a_{i}}\frac{1}{G_{i}(a_{1},a_{2})}\partial a_{i},$
(28)

which results in

$\displaystyle\frac{\partial\tilde{f}^{*}}{\partial
t}+\frac{\partial\tilde{f}^{*}}{\partial\tilde{a}_{1}}=0,\quad\tilde{f}^{*}(t,\tilde{a}_{1},a_{2})=G_{1}(a_{1},a_{2})f(t,a_{1},a_{2}),$
$\displaystyle\frac{\partial\tilde{f}^{**}}{\partial
t}+\frac{\partial\tilde{f}^{**}}{\partial\tilde{a}_{2}}=0,\quad\tilde{f}^{**}(t,a_{1},\tilde{a}_{2})=f^{*}(t+\Delta
t,a_{1},a_{2})\frac{G_{2}(a_{1},a_{2})}{G_{1}(a_{1},a_{2})},$ $\displaystyle
f(t+\Delta t,a_{1},a_{2})=\frac{\hat{f}^{**}(t+\Delta
t,a_{1},a_{2})}{G_{2}(a_{1},a_{2})},$ (29)

where the $\tilde{f}$ indicates that $f$ has been reparameterized in terms of
the transformed variables $\tilde{a}_{i}$.

Five solution strategies are considered alongside the proposed “exact” scheme,

  1. 1.

Apply the upwind scheme to (25) on a uniform mesh (“Con-Uniform,Upwind)

  2. 2.

Apply the WENO scheme to (25) on a uniform mesh (“Expanded-Uniform,WENO”)

  3. 3.

Apply the upwind scheme to (26) on a uniform mesh (“Split-Con-Uniform,Upwind”)

  4. 4.

Apply the upwind scheme to (27) on a uniform mesh (“Split-Trans-
Uniform,Upwind”)

  5. 5.

Apply the upwind scheme to (27) on a non-uniform mesh to locally enforce CFL =
1 for each sub-problem (“Split-Trans-Nonuniform,Upwind)

  6. 6.

Apply the scheme developed in this section (“Split-Exact”).

For the “Con-Uniform,Upwind” scheme, the time step $\Delta t$ is evaluated to
enforce CFL $\leq 1$ for numerical stability,

$\Delta t=\frac{1}{\frac{\max G_{1}(a_{1},a_{2})}{\Delta a_{1}}+\frac{\max
G_{2}(a_{1},a_{2})}{\Delta a_{2}}}.$ (30)

In comparison, the methods using operator splitting on a uniform mesh permit a
larger $\Delta t$,

$\Delta t=\min\\!\left\\{\frac{\Delta a_{1}}{\max
G_{1}(a_{1},a_{2})},\frac{\Delta a_{2}}{\max G_{2}(a_{1},a_{2})}\right\\}.$
(31)

The implementation of the “Split-Trans-Nonuniform,Upwind” scheme is
comparatively involved. One challenge is the mesh construction for each sub-
problem in (27) as each sub-problem solved on its own non-uniform mesh. To
outline the mesh construction process for this scheme, suppose that the
problem domain is specified to be $a_{1}\in[0,L_{1}]$, $a_{2}\in[0,L_{2}]$,
where $L_{1},L_{2}$ are positive constants. For the mesh for the first sub-
problem (in $a_{1}$), we first evaluate grid points in the $a_{2}$ direction
of the mesh at $a_{1}=L_{1}$ as these points are then used as the basis for
stepping backwards to generate the rest of the mesh using the formula,

$a_{1}^{j-1,k}=a_{1}^{j,k}-\Delta t\,G_{1}(a_{1}^{j,k},a_{2}^{k}).$ (32)

This task is underspecified and can be implemented in multiple ways such as
using a linearly spaced array of $a_{2}\in[0,L_{2}]$ at $a_{1}=L_{1}$. We have
elected to use the other growth rate $G_{2}(a_{1},a_{2})$ to construct this
array using the formula

$a_{2}^{k-1}=a_{2}^{k}-\Delta t\,G_{2}(L_{1},a_{2}^{k}).$ (33)

The mesh for the second sub-problem (in $a_{2}$) is similarly created. These
steps can and should be done offline. The mesh construction process for this
scheme involves the formation of jagged arrays (an array of arrays in which
member arrays can have different number of elements) and requires due care. An
example of such a mesh created as part of the scheme can be seen in Figure 1.
The Awkward Array Python library [34] can be used to handle these jagged
arrays.

Both the “Split-Trans-Nonuniform,Upwind” and “Split-Exact” schemes require
interpolation of the solution from one mesh to the other to solve each sub-
problem at each time step. The interpolation step can become computationally
costly and potentially constrain accuracy the solution. The SciPy Python
Library [35] was used to perform the interpolation and Delaunay computation
was performed offline to speed up the interpolation step. Another potential
issue is failure of the interpolation step which can result in the population
of “Not a Number” (NaN) values in the solution. This issue is mitigated
through the use of a low-order interpolation scheme.

Figure 1: Exemplar jagged mesh generated as part of the “Split-Trans-
Nonuniform,Upwind” scheme.

###  2.3 Case 3: PBMs with Time-Dependent Growth Rates

####  2.3.1 PBMs with Growth Rate $G_{i}=G_{i}(t)$

Consider a homogeneous PBM with $\mathbf{G}(t)=(G_{1}(t),G_{2}(t))$ expressed
in conservative form,

$\frac{\partial f}{\partial t}+\frac{\partial(G_{1}(t)f)}{\partial
a_{1}}+\frac{\partial(G_{2}(t)f)}{\partial a_{2}}=0,\quad
f(0,a_{1},a_{2})=f_{0}(a_{1},a_{2}).$ (34)

Bringing $G_{i}(t)$ out of the “spatial” derivatives gives

$\frac{\partial f}{\partial t}+G_{1}(t)\frac{\partial f}{\partial
a_{1}}+G_{2}(t)\frac{\partial f}{\partial a_{2}}=0.$ (35)

The operators $G_{1}(t)\frac{\partial}{\partial a_{1}}$ and
$G_{2}(t)\frac{\partial}{\partial a_{2}}$ commute for sufficiently smooth $f$,

$G_{1}(t)\frac{\partial}{\partial a_{1}}G_{2}(t)\frac{\partial}{\partial
a_{2}}=G_{1}(t)G_{2}(t)\frac{\partial^{2}}{\partial a_{1}\partial
a_{2}}=G_{2}(t)\frac{\partial}{\partial a_{2}}G_{1}(t)\frac{\partial}{\partial
a_{1}}.$ (36)

Therefore, first-order splitting with no splitting error can be applied which
gives

$\displaystyle\frac{\partial f^{*}}{\partial t}+G_{1}(t)\frac{\partial
f^{*}}{\partial a_{1}}=0,\quad f^{*}(t,a_{1},a_{2})=f(t,a_{1},a_{2}),\quad
t\in[t,t+\Delta t],$ $\displaystyle\frac{\partial f^{**}}{\partial
t}+G_{2}(t)\frac{\partial f^{**}}{\partial a_{2}}=0,\quad
f^{**}(t,a_{1},a_{2})=f^{*}(t+\Delta t,a_{1},a_{2}),\quad t\in[t,t+\Delta t],$
$\displaystyle f(t+\Delta t,a_{1},a_{2})=f^{**}(t+\Delta t,a_{1},a_{2}).$ (37)

As shown in previous work [15], each sub-problem can be solved exactly by
employing a variable transformation for $t$,

$\tilde{t}_{i}=\int_{0}^{t}G_{i}(t^{\prime})dt^{\prime}.$ (38)

Introducing this variable transformation for $t$ transforms (37) into

$\displaystyle\frac{\partial\tilde{f}^{*}}{\partial\tilde{t}_{1}}+\frac{\partial\tilde{f}^{*}}{\partial
a_{1}}=0,\quad\tilde{f}^{*}(\tilde{t}_{1},a_{1},a_{2})=f(t,a_{1},a_{2})$
$\displaystyle\frac{\partial\tilde{f}^{**}}{\partial\tilde{t}_{2}}+\frac{\partial\tilde{f}^{**}}{\partial\tilde{a}_{2}}=0,\quad\tilde{f}^{**}(\tilde{t}_{2},a_{1},a_{2})=\tilde{f}^{*}(\tilde{t}_{1}+\Delta\tilde{t}_{1},a_{1},a_{2})$
$\displaystyle f(t+\Delta
t,a_{1},a_{2})=\tilde{f}^{**}(\tilde{t}_{2}+\Delta\tilde{t}_{2},a_{1},a_{2}),$
(39)

where $\tilde{f}$ denotes that $f$ has been reparametrized in terms of
$\tilde{t}_{i}$. Since each sub-problem can be solved exactly and there is no
error arising from operator splitting, this strategy solves the PBM exactly.

####  2.3.2 PBMs with Growth Rate $G_{i}=G_{i,t}(t)G_{i,a_{i}}(a_{i})$

Consider a homogeneous PBM with a separable time- and size-dependent growth
rate $\mathbf{G}(t,\mathbf{a})=(G_{1,t}(t)G_{1,a_{1}}(a_{1}),\allowbreak
G_{2,t}(t)G_{2,a_{2}}(a_{2}))$ in conservative form,

$\frac{\partial f}{\partial
t}+\frac{\partial(G_{1,t}(t)G_{1,a_{1}}(a_{1})f)}{\partial
a_{1}}+\frac{\partial(G_{2,t}(t)G_{2,a_{2}}(a_{2})f)}{\partial a_{2}}=0,\quad
f(0,a_{1},a_{2})=f_{0}(a_{1},a_{2}).$ (40)

Many physical systems can be described by such a separable growth rate as the
growth rate can typically be split into an “environmental” part which is a
function of $t$ and a part that is only a function of the intrinsic variable
[36]. This model can be solved exactly by employing the variable
transformation and splitting strategies explored in previous cases. Factoring
$G_{i,t}(t)$ from the spatial derivatives and multiplying each term in (40) by
$G_{1,a_{1}}(a_{1})G_{2,a_{2}}(a_{2})$ transforms (40) into

$\displaystyle\frac{\partial(G_{1,a_{1}}(a_{1})G_{2,a_{2}}(a_{2})f)}{\partial
t}$
$\displaystyle+G_{1,t}(t)G_{1,a_{1}}(a_{1})\frac{\partial(G_{1,a_{1}}(a_{1})G_{2,a_{2}}(a_{2})f)}{\partial
a_{1}}$
$\displaystyle+G_{2,t}(t)G_{2,a_{2}}(a_{2})\frac{\partial(G_{1,a_{1}}(a_{1})G_{2,a_{2}}(a_{2})f)}{\partial
a_{2}}=0.$ (41)

Defining a new variable $\hat{f}=G_{1,a_{1}}(a_{1})G_{2,a_{2}}(a_{2})f$ and
introducing the variable transforms for $a_{i}$, i.e.,
$\tilde{a}_{i}=\int_{0}^{a_{i}}\frac{1}{G_{i,a_{i}}(a_{i})}da_{i}$ transforms
(41) to

$\frac{\partial\hat{f}}{\partial
t}+G_{1}(t)\frac{\partial\hat{f}}{\partial\tilde{a}_{1}}+G_{2}(t)\frac{\partial\hat{f}}{\partial\tilde{a}_{2}}=0.$
(42)

Reparametrizing $\hat{f}$ in terms of $\tilde{a}_{i}$,

$\tilde{f}(t,\tilde{a}_{1},\tilde{a}_{2})=\hat{f}(t,a_{1},a_{2}),$ (43)

where $\tilde{f}$ denotes that $\hat{f}$ has been reparametrized. This
reparameterization enables (42) to be expressed in a form that can be exactly
solved as shown in the previous case,

$\frac{\partial\tilde{f}}{\partial
t}+G_{1,t}(t)\frac{\partial\tilde{f}}{\partial\tilde{a}_{1}}+G_{2,t}(t)\frac{\partial\tilde{f}}{\partial\tilde{a}_{2}}=0.$
(44)

###  2.4 Case 4: Nonhomogeneous PBMs

####  2.4.1 Nonhomogeneous PBMs with Constant Growth Rates

Consider a PBM with a growth rate $\mathbf{G}=(g_{1},g_{2})$, where $g_{1}$
and $g_{2}$ are positive constants, and a nonhomogeneous term
$h(t,a_{1},a_{2})$,

$\frac{\partial f}{\partial t}+g_{1}\frac{\partial f}{\partial
a_{1}}+g_{2}\frac{\partial f}{\partial a_{2}}=h(t,a_{1},a_{2}),\quad
f(0,a_{1},a_{2})=f_{0}(a_{1},a_{2}).$ (45)

It is not possible to generally transform (45) into a form that can be solved
exactly. However, significant improvement can be achieved by employing
operator splitting and solving each sub-problem while enforcing CFL = 1.
Applying first-order splitting gives

$\displaystyle\frac{\partial f^{*}}{\partial t}+g_{1}\frac{\partial
f^{*}}{\partial a_{1}}=0,\quad f^{*}(t,a_{1},a_{2})=f(t,a_{1},a_{2}),\quad
t\in[t,t+\Delta t]$ $\displaystyle\frac{\partial f^{**}}{\partial
t}+g_{2}\frac{\partial f^{**}}{\partial a_{2}}=0,\quad
f^{**}(t,a_{1},a_{2})=f^{*}(t+\Delta t,a_{1},a_{2}),\quad t\in[t,t+\Delta t]$
$\displaystyle\frac{\partial f^{***}}{\partial t}=h(t,a_{1},a_{2}),\quad
f^{***}(t,a_{1},a_{2})=f^{**}(t+\Delta t,a_{1},a_{2}),\quad t\in[t,t+\Delta
t]$ $\displaystyle f(t+\Delta t,a_{1},a_{2})=f^{***}(t+\Delta t,a_{1},a_{2}).$
(46)

The first two sub-problem can be solved exactly, while the ;ast sub-problem
can be solved efficiently using the forward Euler time stepping scheme. Many
of the cases previously considered can be transformed into a form comparable
to (46). To demonstrate this, consider a nonhomogeneous PBM with a growth rate
$\mathbf{G}=(G_{1,t}(t)G_{1,a_{1}\\!}(a_{1}),G_{2,t}(t)G_{2,a_{2}\\!}(a_{2}))$,

$\frac{\partial f}{\partial
t}+\frac{\partial(G_{1,t}(t)G_{1,a_{1}\\!}(a_{1})f)}{\partial
a_{1}}+\frac{\partial(G_{2,t}(t)G_{2,a_{2}\\!}(a_{2})f)}{\partial
a_{2}}=b(t,a_{1},a_{2}),\quad f(0,a_{1},a_{2})=f_{0}(a_{1},a_{2}).$ (47)

Employing the steps developed in Section 2.3.2 transforms (47) into

$\frac{\partial\tilde{f}}{\partial
t}+G_{1,t}(t)\frac{\partial\tilde{f}}{\partial\tilde{a}_{1}}+G_{2,t}(t)\frac{\partial\tilde{f}}{\partial\tilde{a}_{2}}=\tilde{G}_{1,\tilde{a}_{1}\\!}(\tilde{a}_{1})\tilde{G}_{2,\tilde{a}_{2}\\!}(\tilde{a}_{2})\tilde{b}(t,\tilde{a}_{1},\tilde{a}_{2}),$
(48)

where the superscript $\sim$ over the various functions denotes the
reparameterization in terms of the transformed variables $\tilde{a}_{i}$.
Recognizing that
$\tilde{h}(t,\tilde{a}_{1},\tilde{a}_{2})\equiv\tilde{G}_{1,\tilde{a}_{1}\\!}(\tilde{a}_{1})\tilde{G}_{2,\tilde{a}_{2}\\!}(\tilde{a}_{2})\tilde{b}(t,\tilde{a}_{1},\tilde{a}_{2})$,
and applying first-order splitting with a variable transformation for $t$ to
(48) as in Section 2.3.1 yields

$\displaystyle\frac{\partial\tilde{f}^{*}}{\partial\tilde{t}_{1}}+\frac{\partial\tilde{f}^{*}}{\partial\tilde{a}_{1}}=0,\quad\tilde{f}^{*}(\tilde{t}_{1},\tilde{a}_{1},\tilde{a}_{2})=\tilde{f}(\tilde{t},\tilde{a}_{1},\tilde{a}_{2}),\quad\tilde{t}_{1}\in[\tilde{t}_{1},\tilde{t}_{1}+\Delta\tilde{t}_{1}]$
$\displaystyle\frac{\partial\tilde{f}^{**}}{\partial\tilde{t}_{2}}+\frac{\partial\tilde{f}^{**}}{\partial\tilde{a}_{2}}=\tilde{h}(t,\tilde{a}_{1},\tilde{a}_{2}),\quad
f^{**}(\tilde{t}_{2},\tilde{a}_{1},\tilde{a}_{2})=f^{*}(\tilde{t}_{1}+\Delta\tilde{t}_{1},\tilde{a}_{1},\tilde{a}_{2}),\quad\tilde{t}_{2}\in[\tilde{t}_{2},\tilde{t}_{2}+\Delta\tilde{t}_{2}]$
$\displaystyle f(t+\Delta t,a_{1},a_{2})=f^{**}(t+\Delta t,a_{1},a_{2}).$ (49)

###  2.5 PBMs with a Linear Nonhomogeneous Term

Consider a PBM with a constant growth rate $\mathbf{G}=(g_{1},g_{2})$ and a
linear nonhomogeneous term,

$\frac{\partial f}{\partial t}+g_{1}\frac{\partial f}{\partial
a_{1}}+g_{2}\frac{\partial f}{\partial a_{2}}=-\lambda(t,a_{1},a_{2})f,\quad
f(0,a_{1},a_{2})=f_{0}(a_{1},a_{2}).$ (50)

By employing a variable transform $\hat{f}=\mu f$, where the functional form
of $\mu$ is to be determined, (50) can be transformed into a form that can be
solved exactly,

$\frac{\partial\hat{f}}{\partial t}+g_{1}\frac{\partial\hat{f}}{\partial
a_{1}}+g_{2}\frac{\partial\hat{f}}{\partial a_{2}}=0.$ (51)

This equation can be expanded to yield a PDE for $\mu$,

$\frac{\partial\mu}{\partial t}+g_{1}\frac{\partial\mu}{\partial
a_{1}}+g_{2}\frac{\partial\mu}{\partial a_{2}}=\lambda\mu.$ (52)

An expression for $\mu$ can be derived from the solution of (52).

Table 1: Summary of functional forms of $\mu$ for different forms of $\lambda$ Form of $\lambda$ | Form of $\mu$ | Functional form of $\mu$
---|---|---
Constant | $\mu(a_{1})\lor\mu(a_{2})\lor\mu(t)$ | $e^{\frac{\lambda a_{1}}{g_{1}}}\lor e^{\frac{\lambda a_{2}}{g_{2}}}\lor e^{\lambda t}$
$\lambda(a_{1})$ | $\mu(a_{1})$ | $\exp\\!\left(\frac{1}{g_{1}}\int_{0}^{a_{1}}\lambda(a_{1}^{\prime})da_{1}^{\prime}\right)$
$\lambda(a_{2})$ | $\mu(a_{2})$ | $\exp\\!\left(\frac{1}{g_{2}}\int_{0}^{a_{2}}\lambda(a_{2}^{\prime})da_{2}^{\prime}\right)$
$\lambda(t)$ | $\mu(t)$ | $\exp\\!{\big{(}\\!\int_{0}^{t}\lambda(t^{\prime})dt^{\prime}\big{)}}$
$\lambda(t,a_{1})$ | $\mu(t,a_{1})$ |  $\exp\\!{\left(\\!\int_{0}^{t}\lambda(t^{\prime},g_{1}t^{\prime}+a_{1,0})dt^{\prime}\right)}$, $a_{1,0}=a_{1}-g_{1}t$
$\lambda(t,a_{2})$ | $\mu(t,a_{2})$ |  $\exp\\!{\left(\\!\int_{0}^{t}\lambda(t^{\prime},g_{2}t^{\prime}+a_{2,0})dt^{\prime}\right)}$, $a_{2,0}=a_{2}-g_{2}t$
$\lambda(a_{2},a_{2})$ | $\mu(a_{1},a_{2})$ |  $\exp\\!{\left(\\!\int_{0}^{a_{1}}\frac{1}{g_{1}}\lambda(a_{1}^{\prime},\frac{g_{2}}{g_{1}}a_{1}^{\prime}+a_{2,0})da_{1}^{\prime}\right)}$, $a_{2,0}=a_{2}-\frac{g_{2}}{g_{1}}a_{1}$
$\lambda(t,a_{1},a_{2})$ | $\mu(t,a_{1},a_{2})$ |  $\exp\\!{\left(\\!\int_{0}^{t}\lambda(t^{\prime},g_{1}t^{\prime}+a_{1,0},g_{2}t+a_{2,0})dt^{\prime}\right)}$, $a_{1,0}=a_{1}-g_{1}t$, $a_{2,0}=a_{2}-g_{2}t$

##  3 Results and Discussion

This section is structured as a series of cases which correspond to the
various classes of PBMs explored previously. The error of the various
numerical schemes is compared via the Root Mean Square Error (RMSE) and
Maximum Absolute Error (MAE),

$\text{RMSE}=\sqrt{\frac{\sum_{i=1}^{n}(y_{i}-y_{\text{analytical},i})^{2}}{n}},$
(53) $\text{MAE}=\max_{i}|y_{i}-y_{\text{analytical},i}|.$ (54)

###  3.1 Case 1: PBMs with Constant Growth Rates

Consider the PBM,

$\frac{\partial f}{\partial t}+\frac{\partial f}{\partial
a_{1}}+\frac{\partial f}{\partial a_{2}}=0,\quad
f_{0}(a_{1},a_{2})=50\exp\\!\left(-\frac{(a_{1}-0.4)^{2}}{0.005}-\frac{(a_{2}-0.4)^{2}}{0.005}\right),\quad
a_{i}\in[0,2].$ (55)

A no-flux boundary condition is applied at the top and right ends of the
domain (i.e., $a_{1}=a_{2}=2$), and a modified Dirichlet boundary which
enforces the value of $f$ at the ghost nodes are zero [17] is applied on the
left and bottom ends of the domain, i.e., $a_{1}=a_{2}=0$. The PBM has the
analytical solution,

$f(t,a_{1},a_{2})=f_{0}(a_{1}-t,a_{2}-t)=50\exp\\!\left(-\frac{(a_{1}-0.4-t)^{2}}{0.005}-\frac{(a_{2}-0.4-t)^{2}}{0.005}\right).$
(56)

Exemplar simulation results can be found in Figure 2 and the convergence
analysis can be found in Figure 3. An additional numerical scheme
“Exact,Interpolation”, which adapts the “Exact” scheme by replacing the memory
reallocation step with an interpolation function call to compute the updated
values at $t+\Delta t$, was also implemented to demonstrate that the use of an
interpolation step to compute the values of $f$ at $t+\Delta t$ does not
adversely impact error performance.

(a) Analytical Solution

(b) Exact Scheme

(c) Upwind Scheme

(d) WENO Scheme

Figure 2: Simulation results for Case 1 at $t=1.0$ using the various schemes.
101 grid points are used in both the $a_{1}$ and $a_{2}$ directions for the
Upwind and Exact schemes while 100 cells in both the $a_{1}$ and $a_{2}$
directions for the WENO scheme.

(a) RMSE

(b) MAE

Figure 3: Error analysis for Case 1. The use of dimensional splitting and
employing $\text{CFL}=1$ for each subproblem solves the PBM to machine
precision as expected.

###  3.2 Case 2: PBMs with Growth Rate $G_{i}=G_{i}(a_{i})$

Consider the PBM adapted from [37],

$\frac{\partial f}{\partial t}+\frac{\partial(G_{1}(a_{1})f)}{\partial
a_{1}}+\frac{\partial(G_{2}(a_{2})f)}{\partial a_{2}}=0,\quad
f_{0}(a_{1},a_{2})=50\exp\\!\left(-\frac{(a_{1}-0.4)^{2}}{0.005}-\frac{(a_{2}-0.4)^{2}}{0.005}\right),$
(57)

with,

$G_{1}(a_{1})=0.1+0.05a_{1},\quad G_{2}(a_{2})=0.5+0.25a_{2},\quad
a_{i}\in[0,2].$ (58)

Multiplying each term in (57) by $G_{1}(a_{1})G_{2}(a_{2})$ and defining
$\hat{f}(t,a_{1},a_{2})=G_{1}(a_{1})G_{2}(a_{2})f$ transforms (57) into

$\frac{\partial\hat{f}}{\partial
t}+G_{1}(a_{1})\frac{\partial\hat{f}}{\partial
a_{1}}+G_{2}(a_{2})\frac{\partial\hat{f}}{\partial
a_{2}}=0,\quad\hat{f}(0,a_{1},a_{2})=\hat{f}_{0}(a_{1},a_{2})=G_{1}(a_{1})G_{2}(a_{2})f_{0}(a_{1},a_{2}).$
(59)

This PBM has an analytical solution,

$f(t,a_{1},a_{2})=\frac{\hat{f}_{0}\left((a_{1}+2)e^{-0.05t}-2),(a_{2}+2)e^{-0.25t}-2)\right)}{G_{1}(a_{1})G_{2}(a_{2})}$
(60)

Exemplar numerical solutions to the PBM can be found in Figure 4 and
convergence analysis in Figure 5. To employ the “Exact,Analytical” scheme, the
functions $\tilde{a}_{i}(a_{i})$ given by (18) and their inverse need to be
separately computed analytically and supplied into the scheme. Another variant
of the exact scheme, “Exact,Numerical” is also considered where
$\tilde{a}_{i}(a_{i})$ and its inverse are computed numerically instead. While
the accuracy of the “Exact,Numerical” scheme is constrained by the accuracy of
the quadrature step as can be seen in Figure 5, the scheme is still able to
perform very well and is more user-friendly as it does not require any pre-
computation.

(a) Analytical Solution

(b) Con-Uniform,Upwind

(c) Trans-Uniform,Upwind

(d) Trans-Uniform,WENO

(e) Con-Nonuniform,Upwind

(f) Trans-Nonuniform,Upwind

(g) Exact(Analytical)

(h) Exact(Numerical)

Figure 4: Simulation results for Case 2 at $t=1.0$ using the various schemes.
101 grid points are used in both the $a_{1}$ and $a_{2}$ directions for the
Upwind schemes on a uniform grid and Exact schemes while 100 cells in both the
$a_{1}$ and $a_{2}$ directions for the WENO scheme. The simulations on a
nonuniform grid have 277 and 56 grid points in the $a_{1}$ and $a_{2}$
directions respectively.

(a) RMSE

(b) MAE

Figure 5: Error analysis for Case 2.

###  3.3 Case 3: PBMs with Growth Rate $G_{i}=G_{i}(\mathbf{a})$

Consider the PBM,

$\frac{\partial f}{\partial t}+\frac{\partial(G_{1}(a_{1},a_{2})f)}{\partial
a_{1}}+\frac{\partial(G_{2}(a_{1},a_{2})f)}{\partial a_{2}}=0,\quad
f_{0}(a_{1},a_{2})=50\exp\\!\left(-\frac{(a_{1}-0.4)^{2}}{0.005}-\frac{(a_{2}-0.4)^{2}}{0.005}\right),$
(61)

with,

$G_{1}(a_{1},a_{2})=0.25+0.5(a_{1}+a_{2}),\quad
G_{2}(a_{1},a_{2})=0.5+0.25(a_{1}+a_{2}).$ (62)

The analytical solution to this PBM is

$f(t,a_{1},a_{2})=f_{0}\\!\left(\frac{Bc_{2}-Dc_{1}}{BC-
AD},\frac{Cc_{1}-Ac_{2}}{BC-AD}\right)\exp(-0.75t),$ (63)

where

$\displaystyle c_{1}$ $\displaystyle=3a_{1}+0.75t+2-2e^{0.75t},\quad
c_{2}=3a_{2}-0.75t+1-e^{0.75t},$ $\displaystyle A$
$\displaystyle=1+2e^{0.75t},\quad B=-2+2e^{0.75t},\quad C=-1+e^{0.75t},\quad
D=2+e^{0.75t}.$ (64)

(a) RMSE

(b) MAE

(c) “Split-Exact” RMSE

(d) “Split-Exact” MAE

Figure 6: Error Analysis for Case 3.

(a) Analytical Solution

(b) Con-Uniform,Upwind

(c) Expanded-Uniform,WENO

(d) Split-Con-Uniform, Upwind

(e) Split-Trans-Uniform,Upwind

(f) Split-Trans-Nonuniform,Upwind

(g) Split-Exact

Figure 7: Simulation results for Case 3 at $t=1.0$ using the various schemes.
101 grid points are used in both the $a_{1}$ and $a_{2}$ directions for the
Upwind schemes on a uniform grid and Split-Exact schemes while 100 cells in
both the $a_{1}$ and $a_{2}$ directions for the WENO scheme. The Split-Exact
scheme used a value of $\Delta t=0.1$. The simulation on a nonuniform grid has
25,367 grid points by specifying $\Delta t=0.01$.

None of the schemes employed are able to solve the PBM to a high degree of
accuracy (Figure 6). Even though the “Split-Exact” scheme performs poorer than
the WENO scheme (which is regarded as both mathematically involved and
computationally expensive), the “Split-Exact” scheme enables the use of much
larger $\Delta t$ values as the sub-problems are effectively solved as a
single function call to advance the full $\Delta t$ instead of requiring time-
stepping.

The “Split-Exact scheme” using the method of characteristics can be modified
and enhanced to eliminate the need to interpolate between meshes to solve each
sub-problem and further speed up the solution. However, this enhancement
requires the offline analytical evaluation of two integrals and does not
significantly improve the error performance (see Figure A1). This modification
is discussed in the Appendix (see Section A1). The lack of improvement in the
error performance from using the enhanced scheme indicates that the main
source of error with the “Split-Exact” scheme arises from the operator
splitting step rather than the quadrature and interpolation steps.
Correspondingly, it might be worth exploring more sophisticated splitting
schemes to further improve error performance.

###  3.4 Case 4: PBMs with a Nonhomogeneous Term

Consider the PBM,

$\frac{\partial f}{\partial t}+\frac{\partial f}{\partial
a_{1}}+\frac{\partial f}{\partial a_{2}}=1+a_{1}a_{2},\quad
f_{0}(a_{1},a_{2})=10\exp\\!\left(-\frac{(a_{1}-0.4)^{2}}{0.005}-\frac{(a_{2}-0.4)^{2}}{0.005}\right),$
(65)

with the same boundary conditions as Case 1. An analytical solution cannot be
readily supplied for benchmarking so, to perform the convergence analysis, a
“reference” numerical solution was generated by running the simulation using
the presented scheme on a very fine mesh ($\sim 2.56\times 10^{6}$ grid
points).

(a) Reference Solution

(b) Split Scheme

(c) Upwind Scheme

(d) WENO Scheme

Figure 8: Simulation results for Case 4 at $t=1.0$ using the various schemes.
101 grid points are used in both the $a_{1}$ and $a_{2}$ directions for the
Upwind and Split schemes while 100 cells in both the $a_{1}$ and $a_{2}$
directions for the WENO scheme.

(a) RMSE

(b) MAE

Figure 9: Error analysis for Case 4.

The numerical diffusion is much smaller for the Split and WENO than the Upwind
scheme (Figure 8). At low resolution, the proposed Split scheme has much
higher accuracy than the WENO scheme (Figure 9). Employing operator splitting
and solving the “advection” component of the Split problem exactly results in
much less numerical diffusion on coarser meshes compared to WENO scheme.

###  3.5 Case 5: PBMs with a Linear Nonhomogeneous Term

Consider an example of the 2D von Foerster equation,

$\frac{\partial f}{\partial t}+\frac{\partial f}{\partial
a_{1}}+\frac{\partial f}{\partial a_{2}}=-(a_{1}+a_{2})f,\quad
f_{0}(a_{1},a_{2})=50\exp\\!\left(-\frac{(a_{1}-0.4)^{2}}{0.005}-\frac{(a_{2}-0.4)^{2}}{0.005}\right),$
(66)

with the boundary conditions,

$f(t,a_{1}=0,a_{2})=f(t,a_{1},a_{2}=0)=0.$ (67)

This PBM has the analytical solution,

$f(t,a_{1},a_{2})=\begin{cases}f_{0}(a_{1}-t,a_{2}-t)\exp\\!{\left(-(a_{1}+a_{2})t+t^{2}\right)},&a_{1},a_{2}\geq
t,\\\ 0,&a_{1},a_{2}<t.\end{cases}$ (68)

The simulation and error results (see Figures 8 and 11) are qualitatively
similar to Cases 1 and 2, with the proposed exact scheme solving the PBM to
machine precision, and the upwind and WENO schemes having much lower numerical
accuracy.

(a) Analytical Solution

(b) Exact Scheme

(c) Upwind Scheme

(d) WENO Scheme

Figure 10: Simulation results for Case 5 at $t=1.0$ using the various schemes.
101 grid points are used in both the $a_{1}$ and $a_{2}$ directions for the
Upwind and Exact schemes while 100 cells in both the $a_{1}$ and $a_{2}$
directions for the WENO scheme. The ability of the proposed scheme to solve
the PBM to machine precision is demonstrated.

(a) RMSE

(b) MAE

Figure 11: Error analysis for Case 5.

##  4 Conclusions

With a combination of dimensional splitting, variable transformations, and
operating at the limit of numerical stability, the upwind finite difference
scheme is able to solve many classes of multidimensional PBMs either to
machine precision or with sufficiently high accuracy for most applications.
One of the most significant features of the proposed numerical scheme is its
high computational efficiency, requiring only memory reallocation or, in some
instances, a minimal number of floating point operations and function calls
during the time-stepping. Even in Case 3 in Section 3, where the accuracy of
all of the tested numerical schemes is comparatively poor, the proposed scheme
can be more efficient due to the larger time-steps made possible by the scheme
effectively transforming each sub-problem into a function call. The low
computational cost enables the direct incorporation of the multidimensional
population balance model into on-line optimization-based control, aka model
predictive control.

##  5 Availability of Code

A Python implementation of the various schemes discussed in this work is
available at https://github.com/pavaninguva/PBM_Schemes.

## Acknowledgments

Financial support is acknowledged from the U.S. Food and Drug Administration
(75F40121C00090) and the Agency for Science, Technology and Research (A*STAR),
Singapore.

## References

  * [1] Heiko Briesen  “Simulation of crystal size and shape by means of a reduced two-dimensional population balance model”  In _Chemical Engineering Science_ 61.1, 2006, pp. 104–112  DOI: 10.1016/j.ces.2004.11.062
  * [2] Cai Y. Ma, Xue Z. Wang and Kevin J. Roberts  “Multi-dimensional population balance modeling of the growth of rod-like L-glutamic acid crystals using growth rates estimated from in-process imaging”  In _Advanced Powder Technology_ 18.6, 2007, pp. 707–723  DOI: 10.1163/156855207782514932
  * [3] Yang Zhang, Jing Jing Liu, Jian Wan and Xue Z. Wang  “Two dimensional population balance modelling of crystal growth behaviour under the influence of impurities”  In _Advanced Powder Technology_ 26.2, 2015, pp. 672–678  DOI: 10.1016/j.apt.2015.01.001
  * [4] D. Fysikopoulos, B. Benyahia, A. Borsos, Z.K. Nagy and C.D. Rielly  “A framework for model reliability and estimability analysis of crystallization processes with multi-impurity multi-dimensional population balance models”  In _Computers & Chemical Engineering_ 122, 2019, pp. 275–292  DOI: 10.1016/j.compchemeng.2018.09.007
  * [5] Charles David Immanuel and Francis Joseph Doyle  “Solution technique for a multi-dimensional population balance model describing granulation processes”  In _Powder Technology_ 156.2-3, 2005, pp. 213–225  DOI: 10.1016/j.powtec.2005.04.013
  * [6] Dana Barrasso, Arwa El Hagrasy, James D. Litster and Rohit Ramachandran  “Multi-dimensional population balance model development and validation for a twin screw granulation process”  In _Powder Technology_ 270, 2015, pp. 612–621  DOI: 10.1016/j.powtec.2014.06.035
  * [7] Simon M. Iveson  “Limitations of one-dimensional population balance models of wet granulation processes”  In _Powder Technology_ 124.3, 2002, pp. 219–229  DOI: 10.1016/S0032-5910(02)00026-8
  * [8] Nikolaos V. Mantzaris, Prodromos Daoutidis and Friedrich Srienc  “Numerical solution of multi-variable cell population balance models: I. Finite difference methods”  In _Computers & Chemical Engineering_ 25.11-12, 2001, pp. 1411–1440  DOI: 10.1016/S0098-1354(01)00709-8
  * [9] Robert Dürr, Thomas Müller, Stefanie Duvigneau and Achim Kienle  “An efficient approximate moment method for multi-dimensional population balance models – Application to virus replication in multi-cellular systems”  In _Chemical Engineering Science_ 160, 2017, pp. 321–334  DOI: 10.1016/j.ces.2016.11.015
  * [10] V. Quedeville, H. Ouazaite, B. Polizzi, R.O. Fox, P. Villedieu, P. Fede, F. Létisse and J. Morchain  “A two-dimensional population balance model for cell growth including multiple uptake systems”  In _Chemical Engineering Research and Design_ 132, 2018, pp. 966–981  DOI: 10.1016/j.cherd.2018.02.025
  * [11] Jérôme Morchain, Jean-Christophe Gabelle and Arnaud Cockx  “Coupling of biokinetic and population balance models to account for biological heterogeneity in bioreactors”  In _AIChE Journal_ 59.2, 2013, pp. 369–379  DOI: 10.1002/aic.13820
  * [12] Michael L. Rasche, Mo Jiang and Richard D. Braatz  “Mathematical modeling and optimal design of multi-stage slug-flow crystallization”  In _Computers & Chemical Engineering_ 95, 2016, pp. 240–248  DOI: 10.1016/j.compchemeng.2016.09.010
  * [13] Saeed Shirazian, Hamza Y. Ismail, Mehakpreet Singh, Rahamatullah Shaikh, Denise M. Croker and Gavin M. Walker  “Multi-dimensional population balance modelling of pharmaceutical formulations for continuous twin-screw wet granulation: Determination of liquid distribution”  In _International Journal of Pharmaceutics_ 566, 2019, pp. 352–360  DOI: 10.1016/j.ijpharm.2019.06.001
  * [14] Nicholas J. Mozdzierz, Yongkyu Lee, Moo Sun Hong, Moritz H.P. Benisch, Michael L. Rasche, Uku Erik Tropp, Mo Jiang, Allan S. Myerson and Richard D. Braatz  “Mathematical modeling and experimental validation of continuous slug-flow tubular crystallization with ultrasonication-induced nucleation and spatially varying temperature”  In _Chemical Engineering Research and Design_ 169, 2021, pp. 275–287  DOI: 10.1016/j.cherd.2021.03.026
  * [15] Pavan K. Inguva, Kaylee C. Schickel and Richard D. Braatz  “Efficient numerical schemes for population balance models”  In _Computers & Chemical Engineering_ 162, 2022, pp. 107808  DOI: 10.1016/j.compchemeng.2022.107808
  * [16] David L. Ma, Danesh K. Tafti and Richard D. Braatz  “High-Resolution Simulation of Multidimensional Crystal Growth”  In _Industrial & Engineering Chemistry Research_ 41.25, 2002, pp. 6217–6223  DOI: 10.1021/ie010680u
  * [17] Rudiyanto Gunawan, Irene Fusman and Richard D. Braatz  “High resolution algorithms for multidimensional population balance equations”  In _AIChE Journal_ 50.11, 2004, pp. 2738–2749  DOI: 10.1002/aic.10228
  * [18] Shamsul Qamar and Gerald Warnecke  “Solving population balance equations for two-component aggregation by a finite volume scheme”  In _Chemical Engineering Science_ 62.3, 2007, pp. 679–693  DOI: 10.1016/j.ces.2006.10.001
  * [19] Mark A. Pinto, Charles D. Immanuel and Francis J. Doyle  “A feasible solution technique for higher-dimensional population balance models”  In _Computers & Chemical Engineering_ 31.10, 2007, pp. 1242–1256  DOI: 10.1016/j.compchemeng.2006.10.016
  * [20] Mehakpreet Singh, Randhir Singh, Sukhjit Singh, Gavin Walker and Themis Matsoukas  “Discrete finite volume approach for multidimensional agglomeration population balance equation on unstructured grid”  In _Powder Technology_ 376, 2020, pp. 229–240  DOI: 10.1016/j.powtec.2020.08.022
  * [21] Nikolaos V. Mantzaris, Prodromos Daoutidis and Friedrich Srienc  “Numerical solution of multi-variable cell population balance models. III. Finite element methods”  In _Computers & Chemical Engineering_ 25.11-12, 2001, pp. 1463–1481  DOI: 10.1016/S0098-1354(01)00711-6
  * [22] Sashikumaar Ganesan and Lutz Tobiska  “An operator-splitting finite element method for the efficient parallel solution of multidimensional population balance systems”  In _Chemical Engineering Science_ 69.1, 2012, pp. 59–68  DOI: 10.1016/j.ces.2011.09.031
  * [23] Nikolaos V. Mantzaris, Prodromos Daoutidis and Friedrich Srienc  “Numerical solution of multi-variable cell population balance models. II. Spectral methods”  In _Computers & Chemical Engineering_ 25.11-12, 2001, pp. 1441–1462  DOI: 10.1016/S0098-1354(01)00710-4
  * [24] Aniruddha Majumder, Vinay Kariwala, Santosh Ansumali and Arvind Rajendran  “Lattice Boltzmann method for multi-dimensional population balance models in crystallization”  In _Chemical Engineering Science_ 70, 2012, pp. 121–134  DOI: 10.1016/j.ces.2011.04.041
  * [25] Randall J. LeVeque  “Finite Volume Methods for Hyperbolic Problems”  Cambridge,UK: Cambridge University Press, 2002  DOI: 10.1017/CBO9780511791253
  * [26] David I. Ketcheson, Matteo Parsani and Randall J. LeVeque  “High-Order wave propagation algorithms for hyperbolic systems”  In _SIAM Journal on Scientific Computing_ 35.1, 2013, pp. A351–A377  DOI: 10.1137/110830320
  * [27] David I. Ketcheson, Kyle Mandli, Aron J. Ahmadia, Amal Alghamdi, Manuel Quezada Luna, Matteo Parsani, Matthew G. Knepley and Matthew Emmett  “PyClaw: Accessible, extensible, scalable tools for wave propagation problems”  In _SIAM Journal on Scientific Computing_ 34.4, 2012, pp. C210–C231  DOI: 10.1137/110856976
  * [28] Kyle T. Mandli, Aron J. Ahmadia, Marsha Berger, Donna Calhoun, David L. George, Yiannis Hadjimichael, David I. Ketcheson, Grady I. Lemoine and Randall J. LeVeque  “Clawpack: building an open source ecosystem for solving hyperbolic PDEs”  In _PeerJ Computer Science_ 2, 2016, pp. e68  DOI: 10.7717/peerj-cs.68
  * [29] Romuald Szymkiewicz and Dariusz Ga̧siorowski  “Adaptive method for the solution of 1D and 2D advection–diffusion equations used in environmental engineering”  In _Journal of Hydroinformatics_ 23.6, 2021, pp. 1290–1311  DOI: 10.2166/hydro.2021.062
  * [30] B. Seibold  “Numerical Methods for Partial Differential Equations”  In _18.336_ MIT OpenCourseWare, 2009  URL: https://ocw.mit.edu/
  * [31] Rasool Hosseini and Mehdi Tatari  “Some splitting methods for hyperbolic PDEs”  In _Applied Numerical Mathematics_ 146, 2019, pp. 361–378  DOI: 10.1016/j.apnum.2019.07.005
  * [32] Gilbert Strang  “On the construction and comparison of difference schemes”  In _SIAM Journal on Numerical Analysis_ 5.3, 1968, pp. 506–517  DOI: 10.1137/0705041
  * [33] Raymond L. Speth, William H. Green, Shev MacNamara and Gilbert Strang  “Balanced splitting and rebalanced splitting”  In _SIAM Journal on Numerical Analysis_ 51.6, 2013, pp. 3084–3105  DOI: 10.1137/120878641
  * [34] Jim Pivarski, Peter Elmer and David Lange  “Awkward arrays in Python, C++, and Numba”  In _EPJ Web of Conferences_ 245, 2020, pp. 05023  DOI: 10.1051/epjconf/202024505023
  * [35] Pauli Virtanen et al.  “SciPy 1.0: Fundamental algorithms for scientific computing in Python”  In _Nature Methods_ 17.3, 2020, pp. 261–272  DOI: 10.1038/s41592-019-0686-2
  * [36] H.M. Hulburt and S. Katz  “Some problems in particle technology”  In _Chemical Engineering Science_ 19.8, 1964, pp. 555–574  DOI: 10.1016/0009-2509(64)85047-8
  * [37] Rudiyanto Gunawan, Irene Fusman and Richard D. Braatz  “Parallel high-resolution finite volume simulation of particulate processes”  In _AIChE Journal_ 54.6, 2008, pp. 1449–1458  DOI: 10.1002/aic.11484