# Elzaki Transform Based Accelerated Homotopy Perturbation Method for Multi-
dimensional Smoluchowski’s Coagulation and Coupled Coagulation-fragmentation
Equations

Gourav Arora†, Rajesh Kumar† and Youcef Mammeri††

†Department of Mathematics, Birla Institute of Technology and Science, Pilani,
Rajasthan-333031, India (p20190421@pilani.bits-pilani.ac.in).
∗Corresponding author: Department of Mathematics, Birla Institute of
Technology and Science, Pilani, Rajasthan-333031, India
(rajesh.kumar@pilani.bits-pilani.ac.in).
††Institut Camille Jordan CNRS UMR 5208, Université Jean Monnet, 42100 Saint-
Etienne, France (youcef.mammeri@u-picardie.fr).

> Abstract: This article aims to establish a semi-analytical approach based on
> the homotopy perturbation method (HPM) to find the closed form or
> approximated solutions for the population balance equations such as
> Smoluchowski’s coagulation, fragmentation, coupled coagulation-fragmentation
> and bivariate coagulation equations. An accelerated form of the HPM is
> combined with the Elzaki transformation to improve the accuracy and
> efficiency of the method. One of the significant advantages of the technique
> lies over the classic numerical methods as it allows solving the linear and
> non-linear differential equations without discretization. Further, it has
> benefits over the existing semi-analytical techniques such as Adomian
> decomposition method (ADM), optimized decomposition method (ODM), and
> homotopy analysis method (HAM) in the sense that computation of Adomian
> polynomials and convergence parameters are not required. The novelty of the
> scheme is shown by comparing the numerical findings with the existing
> results obtained via ADM, HPM, HAM and ODM for non-linear coagulation
> equation. This motivates us to extend the scheme for solving the other
> models mentioned above. The supremacy of the proposed scheme is demonstrated
> by taking several numerical examples for each problem. The error between
> exact and series solutions provided in graphs and tables show the accuracy
> and applicability of the method. In addition to this, convergence of the
> series solution is also the key attraction of the work.

Keywords: Population Balance Equation; Aggregation Equation; Semi-analytical
Technique; Elzaki Transformation; Accelerated Homotopy Perturbation Method;
Series Solution; Convergence Analysis.

##  1\. Introduction

Particulate processes have drawn much attention of researchers because of
their technological applications in many engineering and natural science
disciplines, including granulation, crystallization, activated sludge
flocculation, and raindrop generation [1, 2, 3, 4, 5]. The particle size
distribution, which represents the amount of a specific size within the
system, affects the behavior during processing and the final product’s
performance. During processing, distinct mechanisms like nucleation, breakage
(fragmentation), aggregation (coagulation) or growth may occur. Breakage
refers the phenomenon in which a particle divides into two or more particles,
while aggregation refers to two particles merging to form a more extensive
particle. Thus, the total number of particles increases during the breakage
process, whereas it decreases in the aggregation phenomenon as time passes,
but the mass remains conserved in both the situations. The scope of the
article is limited to the pure breakage, aggregation in single and multi-
dimensions as well as coupled aggregation-breakage equations. The mathematical
formulation of pure fragmentation equation [6] is given by

$\displaystyle\dfrac{\partial u(x,t)}{\partial
t}=\int_{x}^{\infty}B(x,y)S(y)u(y,t)dy-S(x)u(x,t),$ (1.1)

and the non-linear Smoluchowski’s coagulation equation in $1-D$ is provided
by, see [7],

$\displaystyle\frac{\partial u(x,t)}{\partial
t}=\frac{1}{2}\int_{0}^{x}K(x-y,y)u(x-y,t)u(y,t)dy-\int_{0}^{\infty}K(x,y)u(x,t)u(y,t)dy,$
(1.2)

with the initial condition

$\displaystyle u(x,0)=f(x).$ (1.3)

Here, $u(x,t)\in(0,\infty)\times[0,T]$ represents the number of particles of
size $x$ at time $t$, $B(x,y)$ gives the breakage function, i.e., the rate at
which the particles of size $y$ break into particles of size $x$ and the rate
at which a particle size $y$ is chosen to break is shown by the selection
function $S(y)$. Further, the term $K(x-y,y)$ denotes the rate at which
particles of sizes $x-y$ and $y$ merge to form a particle of size $x$. In
equations (1.1) and (1.2), the first integral terms provide the birth of a
particle of size $x$ during the process of breakage and aggregation,
respectively, while the second terms in both models indicate the death of
particle size $x$.

Along with the number density $u(x,t)$, some integral properties, such as
moments, grab the attention because of their physical interpretation. The
moments for the number density are defined as

$\displaystyle\mu_{j}(t)=\int_{0}^{\infty}x^{j}u(x,t)dx,\quad j=0,1,2,\cdots.$

The zeroth moment $\mu_{0}(t)$ defines the total number of particles in the
system at time $t$, first moment $\mu_{1}(t)$ gives the total mass in system
and $\mu_{2}(t)$ gives the energy dissipated by the system.

In solid processing, e.g., in foods and pharmaceuticals, product quality is
characterized by multiple particle properties, for example, the volume and
composition of aggregating particles. To model such phenomenon, more then one
dimensional is required. Therefore, in the following, the bivariate case is
considered, i.e., particles (or individual objects) are characterized by two
properties, named $x$ and $y$. The two dimensional aggregation is governed by

$\displaystyle\dfrac{\partial u(x,y,t)}{\partial t}=$
$\displaystyle\frac{1}{2}\int_{0}^{x}\int_{0}^{y}K(x-x^{\prime},y-y^{\prime},x^{\prime},y^{\prime})u(x-x^{\prime},y-y^{\prime},t)u(x^{\prime},y^{\prime},t)dx^{\prime}dy^{\prime}$
$\displaystyle-$
$\displaystyle\int_{0}^{\infty}\int_{0}^{\infty}K(x,x^{\prime},y,y^{\prime})u(x,y,t)u(x^{\prime},y^{\prime},t)dx^{\prime}dy^{\prime},$
(1.4)

with the initial condition

$\displaystyle c(x,y,0)=c_{0}(x,y)\geq 0.$ (1.5)

Due to complexity of these models and unavailability of the analytical
solutions (except some simple cases), several numerical and semi-analytical
techniques are applied to solve these problems approximately. Numerical
schemes to solve breakage equation (1.1) and/or coagulation model (1.2)
includes finite element method [8], quadrature method of moments [9, 10],
finite volume scheme [11, 12, 13, 14], fixed pivot element [15], fast Fourier
transformation method [16], cell average technique [17] and reference therein.
The drawbacks of the schemes are shown in the potential reliance of these
numerical techniques on non-physical assumptions such as discretization,
linearization, sets of basis functions, and many others. Recently, several
authors have developed interest in semi-analytical approaches to overcome
these shortcomings. These series solution techniques offer results without
making such assumptions. Some of the available strategies are Taylor
polynomial and radial basis functions [18], Laplace-variational iteration
method [19], ADM [20], HPM [21], optimal homotopy asymptotic method (OHAM)
[22], HAM [23] and ODM [24]. Interestingly, some of the algorithm provided the
closed form series solutions of coagulation equation (1.2) for the aggregation
kernels

$K(x,y)=1,x+y,xy\text{ and }x^{\frac{2}{3}}+y^{\frac{2}{3}},$

with exponential initial condition ($u(x,0)=e^{-x},e^{-x}/x$), see [20, 21,
23] for more detailed computations. They also dealt with the breakage equation
(1.1) with the breakage rate

$b(x,y)=\frac{\alpha}{y}\left(\frac{x}{y}\right)^{\alpha-2}\quad\forall\quad
1\leq\alpha\leq 2\text{ with selection rate }S(x)=x^{\alpha}$

having the exponential ($e^{-x}$) and mono disperse $(\delta(x-a))$ being the
two different initial conditions. Hammouch and Mekkaoui in [19] developed the
Laplace-variational iteration method for solving the coagulation equation
(1.2) only for two cases of aggregation kernels, constant ($K(x,y)=1$) and
product ($K(x,y)=xy$). Moreover, Hasseine et al. in [25] employed ADM and HPM
to solve the breakage equation for the kernel $B(x,y)=\frac{12(y-x)}{x^{3}}$
with selection function $S(x)=x$. Very recently in [24], ODM is implemented to
solve the coagulation equation using the parameters

$K(x,y)=1,x+y\text{ and }xy\text{ with }u(x,0)=e^{-x}.$

In the literature, it was observed that ADM, HPM, and HAM provide the closed
form solutions, but some drawbacks are observed in these techniques. In [26,
27], it was found that a large number of iterations are required to obtain a
more accurate approximation. When dealing with chaotic systems, Chowdhury and
Hashimstill [28] found that time, time step, and the number of terms must be
handled with extreme caution. Further, Obidat [29] has drawn attention to
various drawbacks of ADM, including its delayed convergence [30] and inability
to handle boundary conditions [31] for solving non-linear PDEs. These
shortcoming were avoided by Obidat in [29]. To overcome these issues,
recently, ODM [24] has been implemented to solve the model, but the accuracy
is still maintained only for a small period of time. Recently, HPM is
accelerated by approximating the nonlinear term and incorporating the Elzaki
transformation for differential equations [32] in order to improve the
accuracy of the truncated solution. Thus, the first aim of this article is to
obtain more accurate solutions to the pure breakage and Smoluchowski’s
coagulation equations by applying the accelerated homotopy perturbation Elzaki
transformation method [32].

For the second task of this work, combined aggregation-breakage equation is
considered which is an intriguing issue for academics. The problem was
resolved using a class of numerical or stochastic methods. Lee and Matsoukas
[33] employed a stochastic process, namely the constant-N Monte Carlo method,
to solve the aggregation with a binary breakage equation. In 2002, Mahoney et
al. used the finite element method for aggregation, growth, and nucleation
equations [34]. Further, number density and moments were computed with the
help of the method of moments by Madras et al. in [35]. The model (1.1) and
(1.2) was also solved by implementing the finite volume method for several
test cases in [36, 37, 38, 39]. Since, numerical schemes have some limitations
and till date, there is no literature on semi-analytical schemes for coupled
aggregation-breakage model, here we implement the AHPETM for solving the
combined equation for two test cases. Moving further, the analytical solutions
for the bivariate aggregation equation are available for limited cases [40,
41, 42, 43]. Several numerical methods, such as moving sectional [44], finite
difference [45], Monte Carlo [46], sectional quadrature [47], dual quadrature
[48], finite volume schemes [49, 50], and many more [51, 52, 53, 54], are
considered to solve the equation. Therefore, our third aim here is to fill
this gap of series solution for finding the approximate results for bivariate
aggregation PBE.

The article is organized as follows: Section 2 discusses a brief outline of
the Elzaki transformation. In Section 3, the general methodology of HPM and
AHPETM are presented. In Section 4, AHPETM is developed for aforementioned
population balance equations. Further, Section 5 gives a detailed convergence
analysis of the proposed iterative scheme. In Section 6, the developed
formulations are adopted to demonstrate solutions for several kernels and the
supremacy of the scheme over HPM, ADM, HAM, and ODM solutions are shown by
means of numerical simulations.

##  2\. Elzaki Transformation and Its Properties

Tarig Elzaki developed the Elzaki transformation in 2011 [55, 56, 57, 58],
which is the modification of the general Laplace and Sumudu transformations to
solve the differential equation in the time domain. In [55, 56], authors show
the efficiency and accuracy of the Elzaki transformation on a large class of
differential and integral equations. To understand the definition of the
transformation, consider a set

$A=\left\\{f(t):\exists M,k_{1},k_{2}>0,|f(t)|<Me^{\frac{|t|}{k_{j}}},\text{
if }t\in(-1)^{j}\times[0,\infty)\right\\}$

then the Elzaki transformation is defined as

$E[f(t)]=T[v]=v\int_{0}^{\infty}f(t)e^{-\frac{t}{v}}dt,\quad t>0,$

and the inverse of Elzaki transformation [59] is defined as

$E^{-1}[T[v]]=\frac{1}{2\pi
i}\int_{0}^{\infty}e^{tv}T\left[\frac{1}{v}\right]vdv.$

Some of the Elzaki transformation for standard functions are listed in TABLE
1.

Table 1. Properties of Elzaki transformation $f(t)$ | $E[f(t)]$
---|---
1 | $v^{2}$
$t^{n}$ | $n!v^{n+2}$
$e^{at}$ | $\dfrac{v^{2}}{1-av}$
$E[f(t)+g(t)]$ | $E[f(t)]+E[g(t)]$
$E[f^{n}(t)]$ | $\quad\quad\frac{T[v]}{v^{n}}-\sum_{k=0}^{n-1}v^{2-n+k}f^{k}(0),\quad n\geq 1$

##  3\. Methodology

In this section, we review the basics of HPM and AHPETM for solving general
differential equations. Then the schemes are applied to solve multi-
dimensional coagulation and coupled coagulation-fragmentation equations.

###  3.1. Review of HPM

Let us consider the general differential equation

$\displaystyle D(c)-h(x)=0,\quad x\in\Omega$ (3.1)

with the boundary conditions

$\displaystyle B\left(c,\frac{\partial c}{\partial
n}\right)=0,r\in\partial\Omega,$ (3.2)

where $D$ and $B$ are the differential and boundary operators, respectively.
One can usually decompose the differential operator into linear ($L$) and non-
linear ($N$) operators, implying that equation (3.1) becomes

$\displaystyle L(c)+N(c)-h(x)=0.$ (3.3)

Now, according to HPM, a homotopy $H:\Omega\times[0,1]\rightarrow\mathbb{R}$
is constructed that satisfies

$\displaystyle H[v(x,p)]=(1-p)[L[v(r,p)]-L[(c_{0})]]+p[D[v(r,p)]-h(x)]=0,$
(3.4)

where $c_{0}$ is the initial guess for the equation (3.1) and $p$ is the
embedding parameter that increases monotonically from 0 to 1. According to the
HPM, we can write the solution of the equation (3.1) in the form of series as

$\displaystyle
v=\sum_{k=0}^{\infty}p^{k}v_{k}=v_{0}+pv_{1}+p^{2}v_{2}+\cdots.$ (3.5)

Substituting equation (3.5) in (3.4) and letting $p\rightarrow 1$, the
solution is obtained as follows

$\displaystyle c=\lim\limits_{p\rightarrow 1}v=\sum_{k=0}^{\infty}v_{k}.$
(3.6)

###  3.2. Accelerated Homotopy Perturbation Elzaki Transformation Method
(AHPETM)

Consider a non-linear differential equation

$\displaystyle\frac{\partial^{n}c}{\partial t^{n}}+L[c(x,t)]+N[c(x,t)]=b(x)$
(3.7)

with the initial conditions $c^{i}(x,0)=g_{i}(x),\quad i=0,1,2,\cdots,n-1$,
where $c^{i}(x,t)$ denotes the $i^{th}$ order derivative of $c(x,t)$ with
respect to $t$. Taking Elzaki transformation and using its properties on
equation (3.7) provide, by following [32],

$\displaystyle
E[c(x,t)]=\sum_{k=0}^{n-1}v^{k+2}c^{k}(x,0)+v^{n}E[b(x)-L[c(x,t)]-N[c(x,t)]].$
(3.8)

Now, applying the homotopy perturbation method to the equation (3.8), we get

$\displaystyle(1-p)(E[c(x,t)]-E[c(x,0)])+p\left(E[c(x,t)]-\sum_{k=0}^{n-1}v^{k+2}c^{k}(x,0)-v^{n}E[b(x)-L[c(x,t)]-N[c(x,t)]]\right)=0.$
(3.9)

Let the unknown function $c(x,t)$ and non-linear operator $N[c(x,t)]$ can be
written in series form as

$\displaystyle c(x,t)$ $\displaystyle=\sum_{n=0}^{\infty}v_{n}p^{n}$ (3.10)

and

$\displaystyle N[c(x,t)]=\sum_{n=0}^{\infty}H_{n}p^{n}$ (3.11)

where $H_{n}$ represents the accelerated He’s polynomial with

$\displaystyle H_{n}(x,t)=N(\sum_{i=0}^{n}v_{i})-\sum_{i=0}^{n-1}H_{i},\text{
for }n\geq 1\text{ and }H_{0}=N(v_{0}).$ (3.12)

Substituting the values of $c(x,t)$ and $N[c(x,t)]$ from the equations (3.10)
and (3.11) into equation (3.9) give

$\displaystyle
E[\sum_{n=0}^{\infty}v_{n}p^{n}]=\sum_{k=0}^{n-1}v^{k+2}c^{k}(x,0)+p\left\\{v^{n}E\left[g(x)-L[\sum_{n=0}^{\infty}v_{n}p^{k}]+\sum_{n=0}^{\infty}H_{n}p^{n}\right]\right\\}.$

Applying inverse Elzaki transformation and comparing the coefficients of
powers of $p$, the components of series solution, i.e., $v_{i}^{\prime}s$ are
given in TABLE 2

Table 2. Components of series solution $v_{0}$ | $c(x,0)$
---|---
$v_{1}$ | $\quad\sum_{k=1}^{n-1}\frac{t^{k}}{k!}c^{k}(x,0)+E^{-1}\\{v^{n}E[b(x)-L[v_{0}]+H_{0}]\\}$
$v_{2}$ | $-E^{-1}\\{v^{n}E[L[v_{1}]+H_{1}]\\}$
$\vdots$ | $\vdots$
$v_{n}$ | $-E^{-1}\\{v^{n}E[L[v_{n-1}]+H_{n-1}]\\}$

and hence the solution of the equation (3.7) is obtained by taking
$p\rightarrow 1$ in the equation (3.10).

##  4\. Development of Mathematical formation using AHPETM

In this AHPETM is extended to solve the Smoluchowski’s coagulation, pure
fragmentation, coupled coagulation-fragmentation and bivariate coagulation
equations.

###  4.1. Smoluchowski’s Coagulation Equation (SCE)

Consider the non-linear aggregation equation (1.2) with initial condition
$u(x,0)=u_{0}(x)$. Applying Elzaki transformation, an integral form is
obtained as

$\displaystyle
E[u(x,t)]=v^{2}u(x,0)+vE\left[\frac{1}{2}\int_{0}^{x}K(x-y,y)u(x-y,t)u(y,t)dy-\int_{0}^{\infty}K(x,y)u(x,t)u(y,t)dy\right].$
(4.1)

In order to apply the scheme, compare equation (4.1) with the transformed
equation (3.8), which provides $L[u]=0,\quad b(x)=0$ and

$\displaystyle
N[u]=-\frac{1}{2}\int_{0}^{x}K(x-y,y)u(x-y,t)u(y,t)dy+\int_{0}^{\infty}K(x,y)u(x,t)u(y,t)dy.$
(4.2)

Now, applying the HPM on equation (4.1) as defined in equation (3.9), we get

$\displaystyle(1-p)(E[c(x,t)]-E[c(x,0)])+p\bigg{(}E[c(x,t)]-v^{2}c(x,0)-vE\bigg{[}\frac{1}{2}\int_{0}^{x}K(x-y,y)u(x-y,t)u(y,t)dy$
$\displaystyle-\int_{0}^{\infty}K(x,y)u(x,t)u(y,t)dy\bigg{]}\bigg{)}=0.$ (4.3)

According to the methodology defined in Section 3.2,
$c(x,t)=\sum_{n=0}^{\infty}v_{n}p^{n}$ and the non-linear operator
$N[u]=\sum_{n=0}^{\infty}H_{n}p^{n}$, where $H_{n}$ for SCE is given by

$\displaystyle
H_{n}=\frac{1}{2}\int_{0}^{x}K(x-y,y)\sum_{i=0}^{n}v_{i}(x-y,t)\sum_{i=0}^{n}v_{i}(y,t)dy-\int_{0}^{\infty}K(x,y)\sum_{i=0}^{n}v_{i}(x,t)\sum_{i=0}^{n}v_{i}(y,t)dy-\sum_{i=0}^{n-1}H_{i},n\geq
1$ (4.4)

with $H_{0}=N[v_{0}].$ Using the above defined decomposition in equation (4.1)
and comparing the powers of $p$, the $n^{th}$ component of the series solution
is

$\displaystyle
v_{n+1}(x,t)=E^{-1}\bigg{\\{}vE\bigg{(}\frac{1}{2}\int_{0}^{x}K(x-y,y)\sum_{i=0}^{n}v_{i}(x-y,t)\sum_{i=0}^{n}v_{i}(y,t)dy$
$\displaystyle-\int_{0}^{\infty}K(x,y)\sum_{i=0}^{n}v_{i}(x,t)\sum_{i=0}^{n}v_{i}(y,t)dy\bigg{)}-\sum_{i=0}^{n}H_{i}\bigg{\\}}$
(4.5)

where $v_{0}(x,t)=u(x,0)$ and hence, the $n$ term truncated series solution is
calculated by

$\displaystyle\Psi_{n}^{CE}(x,t):=\sum_{j=0}^{n}v_{j}(x,t).$ (4.6)

###  4.2. Fragmentation Equation (FE)

Considering the pure fragmentation equation (1.1) and applying Elzaki
transformation, the following integral operator form is achieved

$\displaystyle
E[u(x,t)]=v^{2}u(x,0)+E\left(\int_{x}^{\infty}B(x,y)S(y)u(y,t)dy-S(x)u(x,t)\right).$
(4.7)

Next, equation (4.7) is compared with the equation (3.8) for the
implementation of AHPETM. It is observed that for the case of pure breakage
equation $N[u(x,t)]=b(x)=0$ and

$L[u(x,t)]=-\int_{x}^{\infty}B(x,y)S(y)u(y,t)dy+S(x)u(x,t).$

By following the steps discussed in the previous Section 3.2, a homotopy is
generated as follows

$\displaystyle(1-p)\\{E[u(x,t)]-E[u(x,0)]\\}+p\left(E[u(x,t)]-v^{2}u(x,0)-vE\left[\int_{x}^{\infty}B(x,y)S(y)u(y,t)dy-S(x)u(x,t)\right]\right).$
(4.8)

According to the proposed method, AHPETM introduces the solution of unknown
function $u(x,t)$ in the form of infinite series as
$u(x,t)=\sum_{j=0}^{\infty}v_{j}(x,t)$. Substituting this into equation (4.8)
and comparing the coefficients of the power of $p$, provide the iterations for
the solution as follows

$\displaystyle
v_{n+1}(x,t)=E^{-1}\left\\{vE\left[\int_{x}^{\infty}B(x,y)S(y)v_{n}(x,t)dy-S(x)v_{n}(x,t)\right]\right\\}$
(4.9)

where $v_{0}(x,t)=u(x,0)$ and the $n$ term truncated solution will be provided
as

$\displaystyle\Psi_{n}^{FE}(x,t):=\sum_{j=0}^{n}v_{j}(x,t).$ (4.10)

###  4.3. Coupled Coagulation-fragmentation Equation (CCFE)

The CCFE is governed by

$\displaystyle\frac{\partial u(x,t)}{\partial
t}=\frac{1}{2}\int_{0}^{x}K(x-y,y)u(x-y,t)u(y,t)dy-\int_{0}^{\infty}K(x,y)u(x,t)u(y,t)dy$
$\displaystyle+\int_{x}^{\infty}B(x,y)S(y)u(y,t)dy-S(x)u(x,t).$ (4.11)

Applying Elzaki transformation on both sides leads to

$\displaystyle
E[u(x,t)]=v^{2}u(x,0)+vE\bigg{[}\frac{1}{2}\int_{0}^{x}K(x-y,y)u(x-y,t)u(y,t)dy-\int_{0}^{\infty}K(x,y)u(x,t)u(y,t)dy$
$\displaystyle+\int_{x}^{\infty}B(x,y)S(y)u(y,t)dy-S(x)u(x,t)\bigg{]}.$ (4.12)

For the implementation of AHPETM, expression (4.3) is compared with (3.8) and
the following observations are made

$b(x)=0,\quad L[u]=-\int_{x}^{\infty}B(x,y)S(y)u(y,t)dy+S(x)u(x,t),$

and

$\displaystyle
N[u]=-\frac{1}{2}\int_{0}^{x}K(x-y,y)u(x-y,t)u(y,t)+\int_{0}^{\infty}K(x,y)u(x,t)u(y,t).$

Following the procedure defined in Section 3.2, the iterations to solve the
equation (4.3) are as follows

$\displaystyle v_{n+1}(x,t)=$ $\displaystyle
E^{-1}\bigg{\\{}vE\bigg{(}\frac{1}{2}\int_{0}^{x}K(x-y,y)\sum_{i=0}^{n}v_{i}(x-y,t)\sum_{i=0}^{n}v_{i}(y,t)dy$
$\displaystyle-\int_{0}^{\infty}K(x,y)$
$\displaystyle\sum_{i=0}^{n}v_{i}(x,t)\sum_{i=0}^{n}v_{i}(y,t)dy-\sum_{i=0}^{n}H_{i}+\int_{x}^{\infty}B(x,y)S(y)v_{n}(x,t)dy-S(x)v_{n}(x,t)\bigg{)}\bigg{\\}},$
(4.13)

where $v_{0}(x,t)=u(x,0).$ Let us denote the $n$ term approximated series
solution for CCFE as

$\displaystyle\Psi_{n}^{CCFE}(x,t):=\sum_{j=0}^{n}v_{j}(x,t).$ (4.14)

###  4.4. Bivariate Smoluchowski’s Coagulation Equation (BSCE)

Consider 2D aggregation equation (1) with initial condition
$u(x,y,0)=u_{0}(x,y)$ and applying Elzaki transformation, leads to form as

$\displaystyle E[u(x,y,t)]=$ $\displaystyle
v^{2}u(x,y,0)+vE\bigg{[}\frac{1}{2}\int_{0}^{x}\int_{0}^{y}K(x-x^{\prime},y-y^{\prime},x^{\prime},y^{\prime})u(x-x^{\prime},y-y^{\prime},t)u(x^{\prime},y^{\prime},t)dy^{\prime}dx^{\prime}$
$\displaystyle-\int_{0}^{\infty}\int_{0}^{\infty}K(x,x^{\prime},y,y^{\prime})u(x,y,t)u(x^{\prime},y^{\prime},t)dy^{\prime}dx^{\prime}\bigg{]}.$
(4.15)

In order to apply the AHPETM, equation (4.4) is compared with the transformed
equation (3.8), implying that $L[u]=0,\quad b(x)=0$ and

$\displaystyle N[u]=$
$\displaystyle-\frac{1}{2}\int_{0}^{x}\int_{0}^{y}K(x-x^{\prime},y-y^{\prime},x^{\prime},y^{\prime})u(x-x^{\prime},y-y^{\prime},t)u(x^{\prime},y^{\prime},t)dy^{\prime}dx^{\prime}$
$\displaystyle+$
$\displaystyle\int_{0}^{\infty}\int_{0}^{\infty}K(x,x^{\prime},y,y^{\prime})u(x,y,t)u(x^{\prime},y^{\prime},t)dy^{\prime}dx^{\prime}.$
(4.16)

Thanks to equation (3.9), applying the HPM on equation (4.4) enables us to
have

$\displaystyle(1-p)(E[c(x,t)]-E[c(x,0)])+p\bigg{(}E[c(x,y,t)]-v^{2}c(x,y,0)-vE\bigg{[}\frac{1}{2}\int_{0}^{x}\int_{0}^{y}K(x-x^{\prime},y-y^{\prime},x^{\prime},y^{\prime})$
$\displaystyle
u(x-x^{\prime},y-y^{\prime},t)u(x^{\prime},y^{\prime},t)dy^{\prime}dx^{\prime}-\int_{0}^{\infty}\int_{0}^{\infty}K(x,x^{\prime},y,y^{\prime})u(x,y,t)u(x^{\prime},y^{\prime},t)dy^{\prime}dx^{\prime}\bigg{]}\bigg{)}=0.$
(4.17)

Again, following the idea of Section 3.2,
$u(x,y,t)=\sum_{n=0}^{\infty}v_{n}p^{n}$ and non-linear operator
$N[u]=\sum_{n=0}^{\infty}H_{n}p^{n}$ where $H_{n}$ is being given by

$\displaystyle H_{n}=$
$\displaystyle\frac{1}{2}\int_{0}^{x}K(x-x^{\prime},y-y^{\prime},x^{\prime},y^{\prime})\sum_{i=0}^{n}v_{i}(x-x^{\prime},y-y^{\prime},t)\sum_{i=0}^{n}v_{i}(x^{\prime},y^{\prime},t)dy^{\prime}dx^{\prime}$
$\displaystyle-\int_{0}^{\infty}K(x,x^{\prime},y,y^{\prime})\sum_{i=0}^{n}v_{i}(x,y,t)\sum_{i=0}^{n}v_{i}(x^{\prime},y^{\prime},t)dy^{\prime}dx^{\prime}-\sum_{i=0}^{n-1}H_{i}\text{
with }H_{0}=N[v_{0}].$ (4.18)

Using the above defined decomposition in equation (4.4) and comparing the
powers of $p$, we get the $n^{th}$ component of the series solution as follows

$\displaystyle v_{n+1}(x,y,t)=$ $\displaystyle
E^{-1}\bigg{\\{}vE\bigg{(}\frac{1}{2}\int_{0}^{x}\int_{0}^{y}K(x-x^{\prime},y-y^{\prime},x^{\prime},y^{\prime})\sum_{i=0}^{n}v_{i}(x-x^{\prime},y-y^{\prime},t)\sum_{i=0}^{n}v_{i}(x^{\prime},y^{\prime},t)dy^{\prime}dx^{\prime}$
$\displaystyle-\int_{0}^{\infty}\int_{0}^{\infty}K(x,x^{\prime},y,y^{\prime})\sum_{i=0}^{n}v_{i}(x,y,t)\sum_{i=0}^{n}v_{i}(x^{\prime},y^{\prime},t)dy^{\prime}dx^{\prime}\bigg{)}\bigg{\\}}$
(4.19)

where $v_{0}(x,y,t)=u(x,y,0).$ Let us denote the $n$ term truncated solution
by

$\displaystyle\Psi_{n}^{BSCE}(x,y,t):=\sum_{j=0}^{n}v_{j}(x,y,t).$ (4.20)

##  5\. Convergence Analysis

###  5.1. Smoluchowski’s Coagulation Equation

Consider a Banach space
$\mathbb{X}=\mathbb{C}([0,T]:\mathbb{L}^{1}[0,\infty),\|.\|)$ over the norm
defined as

$\|u\|=\sup_{s\in[0,t_{0}]}\int_{0}^{\infty}|u(x,s)|dx<\infty.$

Let us use equation (4.1) in the operator form as

$\displaystyle u(x,t)=\tilde{\mathcal{N}}[u]$

where

$\displaystyle\tilde{\mathcal{N}}[u]=u(x,0)+E^{-1}\\{vE[N[u]]\\}$ (5.1)

and $N[u]$ is given by

$\displaystyle
N[u]=\frac{1}{2}\int_{0}^{x}K(x-y,y)u(x-y,t)u(y,t)dy-\int_{0}^{\infty}K(x,y)u(x,t)u(y,t)dy.$

######  Theorem 1.

Let us consider the coagulation equation (1.2) with kernel $K(x,y)=1$ for all
$x,y\in(0,\infty)$. If $v_{i}^{s}$ are the components of the series solution
computed using (4.1) and $\Psi_{n}^{CE}$ being the $n$ term truncated solution
provided in equation (4.6), then $\Psi_{n}^{CE}$ converges to the exact
solution $u$ with the error bound

$\|u-\Psi_{m}^{CE}\|\leq\dfrac{\Delta^{m}}{1-\Delta}\|v_{1}\|$

where $\Delta=t_{0}^{2}e^{2t_{0}L}(\|u_{0}\|+2t_{0}L^{2}+2t_{0}L)<1\text{ and
}L=\|u_{0}\|(T+1)$.

###### Proof.

Two separate phases complete the theorem’s proof. The contractive nature of
the non-linear operator $\tilde{\mathcal{N}}$ is initially demonstrated. Then
convergence of the truncated solution towards the exact one is established.
Step 1: As presented in [20], equation (5.1) can be written in the equivalent
form as

$\displaystyle\frac{\partial}{\partial
t}[u(x,t)\exp[H[x,t,u]]]=\frac{1}{2}\exp[H[x,t,u]]\int_{0}^{x}K(x-y,y)u(x-y,t)u(y,t)dy$

where $H[x,t,u]=\int_{0}^{t}\int_{0}^{\infty}K(x,y)u(y,s)dyds$. Thus the
equivalent operator $\tilde{N}$ is given by

$\displaystyle\tilde{N}[u]=u(x,0)\exp[-H(x,t,u)]+\frac{1}{2}\int_{0}^{t}\exp[H(x,s,u)-H(x,t,u)]\int_{0}^{\infty}K(x,y)u(x-y,s)u(y,s)dyds.$

Since $\tilde{N}$ is contractive (Singh et al. established in [20]) and
equivalent to $N[u]$, the non-linear operator $N[u]$ is also contractive,
i.e.,

$\displaystyle\|Nu-Nu^{*}\|\leq\delta\|u-u^{*}\|$ (5.2)

where $\delta:=t_{0}e^{2t_{0}L}(\|u_{0}\|+2t_{0}L^{2}+2t_{0}L)<1$ (for
suitably chosen $t_{0}$) and $L=\|u_{0}\|(T+1)$.
Now, using the definition and basic properties of Elzaki and Laplace
transformations as well as employing (5.2), we get

$\displaystyle\|\tilde{\mathcal{N}}u-\tilde{\mathcal{N}}u^{*}\|$
$\displaystyle=\|E^{-1}\\{vE(N(u))\\}-E^{-1}\\{vE(N(u^{*}))\\}\|$
$\displaystyle=\left\|\frac{1}{2\pi}\int_{0}^{\infty}\left(\frac{1}{v^{2}}\int_{0}^{\infty}(Nu-
Nu^{*})e^{-tv}dt\right)e^{tv}vdv\right\|$
$\displaystyle\leq\frac{1}{2\pi}\int_{0}^{\infty}\left(\frac{1}{v}\int_{0}^{\infty}\delta\|u-u^{*}\|e^{-tv}dt\right)e^{tv}dv$
$\displaystyle=\frac{1}{2\pi}\int_{0}^{\infty}\frac{1}{v}\mathcal{L}(\delta\|u-u^{*}\|)e^{tv}dv$
$\displaystyle=\mathcal{L}^{-1}\left\\{\frac{1}{v^{2}}\mathcal{L}(\delta\|u-u^{*}\|)\right\\}\leq\delta
t_{0}\|u-u^{*}\|\text{ for a suitable }t_{0}.$

Step 2: Now, in this phase, an $n$ term truncated solution is computed using
the iterations defined in (4.1) and then error is estimated. Given that,

$\displaystyle\Psi_{n}^{CE}=$ $\displaystyle\sum_{i=0}^{n}v_{i}(x,t)$
$\displaystyle=$ $\displaystyle
u(x,0)+E^{-1}\\{vE(N(u_{0}))\\}+E^{-1}\\{vE(N(u_{0}+u_{1})-H_{0})\\}+\cdots+E^{-1}\\{vE(N(\sum_{j=0}^{n-1}u_{j})-\sum_{j=0}^{n-2}H_{i})\\}$
$\displaystyle=$ $\displaystyle
u(x,0)+E^{-1}\\{vE(N[v_{0}]+N[v_{0}+v_{1}]+\cdots+N[v_{0}+v_{1}+\cdots+v_{n-1}]-$
$\displaystyle(H_{0}+(H_{0}+H_{1})+\cdots+(H_{0}+H_{1}+\cdots+H_{n-2})))\\}$
$\displaystyle=$ $\displaystyle
u(x,0)+E^{-1}\\{vE(N(\Psi_{n-1}^{CE}))\\}=\mathcal{\tilde{N}}[\Psi_{n-1}^{CE}].$

Using the contractive mapping of $\mathcal{\tilde{N}}$ leads to

$\|\Psi_{n+1}^{CE}-\Psi_{n}^{CE}\|\leq\Delta\|\Psi_{n}^{CE}-\Psi_{n-1}^{CE}\|.$

and thus, we have

$\displaystyle\|\Psi_{n+1}^{CE}-\Psi_{n}^{CE}\|$
$\displaystyle\leq\Delta\|\Psi_{n}^{CE}-\Psi_{n-1}^{CE}\|\leq\Delta^{n}\|\Psi_{1}^{CE}-\Psi_{0}^{CE}\|.$

Using the triangle inequality for all $n,m\in\mathbb{N}$ with $n>m$, we have

$\displaystyle\|\Psi_{n}^{CE}-\Psi_{m}^{CE}\|$
$\displaystyle\leq\|\Psi_{n}^{CE}-\Psi_{n-1}^{CE}\|+\|\Psi_{n-1}^{CE}-\Psi_{n-2}^{CE}\|+\cdots+\|\Psi_{m+1}^{CE}-\Psi_{m}^{CE}\|$
$\displaystyle\leq(\Delta^{n-1}+\Delta^{n-2}+\cdots+\Delta^{m})\|\Psi_{1}^{CE}-\Psi_{0}^{CE}\|$
$\displaystyle=\dfrac{\Delta^{m}(1-\Delta^{n-m})}{1-\Delta}\|u_{1}\|\leq\dfrac{\Delta^{m}}{1-\Delta}\|u_{1}\|,$

which converges to zero as $m\rightarrow\infty$, implies that there exists a
$\Psi$ such that $\lim\limits_{n\rightarrow\infty}\Psi_{n}^{CE}=\Psi.$
Therefore,

$\displaystyle
u(x,t)=\sum_{i=0}^{\infty}v_{i}=\lim\limits_{n\rightarrow\infty}\Psi_{n}^{CE}=\Psi,$

which is the exact solution of the coagulation equation (1.2). The theoretical
error is obtained by fixing $m$ and letting $n\rightarrow\infty$ in the above
formulation. ∎

###  5.2. Pure Breakage Equation

Let $\mathbb{X}=\mathbb{C}([0,T]:\mathbb{L}^{1}[0,\infty),\|.\|])$ be a Banach
space with the norm

$\displaystyle\|u\|=\sup_{t\in[0,t_{0}]}\int_{0}^{\infty}e^{\lambda
x}|u(x,t)|dx,\text{ where }\lambda>0.$ (5.3)

Now, equation (1.1) can be rewritten in the operator form as

$u=\mathcal{\tilde{L}}[u]=u(x,0)+E^{-1}{vE(L[u])}$

with $L[u]$ being the right-hand side of equation (1.1).

######  Theorem 2.

Let $\Psi_{n}^{FE}$ be the $n$ term truncated series solution of the
fragmentation problem defined in equation (1.1). Then $\Psi_{n}^{FE}$
converges to the exact solution and provides the error estimates

$\displaystyle\|u-\Psi_{m}^{FE}\|\leq\dfrac{\vartheta^{m}}{1-\vartheta}\|v_{1}\|,$
(5.4)

if the following conditions hold

  * B1.

$B(x,y)=c\dfrac{x^{r-1}}{y^{r}}$ where $r=1,2,3,\cdots$ and $c$ is a positive
constant satisfying $\int_{0}^{y}xB(x,y)dx=y,$

  * B2.

$S(x)\leq x^{k},\text{ where }k=1,2,3,\cdots,$

  * B3.

$\lambda$ is chosen such that $e^{\lambda y}-1<1,$

  * B4.

$\vartheta:=\dfrac{k!(t_{0})^{2}}{\lambda^{k+1}}<1$ for suitable $t_{0}$.

###### Proof.

Let us begin with the proof that the operator $\mathcal{\tilde{L}}$ is
contractive. In order to do so, we use the fact that the operator $L[u]$ is a
contractive operator under the assumptions mentioned in B1-B3 i.e.,,
$\|L[u]-L[u^{*}]\|\leq\rho\|u-u^{*}\|$ where
$\rho=\dfrac{k!t_{0}}{\lambda^{k+1}}<1$ by following ([20] Theorem 2.1). Now,
thanks to Elzaki and Laplace transformations, one can write

$\displaystyle\|\mathcal{L}[u]-\mathcal{L}[u^{*}]\|$
$\displaystyle=\|E^{-1}\\{vE(L[u])\\}-E^{-1}\\{vE(L[u^{*}])\\}\|$
$\displaystyle=\left\|\frac{1}{2\pi}\int_{0}^{\infty}\left(\frac{1}{v^{2}}\int_{0}^{\infty}(Lu-
Lu^{*})e^{-tv}dt\right)e^{tv}vdv\right\|$
$\displaystyle\leq\frac{1}{2\pi}\int_{0}^{\infty}\left(\frac{1}{v}\int_{0}^{\infty}\rho\|u-u^{*}\|e^{-tv}dt\right)e^{tv}dv$
$\displaystyle=\frac{1}{2\pi}\int_{0}^{\infty}\frac{1}{v}\mathcal{L}(\rho\|u-u^{*}\|)e^{tv}dv$
$\displaystyle=\mathcal{L}^{-1}\left\\{\frac{1}{v}\mathcal{L}(\rho\|u-u^{*}\|)\right\\}\leq\vartheta\|u-u^{*}\|\text{
where }\vartheta=\rho t_{0}.$

We proceed further to obtain the estimate (5.4). By using the iteration
formula (4.9), the $n$-term truncated solution is computed as

$\displaystyle\Psi_{n}^{FE}=$ $\displaystyle
E^{-1}\\{vE[v_{0}]\\}+E^{-1}\\{vE[v_{1}]\\}+\cdots+E^{-1}\\{vE[v_{n-1}]\\}$
$\displaystyle=$ $\displaystyle
E^{-1}\\{vE[v_{0}+v_{1}+\cdots+v_{n-1}]\\}=E^{-1}\\{vE[\Psi_{n-1}^{FE}]\\}.$

Therefore, we have

$\displaystyle\|\Psi_{n+1}^{FE}-\Psi_{n}^{FE}\|\leq\vartheta\|\Psi_{n}^{FE}-\Psi_{n-1}^{FE}\|\leq\vartheta^{n}\|\Psi_{1}^{FE}-\Psi_{0}^{FE}\|.$

The above can be used to establish the following, for all $n,m\in\mathbb{N}$
with $n>m$,

$\displaystyle\|\Psi_{n}^{FE}-\Psi_{m}^{FE}\|$
$\displaystyle\leq\|\Psi_{n}^{FE}-\Psi_{n-1}^{FE}\|+\|\Psi_{n-1}^{FE}-\Psi_{n-2}^{FE}\|+\cdots+\|\Psi_{m+1}^{FE}-\Psi_{m}^{FE}\|$
$\displaystyle\leq(\vartheta^{n-1}+\vartheta^{n-2}+\cdots+\vartheta^{m})\|\Psi_{1}^{FE}-\Psi_{0}^{FE}\|$
$\displaystyle=\dfrac{\vartheta^{m}(1-\vartheta^{n-m})}{1-\vartheta}\|v_{1}\|\leq\dfrac{\vartheta^{m}}{1-\vartheta}\|v_{1}\|.$

Thanks for hypothesis B4, the above tend to zero as $m\rightarrow\infty$ which
means that there exists a $\Psi$ such that
$\lim\limits_{n\rightarrow\infty}\Psi_{n}^{FE}=\Psi$. Thus, we obtain the
exact solution of the breakage equation (1.1) as

$\displaystyle
u(x,t)=\sum_{i=0}^{\infty}v_{i}=\lim\limits_{n\rightarrow\infty}\Psi_{n}^{FE}=\Psi.$

∎

###  5.3. Bivariate Smoluchowski’s Coagulation Equation

Consider a Banach space $X=\mathbb{C}([0,T]:L^{1}[0,\infty)\times
L^{1}[0,\infty),\|.\|)$ with the enduced norm

$\|u\|=\sup_{s\in[0,t_{0}]}\int_{0}^{\infty}\int_{0}^{\infty}|u(x,y,s)|dxdy<\infty.$

To demonstrate the convergence analysis, let us write the operator form of the
equation (4.4) as

$\displaystyle u=\tilde{Q}[u],$ (5.5)

where $\tilde{Q}$ is given by

$\displaystyle\tilde{Q}[u]=u_{0}(x,y)+E^{-1}[vE[Q[u]]]$ (5.6)

and

$\displaystyle Q[u]=$
$\displaystyle-\frac{1}{2}\int_{0}^{x}\int_{0}^{y}K(x-x^{\prime},y-y^{\prime},x^{\prime},y^{\prime})u(x-x^{\prime},y-y^{\prime},t)u(x^{\prime},y^{\prime},t)dy^{\prime}dx^{\prime}$
$\displaystyle+\int_{0}^{\infty}\int_{0}^{\infty}K(x,x^{\prime},y,y^{\prime})u(x,y,t)u(x^{\prime},y^{\prime},t)dy^{\prime}dx^{\prime}.$

The iterative scheme’s convergence concept is splitted into two components,
firstly, we establish that the operator $\tilde{Q}$ is contractive (Theorem 3)
and then proceed further to discuss the worst case upper bound for error
(Theorem 4) below. To show the operator $\tilde{Q}$ contractive, initially we
prove that $Q$ is contractive. To do so, an equivalent form of the equation
(1) is taken as

$\displaystyle\frac{\partial}{\partial
t}[u(x,y,t)\exp[R(x,y,t,c)]]=\frac{1}{2}\exp[R(x,y,t,c)]\int_{0}^{x}\int_{0}^{y}K(x-x^{\prime},x^{\prime},y-y^{\prime},y^{\prime})u(x-x^{\prime},y-y^{\prime},t)u(x^{\prime},y^{\prime},t)dy^{\prime}dx^{\prime},$
(5.7)

where
$R(x,y,t,c)=\int_{0}^{t}\int_{0}^{\infty}\int_{0}^{\infty}K(x,x^{\prime},y,y^{\prime})u(x^{\prime},y^{\prime},t)dx^{\prime}dy^{\prime}dt$.
Thus the equivalent operator $\mathcal{N}$ is given by

$\displaystyle\mathcal{N}[u]=u(x,y,0)\exp[-R(x,y,t,u)]+\frac{1}{2}\int_{0}^{t}\exp[R(x,y,s,u)-R(x,y,t,u)]$
$\displaystyle\int_{0}^{x}\int_{0}^{y}K(x-x^{\prime},x^{\prime},y-y^{\prime},y^{\prime})u(x-x^{\prime}y-y^{\prime},s)u(x^{\prime},y^{\prime},s)dy^{\prime}dx^{\prime}ds.$
(5.8)

Since, $\mathcal{N}$ and $Q$ are equivalent, it is sufficient to show that
$\mathcal{N}$ is contractive.

######  Theorem 3.

The operator $\tilde{Q}$, defined in equation (5.6) is contractive for all
$u,u^{*}\in\mathbb{X}$ if the following conditions

  * •

$K(x,x^{\prime},y,y^{\prime})=1\quad\forall
x,x^{\prime},y,y^{\prime}\in(0,\infty)$ and

  * •

$\delta=2t_{0}^{2}e^{2t_{0}L}(\|u\|+2t_{0}L^{2}+2t_{0}L)<1$ where
$L=\|u_{0}\|(T+1)$ hold.

###### Proof.

Consider $u,u^{*}\in\mathbb{X},$ then

$\displaystyle\mathcal{N}[u]-\mathcal{N}[u^{*}]=$ $\displaystyle
u(x,y,0)\exp[-R(x,y,t,u)]-u^{*}(x,y,0)\exp[-R(x,y,t,u^{*})]$
$\displaystyle+\frac{1}{2}\int_{0}^{t}\exp[R(x,y,s,u)-R(x,y,t,u)]\int_{0}^{x}\int_{0}^{y}u(x-x^{\prime},y-y^{\prime},s)u(x^{\prime},y^{\prime},s)dy^{\prime}dx^{\prime}ds$
$\displaystyle-\frac{1}{2}\int_{0}^{t}\exp[R(x,y,s,c^{*})-R(x,y,t,c^{*})]\int_{0}^{x}\int_{0}^{y}u^{*}(x-x^{\prime}y-y^{\prime},s)u^{*}(x^{\prime},y^{\prime},s)dy^{\prime}dx^{\prime}ds.$

Let us define an another operator

$H[x,y,s,t]=\exp\\{R[x,y,s,u]-R[x,y,t,u]\\}-\exp\\{R[x,y,s,u^{*}]-R[x,y,t,u^{*}]\\}.$

It can be easily proven that

$|H[x,y,s,t]|\leq(t-s)\exp\\{(t-s)B\\}\|u-u^{*}\|\leq B_{1}\|u-u^{*}\|,$

where $B_{1}=te^{tB}$ and $B=\max\\{\|u\|,\|u^{*}\|\\}.$ Further,

$\displaystyle\mathcal{N}[u]-\mathcal{N}[u^{*}]=$ $\displaystyle
u_{0}(x,y)H(x,y,0,t)+\frac{1}{2}\int_{0}^{t}H(x,y,s,t)\int_{0}^{x}\int_{0}^{y}u(x-x^{\prime},y-y^{\prime},s)u(x^{\prime},y^{\prime},s)dy^{\prime}dx^{\prime}ds$
$\displaystyle+\frac{1}{2}\int_{0}^{t}\exp[R(x,y,s,u^{*})-H(x,y,t,u^{*})]$
$\displaystyle\bigg{[}\int_{0}^{x}\int_{0}^{y}u^{*}(x-x^{\prime}y-y^{\prime},s)\\{u(x^{\prime},y^{\prime},s)-u^{*}(x^{\prime},y^{\prime},s)\\}dy^{\prime}dx^{\prime}$
$\displaystyle+\int_{0}^{x}\int_{0}^{y}u(x^{\prime},y^{\prime},s)\\{u(x-x^{\prime},y-y^{\prime},s)-u^{*}(x-x^{\prime},y-y^{\prime},s)\\}dx^{\prime}dy^{\prime}\bigg{]}ds.$
(5.9)

To show the non-linear operator $\mathcal{N}$ contractive, a set $D$ is
defined such as $D=\\{u\in\mathbb{X}:\|u\|\leq 2L\\}.$ Taking norm on both
sides of (5.3) provides

$\displaystyle\|\mathcal{N}[u]-\mathcal{N}[u^{*}]\|$ $\displaystyle\leq
B_{1}\|u-u^{*}\|\|u_{0}\|+B_{1}\|u-u^{*}\|\int_{0}^{t}\left[\frac{1}{2}\|u\|^{2}\right]ds+\int_{0}^{t}B_{1}\left[\frac{1}{2}(\|u\|+\|u^{*}\|)\|u-u^{*}\|\right]ds$
$\displaystyle\leq
B_{1}\left[\|u_{0}\|+\frac{1}{2}t\|u\|^{2}+\frac{1}{2}t(\|u\|+\|u^{*}\|)\right]\|u-u^{*}\|$
$\displaystyle\leq
t_{0}e^{2t_{0}L}\left[\|u_{0}\|+2t_{0}L^{2}+\frac{1}{2}t_{0}(2L+2L)\right]\|u-u^{*}\|$
$\displaystyle=\Delta\|u-u^{*}\|,$

for a suitable choice of $t_{0}$. So, the operator $\mathcal{N}$ is
contractive if
$\Delta=t_{0}e^{2t_{0}L}\left[\|u_{0}\|+2t_{0}L^{2}+2t_{0}L\right]<1$, hence
$Q$ is contractive. Now we are in position to demonstrate that
$\tilde{\mathcal{Q}}$ is contractive. Consider,

$\displaystyle\|\tilde{Q}u-\tilde{Q}u^{*}\|$
$\displaystyle=\|E^{-1}(vE\left[Qu\right])-E^{-1}(vE\left[Qu^{*}\right])\|$
$\displaystyle=\|\frac{1}{2\pi}\int_{0}^{\infty}\left(\frac{1}{v^{2}}\int_{0}^{\infty}(Qv-
Qv^{*})e^{-vt}dt\right)e^{vt}vdv\|$
$\displaystyle\leq\frac{1}{2\pi}\int_{0}^{\infty}\left(\frac{1}{v}\int_{0}^{\infty}\|Qu-
Qu^{*}\|e^{-vt}dt\right)e^{vt}dp$
$\displaystyle\leq\frac{1}{2\pi}\int_{0}^{\infty}\left(\frac{1}{v}\int_{0}^{\infty}\delta\|u-u^{*}\|e^{-vt}dt\right)e^{vt}dv$
$\displaystyle=\frac{1}{2\pi}\int_{0}^{\infty}\frac{1}{v}\mathcal{L}(\delta\|u-u^{*}\|)e^{vt}dv$
$\displaystyle=\mathcal{L}^{-1}\left\\{\frac{1}{v}\mathcal{L}(\delta\|u-u^{*}\|)\right\\}=\delta
t_{0}\|u-u^{*}\|=\Delta\|u-u^{*}\|.$

Which accomplishes the contractive nature of $\tilde{Q}$ by following the
assumption. ∎

######  Theorem 4.

Assuming that the criteria of Theorem (3) hold and $v_{i}^{\prime}s$ are the
elements of the series solution calculated by equation (4.4). Then the series
solution converges to the exact solution with the error bound

$\displaystyle\|u-\Psi_{n}^{FE}\|=\dfrac{\Delta^{n}}{1-\Delta}\|v_{1}\|,$

whenever $\Delta<1$ and $\|v_{1}\|<\infty.$

###### Proof.

The proof is similar to the Theorem 1, hence it is omitted here. ∎

######  Remark 5.1.

It is worth mentioning that the iterations and hence the finite term series
solutions, computed using the HAM [23], HPM [21], ADM [20] and ODM [24] are
identical to the iterations obtained using the AHPETM for the breakage
equation which is linear. As a result, we have omitted the numerical
implementations for pure breakage equation. So the main focus of all the
approaches is on approximating the non-linearity, which has no bearing on the
linearity in the equations.

##  6\. Numerical Results and Discussion

This section verifies numerically the effectiveness of the suggested approach
for coagulation, combined fragmentation-coagulation, and bivariate aggregation
equations. Three physical test cases are considered and results for the number
density and moments are compared with the precise solution as well as
established and recently developed methods (ADM, HPM, HAM, ODM) for SCE. Due
to the improved and significant results noticed in SCE, the numerical
implementation is made for solving the coupled CFE and BSCE. Two test cases of
CFE and one example of BSCE are taken into account to justify the
effectiveness of our scheme.

###  6.1. Smoluchowski’s Coagulation Equation

######  Example 6.1.

Consider the case of constant aggregation kernel $K(x,y)=1$ with the
exponential initial data $u(x,0)=e^{-x}$ and for this the exact solution

$u(x,t)=\frac{4}{(2+t)^{2}}e^{-\frac{2x}{2+t}},$

is discussed in [18].

Employing the equations (4.1) and (4.6), first three components of the series
solutions are given as follows

$v_{0}(x,t)=e^{-x},\quad v_{1}(x,t)=\frac{1}{2}te^{-x}(x-2),$ $\displaystyle
v_{2}(x,t)=t^{3}\left(\frac{x^{3}}{144}-\frac{x^{2}}{12}+\frac{x}{4}-\frac{1}{6}\right)e^{-x}+t^{2}\left(\frac{x^{2}}{8}-\frac{3x}{4}+\frac{3}{4}\right)e^{-x},$
$\displaystyle v_{3}(x,t)=$
$\displaystyle\frac{1}{40642560}t^{3}e^{-x}\bigg{(}t^{4}x^{7}+14t^{3}(7-4t)x^{6}+588(t-2)t^{2}(2t-3)x^{5}-2940t(t(t(4t-21)+36)-24)x^{4}$
$\displaystyle+11760(5(t-4)t((t-3)t+6)+48)x^{3}-35280(t(t(t(4t-35)+120)-240)+192)x^{2}$
$\displaystyle+70560(t(t(t(2t-21)+90)-240)+288)x-10080(t(t(t(4t-49)+252)-840)+1344)\bigg{)}.$

It is essential to mention here that the components $v_{i}$ are quite
complicated and due to the complexity of the terms, it is hard to find a
closed-form solution. Therefore, a three-term truncated solution is
considered. However, thanks to ”MATHEMATICA”, one can compute the higher order
terms using equation (4.1).

(a) AHPETM $(n=3)$ (b) Exact

Figure 1. Number density for AHPETM and exact solutions for Example 6.1

To see the accuracy of our proposed method, the approximated three-term and
exact solutions are plotted in Figures 1(a) and 1(b). One can scrutinize that
the AHPETM solution shows a remarkable agreement with the exact one.

(a) AHPETM error $(n=3)$ (b) HPM error $(n=3)$ (c) ODM error $(n=3)$

Figure 2. AHPETM, HPM/ADM/HAM & ODM errors Table 3. Numerical errors at $t=0.5,1,1.5\text{ and }2$ for $n=3,4,5,6$. $n$ |  $t=0.5$ |  $t=1$ |  $t=1.5$ |  $t=2$
---|---|---|---|---
3 |  0.0014 |  0.0153 |  0.0543 |  0.1239
4 |  1.366$\times 10^{-4}$ |  2.656$\times 10^{-3}$ |  1.294$\times 10^{-2}$ |  3.632$\times 10^{-2}$
5 |  $1.072\times 10^{-5}$ |  3.7972$\times 10^{-4}$ |  2.5718$\times 10^{-3}$ |  9.0682$\times 10^{-3}$
6 |  $7.154\times 10^{-7}$ |  4.6146$\times 10^{-5}$ |  4.3241$\times 10^{-4}$ |  1.8931$\times 10^{-3}$

Further, to see the beauty of our algorithm, errors between exact and AHPETM
solutions are compared with the errors between exact and other well-
established approximated solutions obtained via HPM/ADM/HAM and ODM in Figure
2. It is important to point out here that the HAM [23]/ADM [20]/HPM [21]
provide the same iterations and hence the identical finite term series
solution for the considered model. From Figure 2, it is observed that HAM is
very badly approximated in comparison to ODM and AHPETM further improves the
error of ODM very significantly. Table 3 depicts the numerical errors of
AHPETM at different time levels for $n=3,4,5$ and 6 using the formula provided
in [20]. As one can notice, the inaccuracy grows as time increases for a fixed
number of terms and error decreases when more terms in the approximated
solutions are taken into account.

(a) Number density at $t=2$ (b) Error at $x=5$

Figure 3. Number density and error

The superiority of AHPETM over HAM and ODM is also demonstrated in Figure 3.
Figure 3(a) represents the concentration of particles at time $t=2$ and the
solutions obtained using HAM and ODM blow up where the AHPETM solution matches
well with the exact solution. Results for error from Figure 3(b) indicate that
all the schemes are quite efficient for a short time where as for a
significant time, the errors due to HAM and ODM are relatively very high
compared to AHPETM.

(a) Zeroth moments (b) Second moments

Figure 4. Zeroth and second moments

Further, integral properties associated with number density are plotted in
Figure 4. The zeroth (total number of particles) and second (energy dispersed
by the system) moments are displayed and comparison are made with the precise
moments. In Figure 4(a), AHPETM offers superior approximations in the zeroth
moment while HAM under predicts the result and deviates almost exponentially
from the exact one. ODM shows much better approximation than HAM but still
suffers fluctuations. As shown in Figure 4(b), AHPETM continues to be the best
option as the second moment produced by AHPETM and HAM coincides with the
exact ones but ODM does not offer a decent estimate.

######  Example 6.2.

Let us take aggregation kernel $K(x,y)=x+y$ with the exponential initial
condition $u(x,0)=e^{-x}$. The exact number density is provided in [60] as

$u(x,t)=\frac{e^{\left(e^{-t}-2\right)x-t}I_{1}\left(2\sqrt{1-e^{-t}}x\right)}{\sqrt{1-e^{-t}}x},$

where $I_{1}$ is the Bessel function of the first kind.

Using the equations (4.1) and (4.6), first few components of the series
solution are determined as

$v_{0}(x,t)=e^{-x},\quad
v_{1}(x,t)=\frac{1}{2}te^{-x}\left(x^{2}-2x-2\right),$ $\displaystyle
v_{2}(x,t)=\frac{1}{720}t^{2}e^{-x}\left(tx(x^{5}-10x^{4}-20x^{3}+240x^{2}-120x-240)+60x^{4}-360x^{3}-180x^{2}+1080x+360\right).$

Continuing in a similar fashion, it is easy to compute the higher order
components to find better-approximated results. A four-term truncated solution
is considered here and the results are compared with the HAM and ODM solutions
using the same number of terms.

(a) AHPETM error $(n=4)$ (b) HPM error $(n=4)$ (c) ODM error $(n=4)$

Figure 5. AHPETM, HPM/ADM/HAM & ODM errors

(a) Number density at $t=2$ (b) Error at $x=5$

Figure 6. Number density and error Table 4. Absolute error for $x=5$ at different time levels $t$ |  Exact |  AHPETM |  ODM |  HPM |  AHPETM error |  ODM error |  HPM error
---|---|---|---|---|---|---|---
0.2 |  0.0129 |  0.0129 |  0.0118 |  0.0131 |  2.71288$\times 10^{-5}$ |  1.0352$\times 10^{-3}$ |  1.4248$\times 10^{-4}$
0.4 |  0.0146 |  0.0151 |  0.0053 |  0.0184 |  5.3035$\times 10^{-4}$ |  9.2238$\times 10^{-3}$ |  3.8424$\times 10^{-3}$
0.6 |  0.0138 |  0.0162 |  -0.0177 |  0.03868 |  2.3932$\times 10^{-3}$ |  3.1524$\times 10^{-2}$ |  2.4873$\times 10^{-2}$
0.8 |  0.0121 |  0.0179 |  -0.0601 |  0.1026 |  5.8862$\times 10^{-3}$ |  7.223$\times 10^{-2}$ |  9.060 $\times 10^{-2}$
1.0 |  0.0101 |  0.0203 |  -0.1226 |  0.2523 |  0.0102 |  0.1327 |  0.2422
1.2 |  0.0082 |  0.0220 |  -0.2034 |  0.5428 |  0.0137 |  0.2117 |  0.5345
1.4 |  0.0067 |  0.0197 |  -0.2989 |  1.042 |  0.0131 |  0.3057 |  1.0357
1.6 |  0.00545 |  0.00506 |  -0.4031 |  1.8325 |  0.00038 |  0.4085 |  1.8271

As observed in the previous case, AHPETM is again found to be more accurate
than HAM and ODM, see Figure 5. The error due to AHPETM is not only
significantly smaller than the existing approximated solutions of HAM and ODM
but also close to zero. This is also clear from Figure 6(a) which displays the
error at $t=2$ for all the schemes. Further, Figure 6(b) and Table 4
demonstrate that for a fixed $x$ and time upto $t=1$, HPM and ODM errors are
almost identical and insignificant but as the time increases, the error due to
the first two schemes are still the same but grow almost exponentially while
AHPETM performs almost consistently.

(a) Zeroth moments (b) Second moments

Figure 7. Zeroth and second moments

Moving further, approximated and exact moments are compared for all three
methods in Figure 7. Surprisingly, ODM under predicts and over predicts the
zeroth and second moment, respectively, while HPM and AHPETM gave almost
identical findings and provided an excellent approximations to the exact
zeroth and second moments. We would like to point out that the first moment is
constant in all cases therefore, graphs are omitted.

######  Example 6.3.

Consider the case of product aggregation kernel $K(x,y)=xy$ with the
exponential initial condition $u(x,0)=e^{-x}$ and the precise solution is
provided in [18] as

$u(x,t)=\sum_{k=0}^{\infty}\frac{t^{k}x^{3k}\exp(-(t+1)x)}{(k+1)!\Gamma(2k+2)}.$

Using the recursive scheme defined in equation (4.1), a five term truncated
solution is considered. Due to complexity of the terms, only few given here as

$v_{0}(x,t)=e^{-x},\quad
v_{1}(x,t)=\frac{1}{12}te^{-x}x\left(x^{2}-12\right),$ $\displaystyle
v_{2}(x,t)=\frac{1}{544320}\left(t^{2}e^{-x}x^{2}\left(tx^{7}-144tx^{5}+3024tx^{3}+756x^{4}-45360x^{2}+272160\right)\right).$

(a) AHPETM error $(n=5)$ (b) HPM error $(n=5)$ (c) ODM error $(n=5)$

Figure 8. AHPETM, HPM/ADM/HAM & ODM errors

(a) Number density at $t=2$ (b) Error at $x=5$

Figure 9. Number density and error

(a) Zeroth moments (b) Second moments

Figure 10. Zeroth and second moments

Figure 8 contrasts the error between the exact and truncated solutions
obtained via HAM/HPM/ADM, ODM, and AHPETM. The figure demonstrates that the
AHPETM outperforms both HPM and ODM results. Figure 9 further indicates the
scheme’s superiority as the unexpected behavior of the HPM and ODM solutions
are noticed, where as the AHPETM offers better estimates of the exact
solution. Figure 10 continues by contrasting the analytical moments with the
approximated moments. In situations where the HPM and AHPETM moments are
almost identical and closer to exact moments, as seen in the prior
occurrences, ODM moments explode.

From the above illustrations, it can be seen that in all contexts, AHPETM
performs better than ADM, HPM, HAM, and ODM. Therefore, due to the novelty of
the proposed scheme, we use AHPETM to solve the more complex models such as
coupled aggregation breakage and bivariate aggregation equations.

###  6.2. Coupled Aggregation-Breakage Equation

######  Example 6.4.

Considering the case of constant aggregation rate $(K(x,y)=1)$, binary
breakage $(B(x,y)=2/y)$ with the selection rate $S(x)=\frac{x}{2}$ and for the
initial condition $u(x,0)=4xe^{-2x}$, the exact solution for the problem (4.3)
is provided in [61].

Using the iterations defined in equation (4.3), components $v_{i}^{\prime}s$
of the solutions are computed as follows

$v_{0}(x,t)=4xe^{-2x},\quad
v_{1}(x,t)=\frac{1}{3}te^{-2x}\left(4x^{3}-6x^{2}-6x+3\right),$ $\displaystyle
v_{2}(x,t)=$
$\displaystyle\frac{1}{3780}\bigg{(}t^{2}e^{-2x}\bigg{(}8tx^{7}-56tx^{6}-84tx^{5}+840tx^{4}-420tx^{3}-1260tx^{2}+630tx+504x^{5}$
$\displaystyle-2520x^{4}-1890x^{3}+9450x^{2}+945x-1890\bigg{)}\bigg{)}.$

(a) Number density $(n=4)$ (b) Time distribution

Figure 11. Number density

(a) Truncated error  (b) Error at $t=0.5$ (c) Zeroth moment

Figure 12.  Error and moment

A four-term truncated solution is computed with the aid of ”MATHEMATICA”. At a
specific period $t$, the number density of particles of size $x$ is shown in
Figure 11(a). It is observed from Figures 11(a) and (b) that smaller particles
tend to increase as time goes on, while larger particles start to fragments
into smaller ones. The error between the exact and truncated solutions is
presented in Figure 12(a) and is found to be nearly insignificant. Further,
Figure 12(b) gives the absolute difference between the subsequent components
of the series solution and it is clear that the difference between the second
and third terms nearly vanishes, which serves as the inspiration for the
decision to truncate the solution for three terms. As shown in Figure 12(c),
the truncated solution exhibits steady state behavior for the number of
particles as the zeroth moment is constant. This behavior was also analyzed
analytically in [61], and thus demonstrating the method’s novelty.

######  Example 6.5.

Consider another example of the coupled aggregation-breakage equation (4.3)
with the same parameters as taken in Example 6.4 but with selection rate
$S(x)=2x$ and initial condition $u(x,0)=32xe^{-4x}$. Similar to the previous
case, here as well, steady state behavior of zeroth moment was studied in
[61].

Thanks to the formula (4.3), three terms of the truncated solution are
computed as

$v_{0}(x,t)=32xe^{-4x},\quad
v_{1}(x,t)=\frac{8}{3}te^{-4x}\left(32x^{3}-24x^{2}-12x+3\right),$
$\displaystyle v_{2}(x,t)=$
$\displaystyle\frac{8}{945}t^{2}e^{-4x}\bigg{(}1024tx^{7}-3584tx^{6}-2688tx^{5}+13440tx^{4}-3360tx^{3}-5040tx^{2}+1260tx$
$\displaystyle+8064x^{5}-20160x^{4}-7560x^{3}+18900x^{2}+945x-945\bigg{)}.$

(a) Number density $(n=4)$ (b) Time distribution

Figure 13. Number density

(a) Truncated error  (b) Error at $t=0.5$ (c) Zeroth moment

Figure 14.  Error and moment

Due to the complexity involved in the terms, a four-term truncated solution is
considered. The number density of particle size $x$ in the system is presented
in Figure 13(a). Further, Figure 13(b) presents the concentration of particles
at different time levels, and an increment in smaller particles is
encountered, where as larger particles start breaking as time increases. In
Figures 14(a) and 14(b), the difference between the consecutive terms is
presented, and the error between the second and third terms seems to be
vanishing, which leads us to truncate the solution for three terms. As
expected, in Figure 14(c) AHPETM shows the steady-state nature of zeroth
moment and is exactly matching with the precise total number of particles.

###  6.3. Bivariate Aggregation Equation

######  Example 6.6.

Let us take two-dimensional aggregation equation (1) with the constant
aggregation kernel $K(x,x^{\prime},y,y^{\prime})=1$ and the initial condition
$u(x,y,0)=\frac{16N_{0}xy}{m_{1}^{2}m_{2}^{2}}\exp\\{-\frac{2x}{m_{1}}-\frac{2y}{m_{2}}\\}$
where the parameters and the exact solution are given in Table 5. For more
details, readers may refer to [62].

Table 5. Parameters and exact solution $N_{0}$, $p_{1}$, $p_{2}$ | 1
---|---
$m_{1}$, $m_{2}$ | 0.04
$u(x,t)$ | $\frac{\left(4N_{0}\right)}{\left(m_{1}m_{2}\right)(t+2)^{2}}\left(\left(p_{1}+1\right){}^{p_{1}+1}\left(p_{1}+1\right){}^{p_{2}+1}\right)\exp\left(-\frac{\left(p_{1}+1\right)x}{m_{1}}-\frac{\left(p_{2}+1\right)y}{m_{2}}\right)$
|
$\sum_{k=0}^{\infty}\frac{\left(\frac{t}{t+2}\right)^{k}\left(\left(\left(p_{1}+1\right){}^{p_{1}+1}\right){}^{k}\left(\left(p_{2}+1\right){}^{p_{2}+1}\right){}^{k}\left(\frac{x}{m_{1}}\right){}^{(k+1)\left(p_{1}+1\right)-1}\left(\frac{y}{m_{2}}\right){}^{(k+1)\left(p_{2}+1\right)-1}\right)}{\Gamma\left(\left(p_{1}+1\right)(k+1)\right)\Gamma\left(\left(p_{2}+1\right)(k+1)\right)}$

The first three elements of the series solution are provided below using the
iterations specified in equation (4.4)

$v_{0}(x,y,t)=u_{0}(x,y,0),\quad v_{1}(x,y,t)=5.42535\times
10^{11}txye^{-50x-50y}\left(x^{2}y^{2}-0.1152\times 10^{-4}\right)$
$\displaystyle v_{2}(x,y,t)=$ $\displaystyle 3.10441\times
10^{-10}xye^{-50x-50y}\bigg{(}8.06248\times
10^{27}t^{3}x^{6}y^{6}-9.10222\times 10^{24}t^{3}x^{4}y^{4}$
$\displaystyle+8.73813\times 10^{20}t^{3}x^{2}y^{2}-3.35544\times
10^{15}t^{3}+1.36533\times 10^{25}t^{2}x^{4}y^{4}-2.62144\times
10^{21}t^{2}x^{2}y^{2}$ $\displaystyle+1.50995\times 10^{16}t^{2}+6t\bigg{)}.$

(a) Number density $(n=4)$ (b) Truncated error  (c)
$|\Phi_{2}(x,y,t)-u(x,y,t)|$ (d) $|\Phi_{3}(x,y,t)-u(x,y,t)|$ (e)
$|\Phi_{4}(x,y,t)-u(x,y,t)|$ (f) Moments

Figure 15. Number density, error and moments

Continuing in a similar pattern, a four-term truncated solution is computed
and compared with the exact solution. Figure 15(a) gives the number density at
time $t=0.4$ and it is marked that larger particles almost disappear, and
microscopic particles dominate the system. A minimal error is seen between the
exact and truncated solutions, according to the error curve shown in Figure
15(b). In addition to this, Figures 15(c)- 15(e) present the contour plots of
the errors by taking two, three, and four terms truncated series solutions.
One can observe that as the number of terms increases, the error reduces
significantly. Finally, Figure 15(f) shows that the approximated moments,
namely $\mu_{0,0},\mu_{1,0},\mu_{2,0}$, provide a great agreement with the
corresponding exact moments.

##  7\. Concluding remarks

This study employed AHPETM to solve the fragmentation, multi-dimensional
coagulation, and linked aggregation-fragmentation equations. Due to the
complexity in the models, convergence analysis were discussed for
fragmentation and multi-dimensional aggregation equations considering the
constant kernels. With the help of MATHEMATICA, this article also contained
the detailed numerical investigations for each of the predefined models. It
was observed that, for pure fragmentation equation which is linear, all the
schemes offered the same results. However, for non-linear aggregation
equation, AHPETM significantly outperformed the results of ADM, HAM, HPM and
ODM even after a lengthy period of time. This justified the method’s
reliability and applicability. AHPETM was also designed to solve non-linear
2-D aggregation and combined aggregation-fragmentation equations due to the
accuracy and efficiency observed in the pure aggregation equation and
remarkable results were obtained in each case.

## References

  * [1] H. Briesen, “Simulation of crystal size and shape by means of a reduced two-dimensional population balance model,” _Chemical Engineering Science_ , vol. 61, no. 1, pp. 104–112, 2006.
  * [2] A. Majumder, V. Kariwala, S. Ansumali, and A. Rajendran, “Lattice boltzmann method for population balance equations with simultaneous growth, nucleation, aggregation and breakage,” _Chemical Engineering Science_ , vol. 69, no. 1, pp. 316–328, 2012.
  * [3] I. Nopens, T. Koegst, K. Mahieu, and P. A. Vanrolleghem, “Pbm and activated sludge flocculation: from experimental data to calibrated model,” _AIChE journal_ , vol. 51, no. 5, pp. 1548–1557, 2005.
  * [4] D. Ramkrishna, _Population balances: Theory and pplications to particulate systems in engineering_. Elsevier, 2000\.
  * [5] M. J. Rhodes, _Introduction to particle technology_. John Wiley & Sons, 2008.
  * [6] R. P. Batycky, J. Hanes, R. Langer, and D. A. Edwards, “A theoretical model of erosion and macromolecular drug release from biodegrading microspheres,” _Journal of pharmaceutical sciences_ , vol. 86, no. 12, pp. 1464–1477, 1997\.
  * [7] Z. Yin and H. Liu, “Numerical simulation of nanoparticles diffusion and coagulation in a twin-jet via a temom method,” _International Journal of Numerical Methods for Heat & Fluid Flow_, 2014.
  * [8] Y. Bie, X. Cui, and Z. Li, “A coupling approach of state-based peridynamics with node-based smoothed finite element method,” _Computer Methods in Applied Mechanics and Engineering_ , vol. 331, pp. 675–700, 2018.
  * [9] D. L. Marchisio, R. D. Vigil, and R. O. Fox, “Quadrature method of moments for aggregation–breakage processes,” _Journal of colloid and interface science_ , vol. 258, no. 2, pp. 322–334, 2003.
  * [10] J. Su, Z. Gu, Y. Li, S. Feng, and X. Y. Xu, “Solution of population balance equation using quadrature method of moments with an adjustable factor,” _Chemical Engineering Science_ , vol. 62, no. 21, pp. 5897–5911, 2007.
  * [11] M. Singh, T. Matsoukas, A. B. Albadarin, and G. Walker, “New volume consistent approximation for binary breakage population balance equation and its convergence analysis,” _ESAIM: Mathematical Modelling and Numerical Analysis_ , vol. 53, no. 5, pp. 1695–1713, 2019.
  * [12] M. Singh and G. Walker, “Finite volume approach for fragmentation equation and its mathematical analysis,” _Numerical Algorithms_ , vol. 89, no. 2, pp. 465–486, 2022.
  * [13] F. Filbet and P. Laurençot, “Mass-conserving solutions and non-conservative approximation to the smoluchowski coagulation equation,” _Archiv der Mathematik_ , vol. 83, no. 6, pp. 558–567, 2004.
  * [14] R. Kumar, J. Kumar, and G. Warnecke, “Convergence analysis of a finite volume scheme for solving non-linear aggregation-breakage population balance equations,” _arXiv preprint arXiv:1403.1111_ , 2014.
  * [15] A. K. Giri and E. Hausenblas, “Convergence analysis of sectional methods for solving aggregation population balance equations: The fixed pivot technique,” _Nonlinear Analysis: Real World Applications_ , vol. 14, no. 6, pp. 2068–2090, 2013.
  * [16] R. Ahrens and S. Le Borne, “Fft-based evaluation of multivariate aggregation integrals in population balance equations on uniform tensor grids,” _Journal of Computational and Applied Mathematics_ , vol. 338, pp. 280–297, 2018.
  * [17] J. Kumar, M. Peglow, G. Warnecke, S. Heinrich, and L. Mörl, “Improved accuracy and convergence of discretized population balance for aggregation: The cell average technique,” _Chemical Engineering Science_ , vol. 61, no. 10, pp. 3327–3342, 2006.
  * [18] M. Ranjbar, H. Adibi, and M. Lakestani, “Numerical solution of homogeneous smoluchowski’s coagulation equation,” _International Journal of Computer Mathematics_ , vol. 87, no. 9, pp. 2113–2122, 2010.
  * [19] Z. Hammouch and T. Mekkaoui, “A laplace-variational iteration method for solving the homogeneous smoluchowski coagulation equation,” 2010.
  * [20] R. Singh, J. Saha, and J. Kumar, “Adomian decomposition method for solving fragmentation and aggregation population balance equations,” _Journal of Applied Mathematics and Computing_ , vol. 48, no. 1, pp. 265–292, 2015.
  * [21] G. Kaur, R. Singh, M. Singh, J. Kumar, and T. Matsoukas, “Analytical approach for solving population balances: a homotopy perturbation method,” _Journal of Physics A: Mathematical and Theoretical_ , vol. 52, no. 38, p. 385201, 2019.
  * [22] A. Dutta, Z. Pınar, D. Constales, and T. Öziş, “Population balances involving aggregation and breakage through homotopy approaches,” _International Journal of Chemical Reactor Engineering_ , vol. 16, no. 6, 2018\.
  * [23] G. Kaur, R. Singh, and H. Briesen, “Approximate solutions of aggregation and breakage population balance equations,” _Journal of Mathematical Analysis and Applications_ , vol. 512, no. 2, p. 126166, 2022.
  * [24] S. Kaushik and R. Kumar, “A novel optimized decomposition method for smoluchowski’s aggregation equation,” _Journal of Computational and Applied Mathematics_ , p. 114710, 2022.
  * [25] A. Hasseine, S. Senouci, M. Attarakih, and H.-J. Bart, “Two analytical approaches for solution of population balance equations: Particle breakage process,” _Chemical Engineering & Technology_, vol. 38, no. 9, pp. 1574–1584, 2015.
  * [26] D. Ganji and A. Sadighi, “Application of he’s homotopy-perturbation method to nonlinear coupled systems of reaction-diffusion equations,” _International Journal of Nonlinear Sciences and Numerical Simulation_ , vol. 7, no. 4, pp. 411–418, 2006.
  * [27] M. Dehghan, Y. Rahmani, D. D. Ganji, S. Saedodin, M. S. Valipour, and S. Rashidi, “Convection–radiation heat transfer in solar heat exchangers filled with a porous medium: homotopy perturbation method versus numerical analysis,” _Renewable Energy_ , vol. 74, pp. 448–455, 2015.
  * [28] J.-H. He, “Application of homotopy perturbation method to nonlinear wave equations,” _Chaos, Solitons & Fractals_, vol. 26, no. 3, pp. 695–700, 2005.
  * [29] Z. Odibat, “An optimized decomposition method for nonlinear ordinary and partial differential equations,” _Physica A: Statistical Mechanics and its Applications_ , vol. 541, p. 123323, 2020.
  * [30] Y. Jiao, Y. Yamamoto, C. Dang, and Y. Hao, “An aftertreatment technique for improving the accuracy of adomian’s decomposition method,” _Computers & Mathematics with Applications_, vol. 43, no. 6-7, pp. 783–798, 2002.
  * [31] J. He, “A new approach to nonlinear partial differential equations,” _Communications in Nonlinear Science and Numerical Simulation_ , vol. 2, no. 4, pp. 230–235, 1997.
  * [32] S. Jasrotia and P. Singh, “Accelerated homotopy perturbation elzaki transformation method for solving nonlinear partial differential equations,” in _Journal of Physics: Conference Series_ , vol. 2267, no. 1. IOP Publishing, 2022, p. 012106.
  * [33] K. Lee and T. Matsoukas, “Simultaneous coagulation and break-up using constant-n monte carlo,” _Powder Technology_ , vol. 110, no. 1-2, pp. 82–89, 2000.
  * [34] A. W. Mahoney and D. Ramkrishna, “Efficient solution of population balance equations with discontinuities by finite elements,” _Chemical Engineering Science_ , vol. 57, no. 7, pp. 1107–1119, 2002.
  * [35] G. Madras and B. J. McCoy, “Reversible crystal growth–dissolution and aggregation–breakage: numerical and moment solutions for population balance equations,” _Powder technology_ , vol. 143, pp. 297–307, 2004.
  * [36] J. Kumar, G. Warnecke, M. Peglow, and S. Heinrich, “Comparison of numerical methods for solving population balance equations incorporating aggregation and breakage,” _Powder Technology_ , vol. 189, no. 2, pp. 218–229, 2009\.
  * [37] R. Kumar, J. Kumar, and G. Warnecke, “Moment preserving finite volume schemes for solving population balance equations incorporating aggregation, breakage, growth and source terms,” _Mathematical Models and Methods in Applied Sciences_ , vol. 23, no. 07, pp. 1235–1273, 2013.
  * [38] J. Kumar, J. Saha, and E. Tsotsas, “Development and convergence analysis of a finite volume scheme for solving breakage equation,” _SIAM Journal on Numerical Analysis_ , vol. 53, no. 4, pp. 1672–1689, 2015.
  * [39] J.-P. Bourgade and F. Filbet, “Convergence of a finite volume scheme for coagulation-fragmentation equations,” _Mathematics of Computation_ , vol. 77, no. 262, pp. 851–882, 2008.
  * [40] J. Fernandez-Diaz and G. Gomez-Garcia, “Exact solution of smoluchowski’s continuous multi-component equation with an additive kernel,” _EPL (Europhysics Letters)_ , vol. 78, no. 5, p. 56002, 2007.
  * [41] F. Gelbard and J. H. Seinfeld, “Simulation of multicomponent aerosol dynamics,” _Journal of Colloid and Interface Science_ , vol. 78, no. 2, pp. 485–501, 1980.
  * [42] F. Leyvraz, “Scaling theory and exactly solved models in the kinetics of irreversible aggregation,” _Physics Reports_ , vol. 383, no. 2-3, pp. 95–212, 2003.
  * [43] A. Lushnikov, “From sol to gel exactly,” _Physical Review Letters_ , vol. 93, no. 19, p. 198302, 2004.
  * [44] Y. P. Kim and J. H. Seinfeld, “Simulation of multicomponent aerosol condensation by the moving sectional method,” _Journal of Colloid and Interface Science_ , vol. 135, no. 1, pp. 185–199, 1990.
  * [45] N. V. Mantzaris, P. Daoutidis, and F. Srienc, “Numerical solution of multi-variable cell population balance models: I. finite difference methods,” _Computers & Chemical Engineering_, vol. 25, no. 11-12, pp. 1411–1440, 2001.
  * [46] T. Matsoukas, T. Kim, and K. Lee, “Bicomponent aggregation with composition-dependent rates and the approach to well-mixed state,” _Chemical Engineering Science_ , vol. 64, no. 4, pp. 787–799, 2009.
  * [47] M. Attarakih, M. Jaradat, C. Drumm, H. Bart, S. Tiwari, V. Sharma, J. Kuhnert, and A. Klar, “A multivariate sectional quadrature method of moments for the solution of the population balance equation,” _Computer Aided Chemical Engineering_ , vol. 28, pp. 1551–1556, 2010.
  * [48] J. Favero and P. Lage, “The dual-quadrature method of generalized moments using automatic integration packages,” _Computers & Chemical Engineering_, vol. 38, pp. 1–10, 2012.
  * [49] G. Kaur, J. Kumar, and S. Heinrich, “A weighted finite volume scheme for multivariate aggregation population balance equation,” _Computers & Chemical Engineering_, vol. 101, pp. 1–10, 2017.
  * [50] M. Singh, J. Kumar, A. Bück, and E. Tsotsas, “An improved and efficient finite volume scheme for bivariate aggregation population balance equation,” _Journal of Computational and Applied Mathematics_ , vol. 308, pp. 83–97, 2016.
  * [51] H. M. Vale and T. F. McKenna, “Solution of the population balance equation for two-component aggregation by an extended fixed pivot technique,” _Industrial & Engineering Chemistry Research_, vol. 44, no. 20, pp. 7885–7891, 2005.
  * [52] S. Kumar and D. Ramkrishna, “A general discretization technique for solving population balance equations involving bivariate distributions,” in _AIChE Annual Meeting, Miami Beach, FL, USA, November_ , vol. 12, 1995, p. 17.
  * [53] A. Chaudhury, A. Kapadia, A. V. Prakash, D. Barrasso, and R. Ramachandran, “An extended cell-average technique for a multi-dimensional population balance of granulation describing aggregation and breakage,” _Advanced Powder Technology_ , vol. 24, no. 6, pp. 962–971, 2013.
  * [54] J. Kumar, M. Peglow, G. Warnecke, and S. Heinrich, “The cell average technique for solving multi-dimensional aggregation population balance equations,” _Computers & Chemical Engineering_, vol. 32, no. 8, pp. 1810–1830, 2008\.
  * [55] T. M. Elzaki, “The new integral transform elzaki transform,” _Global Journal of Pure and Applied Mathematics_ , vol. 7, no. 1, pp. 57–64, 2011.
  * [56] ——, “Application of new transform “elzaki transform” to partial differential equations,” _Global Journal of pure and applied Mathematics_ , vol. 7, no. 1, pp. 65–70, 2011.
  * [57] ——, “On the connections between laplace and elzaki transforms,” _Advances in Theoretical and Applied Mathematics_ , vol. 6, no. 1, pp. 1–11, 2011.
  * [58] T. M. Elzaki and S. M. Ezaki, “On the elzaki transform and ordinary differential equation with variable coefficients,” _Advances in Theoretical and Applied Mathematics_ , vol. 6, no. 1, pp. 41–46, 2011.
  * [59] T. M. Elzaki _et al._ , “On the new integral transform”elzaki transform”fundamental properties investigations and applications,” _Global Journal of Mathematical Sciences: Theory and Practical_ , vol. 4, no. 1, pp. 1–13, 2012.
  * [60] W. T. Scott, “Analytic studies of cloud droplet coalescence i,” _Journal of Atmospheric Sciences_ , vol. 25, no. 1, pp. 54–65, 1968.
  * [61] P. Lage, “Comments on the”an analytical solution to the population balance equation with coalescence and breakage-the special case with constant number of particles”by dp patil and jrg andrews [chemical engineering science 53 (3) 599-601],” _Chemical Engineering Science_ , vol. 19, no. 57, pp. 4253–4254, 2002.
  * [62] M. Singh, R. Singh, S. Singh, G. Walker, and T. Matsoukas, “Discrete finite volume approach for multidimensional agglomeration population balance equation on unstructured grid,” _Powder Technology_ , vol. 376, pp. 229–240, 2020.