N. Naghshineh ET AL.

# On the use of asymptotically motivated gauge functions to obtain convergent
series solutions to nonlinear ODEs

NASTARAN NAGHSHINEH W. CADE REINBERGER, NATHANIEL S. BARLOW,∗ MOHAMED A.
SAMAHA  School of Mathematical Sciences, Rochester Institute of Technology,
Rochester, NY, USA  Department of Sciences and Liberal Arts, Rochester
Institute of Technology-Dubai, Dubai 341055, U.A.E.  ∗Corresponding author:
<EMAIL_ADDRESS>School of Mathematical Sciences, Rochester Institute of
Technology, Rochester, NY, USA  School of Mathematical Sciences, Rochester
Institute of Technology, Rochester, NY, USA  Department of Mechanical and
Industrial Engineering,Rochester Institute of Technology-Dubai, Dubai 341055,
U.A.E.  STEVEN J. WEINSTEIN  Department of Chemical Engineering, Rochester
Institute of Technology, Rochester, NY, USA  School of Mathematical Sciences,
Rochester Institute of Technology, Rochester, NY, USA

###### Abstract

We examine the power series solutions of two classical nonlinear ordinary
differential equations of fluid mechanics that are mathematically related by
their large-distance asymptotic behaviors in semi-infinite domains. The first
problem is that of the “Sakiadis” boundary layer over a moving flat wall, for
which no exact analytic solution has been put forward. The second problem is
that of a static air–liquid meniscus with surface tension that intersects a
flat wall at a given contact angle and limits to a flat pool away from the
wall. For the latter problem, the exact analytic solution—given as distance
from the wall as function of meniscus height—has long been known (Batchelor,
1967). Here, we provide an explicit solution as meniscus height vs. distance
from the wall to elucidate structural similarities to the Sakiadis boundary
layer. Although power series solutions are readily obtainable to the governing
nonlinear ordinary differential equations, we show that—in both problems—the
series diverge due to non-physical complex or negative real-valued
singularities. In both cases, these singularities can be moved by expanding in
exponential gauge functions motivated by their respective large distance
asymptotic behaviors to enable series convergence over their full semi-
infinite domains. For the Sakiadis problem, this not only provides a
convergent Taylor series (and conjectured exact) solution to the ODE, but also
a means to evaluate the wall shear parameter (and other properties) to within
any desired precision. Although the nature of nonlinear ODEs precludes general
conclusions, our results indicate that asymptotic behaviors can be useful when
proposing variable transformations to overcome power series divergence.
Sakiadis boundary layer; meniscus; asymptotic expansion; summation of series

##  1 Introduction

Infinite power series solutions to ordinary differential equations (ODEs) are
useful if they are convergent while also satisfying the given constraints over
the domain on which they are defined. It is within this context that we
develop a convergent series solution to a classical nonlinear ODE in fluid
mechanics–the Sakiadis boundary layer flow along a moving wall. To date, this
problem, as well as the related well-known Blasius problem along a stationary
wall, do not have proven exact analytical solutions in the literature,
although approximate analytical solutions have been put forward and are
discussed later in this section. The Sakiadis boundary layer is an important
flow field in configurations where thin liquid films are coated onto moving
substrates (Weinstein & Ruschak,, 2004), and is an essential component of
hydrodynamic assist in high speed curtain coating (Blake et al.,, 1994).

The boundary layer surrounding a flat plate moving through a viscous
incompressible fluid was first examined in the literature by Sakiadis, (1961)
who applied Blasius’s similarity transform to Prandtl’s boundary layer
equations (with appropriate boundary conditions, B.C.s) to arrive at a third-
order nonlinear ODE in $f(\eta)$ refered to here as the ‘Sakiadis Problem’,

$2f^{\prime\prime\prime}+ff^{\prime\prime}=0,~{}~{}~{}0\leq\eta<\infty$ (1a)
$\text{Sakiadis B.C.s: }~{}f(0)=0,~{}f^{\prime}(0)=1,~{}f^{\prime}(\infty)=0.$
(1b)

By contrast, the Blasius problem describing a stationary plate in a moving
fluid is governed by the same operator (1a) but has conditions

$\text{Blasius B.C.s: }~{}f(0)=0,~{}f^{\prime}(0)=0,~{}f^{\prime}(\infty)=1.$
(2)

While both the Sakiadis and Blasius problems can be handled in similar ways
numerically (e.g., shooting, transformation) (Cebeci & Keller,, 1971; Fazio,,
1992; Cortell,, 2010; Eftekhari & Jafari,, 2013; Fazio,, 2015), the difference
in boundary conditions leads to different (approximate) analytical approaches
(Barlow et al., 2017a, ); and, the nonlinear nature of the equations yields
distinctly different solutions. A common measure of the accuracy of any
solution technique applied to either problem is the quantity $\kappa$, defined
as

$\kappa\equiv f^{\prime\prime}(0),$ (3)

which is directly related to the wall shear stress in the boundary-layer flow,
and is typically referred to as the “wall shear” parameter (Cortell,, 2010;
Fazio,, 2015). For a given $\kappa$—which is known numerically and is
determined in this paper algorithmically—an infinite power series solution for
both problems can be obtained through standard means, given as

$f=\sum_{n=0}^{\infty}a_{n}\eta^{n},~{}~{}|\eta|<S,$ (4)

where $S$ is a finite radius of convergence. Singularities that lie outside of
the physical domain at a distance $S$ away from $\eta=0$ are the cause of this
radius. The Blasius problem has 3 well-known singularities lying the same
distance $S$ from the origin (see Boyd, (1999) and historical review therein)
and are reported to lie at values $\eta=S\text{ exp }[i(2j+1)\pi/3]$ ($j$=0,
1, 2) where $S\approx 5.6900380545$. For the Sakiadis problem, $S\approx
4.07217$, arising from singularities lying off the real-line in the left half-
plane (Barlow et al., 2017a, ). Approximate resummations are available, that
bypass the original series’ convergence barrier caused by these singularities
for both the Sakiadis (Barlow et al., 2017a, ) and Blasius (Boyd,, 1997, 1999;
Abbasbandy & Bervillier,, 2011; Barlow et al., 2017a, ) problems.

One approach for avoiding (non-physical) singularities that restrict series
convergence is to re-expand the series by mapping the independent variable
such that the (non-physical) closest singularities no longer affect the
physical domain. The divergent Blasius series has been successfully re-summed
in this way by Boyd, (1999) through re-expanding (4) as

$\text{Blasius:
}f=\sum_{n=0}^{\infty}\tilde{b}_{n}\left[\delta(\eta)\right]^{n},~{}~{}|\delta(\eta)|<\tilde{S},$
(5a) where the expansion variable
$\delta(\eta)=\frac{2\eta^{3}}{S^{3}+\eta^{3}}$ (5b)

is the gauge function by which (5a) may be cast as a formal Taylor series,
specifying how terms are asymptotically ordered (Van Dyke, 1975a, ; Leal,,
1992). In (5b), the three-fold symmetry of the closest singularities (with
modulus $S$) in the complex $\eta$-plane are mapped to infinity, due to the
effect of the above transformation on their orientation. This leads to a new
radius of convergence $\tilde{S}$ that lies outside the original physical
domain, thus creating a convergent series solution for the Blasius problem on
the positive real line. The coefficients $\tilde{b}_{n}$ are obtained by
equating (4) with the expansion of (5b) about $\eta=0$; these may be obtained
recursively (only depending on prior coefficients) due to the gauge function
(5b) and its derivatives equaling $0$ at $\eta=0$. The substitution made in
(5b) may be considered a modified Euler transformation111There are several
definitions of Euler transformations in the literature. The definitions and
interpretations that align most with this work are found in Van Dyke, 1975b ;
Baker, (1990).

While the series (5b) explicitly incorporates singularities in order to bypass
the radius of convergence of the original series (4), another approach is to
consider the other side of the domain ($\eta\to\infty$) as a possible
expansion point. Using the method of dominant balance, the $\eta\to\infty$
behaviors for the Blasius and Sakiadis solutions are given in Barlow et al.,
2017a  as

$\text{Blasius:
}f\sim\eta+B+4Q\frac{\exp[-\eta^{2}/4-B\eta/2]}{(\eta+B)^{2}}\left[1+O\left(\frac{1}{(\eta+B)^{2}}\right)\right],~{}\textrm{as}\
\eta\rightarrow\infty,$ (6) $\text{Sakiadis: }f\sim
C+Ge^{-C\eta/2}+\frac{G^{2}}{4C}e^{-C\eta}+O\left(e^{-3C\eta/2}\right),\
\textrm{as}\ \eta\rightarrow\infty,$ (7)

where $B$, $Q$, $C$, and $G$ are constants arising from integration. We repeat
the above expressions here to highlight a key difference and to motivate the
technique used to examine the Sakiadis problem in this paper. A necessary
condition for an expansion to converge is that it does not introduce new
singularities into the problem that lie in the physical domain. In the Blasius
problem, the constant $B$ takes on a negative value and thus (6) has a
singularity at the positive real value of $-B$. The expansion (7) for the
Sakiadis problem does not have this issue, and it is thus possible that the
expansion converges; this will be discussed in more detail shortly and
revisited throughout the paper.

Barlow et al., 2017a  used the asymptotic form (7) to propose an approximant
to sum the divergent series (4) for the Sakiadis problem as

$\text{Sakiadis:
}f_{A,N}=C+\sum_{n=1}^{N}\mathcal{A}_{n}\left(e^{-C\eta/2}\right)^{n}.$ (8)

In (8), the $\mathcal{A}_{n}$ coefficients are computed by solving an $N\times
N$ system enforcing that the $N$-term Taylor expansion of (8) about $\eta=0$
is exactly the $N^{\mathrm{th}}$ truncation of (4). Unlike the recursive
computation of $\tilde{b}_{n}$ in (5b), all $\mathcal{A}_{n}$ coefficients
change their value as $N$ changes, and thus (8) cannot be considered a formal
power series solution to (1b). Nevertheless, $f_{A,N}$ converges to the
numerical solution of $f$ over the physical domain as $N$ is increased. Note
that, although not remarked on in Barlow et al., 2017a , each
$\mathcal{A}_{n}$ coefficient converges to a specific value as $N$ is
increased. Additionally, (8) was used in Barlow et al., 2017a  to compute
$\kappa$, as well as the asymptotic constants in (7) ($C$ and $G$) to within
12 digits of accuracy (beyond previously reported numerical results) before
hitting a round-off barrier.

Note that, in using approximant (8), Barlow et al., 2017a  did not utilize
knowledge of the singularities that caused the original series (4) to diverge.
However, the convergence of (8) suggests that changing the gauge function from
$\eta$ to $e^{-C\eta/2}$ fortuitously creates a mappping that removes the
influence of the limiting singularities present in (4); this is similar in
behavior to the modified Euler transformation used in (5b) for the Blasius
problem. In the current work, we develop a formal Taylor series solution to
the Sakiadis problem using the exponential gauge transformation mentioned
above, with the goal of providing further explanation for the improved
accuracy of (8).

Interestingly, another classical problem in fluid mechanics that has a similar
mathematical structure is a static air–liquid meniscus (with surface tension)
formed when a flat wall is placed in an infinite horizontal pool. Both the
meniscus and the Sakiadis problem have (a) the same asymptotic decaying
exponential behavior at infinite distance from the wall and (b) a singularity
structure that limits convergence of the standard power series solution, yet
benefits from an exponential transformation. For the meniscus problem, the
exact solution—given as distance from the wall as function of meniscus
height—has long been known Batchelor, (1967). Thus, the emphasis of this 2nd
analysis is not in providing a new solution form, but rather in exploiting the
structural similarities with the Sakiadis problem to better understand the
nature of exponential gauge transformations and their effect on convergence
properties.

The paper is organized as follows. In section 2.1, we review the series and
approximant solutions to the Sakiadis boundary layer problem, expanding upon
the background given above. In section 2.2, an asymptotically motivated gauge
function is used to map the Sakiadis ODE and B.C.s to a domain in which its
series solution converges over the whole un-mapped physical domain. In section
2.3, we provide a procedure for using the analytic solution (from section 2.2)
to compute the asymptotic constants $C$ and $G$, the wall shear parameter
$\kappa$, and locations of the closest singularities (from $\eta$=0) of the
Sakiadis function to arbitrary precision. In section 3, the classical problem
of a meniscus at a flat wall is reviewed and its solution is obtained in
section 3.2 as a divergent power series expansion about the wall location. An
asymptotic expansion away from the wall and a convergent series solution
motivated by this expansion are provided in section 3.3, following a similar
technique to the Sakiadis problem in section 2.2. Remarks on the overall
methodology applied to both problems are provided in section 4.

##  2 The Sakiadis boundary layer

###  2.1 A divergent power series and a convergent approximant

Although the main results of this subsection are found in Barlow et al., 2017a
(and references therein), here we elucidate elements of that work that were
previously not discussed and are relevant to the work herein. The Sakiadis
boundary layer along a moving flat wall is obtained as a solution to the
nonlinear ODE (1b) where $f$ and $\eta$ are the similarity variables, related
to the physical variables in the governing boundary layer equation (Sakiadis,,
1961) through a similarity transform. The general power series solution to the
ODE (1ba) was first developed by Blasius, (1908), and is expressed as

$f(\eta)=\sum_{n=0}^{\infty}a_{n}\eta^{n},~{}~{}|\eta|<S$ (9a) where
$a_{n+3}=\frac{\displaystyle-\sum_{j=0}^{n}(j+1)(j+2)a_{j+2}a_{n-j}}{2(n+1)(n+2)(n+3)},~{}~{}n\geq
0.$ (9b)

For the Sakiadis problem, $a_{0}=0$ and $a_{1}=1$ as prescribed by (1bb), and
$a_{2}=\kappa/2$ as defined in (3). The “wall shear” parameter $\kappa$ is
computed either numerically (Cortell,, 2010) or analytically (Barlow et al.,
2017a, ) in a self-consistent way such that the condition
$f^{\prime}(\infty)=0$ in (1bb) is met. The value of $\kappa$ is given in
Barlow et al., 2017a  as

$\kappa=-0.443748313369\dots,$ (10)

the precision of which is improved upon in section 2.2.

The series (9b) diverges within the physical domain, as indicated in figure 1a
(dashed lines), where it is compared to the numerical solution ($\bullet$’s,
RK4 with $\Delta\eta=10^{-5}$). Using the modified Padé approximants given in
Barlow et al., 2017a , an estimate of the locations of the convergence-
limiting singularities in the complex plane are found222In Barlow et al.,
2017a , only the singularity at $-1.2114+3.8878i$ is reported. Revisiting the
analysis, we find two things: (1) The expected conjugate $-1.2114-3.8878i$ is
also found in the roots of the denominator of the modified Padé and (2) The
two additional digits of $\eta_{s}$ reported in Barlow et al., 2017a  (but not
here) are found to converge to different values as more terms are taken; a
more precise value is given here in section 2.3. to be at
$\eta_{s\pm}\approx-1.2114\pm 3.8878i$, indicating a radius of convergence of
$S=|\eta_{s\pm}|\approx 4.0722$, shown in figure 1a as a vertical line.

(a)  (b)

Figure 1: Analytic solutions of (1b) compared with the numerical solution (RK4
with $\Delta\eta=10^{-5}$) ($\bullet$’s). (a) $N$-term truncations of series
(9b) (dashed curves) and approximant (11) (solid curves). A solid vertical
line indicates the radius of convergence $S$ of (9b). (b) $N$-term truncations
of series (22d) (solid curves) and series (24b) (– - –) of section 2.2.

The series (9b) is analytically continued beyond this radius of convergence in
Barlow et al., 2017a  in an approximate way by constructing an approximant
based on the $\eta\to\infty$ expansion given by (7). As mentioned in section
1, a particular sequence of approximants having the form of the asymptotic
behavior (7) are examined in Barlow et al., 2017a  as

$f_{A,N}(\eta)=C+\sum_{k=1}^{N}\mathcal{A}_{k}\left(e^{-C\eta/2}\right)^{k}$
(11)

where the $\mathcal{A}_{k}\dots\mathcal{A}_{N}$ coefficients are computed such
that the expansion of (11) about $\eta=0$ exactly matches the Taylor series
solution (9b) to order $N$. Taylor expanding each term of (11) to order $N$,
switching the order of summation, and setting the result equal to (9b) leads
to a $N\times N$ Vandermonde system, whose explicit inversion formula is known
(Turner,, 1966). For example, at $N=2$, we have

$\left[\begin{array}[]{cccc}1&&&1\\\
1&&&2\end{array}\right]\left[\begin{array}[]{c}\mathcal{A}_{1}\\\
\mathcal{A}_{2}\end{array}\right]=\left[\begin{array}[]{c}0!~{}a_{0}-C\\\
1!~{}(-2/C)~{}a_{1}\end{array}\right],$ (12)

and for $N$=3 we have

$\left[\begin{array}[]{ccccccc}1&&&1&&&1\\\ 1&&&2&&&3\\\
1&&&4&&&9\end{array}\right]\left[\begin{array}[]{c}\mathcal{A}_{1}\\\
\mathcal{A}_{2}\\\
\mathcal{A}_{3}\end{array}\right]=\left[\begin{array}[]{c}0!~{}a_{0}-C\\\
1!~{}(-2/C)~{}a_{1}\\\ 2!~{}(-2/C)^{2}~{}a_{2}\end{array}\right],$ (13)

and so on. Note that, after inverting the matrices above, $\mathcal{A}_{1}$
and $\mathcal{A}_{2}$ attain different values in (12) than in (13). The
relevance of this issue is discussed in what follows. Also, although the
Vandermonde matrix has an explicit inversion formula, high precision
arithmetic (i.e., beyond double) must be used for large $N$ to avoid round-off
error.

In principle, if $\kappa$ and $C$ are known, (11) provides a complete
approximate solution. A slight modification can also be used to estimate
$\kappa$ and $C$ as well, wherein the $\mathcal{A}_{k}$ coefficients are
chosen to match the Taylor coefficients to order $N$, but $\kappa$ and $C$ are
chosen to make $\mathcal{A}_{N}=\mathcal{A}_{N-1}=0$ (effectively switching
$\mathcal{A}_{N}$ and $\mathcal{A}_{N-1}$ with $\kappa$ and $C$ as unknowns of
the $N\times N$ system). This is precisely what is done in Barlow et al.,
2017a  to compute the value of $\kappa$ given by (10) and to provide a $C$
value of

$C=1.61612544681\dots;$ (14)

the precision of $C$ is improved upon in section 2.2. After using the $\kappa$
and $C$ obtained via the above process, approximant (11) converges uniformly
over the domain as $N$ is increased, as shown in figure 1a (solid curves).
Curves for $f^{\prime}$ and $f^{\prime\prime}$, found analytically from (11),
are also shown to uniformly converge in Barlow et al., 2017a . To revisit the
discussion just above, it is found that, as the approximant (11) itself
converges and the values of $A_{k}$ change with increasing $N$, $A_{1}$ and
$A_{2}$ converge respectively to values

$A_{1}=G=-2.1313459241\dots$ (15)

and $A_{2}=G^{2}/(4C)$, matching the $\eta\to\infty$ expansion (7) obtained
from dominant balance. The other $A_{k}$ coefficients also converge, albeit
slower as $k$ increases, as the portion of the approximant that is used to
match the behavior at $\eta=0$ is pushed to higher-order terms.

Although the approximant (11) is a convergent series, its coefficients depend
on truncation $N$ and so it cannot be considered a formal Taylor series. Thus,
we cannot apply Taylor’s theorem to (11) to relate convergence and singularity
location in the context of switching between a gauge function of $\eta$ in
(9b) to a gauge function of $e^{-C\eta/2}$ in (11). In the next section, we
remove this coefficient dependence on truncation while retaining the
asymptotically motivated gauge function of $e^{-C\eta/2}$; in doing so, we
provide an explanation for why the approximant (11) converges without
explicitly incorporating singular behavior.

###  2.2 Variable transform and convergent power series solution

Motivated by the asymptotic expansion (7) and encouraged by the convergence of
(11), we define the variable transformations

$\displaystyle\omega(\eta)$ $\displaystyle=\exp\left(-\frac{C}{2}\eta\right)$
(16) $\displaystyle f(\eta)\equiv g(\omega(\eta)).$ (17)

Substituting (16) and (17) into (1a), applying the chain rule, and rearranging
terms, the transformed ODE becomes

$\frac{g}{C}\left(\omega\ddot{g}+\dot{g}\right)=\omega^{2}\dddot{g}+3\omega\ddot{g}+\dot{g},$
(18)

where $\dot{g}$ denotes a derivative with respect to $\omega$. The boundary
conditions (1b) become

$g(1)=0,~{}~{}\dot{g}(1)=\frac{-2}{C},~{}~{}\ddot{g}(1)=\frac{2}{C}+\frac{4\kappa}{C^{2}}$
(19)

where the third condition in (1b) here leads to $0=0$ and so the
transformation of condition (3) is written instead for $\ddot{g}(1)$.

In the usual way, we assume a solution to (18) of the form

$g(\omega)=\sum_{n=0}^{\infty}\tilde{a}_{n}\omega^{n},$ (20)

which is readily differentiated term-by-term to compute $\dot{g}$, $\ddot{g}$,
and $\dddot{g}$. After employing Cauchy’s product rule (Churchill,, 1948) to
handle the nonlinear term on the left-hand side of (18), the ODE becomes

$\frac{\omega}{C}\sum_{n=0}^{\infty}\left[\sum_{k=0}^{n}(k+1)(k+2)\tilde{a}_{k+2}\tilde{a}_{n-k}\right]\omega^{n}+\frac{1}{C}\sum_{n=0}^{\infty}\left[\sum_{k=0}^{n}(k+1)\tilde{a}_{k+1}\tilde{a}_{n-k}\right]\omega^{n}\\\
=\omega^{2}\sum_{n=0}^{\infty}(n+1)(n+2)(n+3)\tilde{a}_{n+3}\omega^{n}+3\omega\sum_{n=0}^{\infty}(n+1)(n+2)\tilde{a}_{n+2}\omega^{n}\\\
+\sum_{n=0}^{\infty}(n+1)\tilde{a}_{n+1}\omega^{n}.$ (21)

Equating constant terms on both sides of (21) leads to

$\tilde{a}_{0}=C.$ (22a) Equating $\omega^{1}$ terms on both sides of (21)
leads to $\tilde{a}_{1}=\tilde{a}_{1}$ and to match the notation of the
asymptotic expansion (7), we write $\tilde{a}_{1}=G.$ (22b) Equating
$\omega^{n}$ terms for $n\geq 1$ on both sides of (21) leads to
$\sum_{k=0}^{n-1}(k+1)(k+2)\tilde{a}_{k+1}\tilde{a}_{n-k-1}+\sum_{k=0}^{n}(k+1)\tilde{a}_{k+1}\tilde{a}_{n-k}=C(n+1)^{3}\tilde{a}_{n+1},$
which may be simplified (by appropriate shifts to the $k$ indices) as
$\sum_{k=1}^{n+1}k^{2}\tilde{a}_{k}\tilde{a}_{n-k+1}=C(n+1)^{3}\tilde{a}_{n+1}.$
Solving the above expression for $\tilde{a}_{n+1}$ and using (22a) leads to
the recursion
$\tilde{a}_{n+1}=\frac{1}{Cn(n+1)^{2}}\sum_{k=1}^{n}k^{2}\tilde{a}_{k}\tilde{a}_{n-k+1}.\
n\geq 1.$ (22c) Note that for $n=1$, (22c) leads to
$\tilde{a}_{2}=\displaystyle\frac{G^{2}}{4C},$ which matches the asymptotic
expansion (7). The transformed solution to (1b), written in terms of the
original variables is then
$f(\eta)=\sum_{n=0}^{\infty}\tilde{a}_{n}\left(e^{-C\eta/2}\right)^{n}.$ (22d)

As we can see, (22d) is a Taylor series solution to (1b) in terms of the gauge
function $e^{-C\eta/2}$. Although the implementation of (22d) requires
knowledge of $C$ and $G$, we will show in section 2.3 that these may be
computed (to any desired precision) directly from (22d) via the conditions
(19). Interestingly, (22d) appears to be a recursive way to compute the
asymptotic series (7), originally obtained in Barlow et al., 2017a  via the
method of dominant balance.

(a)  (b)

Figure 2: Absolute error between analytical solutions and the numerical
solution of (1b) shown in figure 1, here plotted versus transformed variable
$\omega$ to encompass the full physical domain $\omega\in[0,1]$. (a) $N$-term
truncations of series (9b) (dashed curves) and approximant (11) (solid
curves). (b) $N$-term truncations of series (22d) (solid curves) and series
(24b) (– - –) of section 2.2.

Using values of $C$ and $G$ determined in section 2.3, convergence of (22d) is
shown as solid lines in figure 1b. An error plot is shown in figure 2b,
indicating convergence to within machine precision over the full domain (shown
versus $\omega\in[0,1]$ for clearer representation of the limit
$\eta\to\infty$). For comparison, error from the divergent Sakiadis series
(9b) and approximant (11) are shown in figure 2a over the same domain.
Although not shown here, error for $f^{\prime}$ and $f^{\prime\prime}$, which
may be analytically obtained from (22d), follow the same error trends as in
figure 2b.

We now provide an explanation of why (22d) (and by extension (11)) converge to
the numerical solution of (1b). As alluded to in section 1, the singularities
that cause the original series (9b) to diverge are mapped via the gauge
function (16) such that their influence lies beyond the physical domain; this
is shown in figure 3 where the complex $\eta$ and $\omega$ planes are
compared. Here, we track the movement of the two singularities $\eta_{s,\pm}$
(shown by $\ast$’s in the figure) closest to $\eta=0$, whose locations are
predicted by the Padé analysis of Barlow et al., 2017a  (see section 2.1) and
later confirmed in section 2.3. In the $\eta$ plane of figure 3, the circle of
convergence for (9b) is drawn based on these two singularities. Note that this
circle intersects the positive real line at the same location as the radius of
convergence of series (9b) drawn in figure 1a, as is expected from Taylor’s
theorem, which only guarantees convergence of Taylor series within such
circles (Watson,, 1914). In the $\omega$ plane of figure 3, the circle of
convergence for (22d) is drawn, centered at $\omega=0$ and with mapped radius
$|e^{-\left(C\eta_{s,\pm}\right)/2}|$; note that this circle now extends
beyond the physical domain itself333To paraphrase Buzz Lightyear: “To
$\eta=\infty$ and beyond”, which explains why (22d) converges over the entire
physical domain.

(a)  (b)

Figure 3: The effect of the gauge transformation (16) on $\eta$-singularities
($\ast$) of the Sakiadis function, showing their placement in the 2nd and 3rd
quadrants of the complex $\eta$ plane (left) and their image on the negative
real $\omega$ line (right). The interior of the dashed circle shown in the
$\eta$ plane delineates the region of convergence of series (9b). The
interiors of the smaller and larger dashed circles shown in the $\omega$ plane
respectively delineate the regions of convergence of (22d) and (24b).

On the mapped $\omega$-plane shown in figure 3, note that, in addition to
drawing the circle of convergence corresponding to expansion (22d), a larger
circle corresponding to the expansion about $\omega=1$ is also drawn, which is
larger because the mapped singularities are farther away. This expansion point
corresponds to $\eta=0$ in the original domain. Using the same procedures
employed above to obtain the series about $\omega=0$, the series about
$\omega=1$ is defined as

$g(\omega)=\sum_{n=0}^{\infty}\hat{a}_{n}(\omega-1)^{n},$ (23)

where the coefficients are given as

$\displaystyle\hat{a}_{n+3}=$
$\displaystyle-\frac{2n+3}{n+3}\hat{a}_{n+2}-\frac{(n+1)^{2}}{(n+2)(n+3)}\hat{a}_{n+1}$
$\displaystyle+\frac{\displaystyle\sum_{k=0}^{n-2}(k+1)(k+2)\hat{a}_{k+2}\hat{a}_{n-k-1}+\sum_{k=0}^{n-1}(k+1)\left[(k+2)\hat{a}_{k+2}+\hat{a}_{k+1}\right]\hat{a}_{n-k}}{C(n+1)(n+2)(n+3)},\
n\geq 0$
$\displaystyle\hat{a}_{0}=0,~{}~{}\hat{a}_{1}=-2/C,~{}~{}\hat{a}_{2}=2\kappa/C^{2}+1/C$
(24a) and, writing (23) in terms of the original variables, we have
$f(\eta)=\sum_{n=0}^{\infty}\hat{a}_{n}\left(e^{-C\eta/2}-1\right)^{n},$ (24b)

which may be considered a series solution to (1b) written in terms of the
gauge function $\left(e^{-C\eta/2}-1\right)$. In figure 2b, the error
associated with (24b) (dashed line) is shown to reduce to machine precision as
additional terms are used in the series. Although it is difficult to identify
an advantage in using (24b) over (22d), as each naturally perform better near
their expansion point as they converge, a “best choice” can be chosen by
examining the rapidity of convergence. This can be determined by considering
the infinity norm of the error (over $\omega\in[0,1]$) versus truncation $N$,
shown in figure 4. In the figure, one can see that, for a fixed number of
terms $N$, the expansion about $\omega=1$ (given by (24b)) is more accurate
than that about $\omega=0$ (given by (22d)). For comparison, the approximant
(11) is also shown in figure 4 and outperforms both (22d) and (24b). The cost
in using (11) is computational, since either inverting an increasingly ill-
conditioned Vandermonde matrix or using a Vandermonde inversion algorithm
requires precision well-beyond double to avoid round-off error.

(a)  Figure 4: The infinity norm of the difference between convergent
analytical solutions and the numerical solution to (1b), taken over the domain
shown in figure 2 and plotted versus series truncation $N$. The legend is the
same as in previous figures, showing (11) series (lowermost solid curve),
series (24b) (middle – - – curve), and series (22d) (uppermost solid curve).
The trends shown here are the same for the norm of the error of $f^{\prime}$
and $f^{\prime\prime}$ using (11), (24b), and (22d).

###  2.3 Computation of $C$, $G$, $\kappa$, and singularities of the Sakiadis
function

Previous estimates for $C$, $G$, and $\kappa$ have been found by quadrature,
shooting, and approximant (11) (Eftekhari & Jafari,, 2013; Cortell,, 2010;
Barlow et al., 2017a, ). Here, we improve upon those estimates and provide
formulae for obtaining them to any desired precision, by considering the
application of conditions (19) to the convergent solution (22d), i.e.,

$\displaystyle g(1)=\sum_{n=0}^{\infty}\tilde{a}_{n}=0$ (25)
$\displaystyle\dot{g}(1)=\sum_{n=1}^{\infty}n\tilde{a}_{n}=\frac{-2}{C}$ (26)
$\displaystyle\ddot{g}(1)=\sum_{n=2}^{\infty}n(n-1)\tilde{a}_{n}=\frac{2}{C}+\frac{4\kappa}{C^{2}}.$
(27)

Noting that (25) and (26) are functions of $C$ and $G$ only (see (22d)), we
have 2 equations and 2 unknowns and one may employ any desired solver (here we
use Newton’s method). After finding $C$ and $G$, one may then use (27) to
explicitly compute $\kappa$. In preparation for Newton’s method, we rewrite
(25) and (26) as

$\displaystyle\phi_{0}(C,G)$
$\displaystyle=\sum_{n=0}^{\infty}\tilde{a}_{n}=0$ (28)
$\displaystyle\phi_{1}(C,G)$
$\displaystyle=\frac{2}{C}+\sum_{n=1}^{\infty}n\tilde{a}_{n}=0,$ (29)

such that the $\nu^{\mathrm{th}}$ Newton iterate becomes

$\left[\begin{array}[]{c}C_{\nu+1}\\\
G_{\nu+1}\end{array}\right]=\left[\begin{array}[]{c}C_{\nu}\\\
G_{\nu}\end{array}\right]-\left[\begin{array}[]{cc}\left(\frac{\partial\phi_{0}}{\partial
C}\right)_{\nu}&\left(\frac{\partial\phi_{0}}{\partial G}\right)_{\nu}\\\
\left(\frac{\partial\phi_{1}}{\partial
C}\right)_{\nu}&\left(\frac{\partial\phi_{1}}{\partial
G}\right)_{\nu}\end{array}\right]^{-1}\left[\begin{array}[]{c}\phi_{0,\nu}\\\
\phi_{1,\nu}\end{array}\right].$ (30a) The partial derivatives in (30a) may be
computed compactly by recognizing the pattern in the series (22d) as
$g=C+G\omega+\frac{1}{4}\frac{G^{2}}{C}\omega^{2}+\frac{5}{72}\frac{G^{3}}{C^{2}}\omega^{3}+\dots+\underbrace{a^{\prime}_{n}G^{n}C^{1-n}}_{\tilde{a}_{n}}\omega^{n}+\dots$
where $a^{\prime}_{n}$ is $\tilde{a}_{n}$ evaluated at $C=G=1$ (i.e., the
coefficient with $C$ and $G$ removed)444Although here we claim this form of
$C$ and $G$ dependence on $\tilde{a}_{n}$ by inspection, this can be verified
by making the substitutions $\omega=G\tilde{\omega}/C$ and $g=C\tilde{g}$ into
(18) and noting the $C$ and $G$ dependence vanishes.. Using the dependence
shown above, the partial derivatives needed for (30a) are found555In cases
where the parameter dependence is more complicated, this can always be done
recursively by directly partially differentiating the recurrence relation. by
first partially differentiating each coefficient with respect to $G$ or $C$,
noting that $a_{n}^{\prime}$ is a constant, and multiplying and dividing by
$G$ or $C$ to make the expression in terms of $\tilde{a}_{n}$. Thus we obtain
$\displaystyle\frac{\partial\phi_{0}}{\partial C}$
$\displaystyle=\frac{1}{C}\sum_{n=0}^{\infty}(1-n)\tilde{a}_{n}$
$\displaystyle\frac{\partial\phi_{0}}{\partial
G}=\frac{1}{G}\sum_{n=0}^{\infty}n\tilde{a}_{n}$
$\displaystyle\frac{\partial\phi_{1}}{\partial C}$
$\displaystyle=\frac{-2}{C^{2}}+\frac{1}{C}\sum_{n=1}^{\infty}n(1-n)\tilde{a}_{n}$
$\displaystyle\frac{\partial\phi_{1}}{\partial
G}=\frac{1}{G}\sum_{n=1}^{\infty}n^{2}\tilde{a}_{n}.$ (30b)

Using the method outlined above, the following improved estimates are obtained
beyond double precision using only a few Newton iterations (if previous
estimates are used as the initial guess) in (30) and a truncation of 50 terms
in the series of (30b):

$C=1.616125446804603717\dots$

$G=-2.1313459240475714821\dots$

Using the above values of $C$ and $G$, we add equations (26) and (27) to each
other, extract $\tilde{A}_{1}=G$ from (26), and compute $\kappa$ as

$\kappa=\frac{C^{2}}{4}\left[G+\sum_{n=2}^{\infty}n^{2}\tilde{a}_{n}\right],$
(31)

where, taking 50 terms in (31), we obtain

$\kappa=-0.443748313368861\dots$

Using the values of $C$ and $G$ given above, we may now explore the
convergence of (22d) with the aim of deducing the precise singularity
locations sketched in figure 3. As shown in figure 5, the radius of
convergence of (22d) can be deduced from the ratio test, and is approximately

$\lim_{N\to\infty}\left|\frac{\tilde{a}_{N}}{\tilde{a}_{N+1}}\right|=2.66149513\dots$

Relating the above to the original variables of the Sakiadis problem via (16),
we have

$\displaystyle\lim_{N\to\infty}\left|\frac{\tilde{a}_{N}}{\tilde{a}_{N+1}}\right|$
$\displaystyle=\left|e^{-C\left(\eta_{s,r}\pm i\eta_{s,i}\right)/2}\right|$
$\displaystyle=e^{-C\eta_{s,r}/2}$ (32)

where $\eta_{s,r}$ and $\eta_{s,i}$ are the respective real and imaginary
parts of the conjugate singularity pair closest to $\eta=0$ in the original
Sakiadis function. While (32) may be used to solve directly for $\eta_{s,r}$,
the imaginary part may be deduced by recognizing that (22d) appears (for all
terms investigated) to be an alternating series and thus the closest
$\omega$-singularity is on the negative real line666This follows from a
variable substitution $\tilde{\omega}=-\omega$ in (20), application of
Pringsheim’s theorem (Flajolet & Sedgewick,, 2009) to the resulting series,
and then mapping back to $\omega$. (indicated in figure 3b) such that

$\arg\left[e^{-C\left(\eta_{s,r}\pm i\eta_{s,i}\right)/2}\right]=\mp
C\eta_{s,i}/2=\mp(2n-1)\pi,~{}~{}n=1,2,3\dots.$ (33)

where, choosing $n=1$ in (33) leads to $\eta_{s,i}$ values consistent with
both the Padé analysis in Barlow et al., 2017a  and the radius of convergence
of (9b) indicated in figure 1 (when combined with $\eta_{s,r}$ above). Hence
we conjecture that $n$ equals 1 in (33) and, consequently, that the closest
singularities to $\eta=0$ in the Sakiadis problem are given, through use of
(32) and (33), by

$\displaystyle\eta_{s,\pm}$
$\displaystyle=\frac{-2}{C}\lim_{N\to\infty}\ln\left|\frac{\tilde{a}_{N}}{\tilde{a}_{N+1}}\right|\pm
i\frac{2\pi}{C}$ $\displaystyle\approx-1.211401\pm 3.887808i.$ (34)

(a)  Figure 5: The limiting behavior of the ratio test for series (22d) versus
truncation $N$

It is worth mentioning that, although we are able to use the convergent series
solution (22d) to deduce the locations of the convergence-limiting
singularities of the divergent series (9b), knowledge of the singularity (34)
does not enter into the formulation of our solution at any point, in contrast
with the modified Euler summation (5b) used for the Blasius problem.
Consequently (and fortunately), the convergence of (22d) does not rely on
retaining “enough” digits of (34).

In section 3, an exponential transformation similar to (16) is used to obtain
a convergent series from a divergent one. However, in contrast with the
Sakiadis problem, the meniscus problem of section 3 has an exact solution that
allows us to anticipate the location convergence-limiting singularity a
priori. Still, we shall show that, as in the problem above, precise knowledge
of this singularity is not required to obtain a convergent expansion.

##  3 Meniscus at a flat wall

###  3.1 Background and formulation

To examine the nature of the exponential gauge function transformation used in
section 2.2, we examine a second classical problem which has similar
structural features. In particular, we consider the well-studied problem of
the shape of a static meniscus rising above an infinite horizontal pool
(Batchelor,, 1967; Probstein,, 2003). The pool of liquid has density, $\rho$,
is subjected to gravity, $g$, and is in contact with air of negligible
density. A flat wall is placed vertically in the pool, and the liquid
intersects the wall with a contact angle $\theta$ (measured through the
liquid) as shown in figure 8b; for purposes of this study, we assume that
$\theta$ lies between $0$ and $\pi/2$. The location of the air–liquid
interface, with surface tension $\sigma$, is parameterized as $y=h(x)$ where
$x$ is the horizontal distance from the wall, and $y=0$ is the undisturbed
location of the interface as $x\to\infty$. For this configuration, the
Young–Laplace equation couples with the hydrostatic field to yield the
following dimensionless equations and boundary conditions (Batchelor,, 1967;
Probstein,, 2003):

$\bar{h}=\frac{\bar{h}^{\prime\prime}}{\left[1+\left(\bar{h}^{\prime}\right)^{2}\right]^{3/2}},$
(35a) $\bar{h}^{\prime}(0)=-\cot\theta,$ (35b) $\bar{h}(\infty)=0.$ (35c)

In (35c), derivatives in $\bar{x}$ are denoted with primes, and the over-bars
denote dimensionless variables defined as

$\bar{h}=\frac{h}{L},~{}\bar{x}=\frac{x}{L},~{}L=\sqrt{\frac{\sigma}{\rho
g}}.$ (36)

In (36), the characteristic length scale $L$ is the well-known capillary
length. Multiplying both sides of (35a) by $\bar{h}^{\prime}$, integrating,
and applying the boundary conditions (35b) and (35c), we obtain

$\bar{h}^{\prime}=-\left[\frac{1}{\left(1-\frac{1}{2}\bar{h}^{2}\right)^{2}}-1\right]^{1/2}$
(37a) with the constraint $\bar{h}(0)=\sqrt{2(1-\sin\theta)},$ (37b)

where (37b) represents the height of the interface at the wall as function of
the contact angle, $\theta\in[0,\pi/2]$. Note that, in (37a),
$\bar{h}^{\prime}\to-\infty$ as $\bar{h}\to\sqrt{2}$.

The exact inverse solution $\bar{x}(\bar{h})$ of (37b) is obtainable via
variable separation and integration (Batchelor,, 1967) and is given as

$\bar{x}=\cosh^{-1}\frac{2}{\bar{h}}-\cosh^{-1}\sqrt{\frac{2}{1-\sin\theta}}+\sqrt{2+2\sin\theta}-\sqrt{4-{\bar{h}^{2}}}.$
(38)

In what follows, we obtain an analytic solution for $\bar{h}(\bar{x})$
directly via series expansion. Our intention in doing so is not to replace
equation (38) in usage, but it is to elucidate structural similarities and
provide insights to the Sakiadis series solution in section (2.2). Equation
(38) is used in what follows to assess the accuracy of the series solutions
provided in sections 3.2 and 3.3.

###  3.2 Divergent power series solution

The standard power series solution of (37b) is found by assuming

$\bar{h}=\sum_{n=0}^{\infty}\alpha_{n}\bar{x}^{n},~{}~{}|\bar{x}|<|\bar{x}_{s}(\theta)|,$
(39a) substituting (39a) into (37b), using Cauchy’s product rule (Churchill,,
1948) to evaluate $\bar{h}^{2}$, and applying JCP Miller’s formula (Henrici,,
1956) for raising a series to a power to obtain the following recursions for
the coefficients:
$\displaystyle\alpha_{n+1}=\frac{-2}{n+1}d_{n},~{}\alpha_{0}=\sqrt{2(1-\sin\theta)}$
(39b) $\displaystyle
d_{n>0}=\displaystyle\frac{1}{n\widetilde{c}_{0}}\sum_{j=1}^{n}\left(\frac{3}{2}j-n\right)\widetilde{c}_{j}d_{n-j},~{}d_{0}=\widetilde{c}_{0}^{\frac{1}{2}},~{}\widetilde{c}_{n>0}=c_{n>0},~{}\widetilde{c}_{0}=c_{0}-\frac{1}{4},$
(39c) $\displaystyle
c_{n>0}=\displaystyle\frac{1}{n\widetilde{b}_{0}}\sum_{j=1}^{n}(-j-n)\widetilde{b}_{j}c_{n-j},~{}c_{0}=(\widetilde{b}_{0})^{-2},~{}\widetilde{b}_{n>0}=b_{n>0},~{}\widetilde{b}_{0}=b_{0}-2,$
(39d) $\displaystyle b_{n}=\displaystyle\sum_{j=0}^{n}\alpha_{j}\alpha_{n-j}.$
(39e)

In (39a), the series is stated to converge within the region
$|\bar{x}|<|\bar{x}_{s}(\theta)|$ with $\bar{x}_{s}(\theta)$ being the as-of-
yet undetermined closest singularity to $\bar{x}=0$. Without formal proof, we
conjecture that $\bar{x}_{s}(\theta)$ is the value that satisfies
$\bar{h}(\bar{x}_{s})=\sqrt{2}$ for a given $\theta$ in the exact solution
(38), causing $\bar{h}^{\prime}\to-\infty$ in (35a). The value of the
conjectured limiting singularity $\bar{x}_{s}(\theta)$ is known exactly in
closed-form from the substitution of $\bar{h}=\sqrt{2}$ into (38) and is given
as

$\bar{x}_{s}(\theta)=\cosh^{-1}\sqrt{2}-\cosh^{-1}\sqrt{\frac{2}{1-\sin\theta}}+\sqrt{2+2\sin\theta}-\sqrt{2},$
(40)

where we note that $\bar{x}_{s}\leq 0$ for $\theta\in[0,\pi/2]$ and, in fact,
$\bar{x}_{s}=0$ when $\theta=0$ (setting the radius of convergence to be 0).
We support the conjecture that $\bar{x}_{s}$ is the closest singularity to
$\bar{x}=0$ (for all $\theta$) in what follows.

Figure (6) shows $N$-term truncations of series (39) (dashed curves), compared
with the exact solution (38) ($\bullet$’s) for $\theta=\pi/4$. As expected the
series (39) agrees for small $\bar{x}$ but ultimately diverges at a finite
radius of convergence (indicated as a solid vertical line in the figure) and
given by $|\bar{x}_{s}(\pi/4)|\approx 0.3$ from (40).

(a)  Figure 6:  The solution to (37b) is shown for a contact angle of
$\theta=\pi/4$. The $N$-term truncations of the divergent series (39) (dashed
curves) and the convergent resummation (49e) (solid curves) are compared
against the exact solution (38) ($\bullet$’s). The solid vertical line shows
the radius of convergence, $|\bar{x}_{s}(\pi/4)|\approx 0.3$, computed from
(40).

Although we have not formally established that $\bar{x}_{s}$ is the closest
singularity to $\bar{x}=0$ in (37b) and thus is responsible for divergence,
evidence to support this conjecture is given by the Domb-Sykes plot (Domb &
Sykes,, 1957) in figure 7, where the magnitude of the ratios of coefficients
of (39) is shown to approach $|\bar{x}_{s}(\pi/4)|$ as $n\to\infty$ (i.e.,
$1/n\to 0$ in the figure). Although only shown in this section for
$\theta=\pi/4$, all permissible $\theta$ values lead to similarly divergent
series, limited by a radius of convergence of $|\bar{x}_{s}(\theta)|$. In
section 3.3, we apply the same type of transformation used in section 2.2 to
overcome this convergence barrier for all contact angles.

(a)  Figure 7: Domb-Sykes plot for (39) with $\theta=\pi/4$, indicating (via
the ratio-test) a radius of convergence of $|\bar{x}_{s}(\pi/4)|\approx 0.3$,
which is consistent with what is observed in (6).

###  3.3 Asymptotic expansion, variable transform, and convergent power series
solution

Now that we have a solution to (37a), given by (39) that—although ultimately
divergent beyond some positive $\bar{x}$-value—provides the correct interface
shape close to the wall, the asymptotic behavior away from the wall is
examined. To this end, the method of dominant balance Bender & Orszag, (1978)
is employed. To meet the constraint (35c), it is assumed that

$\overline{h}(\overline{x})\ll 1,$ (41a) and this also assures that its
derivatives are small in the limit. To lowest order, then, the governing
equation (35a) is approximated as:
$\overline{h}^{\prime\prime}\sim\overline{h},~{}~{}\overline{x}\to\infty.$
(41b)

The asymptotic solution to (41b) is

$\overline{h}\sim C_{0}e^{-\overline{x}},~{}~{}\overline{x}\to\infty,$ (42a)
where $C_{0}$ is an arbitrary constant. To obtain the next correction, we
assume the following expansion $\overline{h}\sim
C_{0}e^{-\overline{x}}+D(\overline{x}),~{}~{}\overline{x}\to\infty,$ (42b)
where $D(\overline{x})\ll C_{0}e^{-\overline{x}},~{}~{}\overline{x}\to\infty.$
(42c) Equation (42b) is substituted in to (35a), subdominant terms are
neglected, and the resulting linear equation is solved in accordance with the
asymptotic relation (42c) to obtain:
$D(\overline{x})\sim\frac{3}{16}C_{0}^{3}e^{-3\overline{x}},~{}~{}\overline{x}\to\infty.$
(42d)

The same process is repeated to generate higher order correction to (42b). We
obtain:

$\overline{h}\sim
C_{0}e^{-\overline{x}}+\frac{3}{16}C_{0}^{3}e^{-3\overline{x}}+\frac{3}{256}C_{0}^{5}e^{-5\overline{x}}+O(e^{-7\overline{x}}),~{}~{}\overline{x}\to\infty.$
(43)

The pattern of exponentials is thus evident in equation (43). We note here
that a similar exponential behavior is also observed in the Sakiadis problem
previously and provides a linkage between the problems in approach and
interpretation to follow.

The asymptotic solution (43) motivates us to transform (37b) to reflect the
exponential pattern of (43). This is achieved by transformations in both the
independent and dependent variables, given respectively as

$U(\bar{x})=e^{-2\bar{x}},$ (44a)
$H\left(U(\bar{x})\right)=\bar{h}(\bar{x})e^{\bar{x}}.$ (44b)

Substituting (44b) into (37b) leads to the transformed ODE:

$\left[H+2U\dot{H}\right]\left[UH^{2}-2\right]=-H\sqrt{4-UH^{2}},~{}~{}$ (45)

with transformed condition

$H(1)=\sqrt{2(1-\sin\theta)},$ (46)

where $\dot{H}$ denotes the derivative with respect to $U$. Although an exact
explicit solution to (45) cannot be found as $H(U)$, an exact implicit
solution in $H$ and $U$ can be found by separating variables, integrating
(45), and applying (46) to arrive at

$\sqrt{2+2\sin\theta}-\sqrt{4-UH^{2}}=\ln\frac{H\left[1+\sqrt{1+\sin\theta}\right]}{\sqrt{1-\sin\theta}\left[2+\sqrt{4-UH^{2}}\right]}.$
(47)

Equation (47) is used to extract the condition

$H(0)=\frac{4\sqrt{1-\sin\theta}~{}e^{-2+\sqrt{2+2\sin\theta}}}{\sqrt{1+\sin\theta}+\sqrt{2}},$
(48)

which implies $H(0)=C_{0}$ in the asymptotic solution (43). This is needed for
the series solution that follows.

Using the same procedure as in sections (2.2) and (3.2), the power series
solution of (45) (with condition (48)) is given by

$H=\sum_{n=0}^{\infty}A_{n}U^{n},$ (49a) where,
$A_{n>0}=\frac{-\displaystyle\sum_{j=0}^{n-1}A_{j}[(1+2j)\mathcal{C}_{n-j}+\mathcal{D}_{n-j}]}{\mathcal{D}_{0}+\mathcal{C}_{0}(1+2n)},~{}~{}A_{0}=\frac{4\sqrt{1-\sin\theta}~{}e^{-2+\sqrt{2+2\sin\theta}}}{\sqrt{1+\sin\theta}+\sqrt{2}},$
(49b)
$\mathcal{D}_{n>0}=\frac{1}{n\widetilde{\mathcal{B}}_{0}}\sum_{j=1}^{n}(\frac{3}{2}j-n)\widetilde{\mathcal{B}}_{j}{\mathcal{D}}_{n-j},~{}~{}\mathcal{D}_{0}=\widetilde{\mathcal{B}}_{0}^{\frac{1}{2}},~{}~{}\mathcal{C}_{n>0}=\mathcal{B}_{n-1},~{}~{}\mathcal{C}_{0}=-2$
(49c)
$\widetilde{\mathcal{B}}_{n>0}=-\mathcal{B}_{n-1},~{}~{}\widetilde{\mathcal{B}}_{0}=4,~{}~{}\mathcal{B}_{n}=\sum_{j=0}^{n}A_{j}A_{n-j}.$
(49d) Transforming back to $\bar{h}(\bar{x})$ space via (44b), our expansion
about $\bar{x}=\infty$ (i.e., $U$=0) is
$\bar{h}=e^{-x}\sum_{n=1}^{\infty}A_{n}\left(e^{-2\bar{x}}\right)^{n},$ (49e)

which, by construction, is consistent with the asymptotic ordering (43) as
$\bar{x}\to\infty$, and shows explicitly the exponential gauge function
$e^{-2\bar{x}}$.

Figure 6 shows $N$-term truncations of (49e) (solid curves) compared with both
the exact solution (38) and divergent series (39). The difference between
(49e) and the exact solution is not discernible on the scale of the figure for
any $N$ shown, and it is noteworthy to mention that the radius of convergence
of the original series (39) has been exceeded. The absolute error between
(49e) and the exact solution (38) is shown in figure 8a, which indicates
convergence (as $N$ increases) for a continuum of angles $\theta$, as
prescribed in figure 8b. Figure 8 is generated for the smallest possible
contact angle ($\theta=0$) and thus, by virtue of (37b) being an autonomous
ODE, contains interface shapes for all contact angles as shifted semi-infinite
domains; this is indicated in the figure. The maximum error occurs at the wall
and is shown versus $N$ for $\theta=0$ in figure 9.

(a)  (b)

Figure 8: (a) Absolute error between $N$-term truncations of series (49e)
(solid lines) and the exact solution to (37b) given by (38) for $\theta=0$.
Labeled vertical lines show the wall at various contact angles specified in
the adjacent figure. For $\theta=\pi/4$, the wall is at point $B$ and the
error between the divergent series (39) and (38) is shown (dashed line). (b)
“Master” solution to (37b) at $\theta=0$, illustrating the $\theta$ values at
points $A$, $B$, $C$, and $D$ for potential wall locations in (a). (a)  Figure
9: Maximum absolute error between $N$-terms truncations of (49e) and (38)
(occurring at $\bar{x}=0$ and $\theta=0$), plotted versus $N$.

As done in (2.2) for the Sakiadis problem, we now provide an explanation of
why the series (49e) converges to the exact solution (38) of the ODE (37b)
describing a meniscus at a flat wall. The answer again lies in the mapping
provided by the gauge function (here, (44a)), as shown in figure 10 where the
complex $\bar{x}$ and $U$ planes are compared. In the $\bar{x}$ plane of
figure 10, circles of convergence for (39) centered around $\bar{x}$=0 are
drawn for various $\theta$ values, based on the conjectured closest
singularity (to $\bar{x}$=0), $\bar{x}_{s}(\theta)$, using (40). For
$\theta=0$, no circle is drawn, since $\bar{x}_{s}(\theta)$ (which sets the
radius) is 0. Note, that these circles of convergence each intersect the
positive real line. In particular, for $\theta=\pi/4$, this intersection
occurs at the same location as the radius of convergence of series (39) drawn
in figure 6; this correspondence holds for all other values of $\theta$, as
expected from Taylor’s theorem. In the $U$ plane of figure 10, the circle of
convergence for (49e) is drawn, centered at $U=0$ and with mapped radius of
$|e^{-2\bar{x}_{s}}|$; note that this circle now extends to the boundary of
the physical domain for $\theta=0$ and extends beyond it for $\theta>0$, which
explains why (49e) converges over the entire physical domain for all values of
$\theta$.

(a)  (b)

Figure 10: The effect of the gauge transformation (44a) on the singularities
$\bar{x}_{s}$ given by (40) for various $\theta$ values, showing their
placement on the negative real $\bar{x}$ line ($\ast$, left) and their image
on the positive real $U$ line ($\circ$, right). The dashed circles shown in
the $\bar{x}$ plane deliniate the regions of convergence of series (39) for
various indicated $\theta$ values. The circles shown in the $U$ plane
delineate the regions of convergence of (49e) for the same $\theta$ values.

The circles of convergence in figure 10 are, of course, still conjectured
because, in constructing the mapping in figure 10, we are assuming that the
singularity $\bar{x}_{s}(\theta)$ (given by (40)) is the closest singularity
to $\bar{x}=0$ and that no other singularities map to $U$ singularities closer
to $U=0$. In addition to evidence given in section 3.2, further evidence that
supports this conjecture is provided in the Domb-Sykes plots in figure 11,
where the radius of convergence deduced from this numerical ratio-test is
consistent with the locations of the mapped singularities in figure 10.

(a)  Figure 11: Domb-Sykes plot for (49e), for $\theta=\pi/4$ (upper curve,
$\circ$) and $\theta=0$ (lower curve, $\bullet$), indicating radii of
convergence consistent with figure 10.

Finally, it is worth noting that—like in the Sakiadis problem—knowledge of
these singularities is not required to implement the exponential gauge
function transformation (44b) that leads to a convergent series solution.

##  4 Conclusions

In this work, we provide convergent power series solutions to the Sakiadis
boundary layer problem and the problem of a meniscus at a flat wall, by means
of transforming the original ODEs in terms of variable substitutions that are
motivated by the asymptotic expansions about $\infty$. In both cases, the
transformations map the dominant convergence-limiting singularities out of the
physical domain; also in both cases, convergence-limiting singularities do not
need to be known a priori but their locations are deduced nonetheless. For the
Sakiadis problem, this provides—in the absence of a formal proof of
exactness—a conjectured exact Taylor series representation of the solution
over the full physical domain. That said, for both the Sakiadis and meniscus
problems, the exponential gauge functions used handle singularities similarly
to achieve these demonstrably convergent series solutions.

Although the nature of nonlinear ODEs precludes general conclusions, our
results indicate that asymptotic behaviors can be useful to motivate gauge
functions to overcome power series divergence. Additionally, the approach used
here supports a growing body of literature (Barlow et al., 2017a, ; Barlow et
al., 2017b, ; Beachley et al.,, 2018; Belden et al.,, 2020; Barlow &
Weinstein,, 2020; Weinstein et al.,, 2020; Ramé et al.,, 2021; Reinberger et
al.,, 2022) underscoring the use of power series solutions as a viable method
for analytically solving nonlinear ODEs.

## References

  * Abbasbandy & Bervillier, (2011) Abbasbandy, S. & Bervillier, C. (2011) Analytic continuation of Taylor series and the boundary value problems of some nonlinear ordinary differential equations. Appl. Math. Comp., 218, 2179–2199.
  * Baker, (1990) Baker, G. A. (1990) Quantitative Theory of Critical Phenomenon, chapter 14: Padé approximants, algebraic aspects. Academic Press.
  * (3) Barlow, N. S., Stanton, C. R., Hill, N., Weinstein, S. J. & Cio, A. G. (2017a) On the summation of divergent, truncated, and underspecified power series via asymptotic approximants. Q. J. Mech. Appl. Math., 70(1), 21–48.
  * Barlow & Weinstein, (2020) Barlow, N. S. & Weinstein, S. J. (2020) Accurate closed-form solution of the SIR epidemic model. Physica D, 408, 132540:1–4.
  * (5) Barlow, N. S., Weinstein, S. J. & Faber, J. A. (2017b) An asymptotically consistent approximant for the equatorial bending angle of light due to Kerr black holes. Class. Quant. Grav., 34(135017), 1–16.
  * Batchelor, (1967) Batchelor, G. K. (1967) An Introduction to Fluid Dynamics, chapter 1: The physical properties of fluids. Cambridge.
  * Beachley et al., (2018) Beachley, R. J., Mistysyn, M., Faber, J. A., Weinstein, S. J. & Barlow, N. S. (2018) Accurate closed-form trajectories of light around a Kerr black hole using asymptotic approximants. Class. Quant. Grav., 35(20), 1–28.
  * Belden et al., (2020) Belden, E. R., Dickman, Z. A., Weinstein, S. J., Archibee, A. D., Burroughs, E. & Barlow, N. S. (2020) Asymptotic approximant for the Falkner-Skan boundary layer equation. Q. J. Mech. Appl. Math., 73(1), 36–50.
  * Bender & Orszag, (1978) Bender, C. M. & Orszag, S. A. (1978) Advanced Mathematical Methods for Scientists and Engineers I: Asymptotic Methods and Perturbation Theory. McGraw-Hill.
  * Blake et al., (1994) Blake, T. D., Clarke, A. & Ruschak, K. J. (1994) Hydrodynamic assist of dynamic wetting. AIChE J., 40(2), 229–242.
  * Blasius, (1908) Blasius, H. (1908) Grenzschichten in Flussigkeiten mit kleiner Reibung. Zeitschrift fur Mathematik und Physik, 56, 1–37.
  * Boyd, (1997) Boyd, J. P. (1997) Padé approximant algorithm for solving nonlinear ordinary differential equation boundary value problems on an unbounded domain. Comp. Phys., 11, 299–303.
  * Boyd, (1999) Boyd, J. P. (1999) The Blasius Function in the Complex Plane. Exper. Math., 8(4), 381–394.
  * Cebeci & Keller, (1971) Cebeci, T. & Keller, H. B. (1971) Shooting and parallel shooting methods for solving the Falkner-Skan boundary-layer equation. J. Comp. Phys., 7, 289–300.
  * Churchill, (1948) Churchill, R. V. (1948) Complex Variables, chapter VI: Power series. McGraw-Hill.
  * Cortell, (2010) Cortell, R. (2010) Numerical Comparisons of Blasius and Sakiadis Flows. MATEMATIKA, 26(2), 187–196.
  * Domb & Sykes, (1957) Domb, C. & Sykes, M. F. (1957) On the Susceptibility of a Ferromagnetic above the Curie Point. Proc. Roy. Soc. London A, 240(1221), 214–228.
  * Eftekhari & Jafari, (2013) Eftekhari, S. A. & Jafari, A. A. (2013) Numerical solution of general boundary layer problems by the method of differential quadrature. Scientia Iranica, 20(4), 1278–1301.
  * Fazio, (1992) Fazio, R. (1992) The Blasius problem formulated as a free boundary problem. Acta Mech., 895, 1–7.
  * Fazio, (2015) Fazio, R. (2015) The iterative transformation method for the Sakiadis problem. Comp. Fluids, 106, 196–200.
  * Flajolet & Sedgewick, (2009) Flajolet, P. & Sedgewick, R. (2009) Analytic Combinatorics, chapter IV.3: Singularities and exponential growth of coefficients. Cambridge University Press.
  * Henrici, (1956) Henrici, P. (1956) Automatic Computations with Power Series. JACM, 3, 10–15.
  * Leal, (1992) Leal, L. G. (1992) Laminar flow and convective transport processes. Scaling Principles and asymptotic analysis, chapter 6.B: Asymptotic Expansions - General Considerations. Butterworth-Heinemann.
  * Probstein, (2003) Probstein, R. F. (2003) Physicochemical Hydrodynamics, chapter 10.3: Coating Flows. Wiley, 2 edition.
  * Ramé et al., (2021) Ramé, E., Weinstein, S. J. & Barlow, N. S. (2021) Free surface shapes in rigid body rotation: Exact solutions, asymptotics and approximants. IMA J. Appl. Math., 86, 1266–1286.
  * Reinberger et al., (2022) Reinberger, W. C., Holland, M. S., Barlow, N. S. & Weinstein, S. J. (2022) On the power series solution to the nonlinear pendulum. Q. J. Mech. Appl. Math, 75(4), 347–369.
  * Sakiadis, (1961) Sakiadis, B. C. (1961) Boundary-layer behavior on continuous solid surfaces: II The boundary layer on a continuous flat surface. AlChE J., 7, 221–225.
  * Turner, (1966) Turner, L. R. (1966) Inverse of the Vandermonde matrix with applications. Technical Note D-3547, NASA.
  * (29) Van Dyke, M. (1975a) Perturbation Methods in Fluid Mechanics, chapter 3.2: Gauge Functions and Order Symbols. Parabolic.
  * (30) Van Dyke, M. (1975b) Perturbation Methods in Fluid Mechanics, chapter 10.8: The Euler Transformation. Parabolic.
  * Watson, (1914) Watson, G. N. (1914) Complex Integration and Cauchy’s Theorem, chapter VII: Expansions in series. Cambridge University Press.
  * Weinstein et al., (2020) Weinstein, S. J., Holland, M. S., Rogers, K. E. & Barlow, N. S. (2020) Analytic solution of the SEIR epidemic model via asymptotic approximant. Physica D, 411, 132633:1–6.
  * Weinstein & Ruschak, (2004) Weinstein, S. J. & Ruschak, K. J. (2004) Coating Flows. Ann. Rev. Fluid Mech., 36, 29–53.


# Transversally exponentially stable Euclidean space extension technique for
discrete time systems

Soham Shanbhag<EMAIL_ADDRESS>School of Electrical Engineering, Korea
Advanced Institute of Science and Technology, Daejeon, Republic of Korea  Dong
Eui Chang Corresponding author<EMAIL_ADDRESS>School of Electrical
Engineering, Korea Advanced Institute of Science and Technology, Daejeon,
Republic of Korea

###### Abstract

We propose a modification technique for discrete time systems for
exponentially fast convergence to compact sets. The extension technique allows
us to use tools defined on Euclidean spaces to systems evolving on manifolds
by modifying the dynamics of the system such that the manifold is an attractor
set. We show the stability properties of this technique using the simulation
of the rigid body rotation system on the unit sphere $S^{3}$. We also show the
improvement afforded due to this technique on a Luenberger like observer
designed for the rigid body rotation system on $S^{3}$.

Keywords: Discrete time systems, transversal extension, embedding, manifolds

##  1 Introduction

With the proliferation of miniaturized microcontrollers which are discrete in
nature, the study of discrete time systems in control theory is gaining
importance. An example of such a system is the rate-integrating gyroscope
measuring a change in rotation $\theta_{k}$ between two instants of time
$t=t_{k}$ to $t=t_{k+1}$, leading to rigid body kinematics of the form

$\displaystyle q(k+1)=q(k)\exp(\theta_{k}/2).$

Moreover, almost all the theory of system identification is in discrete time
due to availability of inputs and outputs at every time step while information
about the derivatives is rarely available. While most of the tools for
studying control systems are developed for continuous time systems, some tools
for linear discrete time systems have been developed, such as the Luenberger
observer[13] and the Kalman filter[9]. Frequency domain analysis of discrete
time systems is also a highly studied topic for linear discrete time systems.

However, such techniques do not extend to nonlinear discrete time systems.
Studies of nonlinear Fourier transforms are sparse, and tools like the Kalman
filter are applied to nonlinear systems[10] by linearizing the nonlinear
system, which leads to divergence of the estimate from the state in case of
highly nonlinear systems. Other filters like the unscented Kalman filter[8]
and the particle filter are computationally expensive. However, some results
on feedback linearization of discrete time systems have been addressed[12, 7]
and an observer has been designed for nonlinear discrete time systems[14].

Moreover, many such systems evolve on manifolds. The study of such systems
evolving on manifolds requires study of properties specific to the manifold,
like a notion of distance on the manifold. An alternative method for such
systems is considering the manifolds as embedded in a higher dimensional
Euclidean space, and considering the systems as evolving on those Euclidean
spaces. Such techniques have been used for various continuous time
applications, such as designing observers for the rigid body rotation and
translation system[16]. These require techniques designed specifically for
such nonlinear systems with the model taken into consideration. Moreover,
using techniques designed for linear discrete time systems for nonlinear
discrete time systems evolving on manifolds embedded in Euclidean spaces may
lead to the state of the system deviating from the manifold.

One of the methods for tackling this problem for continuous time systems is
the was proposed in Feedback Integrator[3], where the author designs an
integration scheme which allows discretization methods to appropriately
estimate the state of the system by adding a transversally stable term to the
system equation. This technique has been further used for various applications
in continuous time systems, like control system design and Kalman filters [1,
15, 6, 2, 11]. The extension technique has also been extended to handle non-
holonomic systems [4, 5].

In this paper, we extend the above technique to discrete time systems. This
enables use of linear discrete time techniques for nonlinear systems while
respecting the underlying manifold. This problem is different from that used
in applications of [3] where discretization has been performed, since here we
study a discrete time system as compared to the discretization of a continuous
time system.

##  2 Main Results

Denote by $\mathbb{R}$ the set of real numbers. The set $f^{-1}(a)$ is defined
as $\\{x\in dom(f)\mid f(x)=a\\}$, where $dom(f)$ is the domain of the map
$f$. The set $f^{-1}(U)$ is then defined as $\cup_{x\in U}f^{-1}(x)$. The norm
of $A\in\mathbb{R}^{n\times m}$ is denoted by $\|A\|=\sqrt{A^{T}A}$. Denote
the ball of radius $r$ around $x\in\mathbb{R}^{n}$ as
$B_{r}(x)=\\{y\in\mathbb{R}^{n}\mid\|y-x\|<r\\}$. Let
$\textrm{dist}\left(A,B\right)$ denote the distance between the sets $A$ and
$B$ defined as $\textrm{dist}\left(A,B\right)=\inf\\{\|x-y\|\mid x\in A,y\in
B\\}$. Denote by $\lambda_{max}(A)$ the largest eigenvalue of $A$.

We now present a modification scheme for discrete time systems which converts
compact sets to local attractors. The following theorem provides the main
result of this paper.

######  Theorem 1.

Consider a discrete time dynamical system on the subset
$U\subset\mathbb{R}^{n}$

$\displaystyle x_{k+1}=f(x_{k},w_{k}),\quad x_{k}\in U,\;w_{k}\in W$

where $f$ is a $C^{1}$ map from $U\times W$ to $U$, where
$W\subset\mathbb{R}^{m}$. We make the following assumptions:

  1. A1

There is a $C^{2}$ function $V:U\to\mathbb{R}_{\geq 0}$ such that $V^{-1}(0)$
is compact and nonempty, and $V\circ f(x,w)=V(x)$ for all $x\in U$ and $w\in
W$.

  2. A2

There exist $b,\epsilon>0$ such that $V^{-1}([0,\epsilon])$ is compact and
$\left\|\nabla V(x)\right\|^{2}\geq bV(x)$ for all $x\in
V^{-1}([0,\epsilon])$.

  3. A3

There is a $\delta>0$ such that the set
$S_{\delta}=\\{x\in\mathbb{R}^{n}\mid\textrm{dist}\left(\\{x\\},V^{-1}([0,\epsilon])\right)<\delta\\}\subset
U$.

Then, the modified dynamical system on $V^{-1}([0,\epsilon])$,

$\displaystyle\tilde{x}_{k+1}=f(\tilde{x}_{k},w_{k})-\alpha\nabla
V(f(\tilde{x}_{k},w_{k})),$ (1)

converges to the set $V^{-1}(0)$ exponentially fast if $\alpha$ is chosen such
that $0<\alpha<\min\left(\frac{2}{d},\frac{\delta}{L}\right)$, where $L$ is
the maximum value of $\|\nabla V(x)\|$ for all $x\in V^{-1}([0,\epsilon])$,
and $d>0$ is chosen such that $d>\lambda_{max}(D^{2}V(x))$ for all $x\in
S_{\delta}$.

###### Proof.

We first show that for all $\tilde{x}_{k}\in V^{-1}([0,\epsilon])$, the state
$\tilde{x}_{k+1}\in S_{\delta}$. Since $V^{-1}([0,\epsilon])$ is compact and
$V$ is a $C^{2}$ function, $\|\nabla V(x)\|<L$ for all $x\in
V^{-1}([0,\epsilon])$. Note that for all $\tilde{x}_{k}\in
V^{-1}([0,\epsilon])$ and $w_{k}\in W$, we have from Assumption A1 that
$f(\tilde{x}_{k},w_{k})\in V^{-1}([0,\epsilon])$. Hence, $\tilde{x}_{k+1}\in
B_{\alpha L}(f(\tilde{x}_{k},w_{k}))\subset S_{\delta}$ since $\alpha
L<\delta$.

From equation (1), we have that

$\displaystyle V(\tilde{x}_{k+1})=V(f(\tilde{x}_{k},w_{k})-\alpha\nabla
V(f(\tilde{x}_{k},w_{k}))).$

For $\tilde{x}_{k}\in V^{-1}([0,\epsilon])$, using the Lagrange form of the
remainder of Taylor series expansion, and noting that due to continuity of
$V$, there exists a $c_{k}$ between the line joining $f(\tilde{x}_{k},w_{k})$
and $\tilde{x}_{k+1}$ due to Assumption A3 such that

$\displaystyle V(\tilde{x}_{k+1})$
$\displaystyle=V(f(\tilde{x}_{k},w_{k}))-\alpha\big{\|}\nabla
V(f(\tilde{x}_{k},w_{k}))\big{\|}^{2}+\frac{\alpha^{2}}{2}\nabla
V^{T}(f(\tilde{x}_{k},w_{k}))D^{2}V(c_{k})\nabla V(f(\tilde{x}_{k},w_{k}))$
$\displaystyle=V(\tilde{x}_{k})-\alpha\big{\|}\nabla
V(f(\tilde{x}_{k},w_{k}))\big{\|}^{2}+\frac{\alpha^{2}}{2}\nabla
V^{T}(f(\tilde{x}_{k},w_{k}))D^{2}V(c_{k})\nabla V(f(\tilde{x}_{k},w_{k}))$
$\displaystyle\leq V(\tilde{x}_{k})-\alpha\big{\|}\nabla
V(f(\tilde{x}_{k},w_{k}))\big{\|}^{2}+\frac{\alpha^{2}}{2}\lambda_{max}(D^{2}V(c_{k}))\big{\|}\nabla
V(f(\tilde{x}_{k},w_{k}))\big{\|}^{2}$

where we have used the property that for any $A\in\mathbb{R}^{n\times n}$,
$x^{T}Ax\leq\lambda_{max}(A)\|x\|^{2}$. Since $\tilde{x}_{k+1}\in B_{\alpha
L}(f(\tilde{x}_{k},w_{k}))$ implies that $c_{k}\in B_{\alpha
L}(f(\tilde{x}_{k},w_{k}))\subset S_{\delta}$, and from the choice of $d$ we
have

$\displaystyle V(\tilde{x}_{k+1})$ $\displaystyle\leq
V(\tilde{x}_{k})-\left(\alpha-\frac{\alpha^{2}}{2}d\right)\big{\|}\nabla
V(f(\tilde{x}_{k},w_{k}))\big{\|}^{2}.$

Since $0<\alpha<2/d$, the value of $V(x_{k})$ is non-increasing in $k$.
Substituting $c:=\alpha-\alpha^{2}d/2>0$, we have

$\displaystyle V(\tilde{x}_{k+1})$ $\displaystyle\leq
V(\tilde{x}_{k})-c\big{\|}\nabla V(f(\tilde{x}_{k},w_{k}))\big{\|}^{2}.$

Since $\tilde{x}_{k}\in V^{-1}([0,\epsilon])$ implies $\tilde{x}_{k+1}\in
S_{\delta}$, we have $V(\tilde{x}_{k+1})\geq 0$. Hence,
$\|V(f(\tilde{x}_{k},w_{k}))\|^{2}\leq\frac{1}{c}V(\tilde{x}_{k})=\frac{1}{c}V(f(\tilde{x}_{k},w_{k}))$.
From Assumption A2, $\|\nabla V(f(\tilde{x}_{k},w_{k}))\|^{2}\geq
bV(f(\tilde{x}_{k},w_{k}))$, hence $b\leq\frac{1}{c}$, or $1-bc\geq 0$. Then

$\displaystyle V(\tilde{x}_{k+1})\leq V(\tilde{x}_{k})-c\|\nabla
V(f(\tilde{x}_{k},w_{k}))\|^{2}\leq
V(\tilde{x}_{k})-bcV(f(\tilde{x}_{k},w_{k}))=V(\tilde{x}_{k})-bcV(\tilde{x}_{k})=(1-bc)V(\tilde{x}_{k}).$

Since $|1-bc|<1$, $V(\tilde{x}_{k})\to 0$ as $t\to\infty$ exponentially fast.
Hence, $V^{-1}(0)$ is exponentially stable with respect to the dynamics in
system (1). ∎

This theorem shows that by a suitable choice of $\alpha$, a discrete time
system can be modified such that a compact set defined as a zero set of a
$C^{2}$ function is a local attractor of the system. Moreover, the convergence
rate is exponential. The above theorem is also applicable when the system is
autonomous, with a similar proof. The following corollary presents the
conditions for the same.

######  Corollary 1.

Consider a discrete time dynamical system on the subset
$U\subset\mathbb{R}^{n}$

$\displaystyle x_{k+1}=f(x_{k}),\quad x_{k}\in U,$

where $f$ is a $C^{1}$ map from $U$ to itself. We make the following
assumptions:

  1. CA1

There is a $C^{2}$ function $V:U\to\mathbb{R}_{\geq 0}$ such that $V^{-1}(0)$
is compact and nonempty, and $V\circ f=V$ on $U$.

  2. CA2

There exist $b,\epsilon>0$ such that $V^{-1}([0,\epsilon])$ is compact and
$\left\|\nabla V(x)\right\|^{2}\geq bV(x)$ for all $x\in
V^{-1}([0,\epsilon])$.

  3. CA3

There is a $\delta>0$ such that the set
$S_{\delta}=\\{x\in\mathbb{R}^{n}\mid\textrm{dist}\left(\\{x\\},V^{-1}([0,\epsilon])\right)<\delta\\}\subset
U$.

Then, the modified dynamical system on $V^{-1}([0,\epsilon])$,

$\displaystyle\tilde{x}_{k+1}=f(\tilde{x}_{k})-\alpha\nabla
V(f(\tilde{x}_{k})),$ (2)

converges to the set $V^{-1}(0)$ exponentially fast for $\alpha$ chosen as
specified in Theorem 1.

###### Proof.

The proof follows from that of Theorem 1 and is omitted here. ∎

The above theorem can be used to extend systems evolving on compact manifolds
to ambient Euclidean spaces using manifold embedding theorems. This also shows
the permissiveness of the assumptions on $V$. Since many systems, like the
rigid body system used in quadcopter and satellite modelling, evolve on the
special orthogonal group, the proposed theorem can be used for a huge class of
problems. Moreover, we here show a simple way to derive the function $V$ from
the manifold. The specified choice of $V$ and $U$ given in the proof of the
theorem enable us to choose $\epsilon$ and $\alpha$ such that the assumptions
of Theorem 1 are satisfied.

######  Theorem 2.

Consider a system evolving on the compact manifold $M:=g^{-1}(0)$ given by the
difference equation

$\displaystyle x_{k+1}=f(x_{k},w_{k}),\quad x_{k}\in M$

where $f$ is a $C^{1}$ map from $M$ to itself, and
$g:\mathbb{R}^{n}\to\mathbb{R}^{k}$ is a $C^{2}$ function. Consider the
modified dynamical system

$\displaystyle\tilde{x}_{k+1}=f(\tilde{x}_{k})-\alpha
Dg(f(\tilde{x}_{k}))^{T}g(\tilde{x}_{k}),\quad\tilde{x}_{k}\in
g^{-1}(B_{\epsilon}(0)),$

where $\epsilon$ is such that $g^{-1}(B_{\epsilon}(0))$ is compact. The
modified dynamical system converges to the manifold $M$ if:

  1. A1’

For some $\varepsilon>\epsilon$ such that $g^{-1}(B_{\varepsilon}(0))$ is
connected, $g(f(x))=g(x)$ for all $x\in g^{-1}(B_{\varepsilon}(0))$.

  2. A2’

For all $x\in g^{-1}(B_{\epsilon}(0))$, there exists a $b>0$ such that
$\sigma_{min}(Dg(x))\geq b$.

  3. A3’

For all $x\in g^{-1}(B_{\varepsilon}(0))$, $\sigma_{max}(Dg(x))$ and
${\lambda_{max}(D^{2}g_{i}(x)),\;i=1,\ldots,k}$ are bounded, where $g_{i}(x)$
is the $i^{th}$ component of $g(x)$.

Here, $\alpha$ is chosen such that

$\displaystyle
0<\alpha<\min\left(\frac{2}{d},\frac{2\delta}{L\epsilon^{2}}\right),$

where $L=\sigma_{max}(Dg(x))$ for all $x\in g^{-1}(B_{\epsilon}(0))$,
$\delta>0$ such that
$\delta<\textrm{dist}\left(g^{-1}(bd(B_{\varepsilon}(0))),g^{-1}(B_{\epsilon}(0))\right)$
and $d>0$ is such that
$d>\lambda_{max}(Dg(x)^{T}Dg(x)+\sum_{i}g_{i}(x)D^{2}g_{i}(x))$ for all $x\in
g^{-1}(B_{\varepsilon}(0))$.

###### Proof.

Choose

$\displaystyle V(x)=\frac{1}{2}\|g(x)\|^{2}\textrm{ and
}U=g^{-1}(B_{\varepsilon}(0))$

such that

$\displaystyle\nabla V(x)=Dg(x)^{T}g(x),\textrm{ and
}D^{2}V(x)=Dg(x)^{T}Dg(x)+\sum_{i}g_{i}(x)D^{2}g_{i}(x).$

Then, our choice of $V$ along with Assumption A1’ satisfies Assumption A1. To
see that Assumption A2 is satisfied, note that

$\displaystyle\|Dg(x)^{T}g(x)\|^{2}\geq\sigma_{min}(Dg(x))^{2}g(x)^{T}g(x)\geq
b^{2}g^{T}(x)g(x)$

from Assumption A2’.

We now show that our system satisfies Assumption A3 of Theorem 1, i.e.
$S_{\delta}\subset U$. For this system, the set $S_{\delta}$ is defined as
$\\{x\in\mathbb{R}^{n}\mid\textrm{dist}\left(\\{x\\},g^{-1}(B_{\epsilon}(0))\right)<\delta\\}$.
From the definition of $S_{\delta}$, if $x\in S_{\delta}$, then $x\in
B_{\delta}(y)$ for some $y\in g^{-1}(B_{\epsilon}(0))$. Note that
$V(y)\leq\frac{1}{2}\epsilon^{2}$. Let us assume that $x\notin U$. Then,
$V(x)>\frac{1}{2}\varepsilon^{2}>\frac{1}{2}\epsilon^{2}$. By the intermediate
value theorem, from the continuity of $g$ and $V$, there exists a $z$ on the
line joining $x$ and $y$ such that $V(z)=\frac{1}{2}\varepsilon^{2}$. This
implies $z\in B_{\delta}(y)\subset S_{\delta}$. However,
$\delta<\textrm{dist}\left(g^{-1}(bd(B_{\varepsilon}(0))),g^{-1}(B_{\epsilon}(0))\right)$.
This is a contradiction. Hence, $x\in U$, i.e. for any $x\in S_{\delta}$,
$x\in U$. This shows that our system satisfies Assumption A3.

In Theorem 1, the closure of the set $S_{\delta}$ is compact, hence the second
derivative of $V$ is bounded. However, since here we do not assume the
boundedness of the set $g^{-1}(B_{\varepsilon}(0))$, we must assume Assumption
A3’ for existence of $d$ since

$\displaystyle\lambda_{max}(Dg(x)^{T}Dg(x)+\sum_{i}g_{i}(x)D^{2}g_{i}(x))$
$\displaystyle\leq\lambda_{max}(Dg(x)^{T}Dg(x))+\sum_{i}g_{i}(x)\lambda_{max}(D^{2}g_{i}(x))$
$\displaystyle=\sigma_{max}(Dg(x))+\sum_{i}g_{i}(x)\lambda_{max}(D^{2}g_{i}(x)),$

which is bounded due to Assumption A3’ and boundedness of $g(x)$. To calculate
the bounds of $\alpha$, note that

$\displaystyle\frac{1}{2}L\epsilon^{2}$
$\displaystyle\geq\sigma_{max}(Dg(x))\|g(x)\|^{2}\geq\|Dg(x)^{T}g(x)\|$

for all $x\in g^{-1}(B_{\epsilon}(0))$. ∎

The manifold extension formulation is beneficial when tools developed for
Euclidean spaces are to be used for systems on manifolds. Usually, study of
systems evolving on manifolds require the notion of distance defined on the
manifold, and an “error term” to be suitably defined using this distance. An
alternative method would be embedding the manifold into ambient Euclidean
spaces, and using the subtraction operator defined on the Euclidean space.
However, when using this method, the manifold may not be forward invariant,
i.e. the state of the system may diverge from the manifold. A transversally
exponentially stable term will lead to the system state to converge to the
manifold, and thus lead to improved results. An example of this is given in
subsection 3.2, where we design a Luenberger observer for a system evolving on
$S^{3}$ and compare two observers, one with the extension term and one without
the extension term. We consider randomly chosen initial conditions and show
that in all cases, the observer with the transversally stable extension term
estimates the state of the system much better as compared to the one without
the extension term.

The manifold formulation also shows that although the restrictions on the
choice of $V$ exist, they are permissive enough that for many systems, a
suitable $V$ can be chosen. Moreover, the choice of $V$ may be straightforward
in a lot of cases, as shown by the choice in Theorem 2.

##  3 Simulation

We now present simulations for showing the effectiveness of this technique. To
show convergence properties, we perform Monte Carlo simulations. The
simulation is that of a rigid body rotation system where the state is defined
using quaternions. This system evolves on $S^{3}$, the unit sphere in the
quaternion space. We show the convergence of the system state of the
quaternion to $S^{3}$. Lastly, to show implementation in real world systems,
we show the effect the addition of the extension term leads to in the case of
a Luenberger like observer.

###  3.1 Simulation on the unit sphere in quaternion space

The algebra of quaternions is denoted by $\mathbb{H}$, and let $S^{3}$ denotes
the set of unit quaternions in $\mathbb{H}$. A vector quaternion is a
quaternion whose real part is $0$. The quaternion $[1,0,0,0]$ is denoted by
$1$. The conjugate quaternion of $q$ is denoted by $\bar{q}$. The norm of the
quaternion $q$ is denoted by $\|q\|=\sqrt{\bar{q}q}$.

We consider the system equations

$\displaystyle q_{k+1}=q_{k}\exp(\Omega_{k}/2),$ (3)

where $q_{k}\in S^{3}$ denotes the rotational state of the body, and
$\Omega_{k}\in\mathbb{H}$ is a vector quaternion denoting the angular rotation
of the body.

Define

$\displaystyle V(x)=\|\bar{x}x-1\|^{2},\quad x\in\mathbb{H},\textrm{ and
}U=\mathbb{H}.$

Using the above choice of $V$, we get that

$\displaystyle\nabla V(x)=4x(\bar{x}x-1)\textrm{ and
}\lambda_{max}(D^{2}V(x))=12x^{T}x.$

Note that the assumptions of Theorem 1 hold here, since:

  * •

$V(f(x_{k},\Omega_{k}))=\|\exp(-\Omega_{k}/2)\bar{x}x\exp(\Omega_{k}/2)-1\|^{2}=\|\bar{x}x-1\|^{2}=V(x)$,
and $V^{-1}(0)=\\{x\mid\bar{x}x=1\\}$ is compact and nonempty.

  * •

For any $\epsilon<1$, $V^{-1}([0,\epsilon])=\\{x\mid
1-\sqrt{\epsilon}\leq\bar{x}x\leq 1+\sqrt{\epsilon}\\}$ is compact, and
choosing $b=16(1-\sqrt{\epsilon})$,

$\displaystyle\|\nabla V(x)\|^{2}=16\|x\|^{2}V(x)\geq bV(x).$

Then, the system equations are modified using the proposed extension technique
as

$\displaystyle\tilde{q}_{k+1}=\tilde{q}_{k}\exp(\Omega_{k}/2)-4\alpha\tilde{q}_{k}(\bar{\tilde{q}}_{k}\tilde{q}_{k}-1)\exp(\Omega_{k}/2)$
(4)

where $\tilde{q}_{k}\in\mathbb{H}$ is the state of the quaternion.

For the simulation, consider $\epsilon=0.5$, $\delta=0.059$, such that
$L=3.69$ and $d=15.87$. Then $\alpha<0.016$. Specifically, we choose
$\alpha=0.01$. The angular rotation of the system is chosen uniformly between
$(0,0,0,0)$ and $(0,0.5,0.5,0.5)$. For showing the convergence of the state of
the system to the manifold $S^{3}$, we perform simulation with the initial
condition given by $[w,0,0,0]$, where $w$ is a random number chosen uniformly
between $\sqrt{1-\sqrt{\epsilon}}$ and $\sqrt{1+\sqrt{\epsilon}}$. We choose
1000 such initial conditions and run each simulation for 100 epochs.

$0$$0.1$$0.2$$0.3$$0.4$$0.5$$0.6$$0.7$$0.8$$0$$10$$20$$30$$40$$50$$60$$70$$80$$90$$100$

$\textrm{dist}\left(q,S^{3}\right)$

EpochsStatesMean State

Figure 1: Convergence of the state of the quaternion system to the manifold
$S^{3}$

$0.1$$0.2$$0.3$$0.4$$0.5$$0.6$$0.7$$0.8$$0.9$$1$$1.1$$0$$10$$20$$30$$40$$50$$60$$70$$80$$90$$100$

$\|q-\hat{q}\|$

Epochs$w$-observer$wo$-observer

(a)

$0$$0.5$$1$$1.5$$2$$2.5$$3$$3.5$$4$$4.5$$5$$0$$10$$20$$30$$40$$50$$60$$70$$80$$90$$100$

$\|q-\hat{q}\|$

Epochs$w$-observer$wo$-observer

(b)

$0$$0.05$$0.1$$0.15$$0.2$$0.25$$0.3$$0.35$$0.4$$0.45$$0$$10$$20$$30$$40$$50$$60$$70$$80$$90$$100$

$\textrm{dist}\left(\hat{q},S^{3}\right)$

Epochs$w$-observer$wo$-observer

(c)

$0$$5$$10$$15$$20$$25$$30$$0$$10$$20$$30$$40$$50$$60$$70$$80$$90$$100$

$\textrm{dist}\left(\hat{q},S^{3}\right)$

Epochs$w$-observer$wo$-observer

(d)

Figure 2: Convergence of the state of the quaternion system to the manifold
$S^{3}$

The distance of the state from the manifold $S^{3}$ can be seen in Figure 1.
The black line shows the mean of all the states of the system. It can be seen
from the simulation that the mean of the states of the system converges to the
underlying manifold exponentially fast. This can be seen for all the states
starting in $V^{-1}([0,\epsilon])$. The gray shaded region shows the spread of
the states, which is highest at the initial epoch due to the randomly chosen
initial states. The spread also exponentially descreases, implying that all
the systems eventually converge to the manifold.

###  3.2 Simulation of observer on $S^{3}$

We show the effect of using such an additive term in designing observers for
systems. A formal proof and the study of the convergence properties along with
comparisons would be a complete work on its own, hence we only show Monte
Carlo simulations here. We consider the system defined in subsection 3.1.

We first note that to design an observer for this system on the manifold, a
notion of distance on the manifold between two elements of the manifold is to
be defined. An observer then needs to be designed such that this distance
tends to $0$ as $t\to\infty$. We sidestep this requirement by considering the
manifold to be embedded in the ambient Euclidean space, and using the
Euclidean distance as the distance between two point on the manifold. However,
this formulation will push the state of the system outside the manifold,
whereas a desirable property is to restrict the state to the manifold. To
remedy this, the system (3) is extended to the ambient Euclidean space as
shown in equation (4). To show the improvement due to the extension term, the
performance of an observer designed using the modified system (4) is compared
with an observer designed using the system (3) embedded into $\mathbb{H}$
without the extension term.

Motivated by a Luenberger like observer system, the observer equations are
considered as

$\displaystyle\hat{q}_{w,k+1}$
$\displaystyle=\hat{q}_{w,k}\exp(\Omega_{k}/2)+(q_{k}-\hat{q}_{w,k})L-4\alpha\hat{q}_{w,k}(\bar{\hat{q}}_{w,k}\hat{q}_{w,k}-1)\exp(\Omega_{k}/2),$
$\displaystyle\hat{q}_{wo,k+1}$
$\displaystyle=\hat{q}_{wo,k}\exp(\Omega_{k}/2)+(q_{k}-\hat{q}_{wo,k})L,$

where the subscript $w$ denotes the observer with the transversally stable
term and the subscript $wo$ denotes the observer without this term. Here
$q_{k}$ denotes the measured state of the system evolving using the difference
equation (3). For the sake of brevity, we use the term $w$-observer and
$wo$-observer for systems containing the transversally stable term and for
those without the transversally stable term.

The initial condition for the system state is chosen as a random state in
$S^{3}$. The initial condition for the observer states are chosen as $1$. We
run 1000 such simulations with each simulation run for 100 epochs. The state
of the system at epoch $k$ is measured as the true state of the system
multiplied by $\exp(\nu_{k})$ where $\nu_{k}$ is Gaussian noise with mean $0$
and standard deviation $0.1$. The angular velocity of the system is chosen as
$(0,\omega_{x},\omega_{y},\omega_{z})$ such that
$0<\omega_{x},\omega_{y},\omega_{z}<1$ with the time step between two epochs
as $10$, such that $\Omega_{k}$ is of the form
$(0,\theta_{x},\theta_{y},\theta_{z})$ such that
$0<\theta_{x},\theta_{y},\theta_{z}<10$. We choose
$L=(0.11,-0.18,-0.18,-0.18)$ as the mean of $\exp(\Omega_{k}/2)$ since there
exists no choice of $L$ such that $\|\exp(\Omega_{k}/2)-L\|<1$ for all
$\Omega_{k}$. The constants $\epsilon,\delta$ and $\alpha$ are chosen as in
subsection 3.1.

The convergence of the observers to the system state is shown in Figure 2.
Here, mean state denotes the mean of all the states of the system. The maximum
error in the observer states, minimum error in the observer states (both
calculated for all $k>25$ and among all simulations) and the final error mean
and standard deviation are shown in Table 1. We note that while both the
observers converge to the system state, the mean of the error of the
$w$-observer is less than that of the $wo$-observer. This difference in mean
is significant in case of the distance from the $S^{3}$ manifold. If the
observer state is intended to be restricted to $S^{3}$, by projecting the
unrestricted observer state to $S^{3}$, the error in the state of the observer
from the manifold also contributes to the error in the observer. Moreover, the
maximum error in case of the $w$-observer is also considerably less than that
of the $wo$-observer, implying that the variation of the error in case of the
$w$-observer is also lower than that of the $wo$-observer. This difference in
variation can be seen in the standard deviation of the observers, where the
standard deviation of the distance of the $wo$-observer state from the
manifold is 1.5 times that of the $w$-observer. This simulation shows the
benefits of extending the system to the ambient space using a transversally
stable term before designing the observer system.

Table 1:  Error in the observer states for Luenberger observer system. Observer | $\|q-\hat{q}\|$ | $\textrm{dist}\left(q,S^{3}\right)$
---|---|---
Mean | Std | Max | Min | Mean | Std | Max | Min
$w$-observer | 0.138 | 0.07 | 1.816 | 0.007 | 0.084 | 0.072 | 2.182 | 0.
$wo$-observer | 0.148 | 0.076 | 4.747 | 0.007 | 0.112 | 0.101 | 25.121 | 0.

##  4 Conclusion

In this paper, we extend the ambient space transversally stable extension
technique[3] to discrete time systems. The use of transversally stable
extension has been used to improve existing results in control theory in the
case of continuous time systems. We expect that the provided formulation in
Theorem 1 enables such improvement in the case of discrete time systems. An
example of improvements afforded by such a formulation has been shown in
subsection 3.2. A quantitative analysis of such improvement for different
problems in control theory and their comparison with corresponding manifold
respecting formulations would be a future area of research.

## References

  * [1] Fanchen Bu and Dong Eui Chang.  Feedback gradient descent: Efficient and stable optimization with orthogonality for DNNs.  Proceedings of the AAAI Conference on Artificial Intelligence, 36(6):6106–6114, June 2022.
  * [2] Dong Eui Chang.  On controller design for systems on manifolds in Euclidean space.  International Journal of Robust and Nonlinear Control, 28(16):4981–4998, 2018.
  * [3] Dong Eui Chang, Fernando Jiménez, and Matthew Perlmutter.  Feedback integrators.  Journal of Nonlinear Science, 26(6):1693–1721, December 2016.
  * [4] Dong Eui Chang and Matthew Perlmutter.  Feedback integrators for nonholonomic mechanical systems.  Journal of Nonlinear Science, 29(3):1165–1204, June 2019.
  * [5] Dong Eui Chang, Matthew Perlmutter, and Joris Vankerschaver.  Feedback integrators for mechanical systems with holonomic constraints.  Sensors, 22(17):6487, January 2022.
  * [6] Dong Eui Chang, Karmvir Singh Phogat, and Jongeun Choi.  Model predictive tracking control for invariant systems on matrix Lie groups via stable embedding into Euclidean spaces.  IEEE Transactions on Automatic Control, 65(7):3191–3198, July 2020\.
  * [7] T. Fliegner, Ü. Kotta, and H. Nijmeijer.  Solvability and Right-Inversion of Implicit Nonlinear Discrete-Time Systems.  SIAM Journal on Control and Optimization, 34(6):2092–2115, November 1996.
  * [8] S.J. Julier and J.K. Uhlmann.  Unscented filtering and nonlinear estimation.  Proceedings of the IEEE, 92(3):401–422, March 2004.
  * [9] R. E. Kalman.  A New Approach to Linear Filtering and Prediction Problems.  Journal of Basic Engineering, 82(1):35–45, March 1960.
  * [10] Rudolf Emil Kalman and Richard S. Bucy.  New results in linear filtering and prediction theory.  Journal of Basic Engineering, 83(1):95–108, March 1961.
  * [11] Tejaswi K.C. and Sukumar Srikant.  Attitude control via a feedback integrator based observer.  Automatica, 151:110882, May 2023.
  * [12] Ülle Kotta.  Inversion Method in the Discrete-time Nonlinear Control Systems Synthesis Problems, volume 205 of Lecture Notes in Control and Information Sciences.  Springer, Berlin, Heidelberg, 1995.
  * [13] D. Luenberger.  An introduction to observers.  IEEE Transactions on Automatic Control, 16(6):596–602, December 1971.
  * [14] P.E. Moraal and J.W. Grizzle.  Observer design for nonlinear systems with discrete-time measurements.  IEEE Transactions on Automatic Control, 40(3):395–404, March 1995\.
  * [15] Jae Hyeon Park, Karmvir Singh Phogat, Whimin Kim, and Dong Eui Chang.  Transversely stable extended Kalman filters for systems on manifolds in Euclidean spaces.  Journal of Dynamic Systems, Measurement, and Control, 143(6), February 2021.
  * [16] Soham Shanbhag and Dong Eui Chang.  Globally exponentially convergent observer for the rigid body system on SE(3).  In 2022 IEEE 61st Conference on Decision and Control (CDC), pages 1257–1262. 2022 IEEE 61st Conference on Decision and Control (CDC), December 2022.