# BLUES iteration applied to nonlinear ordinary differential equations for
wave propagation and heat transfer

Jonas Berx  Institute for Theoretical Physics, KU Leuven, B-3001 Leuven,
Belgium  Joseph O. Indekeu  Institute for Theoretical Physics, KU Leuven,
B-3001 Leuven, Belgium

(August 27, 2024)

###### Abstract

The iteration sequence based on the BLUES (Beyond Linear Use of Equation
Superposition) function method for calculating analytic approximants to
solutions of nonlinear ordinary differential equations with sources is
elaborated upon. Diverse problems in physics are studied and approximate
analytic solutions are found. We first treat a damped driven nonlinear
oscillator and show that the method can correctly reproduce oscillatory
behaviour. Next, a fractional differential equation describing heat transfer
in a semi-infinite rod with Stefan-Boltzmann cooling is handled. In this case,
a detailed comparison is made with the Adomian decomposition method, the
outcome of which is favourable for the BLUES method. As a final problem, the
Fisher equation from population biology is dealt with. For all cases, it is
shown that the solutions converge exponentially fast to the numerically exact
solution, either globally or, for the Fisher problem, locally.

##  I Introduction

Differential equations (DEs) are the cornerstones of modern physics and
mathematics. While linear DEs are quite well understood, nonlinear DEs remain
elusive objects. Notwithstanding this, most physical systems exhibit some form
of nonlinearity. It is a challenge to obtain exact analytical solutions
whenever possible. Some analytical or semi-analytical techniques such as the
Adomian decomposition method (ADM), the homotopy analysis method (HAM) or
perturbative techniques such as the soliton perturbation theory have been
proposed and proved to be extremely useful Adomian (1994); Liao (2012);
Karpman and Maslow (1977); Keener and McLaughlin (1977). When an external
force is also present in the form of an inhomogeneous source or sink Maeda and
Colonius (2017); Ermakov and Stepanyants (2019), one often turns to numerical
solutions. These solutions do not offer the richness of information that
analytical solutions can provide. While linear DEs with sources can be treated
by the theory of Green functions and the principle of superposition, nonlinear
DEs violate the superposition principle and a solution cannot be constructed
in this way. In this paper, following up on our Letter Berx and Indekeu
(2019), we investigate further how the practice of Green functions can be
usefully extended to nonlinear DEs with an inhomogeneous source or sink,
effectively using the superposition principle beyond the linear domain.

Relative to previous works on this method our contribution is situated as
follows. In Indekeu and Smets (2017) the usefulness was demonstrated of simple
exponential tail solutions of nonlinear reaction-diffusion-convection DEs
describing traveling wave fronts. These simple solutions are exact provided a
co-moving Dirac delta function source is added to the DE. In Indekeu and
Müller-Nedebock (2018) it was observed that for some problems one and the same
simple exponential tail solution simultaneously solves the nonlinear DE with a
Dirac delta source as well as a related linear DE with the same Dirac delta
source. This observation led to the idea to formulate an analytic method that
uses the concept of Green function, and its convolution with an arbitrary
source, beyond the linear domain. This approach was named BLUES (Beyond Linear
Use of Equation Superposition) and it was suggested that the method can be
useful as a perturbation expansion, the small parameter being the ratio of the
width of the source to the decay length of the tail. Subsequently, in Berx and
Indekeu (2019) quantitative calculations were performed, demonstrating that
the method entails an analytic iteration procedure which converges
exponentially rapidly for a variety of problems, and which is non-
perturbative. There is no need for a small parameter. Furthermore, the method
was also used the other way around, starting from a linear DE and freely
adding a nonlinearity to it, instead of starting from a nonlinear DE and
looking for a related linear DE. Applications were presented to solitary wave
solutions of the Camassa-Holm equation and traveling wavefront solutions of
the Burgers equation. Our contribution now is to extend substantially the
applicability of the method and demonstrate its accuracy for problems
involving oscillatory waves, phenomena described by fractional differential
equations (FDEs), and systems with nonlinear growth. We present a detailed
comparison with an alternative method, the ADM. In the new applications we
show that the convergence of the method is sometimes local instead of the
global convergence illustrated in Berx and Indekeu (2019). The applications
until now were for problems that can be reduced to an ordinary DE with a co-
moving source. This limitation is henceforth removed.

This paper is organized as follows. In Section II we recall briefly the
analytic iteration procedure in order to make the paper sufficiently self-
contained. In Section III, we study the damped nonlinear oscillator and show
that the method allows one to approximate correctly the decaying oscillations.
In Section IV, a model for heat transfer in a semi-infinite rod with a
fractional-order derivative is treated, and employed to show that the BLUES
function method can be extended to the arena of fractional differential
equations. Also in this Section we compare the performance of the method to
that of the ADM. Finally, in Section V, a wave solution of a nonlinear
extension of the heat equation involving reaction, i.e., the Fisher equation,
is treated. We close the paper with our conclusions and outlook.

##  II The BLUES iteration method

Here we recapitulate concisely the BLUES function concept Indekeu and Müller-
Nedebock (2018) from the viewpoint of the efficient iteration procedure that
was developed from it in Berx and Indekeu (2019). One starts from a linear
ordinary DE which can be written as an operator $\mathcal{L}_{z}$ acting on a
function $U(z)$ and assumes that one knows the (piecewise analytic) Green
function $U(z)=G(z)$ which solves

$\mathcal{L}_{z}U(z)=\delta(z),$ (1)

with suitable boundary conditions. The Dirac delta source compensates a
possible discontinuity of the derivative of order $n-1$ at $z=0$ in the case
of an $n$-th order DE. Next one considers the linear DE with an arbitrary
source $\psi(z)$. The solution is then the convolution product $G\ast\psi$.

One may add a nonlinearity rather freely, but so that the same boundary
conditions are respected, and arrive at the nonlinear ordinary DE

$\mathcal{N}_{z}U(z)=\psi(z),$ (2)

A solution is then proposed in the form $U(z)=B\ast\phi$. The function $B(z)$
is called BLUES function and it is taken to be the Green function of the
linear DE, so $B(z)=G(z)$. The challenge is to calculate the associated source
$\phi(z)$ for the given source $\psi(z)$, knowing that $B\ast\psi$ solves the
linear DE with source $\psi(z)$. For achieving this one defines a residual
operator $\mathcal{R}_{z}\equiv\mathcal{L}_{z}-\mathcal{N}_{z}$ and makes use
of the implicit identity

$\phi(z)=\psi(z)+\mathcal{R}_{z}(B\ast\phi)(z).$ (3)

To obtain the solution to the nonlinear DE (2), equation (3) can be iterated
in order to calculate an approximation in the form of a sequence in powers of
the residual $\mathcal{R}_{z}$. To zeroth order, the sources $\phi^{(0)}(z)$
and $\psi(z)$ are identical and the approximation is the convolution product

$U_{\psi}^{(0)}(z)=(B\ast\phi^{(0)})(z)=(B\ast\psi)(z).$ (4)

To $n$th order ($n\geq 1$), the approximate solution can be found by iterating
(3) and taking the convolution product with $B(z)$ Berx and Indekeu (2019)

$U_{\psi}^{(n)}=(B\ast\phi^{(n)})(z)=U_{\psi}^{(0)}(z)+\left(B\ast\mathcal{R}_{z}U_{\psi}^{(n-1)}\right)(z).$
(5)

We now proceed to applications beyond what was presented in Berx and Indekeu
(2019).

##  III The damped nonlinear oscillator

We start from a general linear wave equation in one dimension with a co-moving
Dirac delta source with reduced amplitude $s$

$u_{tt}-u_{xx}+\gamma u_{x}+u=s\,\delta(x-ct),$ (6)

in which the displacement $u$, time $t$, space $x$, amplitude $\gamma$ and
velocity $c$ are also reduced so as to be dimensionless, and look for
traveling wave solutions by transforming to the coordinate $z=x-ct$ and
restricting $U(z)\equiv u(x,t)$,

$\alpha U_{zz}+\gamma U_{z}+U=s\,\delta(z),$ (7)

where $\alpha=c^{2}-1$. Derivatives (with respect to $t$, $x$ or $z$) have
been denoted by subscripts.

This DE is solved by the Green function, to be used as BLUES function,

$B(z)\equiv\begin{cases}0,&z<0\\\ \sin{\left(\frac{\lambda
z}{2\alpha}\right)}\mathrm{e}^{-\frac{\gamma z}{2\alpha}},&z\geq
0,\end{cases}$ (8)

with $\lambda=\sqrt{4\alpha-\gamma^{2}}$ and source amplitude $s=\lambda/2$.
Now an arbitrary nonlinear term can be added, which we choose to be the cubic-
quintic function $\beta U^{3}+\xi U^{5}$, where $\beta$ and $\xi$ are tuneable
parameters. Altogether the nonlinear wave equation with an arbitrary source is

$\alpha U_{zz}+\gamma U_{z}+U+\beta U^{3}+\xi U^{5}=s\psi(z),$ (9)

where again the amplitude is $s=\lambda/2$ and the 1-norm (i.e., the integral
over $z$) of the source $\psi(z)$ is unity. This DE is a basic model for a
myriad of physical systems. When $\beta$ and $\xi$ are chosen to be $-1/3!$
and $1/5!$ respectively, the terms $U+\beta U^{3}+\xi U^{5}$ can be
interpreted as the first three nonzero terms in the sine Taylor series.
Including higher-order terms in the series, one can construct the nonlinear DE

$\alpha U_{zz}+\gamma U_{z}+\sin{U}=s\psi(z),$ (10)

which is an equation for the damped and driven Sine-Gordon model Ekomasov  _et
al._ (2018), often used to describe the dynamics of Josephson junctions in
superconductors Gul  _et al._ (2018); Starodub and Zolotaryuk (2019). Another
important application of equation (9) is the cubic-quintic Duffing oscillator,
which is used to describe damped harmonic motion in a nontrivial potential and
has become a paradigm for the study of chaos. For computational convenience we
will only include the terms in the sine Taylor series up to and including
third order, so we will choose $\xi=0$.

Following the procedure outlined in Section II, we extract the residual
operator and identify its action on $U$,

$\mathcal{R}_{z}U\equiv-\frac{\beta}{s}U^{3}$ (11)

and use this to calculate higher-order approximants using $B(z)$ in the
iteration sequence (5). Note that $\mathcal{R}_{z}B(z)\neq 0$. For the source
$\psi$ we can choose, e.g., an exponential corner source which possesses a
tunable dimensionless decay length $K$,

$\psi(z)=\frac{\mathrm{e}^{-|z|/K}}{2K}$ (12)

The choice of the source (12) is not limited to functions with non-zero
1-norm; we have also performed calculations for sources with 1-norm zero and
found that this has no significant impact on the convergence or accuracy of
the method.

The zeroth-order approximant for the solution of (9) can be calculated by
performing the convolution integral (4) with BLUES function (8) and normalized
exponential corner source (12)

$U^{(0)}_{\psi}(z)=\begin{cases}\frac{K\lambda}{4C_{+}}\mathrm{e}^{z/K},&z<0\\\
\left[A\sin{\left(\frac{\lambda
z}{2\alpha}\right)}-B\lambda\cos{\left(\frac{\lambda
z}{2\alpha}\right)}\right]\frac{\mathrm{e}^{-\frac{\gamma
z}{2\alpha}}}{2}+\frac{K\lambda}{4C_{-}}\mathrm{e}^{-z/K},&z\geq
0,\end{cases}$ (13)

where the constants $A,B,C_{\pm}$ are introduced to simplify notation. They
are given by combinations of $\alpha$, $\gamma$ and $K$:

$\begin{split}A&=\frac{2\alpha^{2}-K^{2}\gamma^{2}+2\alpha
K^{2}}{\alpha^{2}-K^{2}\gamma^{2}+2\alpha K^{2}+K^{4}}\\\
B&=\frac{K^{2}\gamma}{\alpha^{2}-K^{2}\gamma^{2}+2\alpha K^{2}+K^{4}}\\\
C_{\pm}&=\alpha\pm K\gamma+K^{2}.\end{split}$ (14)

Higher orders can in principle be calculated using equation (5) but in
practice this is not feasible without the aid of mathematical software.

Figure 1:  Traveling wave solution to the nonlinear wave equation (9) with an
exponential corner source (12). The numerical solution $U_{\rm num}$ (red,
full line), the zeroth-order $U_{\psi}^{(0)}$ (black, dashed line), the first-
order approximant $U_{\psi}^{(1)}$ (black, wider spaced dashed line) and the
second-order approximant $U_{\psi}^{(2)}$ (black, dot-dashed line, calculated
for $z<9$) are compared. The BLUES function (gray, solid line) is also shown.
Parameter values are $\alpha=3$, $\gamma=1$, $\beta=1$, $\xi=0$ and $K=1/5$.
Figure 2: Wavefront values $U_{\psi}^{(n)}(z=4)$ versus order $n$ for the
nonlinear wave equation (9). The numerically exact value (red, dashed line) is
also shown. Inset: A $\mathop{\log_{10}}$ semi-log plot of the increments
$|\Delta U_{\psi}^{(n,num)}(4)|\equiv|U_{\psi}^{(n)}(4)-U^{(num)}(4)|$ of the
approximants versus $n$, and a linear fit. Parameter values are $\alpha=3$,
$\gamma=1$, $\beta=1$, $\xi=0$ and $K=1/5$. Figure 3: Residual function
$\mathcal{R}_{z}U_{\psi}^{(n)}$ of the approximants ($n=0,1,2$) for the
nonlinear wave equation (9) with exponential corner source (12). The residual
operator applied to the BLUES function and to the numerical solution (red,
full line) are also shown. Parameter values are $\alpha=3$, $\gamma=1$,
$\beta=1$, $\xi=0$ and $K=1/5$.

In Fig. 1, a comparison between the numerically exact solution and the
zeroth-, first- and second-order approximants is made. Note that the zeroth-
order approximant almost coincides with the BLUES function (also shown), for
$z>0$, while the second-order approximant is on top of the numerically exact
solution. Next, we analyse the convergence of the iteration sequence to the
exact solution. In Fig. 2 the values of the approximants at $z=4$ for orders
$n=0,1$ and 2 are shown. The inset shows that the convergence to the exact
solution is exponentially fast.

Special attention is given to the calculation of the residual, at iteration
$n$, which provides another means of monitoring the convergence of the method.
The residual in the BLUES function and the residual functions for $n=0,1$ and
$2$ are shown in Fig. 3, together with the residual operator (11) applied to
the numerically exact solution (red, full line). Note that the residual
functions are negative in the domain around the global maximum of $U(z)$,
where $U>0$, that they are zero (with a cubic dependence on $z$) where the
curves of Fig. 1 change sign, and that they are positive in the domain around
their first minimum, where $U<0$. This is obvious in view of the simple form
(11).

In the limit $K\rightarrow 0$, the chosen source converges to the Dirac delta
source used to calculate the BLUES function (8). Because the BLUES function
does not solve (9), also not in this limit, the iteration sequence converges
(exponentially fast) to a nontrivial function. One can calculate the first two
terms by setting $\psi(z)=\delta(z)$ in equations (4) and (5). To zeroth
order, the convolution is identical to the BLUES function. To first order in
the iteration, the solution with a Dirac delta function source is given by

$U_{\delta}^{(1)}(z)=B(z)-\frac{\beta}{s}\int\limits_{-\infty}^{\infty}B(z-z_{0})B(z_{0})^{3}\mathrm{d}z_{0}.$
(15)

Note that this does not provide the first-order correction in a series
expansion in the parameter $\beta$ because higher-order iterations generate
additional contributions linear in $\beta$.

##  IV Fractional heat transfer equation

In this section we consider the following nonlinear ordinary fractional
differential equation (FDE) for $U(t)$ defined on the semi-infinite real line
$t\in[0,\infty)$ with differential order $0<\alpha\leq 1$ and exponent $n\geq
1$ and with initial condition $U(t=0)=C_{0}$, with $C_{0}\geq 0$.

$\mathcal{N}_{t}U=D^{\alpha}_{t}U+U^{n}=\psi(t)\,,$ (16)

where $D^{\alpha}_{t}$ is the Riemann- Liouville fractional derivative defined
as follows, for $\alpha>0$ and $t>0$,

$D^{\alpha}_{t}f(t)\equiv\begin{dcases}\frac{1}{\Gamma(m-\alpha)}\frac{d^{m}}{dt^{m}}\int_{0}^{t}\frac{f(\tau)}{(t-\tau)^{\alpha+1-m}}\mathrm{d}\tau,&m-1<\alpha<m\in\mathbb{N}\\\
\frac{d^{m}}{dt^{m}}f(t),&\alpha=m\in\mathbb{N}\end{dcases}$ (17)

where $\Gamma(.)$ is the gamma function. This equation has previously been
studied in the context of nonlinear heat transfer for the case that
$\alpha=1/2$, $C_{0}=0$ and $n=4$ (Stefan-Boltzmann cooling) Duan  _et al._
(2012); Wazwaz and Khuri (1996). The calculations that follow are valid for
all values of $0<\alpha\leq 1$ and $n\geq 1$. It has been shown that if
$\psi(t)$ is a piecewise continuous bounded function, equation (16) is
guaranteed to have a unique solution. If $\psi(t)$ is nondecreasing in an
interval $0<t<s$, $s\in(0,\infty)$ then the solution is also nondecreasing in
that interval Padmavally (1958); Keller and Olmstead (1972). Note that the
differential order can in principle be higher than $\alpha=1$. One can then
separate the order $\alpha=m+\beta$ in an integer part $m\in\mathbb{N}$
corresponding to a regular integer-order differential operator, and a
fractional part $0<\beta\leq 1$ which again corresponds to a fractional
differential operator. The FDE (16) should consequently be supplemented with
additional boundary (or initial) conditions up to a number $m+1$. In the
remainder of this work, we will assume $0<\alpha\leq 1$.

Following the steps in the BLUES procedure outlined in Section II, we can
simply start from the linear DE obtained by dropping the nonlinear $U^{n}(t)$
term. We write the linear FDE in operator form

$\mathcal{L}_{t}U=D^{\alpha}_{t}U=\psi(t)\,,$ (18)

with arbitrary source $\psi(t)$ and the initial condition chosen to be
$U(0)=0$. The Green function for (18) can now be calculated by considering a
Dirac delta function source instead of $\psi(t)$

$D^{\alpha}_{t}G(t,t^{\prime})=\delta(t-t^{\prime})\,,$ (19)

where $t-t^{\prime}\geq 0$ because the problem is formulated on the semi-
infinite real line. This Green function is readily calculated to be

$G(t,t^{\prime})=\frac{(t-t^{\prime})^{\alpha-1}}{\Gamma(\alpha)}$ (20)

Consequently the solution of the linear FDE (18) is the convolution integral
of the Green function (20) and the source $\psi(t)$, which we will choose from
now on to be the constant function $\psi(t)=1$ for $t\geq 0$, as was done in
references Wazwaz and Khuri (1996); Duan  _et al._ (2012).

It is worth emphasizing that this application to heat transfer is
fundamentally different from the other ones considered in this paper as well
as in Berx and Indekeu (2019) in that the source is not assumed to be
originating from a disturbance that is “co-moving” with the solution. The
variable here is time and not a co-moving coordinate, and the source arises as
a natural physical ingredient of the problem. Therefore this example
constitutes a non-trivial extension of the domain of applicability of the
method not only in the type of DE (from DE to FDE) but also in the character
and interpretation of the source term in the DE.

We obtain

$U^{(0)}(t)=\int\limits_{0}^{t}G(t,t^{\prime})\psi(t^{\prime})\mathrm{d}t^{\prime}=\frac{1}{\Gamma(\alpha)}\int\limits_{0}^{t}(t-t^{\prime})^{\alpha-1}\mathrm{d}t^{\prime}=\frac{t^{\alpha}}{\alpha\Gamma(\alpha)}$
(21)

The residual operator is the difference between the operators of the linear
FDE (18) and the nonlinear FDE (16) and is defined by the action on $U(t)$,
i.e.,

$\mathcal{R}_{t}U=\mathcal{L}_{t}U-\mathcal{N}_{t}U=-U^{n}$ (22)

Now the $p-$th order approximant to equation (16) can be calculated by using
the BLUES iteration sequence (5)

$U^{(p)}(t)=U^{(0)}(t)+\int\limits_{0}^{t}G(t,t^{\prime})\mathcal{R}_{t^{\prime}}U^{(p-1)}(t^{\prime})\mathrm{d}t^{\prime}$
(23)

The first-order approximant to the nonlinear problem can easily be calculated
using (20), (21) and the iteration sequence definition (23), with the choice
$n=4$,

$\begin{split}U^{(1)}(t)&=U^{(0)}(t)-\int\limits_{0}^{t}\frac{(t-t^{\prime})^{\alpha-1}}{\Gamma(\alpha)}\frac{t^{\prime
4\alpha}}{\alpha^{4}\Gamma^{4}(\alpha)}\mathrm{d}t^{\prime}\\\
&=\frac{1}{\alpha\Gamma(\alpha)}t^{\alpha}-\frac{\Gamma(1+4\alpha)}{\alpha^{4}\Gamma^{4}(\alpha)\Gamma(1+5\alpha)}t^{5\alpha}\end{split}$
(24)

One can now iterate (23) to generate higher-order approximants to the solution
of the nonlinear FDE (16). In Fig. 4 and Fig. 5, the approximants for
different values of $\alpha$ are compared with the numerically exact solution.

(a)

(b)

Figure 4: Solution to the nonlinear FDE (16) with constant source $\psi(t)=1$
and fractional order (a) $\alpha=1/4$ and (b) $\alpha=1/2$. The numerical
solution (red, full line) is compared with the approximants up to fourth order
(black/gray lines).

(a)

(b)

Figure 5: Solution to the nonlinear FDE (16) with constant source $\psi(t)=1$
and fractional order (a) $\alpha=3/4$ and (b) $\alpha=1$. The numerical
solution (red, full line) is compared with the approximants up to fourth order
(black/gray lines).

For the choice $\alpha=1/2$, equation (16) is associated with the heat
transfer equations for a semi-infinite solid Keller and Olmstead (1972) with
external heating $\psi(t)$ and either linear Newton (for $n=1$) or nonlinear
Stefan-Boltzmann cooling (for $n=4$). In Keller and Olmstead (1972), it was
shown that for $n\geq 3$, some of the energy entering the solid will remain,
while for $n\leq 2$, all energy is eventually radiated away. For the remainder
of this work, we will use nonlinear Stefan-Boltzmann cooling, $n=4$. The
zeroth-, first-, and second-order approximants for $\alpha=1/2$ and $n=4$ are,
respectively, given by

$\begin{split}U^{(0)}(t)&=2\left(\frac{t}{\pi}\right)^{1/2}\\\
U^{(1)}(t)&=2\left(\frac{t}{\pi}\right)^{1/2}-\frac{256}{15}\left(\frac{t}{\pi}\right)^{5/2}\\\
U^{(2)}(t)&=2\left(\frac{t}{\pi}\right)^{1/2}-\frac{256}{15}\left(\frac{t}{\pi}\right)^{5/2}+\frac{2097152}{4725}\left(\frac{t}{\pi}\right)^{9/2}-\frac{1073741824}{225225}\left(\frac{t}{\pi}\right)^{13/2}\\\
&+\frac{8796093022208}{369208125}\left(\frac{t}{\pi}\right)^{17/2}-\frac{2251799813685248}{49104680625}\left(\frac{t}{\pi}\right)^{21/2}\end{split}$
(25)

The approximants obtained by the BLUES function method can be compared to
those obtained with the Adomian decomposition method (ADM) by making use of
the following recursion relation Duan  _et al._ (2012) for the coefficients
$a_{n}$ of the solution series of the ADM,

$U_{\rm ADM}(t)=\sum_{n=0}^{\infty}a_{n}\,t^{\alpha n},$ (26)

with

$\begin{split}a_{0}&=0\\\ a_{1}&=\frac{1}{\Gamma(\alpha+1)}\\\
a_{n+1}&=-\frac{\Gamma(n\alpha+1)}{\Gamma(n\alpha+\alpha+1)}A_{n}\end{split}$
(27)

where the $A_{n}$, $n\geq 1$ are defined as follows

$A_{n}=\sum_{l=0}^{n}\sum_{s=0}^{l}\sum_{k=0}^{s}a_{n-l}a_{l-s}a_{s-k}a_{k}$
(28)

From equations (27) and (28) it can be deduced that for $\alpha=1/2$ one has
$a_{n}=0$ for $n\neq 4k+1$, with $k=0,1,2,...$. We will henceforth define the
$p$-th order ADM approximant as the truncated series

$U_{\rm ADM}^{(p)}(t)\equiv\sum_{n=0}^{p}a_{n}\,t^{\alpha n},$ (29)

A comparison is now made between BLUES and ADM. In Fig.6 the 21st-order
approximant for the ADM (containing five nontrivial exact terms) is shown,
together with the 4th approximant for the BLUES function method (containing
many more terms, but also five nontrivial exact ones – see further) and the
numerically exact solution.

Figure 6: Comparison between the 21st-order approximant of the ADM (dotted
line) and the 4th iteration of the BLUES function method (dashed line) for
$\alpha=1/2$ and $n=4$. The numerically exact solution (red, full line) is
also shown.

The number of nonzero terms $g(i,n)$ for the $i$-th approximant generated by
the BLUES function method can be calculated for a general nonlinearity
exponent $n\in\mathbb{N}$, $n\geq 1$, i.e.,

$g(i,n)=\begin{cases}\frac{1}{n-1}\left(n-2+n^{i}\right)&n\geq 2\\\
i+1&n=1\end{cases}$ (30)

so for $n=4$, the 4th approximant already contains $g(4,4)=86$ nonzero terms.
Note that the number of terms in BLUES increases exponentially with the order
of iteration. In comparison, the ADM generates a series with a number of terms
which grows linearly with the order of the approximation. Note that the ADM
generates the exact coefficients in a series expansion of the solution while
the BLUES function method does not. In contrast, BLUES generates many more
terms, of which only the lower-order ones are exact. The higher-order
coefficients have not yet settled or converged to their exact value. We find
empirically that in the BLUES function method in iteration $p$ the first $p+1$
coefficients are exact, which is a linear progression like in ADM. For
example, in (25) the expressions for $U^{(0)}(t)$ and $U^{(1)}(t)$ contain the
exact coefficients, while only the first three terms in $U^{(2)}(t)$ are
exact. Although the higher-order terms are not exact, it appears that their
presence (in BLUES) leads to a more accurate approximant than their absence
(in ADM), when the comparison is made with equal numbers of exact terms.

The BLUES function method generates, in each iteration, a (huge) number of
scout terms that probe the emerging series expansion and gradually gain
precision. In this respect, the BLUES function method is reminiscent in spirit
of a Padé approximation applied to a series expansion. The coefficients
$a_{n}$ in the BLUES function method saturate roughly linearly with increasing
order of iteration, as can be seen, e.g., when keeping track of the
coefficient $a_{17}$ of the $t^{17/2}$ term. This coefficient is first
generated in the second approximant $U^{(2)}$ with a provisional value of
$a_{17}=1.41659$. In the third approximant $U^{(3)}$ the value of this
coefficient increases to $a_{17}=19.8026$ and settles in the fourth
approximant $U^{(4)}$, attaining its exact value $a_{17}=30.8436$, which is
the same value as is found with the ADM. Similar observations can be made for
the coefficient $a_{21}$ of the $t^{21/2}$ term. This term is first generated
in the second approximant with a value of $a_{21}=-0.27627$ and increases (in
absolute value) to $a_{21}=-40.9762$ in the third iteration, then jumps to
$a_{21}=-99.0372$ in the fourth iteration and settles at $a_{21}=-118.387$ in
the fifth iteration, which is the exact result as found by the ADM.

The higher-order residual functions are shown in Fig.7 together with the
residual operator (22) applied to the numerically exact solution (red, full
line). Note that for this model the residuals are not localized, in contrast
to all previous examples. The approximants to the solution of equation (16) do
not converge to the correct numerical value at $t\rightarrow\infty$ owing to
the divergence of the residual in every order of iteration. While the
approximants diverge for larger values of $t$, a radius of convergence can
still be identified.

Figure 7: Residual $\mathcal{R}_{t}U(t)=-U^{4}(t)$ for $\alpha=1/2$ for
different orders of iteration. The numerically exact residual
$\mathcal{R}_{t}U_{num}$ is also shown (red, full line).

##  V Fisher equation

As a starting point for our final example, consider the diffusion equation
which describes the propagation of a density $u(x,t)$

$u_{t}-\nu u_{xx}=0,$ (31)

with $\nu$ the (dimensionless) diffusion coefficient. Adopting a traveling-
wave Ansatz $z=x-ct$, the diffusion equation becomes a linear ordinary DE. We
add a co-moving Dirac delta source,

$\mathcal{L}_{z}U\equiv-U_{z}-kU_{zz}=\delta(z),$ (32)

where $k$ is a dimensionless constant. We consider the wavefront boundary
conditions $U_{z}(z\rightarrow-\infty)=0$ (and $U(z\rightarrow-\infty)>0$) and
$U(z\rightarrow\infty)=0$. The exact solution (in every point including $z=0$)
is the piecewise analytic exponential tail,

$B(z)=\begin{cases}1,&z<0\\\ \mathrm{e}^{-z/k},&z\geq 0\end{cases}$ (33)

and the wavefront velocity is $c(k)=1/k$. We now add a reaction-type
nonlinearity (growth term) and a source $\psi(z)$ to obtain the forced Fisher
equation Fisher (1937) in co-moving coordinates, i.e.,

$\mathcal{N}_{z}U=-U_{z}-kU_{zz}-kU(1-U)=\psi(z)$ (34)

with boundary conditions $U(z\rightarrow\infty)\rightarrow 0$ and
$U(z\rightarrow-\infty)\rightarrow 1$. Equation (34) governs the dimensionless
density of some bio-chemical substance or biological population experiencing
diffusion and growth. The limit at negative infinity signifies the saturation
of the density at the normalized value 1. The residual operator is now acting
as follows

$\mathcal{R}_{z}U=kU(1-U)$ (35)

If we now choose $\psi(z)$ to be the exponential corner source (12), the
zeroth-order approximant to the nonlinear DE (34) is the same as was
calculated for the Burgers equation in Berx and Indekeu (2019), and is
repeated here in appendix A. Higher orders can easily be calculated by
iteration but will not be given here. We will, however, present the calculated
first-order approximant for $k\neq K$ in the appendix. Note that the first-
order approximant approaches a constant which is not unity at negative
infinity for $z$ and consequently does not obey the boundary condition. One
can calculate the non-trivial constant by considering the limit of the first-
order approximant at negative infinite $z$, i.e.,

$\begin{split}U_{c}&\equiv\lim_{z\rightarrow-\infty}U^{(1)}(z)=1+\lim_{z\rightarrow-\infty}\int\limits_{-\infty}^{\infty}B(z-z^{\prime})\mathcal{R}_{z^{\prime}}U^{(0)}_{\psi}(z^{\prime})\mathrm{d}z^{\prime}\\\
&=1+\int\limits_{-\infty}^{\infty}\mathcal{R}_{z^{\prime}}U^{(0)}_{\psi}(z^{\prime})\mathrm{d}z^{\prime}\\\
&=1+\frac{k\left(2k^{3}+4k^{2}K+6kK^{2}+3K^{3}\right)}{4(k+K)^{2}}\end{split}$
(36)

In Fig.8(a), the zeroth- and first-order approximants are shown together with
the numerically exact solution and the BLUES function (33). In Fig.8(b), a
zoomed-in representation of the shoulder of the wavefront is shown. The
numerical solution is compared with approximants up to fourth order. Note that
while all approximants obey the boundary condition
$U(z\rightarrow\infty)\rightarrow 0$, only the zeroth-order approximant
approaches unity for $z\rightarrow-\infty$. This is a consequence of the lack
of localization of the residual for higher orders. This is shown in Fig.9. The
numerical residual function $\mathcal{R}_{z}U_{num}$ and the zeroth-order
residual function are localized, but higher-order residual functions are not
anymore. This corresponds to a divergence of the approximants of higher
orders.

(a)

(b)

Figure 8: (a) Traveling wavefront solution to the nonlinear Fisher DE (34)
with an exponential corner source (12). The numerical solution $U_{\rm num}$
(red, full line), the zeroth-order $U_{\psi}^{(0)}$ (black, dashed line) and
the first-order approximant $U_{\psi}^{(1)}$ (black, wider spaced dashed line)
are compared. The BLUES function (gray) is also shown. (b) A zoomed-in view
around the shoulder of the wavefront. The approximants are shown up to and
including 4th order. The latter (gray, full line) lies just below the
numerical one (red, full line). Parameter values are $k=1/3$ and $K/k=1/2$.
Figure 9: Residual function $\mathcal{R}_{z}U_{\psi}^{(n)}$ of the
approximants of order $n=0,1,2,3$ for the nonlinear Fisher DE with exponential
corner source (12). The residual operator applied to the numerical solution
(red, full line) is also shown. Parameter values are $k=1/3$ and $K/k=1/2$.

The local convergence of the approximants can nevertheless be assessed by
studying the value of the approximants for a fixed value of $z$. In view of
the divergence of the approximants for $z\rightarrow-\infty$, one has to be
careful in choosing the value of $z$. In this case, we have opted for $z=-1$,
which can be seen to lie within a reasonable region of convergence. The
results are shown in Fig.10. Note that within the region of convergence, the
approximants converge exponentially fast to the numerically exact solution.

Figure 10: Wavefront values $U_{\psi}^{(n)}(-1)$ at $z=-1$ versus order $n$
for the Fisher equation (34). The numerically exact value (red, dashed line)
is also shown. Inset: A $\mathop{\log_{10}}$ semi-log plot of the increments
$|\Delta U_{\psi}^{(n,n-1)}(-1)|$ of the approximants versus $n$, and a linear
fit. These increments are defined as the difference of the $n$-th and
$(n-1)$-st approximant. Parameter values are $k=1/3$ and $K/k=1/2$.

##  VI Conclusions and outlook

In this paper we have extended a useful and accurate approach for solving
nonlinear ordinary DEs with sources, based on the BLUES function method
introduced in Indekeu and Müller-Nedebock (2018) and first applied
quantitatively in Berx and Indekeu (2019). We have shown that the method can
be applied to obtain useful approximations for the oscillating wave solution
of the damped nonlinear oscillator and for traveling wavefront solutions of
the Fisher equation. It can also be extended to the arena of fractional DEs,
as was shown for a heat transfer problem. In problems of wave propagation,
typically a partial DE is transformed to an ordinary DE by means of a
traveling wave Ansatz and the source must be assumed to be co-moving. This
limitation has, however, been removed by considering a broader class of
physical problems, as we have illustrated in the case of the nonlinear
ordinary heat transfer DE.

In the application, in Section IV, to the heat transfer problem we have made a
detailed comparison of our iteration procedure with the Adomian decomposition
method. It is found that both methods generate equal numbers of exact
coefficients in a given order of the series expansion (ADM) or of the
iteration (BLUES). However, the BLUES method generates an exponentially large
number of approximate coefficients of higher-order terms that are not present
in the ADM. The presence of these higher-order terms appears to accelerate the
convergence of the approximation.

One can think of further applications of the method. For an overview of
pertinent types of DEs for which this method could be tried, see for example
Wazwaz (2009). We have in mind, for example, systems described by the
nonlinear Schrödinger (or Gross-Pitaevskii) equation. This DE has been studied
intensively in the past decades and is of general interest in the research on
Bose-Einstein condensates, nonlinear optics and semiconductor physics Yang
_et al._ (2015); Yan  _et al._ (2011); Kippenberg  _et al._ (2018); Ma  _et
al._ (2017). In the field of fluid mechanics, when studying solitary waves on
shallow water, it would be interesting to see whether the BLUES function
method can be used to calculate physically relevant profiles of tidal bores in
the context of the recent minimal analytic model for this problem Berry (2018,
2019). Other opportunities arise when one considers different kinds of DEs,
e.g., nonlinear partial DEs, either deterministic or stochastic Kardar  _et
al._ (1986); Takeuchi (2018) and systems of coupled DEs Horowitz and Kardar
(2019).

Finally, we announce that the method is also applicable to nonlinear partial
differential equations (PDEs), for example in time $t$ and space $x$, which
cannot be reduced to ordinary DEs. For PDEs the initial condition can be
employed as the source term in the method. This is the subject of a
substantial further development on which we will report in a future
communication.

##  Appendix A First-order approximant for the Fisher equation

The zeroth-order approximant to the solution of the nonlinear Fisher equation
(34) with exponential corner source (12) is

$U_{\psi}^{(0)}(z)\equiv(B\ast\psi)(z)=\frac{1}{2}\begin{cases}2-\frac{K}{K+k}\mathrm{e}^{z/K},&z<0\\\
\frac{K}{K-k}\mathrm{e}^{-z/K}-\frac{2k^{2}}{K^{2}-k^{2}}\mathrm{e}^{-z/k},&z\geq
0,\\\ \end{cases}$ (37)

while for $K=k$ the convolution product results in

$U_{\psi}^{(0)}(z)\equiv(B\ast\psi)(z)=\begin{cases}1-\frac{\mathrm{e}^{z/k}}{4},&z<0\\\
\left(\frac{3}{4}+\frac{z}{2k}\right)\mathrm{e}^{-z/k},&z\geq 0.\\\
\end{cases}$ (38)

To calculate the first-order approximation to the solution of the Fisher
equation (34), the residual operator (35) is applied to the zeroth-order
approximant in accordance with the iteration equation (5)

$\mathcal{R}_{z}U_{\psi}^{(0)}(z)=-\frac{k}{4}\begin{cases}\frac{Ke^{z/K}\left(Ke^{z/K}-2(k+K)\right)}{(k+K)^{2}},&z<0,\\\
\frac{e^{-2z\left(\frac{1}{k}+\frac{1}{K}\right)}\left(K(k+K)e^{z/k}-2k^{2}e^{z/K}\right)\left((k+K)e^{z/k}\left(K-2(K-k)e^{z/K}\right)-2k^{2}e^{z/K}\right)}{\left(k^{2}-K^{2}\right)^{2}},&z\geq
0\\\ \end{cases}$ (39)

Now a convolution product of the BLUES function (33) with the residual (39) is
calculated. The increment between the zeroth-order and first-order
approximants $\Delta U_{\psi}^{(1,0)}=U_{\psi}^{(1)}-U_{\psi}^{(0)}$ is
presented here, which is the convolution product $\Delta
U_{\psi}^{(1,0)}(z)=(B\ast\mathcal{R}_{z}U_{\psi}^{(0)})(z)$

$\Delta
U_{\psi}^{(1,0)}(z)=\begin{cases}\frac{k}{8\alpha^{2}}\left(\frac{K^{4}}{2k+K}\mathrm{e}^{2z/K}-4K^{3}\mathrm{e}^{z/K}+2(2k^{3}+4k^{2}K+6kK^{2}+3K^{3})\right),&z<0,\\\
\frac{k}{8\alpha^{2}\beta^{2}}\left(4k^{5}\mathrm{e}^{-2z/k}+\frac{8k^{2}K^{3}(5k^{2}+K^{2})}{\gamma}\mathrm{e}^{-z/k}-8k^{2}K^{3}\mathrm{e}^{-\alpha
z/kK}\right.&\\\
\left.+4K^{3}\alpha^{2}\mathrm{e}^{-z/K}+\frac{K^{4}\alpha^{2}}{2k-K}\mathrm{e}^{-2z/K}\right),&z\geq
0,\\\ \end{cases}$ (40)

where $\alpha=k+K$, $\beta=k-K$ and $\gamma=K^{2}-4k^{2}$. The cases for which
$K=k$ or $K=2k$ must be treated separately but these (elementary) calculations
will not be performed here. Note that the first-order approximant
$U_{\psi}^{(1)}(z)$ approaches the constant value $U_{c}$ for
$z\rightarrow-\infty$, which has already been presented in equation (36).

## References

  * Adomian (1994) G. Adomian,  _Solving Frontier Problems of Physics: The Decomposition Method_, 1st ed., Fundamental Theories of Physics, Vol. 60 (Springer Netherlands, 1994).
  * Liao (2012) S. Liao,  _Homotopy Analysis Method in Nonlinear Differential Equations_, 1st ed. (Springer-Verlag Berlin Heidelberg, 2012).
  * Karpman and Maslow (1977) V. I. Karpman and E. M. Maslow, “Perturbation theory for solitons,” J. Exp. Theor. Phys. 73, 537–559 (1977).
  * Keener and McLaughlin (1977) J. P. Keener and D. W. McLaughlin, “Solitons under perturbations,” Phys. Rev. A 16, 777–790 (1977).
  * Maeda and Colonius (2017) K. Maeda and T. Colonius, “A source term approach for generation of one-way acoustic waves in the Euler and Navier-Stokes equations,” Wave Motion 75, 36 (2017).
  * Ermakov and Stepanyants (2019) A. Ermakov and Y. Stepanyants, “Soliton interaction with external forcing within the Korteweg-de Vries equation,” Chaos 29, 013117 (2019).
  * Berx and Indekeu (2019) J. Berx and J. O. Indekeu, “Analytic iteration procedure for solitons and traveling wavefronts with sources,” J. Phys. A-Math. Theor. 52, 38LT01 (2019).
  * Indekeu and Smets (2017) J. O. Indekeu and R. Smets, “Traveling wavefront solutions to nonlinear reaction-diffusion-convection equations,” J. Phys. A-Math. Theor. 50, 315601 (2017).
  * Indekeu and Müller-Nedebock (2018) J. O. Indekeu and K. K. Müller-Nedebock, “BLUES function method in computational physics,” J. Phys. A-Math. Theor. 51, 165201 (2018).
  * Ekomasov  _et al._ (2018) E. G. Ekomasov, A. M. Gumerov, R. V. Kudryavtsev, S. V. Dmitriev, and V. N. Nazarov, “Multisoliton dynamics in the Sine-Gordon model with two point impurities,” Braz. J. Phys. 48, 576 (2018).
  * Gul  _et al._ (2018) Z. Gul, A. Ali, and A. Ullah, “Localized modes in parametrically driven long Josephson junctions with a double-well potential,” J. Phys. A-Math. Theor. 52, 015203 (2018).
  * Starodub and Zolotaryuk (2019) I. O. Starodub and Y. Zolotaryuk, “Fluxon interaction with the finite-size dipole impurity,” Phys. Lett. A 383, 1419 – 1426 (2019).
  * Duan  _et al._ (2012) J.-S. Duan, R. Rach, D. Baleanu, and A. Wazwaz, “A review of the Adomian decomposition method and its applications to fractional differential equations,” Commun. Frac. Calc. 3, 73–99 (2012).
  * Wazwaz and Khuri (1996) A. Wazwaz and S. Khuri, “A reliable technique for solving the weakly singular second-kind Volterra-type integral equations,” Appl. Math. Comput. 80, 287 – 299 (1996).
  * Padmavally (1958) K. Padmavally, “On a non-linear integral equation,” J. Math. Mech. 7, 533–555 (1958).
  * Keller and Olmstead (1972) J. Keller and W. Olmstead, “Temperature of a nonlinearly radiating semi-infinite solid,” Q. Appl. Math. 29, 559–566 (1972).
  * Fisher (1937) R. Fisher, “The wave of advance of advantageous genes,” Ann. Eugen. 7, 355–369 (1937).
  * Wazwaz (2009) A.-M. Wazwaz,  _Partial Differential Equations and Solitary Waves Theory_, Nonlinear Physical Science (Springer Berlin Heidelberg, Berlin, Heidelberg, 2009).
  * Yang  _et al._ (2015) Y. Yang, Z. Yan, and D. Mihalache, “Controlling temporal solitary waves in the generalized inhomogeneous coupled nonlinear Schrödinger equations with varying source terms,” J. Math. Phys. 56, 053508 (2015).
  * Yan  _et al._ (2011) Z. Yan, X.-F. Zhang, and W. M. Liu, “Nonautonomous matter waves in a waveguide,” Phys. Rev. A 84, 023627 (2011).
  * Kippenberg  _et al._ (2018) T. J. Kippenberg, A. L. Gaeta, M. Lipson, and M. L. Gorodetsky, “Dissipative Kerr solitons in optical microresonators,” Science 361 (2018).
  * Ma  _et al._ (2017) X. Ma, O. A. Egorov, and S. Schumacher, “Creation and manipulation of stable dark solitons and vortices in microcavity polariton condensates,” Phys. Rev. Lett. 118, 157401 (2017).
  * Berry (2018) M. V. Berry, “Minimal analytical model for undular tidal bore profile; quantum and Hawking effect analogies,” New J. Phys. 20, 053066 (2018).
  * Berry (2019) M. V. Berry, “Minimal model for tidal bore revisited,” New J. Phys 21, 073021 (2019).
  * Kardar  _et al._ (1986) M. Kardar, G. Parisi, and Y.-C. Zhang, “Dynamic scaling of growing interfaces,” Phys. Rev. Lett. 56, 889–892 (1986).
  * Takeuchi (2018) K. A. Takeuchi, “An appetizer to modern developments on the Kardar-Parisi-Zhang universality class,” Physica A 504, 77 (2018).
  * Horowitz and Kardar (2019) J. M. Horowitz and M. Kardar, “Bacterial range expansions on a growing front: Roughness, fixation, and directed percolation,” Phys. Rev. E 99, 042134 (2019).


# Latency optimized Deep Neural Networks (DNNs): An Artificial Intelligence
approach at the Edge using Multiprocessor System on Chip (MPSoC)

Seyed Nima Omidsajedi2, Rekha Reddy2, Jianming Yi2, Jan Herbst2, Christoph
Lipps2 and Hans Dieter Schotten12  2Intelligent Networks Research Group,
German Research Center for Artificial Intelligence
D-67663 Kaiserslautern, Email<EMAIL_ADDRESS>1Institute for
Wireless Communication and Navigation, University of Kaiserslautern
D-67663 Kaiserslautern, mail<EMAIL_ADDRESS>
###### Abstract

Almost in every heavily computation-dependent application, from 6G
communication systems to autonomous driving platforms, a large portion of
computing should be near to the client side. Edge computing (AI at Edge) in
mobile devices is one of the optimized approaches for addressing this
requirement.

Therefore, in this work, the possibilities and challenges of implementing a
low-latency and power-optimized smart mobile system are examined. Utilizing
Field Programmable Gate Array (FPGA) based solutions at the edge will lead to
bandwidth-optimized designs and as a consequence can boost the computational
effectiveness at a system-level deadline. Moreover, various performance
aspects and implementation feasibilities of Neural Networks (NNs) on both
embedded FPGA edge devices (using Xilinx Multiprocessor System on Chip
(MPSoC)) and Cloud are discussed throughout this research. The main goal of
this work is to demonstrate a hybrid system that uses the deep learning
programmable engine developed by Xilinx Inc. as the main component of the
hardware accelerator. Then based on this design, an efficient system for
mobile edge computing is represented by utilizing an embedded solution.

###### Index Terms:

Edge computing, Multiprocessor System on Chip (MPSoC), Edge Artificial
Intelligence (Edge AI), Deep Neural Networks (DNN), Deep learning Processing
Unit (DPU).

##  I Introduction: Edge computing platforms

With new emerging technologies in the field of communication and mobile
systems, demands for smart and secure platforms are arising in real-world
situations. Moreover, increasing the complexity of target applications with
exclusive purposes leads to a wide spectrum of research focusing on
constructing smart systems employing Artificial Intelligence (AI) solutions.
For large-scale application models, the need for parallel computational
capabilities with higher performance raises a lot of design challenges. For
instance, Raina, Madhavan, and Ng [1] discuss the implementation of
unsupervised learning methods over the modern graphics processor which scales
70 times faster than a dual-core Central Processing Unit (CPU). The
development of Cloud computing technologies has enabled the provision of
higher computational capabilities. Thus, the custom methodology for creating
powerful AI systems is using Cloud services due to their extensive computation
power. However, the traditional Cloud solutions cannot provide the desired
efficiency for daily increasing mobile systems [2][3]. As a result, a new
computing architecture based on end devices (also known as IoT devices), edge
devices, and Cloud has already been introduced and is used to improve the
performance of systems that are extremely computationally dependent [4].

Cloud servers could be used as a very powerful computing tool, particularly
for running AI inference systems; nevertheless, due to the significant data
transmission between the end device and the Cloud, the issue of reliability
arises. Moreover, the available bandwidth in a network is restricted,
therefore expanding the size of the network may result in further bandwidth
reductions concerning the increased data path. In addition, the data transfer
within a network can compromise the data security and as a result, endanger
the system security [5].

A commonly used choice as a tool for edge computing is utilizing Graphics
Processing Unit (GPU) systems; this is due to the strong computation
capabilities and also the high memory bandwidth of GPU devices. However, these
systems have a lot of issues for battery-dependent systems, especially for
mobile applications. Two of the most important design parameters for every
edge device are i) the energy efficiency and ii) the security of the design.
GPU systems are designed for graphical purposes and would not provide a fully
reliable smart system. One considerable solution for implementing bandwidth-
optimized, energy-efficient, and secure devices is using ARM-FPGA hybrid
systems. These systems consist of ARM processors as the heart of software
processing which act as the Central Control Unit (CCU) alongside a dedicated
hardware accelerator on FPGA blocks for boosting the desired Machine Learning
(ML) algorithm using the concept of parallelism. In 1992, [6] demonstrated
that a large neural network application with a single ANNA chip ( Analog
Neural Network Arithmetic and logic unit (ANNA)) enables it to operate with
the speed advantage of 50 to 500 over conventional hardware evaluated with
Floating-Point (FP) precision.

In the current paper, an FPGA based edge solution is presented in comparison
to the GPU methods over the edge. The remainder of this work is organized as
follows. Section II provides the current state of the art for the embedded
FPGA’s using MPSoCs. The detailed structure of our proposed system with a
dedicated hardware accelerator is presented in Section III. The implementation
comparison of embedded FPGA and GPU (Edge vs Cloud) are discussed in Section
IV. An outline of future work and this work is concluded in Section V.

##  II Related works: Using Embedded FPGAs as Edge computing systems

First ideas of reducing network traffic by the use of edge-computation are
provided by [4]. Therefore they tackled the problem of crowdsourced
applications which are geo-distributed globally with highly heterogeneous
crowdsourced data with the use of edge computing. Especially they want to
combat the challenges in crowdsourced deep learning applications. They claim
to be able to reduce the network traffic by 80 % and running time by 69 % in
comparison to state-of-the-art cloud solutions. One of the big advantages is
the improvement of network latency between end-user devices and edge servers
compared to cloud servers. [7] and [8] try to solve the problem of high
computation-intensive Deep Neural Network (DNN) based tasks on mobile devices
with the development of a specific Framework, which claims to use a device-
edge-synergy. Therefore they use algorithms to adaptively partition
computations between device and edge. Furthermore, they are reducing computing
latency via early exiting inference provided by the use of BranchyNet [9] at
the edge.

[10] examines latency and power optimizations of heavily computation-dependent
applications on the client side with the use of FPGAs. They investigate the
efficient hardware implementation of Cellular Neural Networks (CeNNs) in FPGA
with the help of different optimization algorithms for compressing image
processing techniques. By using a so-called incremental quantization which
includes parameter partition, parameter quantization, and re-training, they
can quantize the numbers in CeNN templates to the power of two, so they are
able to use the full potential of embedded FPGAs. This solves the problem of
the very high need for multiplications needed in CeNNs, which leads to a
bottleneck in an FPGA because of the limited quantity of embedded multipliers
one contains. By reducing the power of numbers to two the multiplications can
be solved with the use of logic shifts, which can be achieved with logic
elements and registers. In the represented work CNNs are used instead of
CeNNs, moreover, it will be the Xilinx deep learning IP (Xilinx DPU) which
delivers sufficient speed up for inference tasks.

Also [11] worked on more efficient ways of implementing efficient DNNs on
FPGAs. Therefore they developed a Framework named DNNWEAVER which can
automatically generate synthesizable accelerators for a given pair of DNN and
FPGA. With performance tests, they compared three different FPGAs (Xilinx
Zynq, Altera Stratix V, Altera Arria 10) against different many-core GPUs.
They were able to achieve higher Performance-per-Watt against all tested many-
core GPUs for the Xilinx Zynq and Altera Arria 10 FPGAs by using the framework
synthesized DNNs.

In [12] the authors went one step further and claim to create the first hybrid
system of GPU-FPGAs for training a network. Therefore they propose a new
framework with a focus on effective DNN training on GPU-FPGA-hybrids. They
used different energy-efficient schemes for DNN training procedures to execute
individual training operations on either high-performance (GPUs) or power-
efficient hardware (FPGAs). With that, they claim to be able to reduce the
average power consumption by 44.3 %. All in all, it is to say the higher
adaptability of FPGAs and the higher more advanced processing power of
embedded FPGAs compared to normal GPUs result in lower cost and relatively
lower power consumption which make them a serious threat to pure GPU solutions
still used in many classical cloud-based-solutions so far.

##  III Proposed system: Hardware accelerator on embedded FPGA for DNNs

This Section describes the proposed design with a dedicated hardware
accelerator on FPGA using Programmable Logic (PL) system of Xilinx MPSoC
evaluation board (Xilinx ZCU102). The heart of this accelerator includes
Xilinx Deep Learning Processing Unit (DPU) [13] which is responsible for
computing given inference tasks for a DNN system. The main functionality of
this system is to calculate massive series of Multiply–Accumulate operations
(MAC) on network parameters of DNNs. On the other hand, the initialization
task, as well as data coordination of a DNN, will be done using ARM processors
in the Processing System (PS) of the MPSoC.

The overall workflow for preparing and implementing a Neural Network (NN) on a
target embedded FPGA board (Xilinx SoC/MPSoC) is depicted in Figure 1. This
diagram consists of three sub-workflows. The most left path (depicted in dark
gray) explains the required steps for preparing a dedicated hardware
accelerator on FPGA logic gates using the Xilinx DPU soft IP. This path also
includes all the desired pre-processing at the hardware level. Then using the
implemented design in the PL side, computational processes of a NN, which are
mainly multiplication and addition functions, are performed on the FPGA side.
The design strategy behind the DPU IP is utilizing Digital Signal Processor
(DSP) tiles for conducting MAC (Multiply-Accumulate) operations.

Figure 1: Overall workflow for implementing NNs on the target SoC/MPSoC board

The middle sub-workflow (depicted in light grey) describes the required steps
for converting the desired NN from its initial form to an equivalent FPGA
compatible form. This step involves model quantization (converting floating-
based weights and biases to equivalent integer numbers). The quantization step
is very important not only because it reduces the required storage memory but
also optimizes the required bandwidth for the data transmission between the
off-chip memory (DDR memory connected to the PS side) and the hardware
accelerator on the FPGA side. After the quantization step, the graph-based NN
must be converted to a binary form which is executable on logic gates of the
PL subsystem. In this step, the pre-quantized model is fed into the Xilinx
compiler for generating the xmodel file of the initial NN.

For controlling the data transmission among various parts of the SoC/MPSoC as
well as for servicing the interrupts, a user-defined Linux application should
be used. This part is shown in the most right sub-workflow (depicted in
white). Based on the requirements of the project, there may be a need for
cross-compiling the given application (C++ application).

Finally, by booting the embedded Linux (PetaLinux) with a given NN on the
target board, further experiments are possible. In this case, the classic
ResNet50 [14] DNN is used which is a pre-trained classification network on
$224\times 224\times 3$ ($224\times 224$ pixels with 3 color channels) images
from ImageNet dataset [15]. The current version of the ResNet50 network can
classify 1000 different classes based on given input images.

The overall design scheme for the hardware accelerator used in this paper is
depicted in Figure 2. This design includes DPU core(s) used for the inference
computation of the target NN. The DPU core is responsible for the computation
of pixels in the feature map of input images. The number of Processing
Elements (PEs) inside the convolution engine is equivalent to the pixel
parallelism of the feature map. The implemented NN parameters are stored in
off-chip memory (DDR4 memory) and will be transmitted to the PL system with
the Full Power Domain (FPD) AXI connections. For boosting the performance of
the accelerator, connections between hardware logic in FPGA and off-chip
memory are built through a Direct Memory Access (DMA) system. In this case, PL
can transmit data directly to or from external memory. The ARM processors in
the APU subsystem are responsible for data coordination as well as required
initialization for given tasks to DPU core(s).

Figure 2: The hardware accelerator design for implementing a DNN on the target
MPSoC

##  IV Performance results of DNN: On Edge versus on Cloud

The DNN implementations over the Edge and Cloud are presented in this Section.
Moreover, the detailed description of the environment setup for used Edge and
Cloud components, FPGA and GPU features, and the model used for the execution
of the DNN are mentioned in related parts.

###  IV-A Deep Neural Network Implementation over Edge

The implementation details of a DNN on an edge-based device are described in
this Subsection.

####  IV-A1 Physical Infrastructure

The edge device used in this set-up is the ZCU102 evaluation board which is an
embedded FPGA from Xilinx Inc. This board is considered as a MPSoC with the
Zynq UltraScale+ architecture. However, the overall workflow presented in this
paper can be replicated in any Zynq and Zynq UltraScale+ based evaluation or
custom boards (for example, Zynq-7000 ARM-FPGA SoC boards). The main reason
for selecting a Zynq UltraScale+ system is the scalable capability of MPSoC
family boards, which is an important factor for creating smart networks.
Moreover, ZCU102 provides more logic resources in FPGA (including DSP blocks,
which are the main computational units) for performing required computations
of the inference graph of DNNs. Another reason for selecting the Zynq
UltraScale+ over the Zynq architecture is the possibility of implementing the
SoftMax activation function at the hardware level. The hardware implementation
of this activation function can boost the performance of the system much more
than the equivalent function at the software level. The ZCU102 consists of ARM
Cortex processors in the Application Processing Unit (APU) and Real-time
Processing Unit (RPU) sub-systems which are implemented in the PS. Moreover, a
mobile GPU is encompassed in the PS side of the board. The system
specifications of Xilinx ZCU102 MPSoC are described in Table I.

TABLE I: Embedded FPGA Hardware Specification Processing System (PS): |
---|---
APU | Quad-core ARM Cortex-A53
RPU | Dual-core ARM Cortex-R5
GPU | Mali-400
Programmable Logic (PL): |
System Logic Cells | 600 K
Block Ram (BRAM) | 32.1 MB
DSP Slices | 2520
Memory types: |
PS-Side DDR4 | 4GB 64-bit
PL-Side DDR4 | 512MB 16-bit

####  IV-A2 DNN Implementation over Edge device

The implementation results for ResNet50 DNN on the Xilinx ZCU102 MPSoC board
are shown in Table II. In this work, several test scenarios are implemented on
the target board using different sizes and numbers of DPU cores. The achieved
throughput and power consumption of the system depend on the actual
requirements of a project; however, in these test scenarios assumed a power
budget of around 10 to 20 Watts is available and also the target throughput of
the system is around 100 to 200 Images per Second (Img/S). Achieving values
higher than the proposed range is possible by utilizing more logic resources
on FPGA, or using more expensive and advanced embedded FPGAs, or providing a
bigger power supply with a proper cooling system. One important design
parameter for these experiments is the Convolution architectures of the DPU
core. According to the Xilinx documentation [13], this design parameter
determines the number of operations the DPU core can perform in one Clock
Cycle (CC). For example, 1/B4096 architecture equals 1 DPU core capable of
performing 4096 operations in one CC. The bigger the Convolution architectures
means higher the logic resources used for creating that DPU core. The maximum
number of implementable DPUs is limited due to the high DSP block usage in the
internal structure of this AI unit. The current version of this soft IP core
supports up to 4 homogeneous cores. Considering the results from these
experiments, the direct relationship between the power consumption and the DPU
core size is considerable. The number of test images used for each scenario is
equal to 2000 images.

TABLE II: Power, Latency and Throughput of the Edge device DPU Cores/Arch | Power(W) | Latency(S) | Throughput(Img/S)
---|---|---|---
1/B4096 | 9.97 | 22.59 | 88.5
2/B4096 | 16.29 | 12.40 | 161.2
3/B4096 | 22.86 | 9.73 | 205.4
4/B2304 | 21.98 | 11.68 | 171.1

###  IV-B Deep Neural Network Implementation over Cloud

A detailed description of DNN implementation constructed on a Cloud-based GPU
cluster infrastructure is discussed.

####  IV-B1 Physical Infrastructure

The experimental setup consists of a Cloud Computing cluster with GPU servers
namely GTX1080Ti and RTX A6000. Table III provides the overall GPU
specifications used in the implementation. Slurm [16] is used for managing
compute resources, scheduling jobs, and execution on worker nodes, it provides
fault-tolerant and highly scalable cluster management.

TABLE III: GPU Cluster Hardware Specification Specifications: | GTX1080Ti | RTX A6000
---|---|---
Architecture | Pascal | Ampere
GPU Memory(GB) | 11 | 48
GPU per node | 8 | 8
CPU per GPU | 5-9 | 12

####  IV-B2 DNN Implementation over GPU Cluster

With the deployed cluster environment, DNN is constructed using a
containerized environment using enroot, administrated, and monitored through
the Slurm cluster. The framework is built using Keras with TensorFlow [17]
backend to run the image identifier model using ResNet50. The ResNet50 model
employed in this study is the same pre-trained model used in the Edge
implementation. GPUs were subjected to iterative image prediction analysis and
throughput as well as the total power consumption are recorded in Table IV.
The implementation is tested over 2 different GPU servers working with GPU in
sets of 2 and 4 along with 2 CPUs. The CPU used alongside the GTX1080Ti GPU is
the Intel(R) Xeon(R) CPU E5-2630 v4 and the CPU used with the RTXA6000 GPU is
the AMD(R) EPYC 7F72 24-Core Processor. The overall latency is the total time
taken during the resource allocation, containerized environment creation,
initialization, and final image prediction. GPUs are allocated with a batch
size of 8, to obtain a minimalistic comparison to the edge device. With the
increased batch size, high throughput can be achieved with higher GPU
utilization by compromising to operate with huge power. Thus, the experiment
performed uses lower GPU capacity, providing a lower throughput. The total
number of test images for each experiment is equal to 2000.

Firstly, the system is tested with resources of 2 GPUs and 2 CPUs are
allocated and tested on servers RTXA6000 and GTX1080Ti, and throughput
obtained are 180 and 175 with the latency of 65s and 46s. Then, the allocated
number of GPUs was increased to 4, and throughputs of 187 and 137 were
obtained. One important fact about the power consumption of these systems is
the GPU utilization, the test scenario with 2/RTXA6000 has higher utilization
than 4/RTXA6000 and this is the reason for the higher power consumption.

TABLE IV: Power, Latency and Throughput of the GPU Cluster Num/GPU | Power(W) | Latency(S) | Throughput(Img/S)
---|---|---|---
2/RTXA6000 | 215.8 | 65 | 180
4/RTXA6000 | 174.4 | 83 | 187
2/GTX1080Ti | 95.9 | 46 | 175
4/GTX1080Ti | 138.4 | 68 | 137

###  IV-C Comparison results of DNN Implementation over Edge versus GPU
Cluster

Considering the obtained results from implementing our target DNN (Classic
ResNet50) over Edge (embedded FPGA device) versus GPU Cluster (as the Cloud
tool) shows drastically lower latency using Edge solution. Figure 3 shows the
latency comparison of these two methodologies. The main reason for lower
latency using embedded FPGAs is the negligible additional overhead for data
transmission and initialization in the target board. Therefore, in an Edge
device, the latency value equals the time it takes to compute the input data
on the board (computation latency). However, in a Cloud implementation, there
are other latency sources for data transmission over the Cloud, creation of
new instances, and initialization of values. In this work, all the additional
latency sources in the Cloud are called ”Transmission and initialization
latency”. Our experiments show high flexibility and low latency of Edge
devices using hardware accelerators built on FPGA logical units. Figure 4
shows the detailed latency sources using the Cloud tool.

Figure 3: Latency comparison between the Edge and Cloud implementations Figure
4: Latency sources in Cloud implementations

Moreover, using embedded FPGAs is more energy-efficient than GPUs and this
design parameter is a decisive point for mobile projects with limited power
sources. The noticeable point is, GPU clusters are initially able to achieve
much higher throughput than our test scenarios, but this is equal to a much
higher power consumption of the system. Increasing the batch size leads to
higher GPU utilization, and this will result in higher throughput. For
instance, in our experiment with increasing the batch size from 8 to 200, the
throughput will jump from 137 Img/S to 751.49 Img/S in the 4/GTX1080Ti;
however, the total power consumption changes from 138.4 Watts to 747 Watts.
With replicating the same experiment using the 4/RTXA6000 scheme, the
throughput jumps from 187 Img/S to 1886.13 Img/S, and at the same time overall
consumed power changes from 174.4 Watts to 899 Watts. Considering the power
budget as well as the practically required throughput of the DNN, our
experiments show better energy efficiency in all the test cases.

##  V Conclusion & Future Work

In this paper, we proposed a Deep Neural Network implementation strategy on
Edge devices using aMultiprocessor System on Chip. Our experiments show that
using Edge devices for Artificial Intelligence inferences has superiority in
comparison to the Cloud implementation of the same network not only in the
latency optimization but also in the energy efficiency of the system. One of
the main components used in the hardware accelerator of this work is the
Xilinx dedicated soft IP core for deep learning purposes (Xilinx DPU). The
presented system can boost the computation process of DNNs using Digital
Signal Processor tiles and logic gates on the FPGA side as well as using ARM
processing for controlling and coordination purposes. In future research, the
possibility of implementing multiple DNNs on MPSoCs will be examined to
observe the effect of co-implementation of DNNs. Moreover, another possibility
for future work is using the Xilinx Versal boards [18] as the next generation
of AI cores with dedicated AI engines instead of soft IP cores.

## References

  * [1] R. Raina, A. Madhavan, and A. Y. Ng, “Large-scale deep unsupervised learning using graphics processors,” in _Proceedings of the 26th annual international conference on machine learning_ , 2009, pp. 873–880.
  * [2] “Imt traffic estimates for the years 2020 to 2030,” International Telecommunication Union, M Series Mobile, radiodetermination, amateur and related satellite services Report ITU-R M.2370-0, 2015.
  * [3] W. Jiang, B. Han, M. A. Habibi, and H. D. Schotten, “The road towards 6g: A comprehensive survey,” _IEEE Open Journal of the Communications Society_ , vol. 2, pp. 334–366, 2021.
  * [4] Y. Huang, X. Ma, X. Fan, J. Liu, and W. Gong, “When deep learning meets edge computing,” in _2017 IEEE 25th international conference on network protocols (ICNP)_. IEEE, 2017, pp. 1–2.
  * [5] C. Lipps, P. Ahr, M. Strufe, and H. D. Schotten, “The physec thing: About trust and security in industrial iot systems,” _Journal of Information Warfare (JIW)_ , vol. 19, no. 3, pp. 35–49, 2020, iSSN Print: 1445-3312; ISSN Online: 1445-3347.
  * [6] E. Säckinger, B. E. Boser, J. M. Bromley, Y. LeCun, and L. D. Jackel, “Application of the anna neural network chip to high-speed character recognition,” _IEEE Transactions on Neural Networks_ , vol. 3, no. 3, pp. 498–505, 1992.
  * [7] E. Li, Z. Zhou, and X. Chen, “Edge intelligence: On-demand deep learning model co-inference with device-edge synergy,” _MECOMM 2018 - Proceedings of the 2018 Workshop on Mobile Edge Communications, Part of SIGCOMM 2018_ , pp. 31–36, 2018.
  * [8] E. Li, L. Zeng, Z. Zhou, and X. Chen, “Edge ai: On-demand accelerating deep neural network inference via edge computing,” _IEEE Transactions on Wireless Communications_ , vol. 19, no. 1, pp. 447–457, 2020.
  * [9] S. Teerapittayanon, B. McDanel, and H.-T. Kung, “Branchynet: Fast inference via early exiting from deep neural networks,” in _2016 23rd International Conference on Pattern Recognition (ICPR)_. IEEE, 2016, pp. 2464–2469.
  * [10] X. Xu, Q. Lu, T. Wang, Y. Hu, C. Zhuo, J. Liu, and Y. Shi, “Efficient hardware implementation of cellular neural networks with incremental quantization and early exit,” _ACM Journal on Emerging Technologies in Computing Systems (JETC)_ , vol. 14, no. 4, pp. 1–20, 2018.
  * [11] H. Sharma, J. Park, D. Mahajan, E. Amaro, J. K. Kim, C. Shao, A. Mishra, and H. Esmaeilzadeh, “From high-level deep neural models to fpgas,” in _2016 49th Annual IEEE/ACM International Symposium on Microarchitecture (MICRO)_. IEEE, oct 2016, pp. 1–12. [Online]. Available: http://ieeexplore.ieee.org/document/7783720/
  * [12] X. He, J. Liu, Z. Xie, H. Chen, G. Chen, W. Zhang, and D. Li, “Enabling energy-efficient dnn training on hybrid gpu-fpga accelerators,” in _Proceedings of the ACM International Conference on Supercomputing_ , 2021, pp. 227–241.
  * [13] _Zynq DPU v3.3: Product Guide_ , Xilinx, Inc.
  * [14] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2016, pp. 770–778.
  * [15] J. Deng, W. Dong, R. Socher, L.-J. Li, K. Li, and L. Fei-Fei, “Imagenet: A large-scale hierarchical image database,” in _2009 IEEE conference on computer vision and pattern recognition_. Ieee, 2009, pp. 248–255.
  * [16] “Slurm workload manager version 20.02.7.” [Online]. Available: https://slurm.schedmd.com/
  * [17] “Tensorflow,” https://www.tensorflow.org/, 2015.
  * [18] K. Vissers, “Versal: The xilinx adaptive compute acceleration platform (acap),” in _Proceedings of the 2019 ACM/SIGDA International Symposium on Field-Programmable Gate Arrays_ , 2019, pp. 83–83.

  *[FPGA]: Field Programmable Gate Array
  *[MPSoC]:  Multiprocessor System on Chip
  *[AI]: Artificial Intelligence
  *[CPU]: Central Processing Unit
  *[GPU]: Graphics Processing Unit
  *[CCU]: Central Control Unit
  *[ML]: Machine Learning
  *[ANNA]: Analog Neural Network Arithmetic and logic unit
  *[MPSoCs]:  Multiprocessor System on Chip
  *[DNN]: Deep Neural Network
  *[FPGAs]: Field Programmable Gate Array
  *[DNNs]: Deep Neural Network
  *[PL]: Programmable Logic
  *[DPU]: Deep Learning Processing Unit
  *[PS]: Processing System
  *[NN]: Neural Network
  *[DSP]: Digital Signal Processor