# On the summation of divergent, truncated, and underspecified power series
via asymptotic approximants

Nathaniel S. Barlow1, Christopher R. Stanton1, Nicole Hill2,
Steven J. Weinstein2, and Allyssa G. Cio3.
1School of Mathematical Sciences, 2Department of Chemical Engineering,
3Department of Industrial Engineering
Rochester Institute of Technology, Rochester, NY 14623, USA

###### Abstract

A compact and accurate solution method is provided for problems whose infinite
power series solution diverges and/or whose series coefficients are only known
up to a finite order. The method only requires that either the power series
solution or some truncation of the power series solution be available and that
some asymptotic behavior of the solution is known away from the series’
expansion point. Here, we formalize the method of asymptotic approximants that
has found recent success in its application to thermodynamic virial series
where only a few to (at most) a dozen series coefficients are typically known.
We demonstrate how asymptotic approximants may be constructed using simple
recurrence relations, obtained through the use of a few known rules of series
manipulation. The result is an approximant that bridges two asymptotic regions
of the unknown exact solution, while maintaining accuracy in-between. A
general algorithm is provided to construct such approximants. To demonstrate
the versatility of the method, approximants are constructed for three
nonlinear problems relevant to mathematical physics: the Sakiadis boundary
layer, the Blasius boundary layer, and the Flierl-Petviashvili monopole. The
power series solution to each of these problems is underspecified since, in
the absence of numerical simulation, one lower-order coefficient is not known;
consequently, higher-order coefficients that depend recursively on this
coefficient are also unknown. The constructed approximants are capable of
predicting this unknown coefficient as well as other important properties
inherent to each problem. The approximants lead to new benchmark values for
the Sakiadis boundary layer and agree with recent numerical values for
properties of the Blasius boundary layer and Flierl-Petviashvili monopole.

This article has been accepted for publication in the Quarterly Journal of
Mechanics and Applied Mathematics. The final published version can be found at
https://academic.oup.com/qjmam/article/2895119/On-the-Summation-of-Divergent-
Truncated-and

##  1 Introduction

Power series arise in virtually all applications of mathematical physics. The
utility of such series is evident in the construction of approximations (e.g.
finite differences) or as a rigorously determined solution to a problem.
Limitations generally inherent to power series solutions often inhibit their
direct use, and they are more often useful in the implementation of numerical
schemes. For instance, a Taylor series representation of an unknown function
may not converge, as it may have a finite radius of convergence arising from
singularities (often complex) in the function it represents [1]. Even when
singularities are not a concern, higher-order terms of the series may be
exceedingly difficult to compute [2], which is especially problematic if the
series converges slowly.

Several techniques have been put forward to analytically continue and to
accelerate convergence of divergent or slowly converging series; see, for
example [3] (ch. 8), [4] (chs. 14, 19, 20), [5], [6], [7], [8], [9], [10]. One
of the more well-known techniques is the Padé approximant method [11], which
has similarities with the method of asymptotic approximants described herein.
A Padé approximant (referred to as a Padé) is the quotient of two polynomials,
with $N$ total polynomial coefficients (distributed between the numerator and
denominator), chosen such that the Taylor expansion of the Padé reproduces the
power series of interest up to $N$ terms. This involves solving an algebraic
system with $N$ unknowns (the Padé coefficients). After its construction, one
may use the Padé in lieu of the original series. The $N^{\text{th}}$-order
Padé and $N^{\text{th}}$-order truncated original series will follow one
another near the expansion point. Ideally, the Padé will then continue on
beyond the radius of convergence (if one is present) and in general represent
the actual solution better than the truncated series. If chosen judiciously, a
Padé sequence (specified by a fixed difference in order between the numerator
and denominator) may converge rapidly to the correct solution as $N$ is
increased, whereas the original series may converge slowly or diverge. In the
case of a divergent series, the Padé approximates the singularity that is
presumed to be responsible for divergence [12], and in doing so enables an
accurate summation of the truncated series.

One drawback of the Padé method is that it is not always clear beforehand
which sequence is best suited for a given problem [12]. If one has a power
series alone, it is difficult to choose a correct Padé form (see [13, 14, 15]
for example). However, if the power series arises from a physical problem, it
is likely that some additional conditions or context may be gleaned. These
conditions may be conveniently available, or may need to be independently
derived. For instance, if one has a power series representation of $f(x)$ and
it is known that $f$ asymptotically approaches a constant value as
$x\to\infty$, a “symmetric” Padé sequence of equal order in the numerator and
denominator will contain Padés that all preserve this asymptotic condition,
and this sequence will uniformly converge towards the correct solution as $N$
is increased. Baker and Gammel [11] recognized this important result and went
further to state: If a series for $f(x)$ is known and $f\sim Cx^{p}$ as
$x\to\infty$ ($p$ being an integer), the uniformly convergent Padé sequence is
the one with a difference of $p$ between the numerator and denominator order,
hence preserving the $x\to\infty$ behavior as $N$ is increased.

Baker $\&$ Gammel’s statement may be extended to asymptotic behaviors beyond
integer power laws using approximants other than Padés [16, 17, 18]. An
example of such an approximant is found in the review by Frost and Harper
[10], used to find an approximate solution for the drag coefficient on a
sphere in fluid-filled tube; the approximant they used incorporates the
correct non-integer power-law asymptotic behavior. The present work may be
considered an extension of the perspective given in [10]\- namely the
construction of approximants that are asymptotically consistent with known
behavior in the vicinity of the domain boundaries while maintaining accuracy
in-between. In this paper, we aim to formalize this approach by defining
asymptotic approximants in general:

###### Definition 1.1

Given a power series representation of some function $f(x)$:

$f=\sum_{n=0}^{\infty}a_{n}(x-x_{0})^{n},$ (1)

and an asymptotic behavior

$f\sim Cf_{a}(x)~{}\text{as}~{}x\to x_{a},$ (2)

where $C$ is a constant, an asymptotic approximant is any function $f_{A}(x)$
that may be expressed analytically in closed form and that satisfies the
following three properties:

  1. 1.

The $N$-term Taylor expansion of $f_{A}$ about $x_{0}$ is identical to the
$N$-term truncation of (1).

  2. 2.

$\displaystyle{\lim_{x\to x_{a}}(f_{A}/f_{a})=}$ constant for any $N$.

  3. 3.

The sequence of approximants converges for increasing $N$.

Choosing an approximant that satisfies the above definition will lead to a
uniformly convergent sequence as $N$ is increased that preserves the correct
asymptotic behavior. Note that in the above definition, $f_{a}$ need not be
exact; in fact, only the leading order is typically needed to construct an
adequate approximant. We seek approximants whose unknown coefficients can be
generated with ease. We do not wish for the construction of approximants to
undergo more floating-point operations than performed in an accurate numerical
solution. Most of the approximant coefficients calculated here may be obtained
recursively using a few simple series relations listed in Appendix A. This is
in contrast with Padés and their extensions, which often require the inversion
of a linear system [3, 4, 10]. As shall be apparent, a theme of this work is
to find/impose the simplest asymptotic approximant form possible while
maintaining the desired accuracy and precision.

Approximants are closed form functions that may be analytically differentiated
to obtain auxiliary properties without the loss of accuracy that occurs with
discretized solutions. Also, physically relevant properties may be cast as
unknowns within an approximant, and the approximant can be used as a predictor
for such properties [19, 16, 17, 18]. This feature of approximants provides a
significant problem-solving advantage.

The paper is organized as follows. In Section 2, we review recent applications
of asymptotic approximants towards the truncated virial series of
thermodynamics. The analysis is recast here as part of a new unified framework
that only requires the use of a few simple series relations. In Section 3, an
algorithm is provided to construct asymptotic approximants for problems in
general. The method is then applied to three nonlinear boundary value problems
in Section 4. In Section 4.1, the Sakiadis boundary layer problem is solved
using two types of asymptotic approximants - one which is both accurate and
simple in form, and another which matches its higher-order asymptotic behavior
and is capable of predicting the wall-shear coefficient and other properties
to a higher precision than reported before. In Section 4.2, the Blasius
boundary layer problem is solved using a simple asymptotic approximant, which
is also capable of predicting important quantities such as the wall-shear
coefficient. In Section 4.3, an asymptotic approximant is used to solve the
Flierl-Petviashvili monopole problem and predict relevant properties, which
are shown to agree with newly generated numerical results. Key findings are
summarized in Section 5.

##  2 Relevant background and results from virial-based approximants

In this section, we review the recent application of asymptotic approximants
towards problems in thermodynamics. In doing the review, simplifications and
details not previously disclosed are elucidated here and also incorporated
into the algorithm of Section 3.

The virial series is an equation of state formulated as a series expansion
about the ideal gas limit [20] ($\rho\to 0$), and is expressed as

$P=kT\sum_{n=1}^{N}B_{n}(T)\rho^{n},\hskip 7.22743ptB_{1}=1$ (3)

where $P$ is the pressure, $\rho$ is the number density, $k$ is the Boltzmann
constant, and $T$ is the temperature. The virial coefficients $B_{n}$ are
functions of $T$ as well as other physical parameters describing a given
intermolecular potential. Specifically, the $n^{th}$ virial coefficient is
defined as an integral over the position of $n$ molecules [21]. In the absence
of an exact equation of state, the virial series represents fluid properties
at low density more accurately than any other theoretical or computational
method, since it is effectively the Taylor series of the exact solution about
$\rho$=0. A barrier to its usage is that the number of integrals appearing for
each viral coefficient increases rapidly and nonlinearly with the order of the
coefficient. Consequently, significant effort is spent to develop efficient
algorithms to compute higher-order coefficients [22, 23]. As an indication of
the difficulty involved, only twelve virial coefficients are currently known
(to within some precision) for the hard-sphere model fluid [23]. For more
realistic model fluids, significantly fewer coefficients are known.

Even with several coefficients, the virial series often converges only in a
small region near $\rho=0$. Padé approximants have been applied to extend this
region, with a common example being the Carnahan-Starling equation of state
for hard-sphere fluids [24]. With only a finite number of terms available, the
choice of approximant (be it Padé or not), is important as convergence may be
accelerated by choosing an approximant with consistent asymptotic behavior
following Definition 1.1. In the following subsections, we review the recent
application of asymptotic approximants towards divergent and slowly converging
truncated virial series to construct accurate isotherms and predict critical
properties.

###  2.1 Soft-sphere fluids

To our knowledge, the first application of a non-Padé asymptotic approximant
towards virial series is given in [16] and is used to describe isotherms for
soft-sphere fluids. This problem highlights the key steps in formulating an
asymptotic approximant used later in this paper, and we thus review the
technique carefully. Soft-sphere fluids [25] are defined by a molecular
potential that behaves like $r^{-h}$, where $r$ is the distance between
molecules and $h$ is a “hardness” parameter. In the limit as $h\to\infty$, the
soft-sphere molecular potential limits to that of the hard-sphere. For soft
spheres, it is typical to rewrite (3) in terms of reduced virial coefficients
$\bar{B}_{n}$ and reduced density $\tilde{\rho}$ using the scalings provided
in [16]. Equation (3) becomes

$Z\equiv\frac{P}{\rho
kT}=1+\sum_{n=2}^{N}\bar{B}_{n}(h)\tilde{\rho}^{~{}n-1},$ (4)

where $Z$ is the compressibility factor and the coefficients $\bar{B}_{n}$ are
now functions of the hardness $h$. The above series diverges (see Fig. 1 for
$h=4$), with a radius of convergence that decreases with decreasing $h$ [16].
While (4) accurately describes the low-density behavior, an approximant can be
constructed by bridging (4) with the large density limit for soft-spheres:

$Z\sim C\tilde{\rho}^{~{}h/3}\text{ as }\tilde{\rho}\to\infty.$ (5)

Note that if $h/3$ is not an integer, it is impossible to capture the above
behavior with a standard Padé. However, an asymptotic approximant such as

$Z_{A}=\left(1+\sum_{n=1}^{N-1}A_{n}(h)\tilde{\rho}^{~{}n}\right)^{\frac{h/3}{N-1}}$
(6a) limits to (5) for any $h$ and allows for an easy evaluation of its
coefficients, such that the Taylor expansion of (6a) about $\tilde{\rho}=0$ is
exactly (4). An expression for the unknown coefficients $A_{n}$ is found by
equating the original series (4) with the approximant form (6a) and then
solving for the series in $A_{n}$:
$\left(1+\sum_{n=2}^{N}\bar{B}_{n}(h)\tilde{\rho}^{~{}n-1}\right)^{\frac{N-1}{h/3}}=1+\sum_{n=1}^{N-1}A_{n}(h)\tilde{\rho}^{~{}n}.$
It is now clear that the $A_{n}$ coefficients are the coefficients of the
Taylor expansion of the left-hand side of the above expression about
$\tilde{\rho}=0$. Such coefficients are easily obtained using J. C. P.
Miller’s formula [26] for recursively evaluating the expansion of a series
raised to any power; see equation (29) in Appendix A. The recursion for the
coefficients is
$A_{n>0}=\frac{1}{n}\sum_{j=1}^{n}\left[\frac{j(N-1)}{h/3}-n+j\right]\bar{B}_{j+1}A_{n-j},~{}~{}A_{0}=1.$
(6b)

Together, (6a) and (6b) provide an asymptotic approximant which converges for
all $\tilde{\rho}$ in the fluid regime, as shown in [16] and Fig. 1 (for
$h$=4). Note that (6b) reduces to the formulae for $A_{2}$ through $A_{10}$
given in the appendix of [16], and can be used to generate higher-order
approximants as higher-order virial coefficients become available. The same
procedure used above will be applied in Section 4 to obtain asymptotic
approximants for nonlinear boundary value problems of mathematical physics.

Figure 1: Compressibility factor $Z$ versus reduced density $\tilde{\rho}$ for
the $r^{-4}$ soft-sphere fluid. A comparison is shown between the $N$-term
virial series (3) ($\\---$) labeled as S$N$, the corresponding approximant (6)
($\\--$) labeled as A$N$, and data from molecular simulation [27] ($\bullet$).
The plot spans the entire fluid regime, as the last data point is where
freezing occurs. Virial coefficients used to generate the curves are taken
from [16].

###  2.2 Prediction of Critical Properties

An important property of asymptotic approximants is their ability to predict
unknown quantities. Padé and other approximants have long been used to predict
critical properties of lattice models [28, 29, 6, 7]. With the advent of
asymptotic approximants, this approach has recently yielded accurate
predictions of critical properties of model fluids [17, 18]. The prediction
begins by constructing the approximant in the same manner as the previous
section.

The virial series describing $P(\rho)$ along the critical isotherm (i.e. using
coefficients evaluated at the critical temperature $T=T_{c}$) captures low
density behavior but diverges due to a branch point singularity at the
critical density $\rho_{c}$ [30]. The behavior of fluids near the critical
point obeys a universal scaling law [31]

$\left(\frac{P}{P_{c}}-1\right)\sim\pm
D\left|\frac{\rho}{\rho_{c}}-1\right|^{\delta},~{}T=T_{c},~{}\rho\rightarrow\rho_{c},$
(7)

where $P_{c}$ is the critical pressure, $D$ is a critical amplitude, and (for
real fluids) $\delta$ is a positive non-integer critical exponent. The
knowledge of a low density series as well as an asymptotic behavior at higher
density are the necessary ingredients for the application of an asymptotic
approximant. In particular, the series (3) is accurate in a region where (7)
does not apply and inaccurate in the region where (7) is accurate. An
approximant is provided in [17] that captures both regions (low density and
the critical region); it is given by

$P_{A}=P_{c}-\sum_{n=0}^{N}A_{n}\rho^{n}\left(1-\frac{\rho}{\rho_{c}}\right)^{\delta}.$
(8a) Following the same approach of the previous subsection, one may equate
(8a) with the virial series (3) evaluated at $T=T_{c}$ and solve for the
series in $A_{n}$:
$\left(P_{c}-kT_{c}\sum_{n=1}^{N}B_{n}(T_{c})\rho^{n}\right)\left(1-\frac{\rho}{\rho_{c}}\right)^{-\delta}=\sum_{n=0}^{N}A_{n}\rho^{n}$
where it becomes clear that the $A_{n}$ coefficients are, in fact, those in
the Taylor expansion about $\rho=0$ of the left-hand side of the above
expression. The $A_{n}$ coefficients are obtained by replacing
$(1-\rho/\rho_{c})^{-\delta}$ with its expansion about $\rho=0$ and then
taking the Cauchy product (equation (27) in Appendix A) of the two series in
the left-hand side of the above expression. This leads to Equation 4b of [17]:
$A_{n>0}=\frac{P_{c}\Gamma(\delta+n)}{n!\rho_{c}^{n}\Gamma(\delta)}-\frac{kT_{c}}{\Gamma(\delta)}\sum_{j=0}^{n-1}\frac{B_{n-j}\Gamma(\delta+j)}{\rho_{c}^{j}~{}j!},~{}A_{0}=P_{c}.$
(8b) If $\rho_{c}$, $T_{c}$, and $P_{c}$ are known, one then inserts (8b) into
(8a) to construct the approximant, which captures both the low-density and
critical region, as shown in [17]. However, it is typically the case that
critical properties are poorly known and difficult to compute through
molecular simulation. Fortunately, the approximant may be used to predict
these quantities. For instance, if $\rho_{c}$ is unknown and all other
parameters are known, one may sacrifice the highest unknown coefficient
$A_{N}$ and instead predict $\rho_{c}$. This is equivalent to setting
$A_{N}$=0 in (8b), leading to Equation 7 of [17],
$P_{c}-\frac{kT_{c}N!}{\Gamma(\delta+N)}\sum_{j=1}^{N}\frac{\Gamma(\delta+N-j)}{(N-j)!}B_{j}\rho_{c}^{j}=0$
(8c)

which is a polynomial in $\rho_{c}$. Although this predicts multiple roots for
$\rho_{c}$, complex or negative values may be eliminated. As one increases $N$
and tracks the roots, typically one set of positive roots will converge to a
limit point more rapidly than others, as shown in Fig. 2a for the Lennard-
Jones model fluid. Taking the fastest converging $\rho_{c}$ sequence in the
figure and substituting it into (8b) and (8a) for each $N$ leads to a set of
approximant isotherms, shown in Fig. 2b. Note in the figure that the
approximant converges more rapidly and uniformly than the virial series as $N$
is increased.

(a)  (b)
---
Figure 2: (a) Predictions of the critical density from (8c) for the Lennard-
Jones fluid. A dashed line indicates the value predicted from Monte-Carlo
simulations [32]. (b) Critical isotherm of the Lennard-Jones fluid. A
comparison is shown between the $N$-term virial series (3) ($\\---$) labeled
as S$N$ and the corresponding approximant (11) ($\\--$) labeled as A$N$. The
most rapidly converging branch of $\rho_{c}$ roots (lowermost points) in (a)
is used as an input to the approximants shown in (b). Virial coefficients used
to generate the plots are taken from [17].

Now if both $\rho_{c}$ and $P_{c}$ are unknown, one may sacrifice the two
highest coefficients (i.e. $A_{N}=A_{N-1}=0$), leading to a system of 2
nonlinear equations to solve for $\rho_{c}$ and $P_{c}$. Recently, this
approach has been used to predict properties and construct approximants that
describe the entire critical region [18] \- not just the critical isotherm.
Results in Section 4 demonstrate that this method is also useful in the
prediction of quantities important to nonlinear boundary value problems.

As mentioned in the introduction, the virtue of approximants such as (6) and
(8) is the ease with which accurate derivative and integral properties are
generated. For example, the more comprehensive critical approximant given in
[18] for $P(\rho,T)$ can be integrated with respect to $\rho$ to accurately
calculate the Helmholtz free energy, which may then be differentiated with
respect to $T$ to find the internal energy and specific heat. Even if these
quantities are extracted numerically after an approximant is constructed, they
may be resolved on a grid of arbitrary resolution with minimal computational
expense (only memory usage), since the approximant is an analytic expression.

##  3 General procedure for constructing an asymptotic approximant

With the aim of formalizing the method of asymptotic approximants introduced
in Section 2, we provide the following algorithm. This procedure is
subsequently applied to power series solutions of nonlinear boundary-value
problems in Section 4.

###### Algorithm 3.1

Steps for constructing an asymptotic approximant for a problem whose analytic
solution $f(x)$ is elusive:

  1. i.

Find a power series solution of the given problem and truncate to $N$ terms:

$f_{S}=\sum_{n=0}^{N}a_{n}(x-x_{0})^{n}.$

Define $m$ as the number of unknowns embedded in the coefficients $a_{n}$;
this value is used in step v below. Note that step i is independent of step ii
and thus their order may be interchanged.

  2. ii.

Seek the asymptotic behavior $f\sim Cf_{a}$ as $x\to x_{a}$, where $x_{a}\neq
x_{0}$ and $C$ is a content. This may come from either an asymptotic expansion
about $x_{a}$, a boundary condition, or some independently known limiting
behavior of $f$. Define $p$ as the number of unknowns associated with the
asymptotic behavior; this value is used in step v below.

  3. iii.

Create an approximant function $f_{A}$ with $N+1$ unknowns $A_{0}\dots A_{N}$.
The function must (a) match the correct asymptotic behavior $f\sim Cf_{a}$ as
$x\to x_{a}$ to at least zeroth order and (b) be generalized to handle
arbitrary $N$; a function of an “approximant series” is well suited for this:

$f_{A}=f_{A}\left(\sum_{n=0}^{N}A_{n}[g(x)]^{n}\right).$

This is the creative step. The goal is to capture known features of $f$ while
leaving flexibility to continuously merge two (potentially disparate) regions.

  4. iv.

Solve for $A_{0}\dots A_{N}$ by imposing the condition that the $N^{th}$-order
Taylor expansion of $f_{A}$ about $x_{0}$ is equal to $f_{S}$. If $g(x)=x$,
this can be accomplished by setting $f_{S}=f_{A}$, isolating the approximant
series within $f_{A}$ to one side and then Taylor expanding the other side
using the relations in Appendix A. If $g(x)\neq x$, a matrix inversion may be
required. In either case, the resultant expressions for $A_{0}\dots A_{N}$
will be functions of the $m+p$ unknowns from steps i and ii.

  5. v.

Solve the algebraic system $A_{N}=0,~{}A_{N-1}=0,~{}\dots~{}A_{N-m-p+1}=0$ for
the $m+p$ unknowns from steps i and ii. Choose the most rapidly converging and
physically plausible roots (for example, the physical quantity the root
represents may be real and positive). In this step, unknown coefficients in
the approximant series are sacrificed in order to solve for other unknowns
embedded in these coefficients, while preserving the correct number of degrees
of freedom.

  6. vi.

Evaluate all $A_{n}$ coefficients at the values obtained in step v and
substitute these into the approximant $f_{A}$ from step iii.

  7. vii.

Ensure that $f_{A}$ converges as $N$ is increased. As an additional metric of
the approximant’s ability to approach the correct asymptotic behavior, one may
plot $f_{A}/f_{a}$ and confirm that, for increasing $N$, this ratio converges
to a constant as $x\to x_{a}$.

  8. viii.

Consistency Checks: Expand $f_{A}$ about $x_{0}$ to $N^{th}$ order; if the
steps were executed correctly this will match the first $N$ terms of $f_{S}$.
Take the limit of $f_{A}$ as $x\to x_{a}$ and confirm that it agrees with
$f_{a}$ to zeroth order. Finally, if available, verify that the approximant
agrees with a full numerical solution of the problem.

##  4 Application to nonlinear boundary value problems

Asympotic approximants are now applied to three nonlinear boundary value
problems. In Section 4.1, we find approximant solutions to the Sakiadis
boundary layer problem, where we are able to take advantage of the full
asymptotic behavior, enabling us to obtain new benchmark values. In Section
4.2, an approximant is constructed for the Blasius boundary layer problem
using the leading order asymptotic behavior; this approximant confirms
existing numerically-obtained benchmark values. In Section 4.3, an approximant
for the Flierl-Petviashvili monopole problem is constructed using the leading
order asymptotic behavior. In this application, newly obtained numerical
results are confirmed and a new asymptotic constant is explored.

###  4.1 The Sakiadis problem

The Sakiadis problem describes steady-state developing flow created by a
moving plate in an otherwise stagnant fluid [33]. It shares the same boundary
layer governing equation as that of the Blasius problem (discussed in Section
4.2), but with different boundary conditions. As shown in [33], a similarity
transform of these equations leads to the following nonlinear boundary value
problem in $f(\eta)$:

$2f^{\prime\prime\prime}+ff^{\prime\prime}=0$ (9a)
$f(0)=0,~{}f^{\prime}(0)=1,~{}f^{\prime}(\infty)=0.$ (9b)

As shown by Blasius [34], substituting the solution form

$f=\sum\limits_{n=0}^{\infty}a_{n}\eta^{n}$ (10a) into (9a) leads to the
following recursion for the series coefficients
$a_{n+3}=\frac{-\sum\limits_{j=0}^{n}(j+1)(j+2)~{}a_{j+2}~{}a_{n-j}}{2(n+1)(n+2)(n+3)}.$
(10b)

The recursion above requires knowledge of the first three coefficients in
order to generate the full series. The first two coefficients $a_{0}=f(0)$ and
$a_{1}=f^{\prime}(0)$ are given by the first two boundary conditions in (9b).
The third coefficient is prescribed by the “wall shear” parameter

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

and is typically obtained from a numerical solution of (9b) (see [35, 36, 37]
for example). Recently, a semi-analytical fixed point method was used to
estimate $\kappa$; see [38] for details and a comprehensive review and
comparison of previous numerical predictions. Here, asymptotic approximants
will be used to predict $\kappa$.

From the infinite condition in (9b), it follows that the solution limits to a
constant $C$ at large $\eta$:

$\lim_{\eta\to\infty}f\equiv C.$ (11)

Higher order asymptotic methods provided in Appendix B.1 enable the
identification of an additional asymptotic constant, $G$, as:

$\lim_{\eta\to\infty}\left[(f-C)e^{\eta C/2}\right]\equiv G.$ (12)

The constants $G$, $C$, and $\kappa$ will be predicted using asymptotic
approximants in the following subsections.

####  4.1.1 Simple asymptotic approximant

First, we consider a simple approximant of the form

$f_{A}=C-C\left(1+\sum_{n=1}^{N}A_{n}~{}\eta^{n}\right)^{-1},$ (13a) which
automatically satisfies both (9b) and (11). The assumed form (13a) has the
appearance of a Padé, but in fact it is not. Since $C$ is an unknown to be
predicted by the approximant, the above cannot be considered a Padé for the
quantity $f-C$. Also, if the terms in (13a) were combined through a common
denominator, a specific relationship between the numerator and denominator
coefficients would be obtained. Thus (13a) cannot be considered a Padé for
$f$. We shall see below that the approximant given by (13a) contains
coefficients which are, in fact, easier to compute than those of a standard
Padé.

The coefficients $A_{0}\dots A_{N}$, $C$, and $\kappa$ are calculated such
that the $N$-term Taylor expansion of (13a) about $\eta=0$ is exactly equal to
the $N$-term truncation of (10b). Equating (13a) with (10b) and re-arranging
leads to the following

$\left[1-\frac{1}{C}\sum\limits_{n=0}^{N}a_{n}\eta^{n}\right]^{-1}=1+\sum_{n=1}^{N}A_{n}~{}\eta^{n}.$

Expanding the left-hand side above (utilizing equation (28) in Appendix A)
provides the following recursion for the coefficients, which are now functions
of $\kappa$ and $C$:

$A_{n>0}=\frac{1}{C}\sum_{j=1}^{n}a_{j}~{}A_{n-j},~{}~{}A_{0}=1.$ (13b)

We now sacrifice the coefficients $A_{N}$ and $A_{N-1}$ to simultaneously
predict $\kappa$ and $C$. Setting $A_{N}=A_{N-1}=0$ in (13b), we arrive at

$\displaystyle\sum_{j=1}^{N}a_{j}A_{N-j}=0,~{}~{}~{}\sum_{j=1}^{N-1}a_{j}A_{N-1-j}=0.$
(13c)

Once the preceding coefficients are obtained from (13b) and substituted into
the above, (13c) becomes a system of nonlinear equations in $\kappa$ and $C$,
and may be solved algebraically. The most rapidly converging $\kappa$ and $C$
roots of (13c) are given in Table 1 up to an optimal asymptotic truncation of
$N=11$. That is, for $N>11$, $|\kappa_{N+1}-\kappa_{N}|$ no longer
consistently decreases as $N$ is increased.

We also record the value $S$ as the magnitude of the singularity, $\eta_{s}$,
closest to $\eta$=0 in the complex $\eta$ plane for the Sakiadis problem, as
predicted by (13) and listed in Table (1). This value limits the radius of
convergence of the series (10b).

Table 1: Predictions from approximants for the Sakiadis problem. For comparison, the numerical values of $\kappa$ and $C$ in the bottom row were computed using the Chebfun [39] package to solve (9b) via rectangular collocation with $\eta=\infty$ replaced with a finite surrogate $\eta_{\infty}$ given in the table. Note that the recent $\kappa$ estimates given in [40, 38] (reported to 7 digits) agree with the first 7 digits of our numerical predictions and predictions from (15). The numerical values of $S$ were obtained by inserting the numerical $\kappa$ and $C$ values into (13), locating the singularity of smallest magnitude, and systematically increasing the approximant order (up to $N=200$), where the value of $S$ has converged to within the reported digits. $N$ |  | $\kappa$
---
from (13)
| $C$
---
from (13)
| $S$
---
from (13)
| $\kappa$
---
from (15)
| $C$
---
from (15)
| $G$
---
from (15)
5 | -0.3879 | 2.1607 | 3.76 | -0.464241808775 | 1.57791591603 | -1.9379232426
7 | -0.4341 | 1.7275 | 3.98 | -0.445957694266 | 1.61001774844 | -2.0928590899
9 | -0.4421 | 1.6418 | 4.05 | -0.443970941603 | 1.61533742596 | -2.1255327825
11 | -0.4430 | 1.6284 | 4.07 | -0.443769095918 | 1.61603802403 | -2.1306162759
13 | $\\--$ | $\\--$ | $\\--$ | -0.443750163304 | 1.61611659263 | -2.1312640887
15 | $\\--$ | $\\--$ | $\\--$ | -0.443748473247 | 1.61612459984 | -2.1313373822
17 | $\\--$ | $\\--$ | $\\--$ | -0.443748326909 | 1.61612536888 | -2.1313450758
19 | $\\--$ | $\\--$ | $\\--$ | -0.443748314499 | 1.61612543984 | -2.1313458428
21 | $\\--$ | $\\--$ | $\\--$ | -0.443748313462 | 1.61612544619 | -2.1313459165
23 | $\\--$ | $\\--$ | $\\--$ | -0.443748313376 | 1.61612544675 | -2.1313459234
25 | $\\--$ | $\\--$ | $\\--$ | -0.443748313368 | 1.61612544681 | -2.1313459241
26 | $\\--$ | $\\--$ | $\\--$ | -0.443748313384 | 1.61612544669 | -2.1313459226
27 | $\\--$ | $\\--$ | $\\--$ | -0.443748313370 | 1.61612544680 | -2.1313459239
28 | $\\--$ | $\\--$ | $\\--$ | -0.443748313369 | 1.61612544680 | -2.1313459240
29 | $\\--$ | $\\--$ | $\\--$ | -0.443748313369 | 1.61612544681 | -2.1313459241
30 | $\\--$ | $\\--$ | $\\--$ | -0.443748313369 | 1.61612544681 | -2.1313459241
numerical | $\kappa$ | $C$ |  | $S$
---
from (13)
$\eta_{\infty}$ |  |
shooting [35] | -0.44374733 | 1.61611200 | 4.07217 | 20 |  |
quadrature [36] | -0.44374831 | 1.61612518 | 4.07217 | 23.20512 |  |
collocation | -0.44374831 | 1.61612545 | 4.07217 | 30 |  |

Lastly, $C$ and $\kappa$ are substituted into (13b) for all non-zero
coefficients ($A_{1}$ through $A_{N-2}$) to construct the approximant series
in (13a). The approximant (13) is compared with the series and numerical
solution of (9b) in Fig. 3. Although we are using a simple approximant form
that enables correspondingly straightforward coefficient generation, the
approximant converges surprisingly well, as can be seen in Fig. 3a. In Fig.
3b, we confirm that $f^{\prime}$ and $f^{\prime\prime}$ are accurately
obtained from the approximant, once $f$ has reasonably converged (here, for
$N=11$). Unfortunately, our estimates for $\kappa$ and $C$ cease to converge
for $N>11$, which limits our ability to predict the asymptotic constant $G$
given by (12). This can been seen in a plot of the effective $G$

$G_{\text{eff}}=(f-C)e^{\eta C/2},~{}G_{\text{eff}}\to
G~{}\text{as}~{}\eta\to\infty$ (14)

versus $\eta$, which, when computed by the approximant (13), only converges in
a small region near $\eta=0$ as shown in Fig. 3b. The reason for this poor
convergence behavior is that the approximant form (13) does not incorporate
the asymptotic correction behavior that leads to (12), and thus
$G_{\text{eff}}$ cannot converge uniformly.

(a)  (b)
---
Figure 3: (a) Comparison between the $N$-term series solution (10b) ($\\---$)
labeled as S$N$, the corresponding approximant (13) ($\\--$) labeled as A$N$,
and the numerical solution [35] ($\bullet$) to the Sakiadis problem (9b). (b)
Derivatives of approximant (13) with $N$=11 and $G_{\text{eff}}$ given by (14)
compared with the numerical solution [35] ($\bullet$).

####  4.1.2 Asymptotic approximant with exponential corrections

We now consider an approximant,

$f_{A}=C+\sum_{n=1}^{N}A_{n}e^{-nC\eta/2},$ (15a) that not only matches the
zeroth order infinite behavior, but also the form of the exponential
corrections from the asymptotic expansion of $f$ as $\eta\to\infty$ provided
in Appendix B.1. Note that $A_{1}$ in the series above is the unknown
asymptotic constant $G$ defined in (12). After expanding (15a) about $\eta=0$
and matching like terms with those of the series solution (10b), we arrive at
a system of $N$ linear equations,

$\left[\begin{array}[]{ccccc}1^{0}&2^{0}&3^{0}&\cdots&N^{0}\\\
1^{1}&2^{1}&3^{1}&\cdots&N^{1}\\\ 1^{2}&2^{2}&3^{2}&\cdots&N^{2}\\\
\vdots&\vdots&\vdots&\vdots&\vdots\\\
1^{N-1}&2^{N-1}&3^{N-1}&\cdots&N^{N-1}\end{array}\right]\left[\begin{array}[]{c}A_{1}\\\
A_{2}\\\ A_{3}\\\ \vdots\\\
A_{N}\end{array}\right]=\left[\begin{array}[]{c}0!~{}a_{0}-C\\\
1!~{}(-2/C)~{}a_{1}\\\ 2!~{}(-2/C)^{2}~{}a_{2}\\\ \vdots\\\
(N-1)!~{}(-2/C)^{N-1}~{}a_{N-1}\end{array}\right],$ (15b)

where the coefficient matrix is a Vandermonde matrix with a known explicit
formula for its inverse [41]. Once the system (15b) is solved, the $A_{n}$
coefficients become functions of the unknowns $C$ and $\kappa$. Like the
previous approximant, we sacrifice the last two coefficients and solve the
nonlinear equations

$A_{N}(\kappa,C)=0,~{}A_{N-1}(\kappa,C)=0$ (15c)

to obtain values for $C$ and $\kappa$, provided in the 5${}^{\text{th}}$ and
6${}^{\text{th}}$ columns of Table 1. In the table, we report values up to
$N=30$, after which the 13${}^{\text{th}}$ digit oscillates without converging
to a discernible limit point. The converged values

$\kappa=-0.443748313369$

and

$C=1.61612544681$

obtained from approximant (15) agree with numerical predictions, and provide 5
additional digits of precision. Also, for the first time, a value of the
asymptotic constant

$G=-2.1313459241$

is reported; convergence to this value is shown in the rightmost column of
Table 1, obtained by substituting the predictions for $C$ and $\kappa$ into
the expression for $A_{1}$. Once $C$ and $\kappa$ are substituted into the
expressions (arising from the solution to (15b)) for all non-zero coefficients
($A_{1}$ through $A_{N-2}$), the approximant series in (15a) may be
constructed.

The approximant (15) is compared with the series and numerical solution of
(9b) in Fig. 4. The new approximant converges more rapidly than the simpler
approximant (see Fig. 3 for comparison). Also, the new approximant captures
the exponential corrections inasmuch as $G_{\text{eff}}$ converges into the
large $\eta$ asymptotic regime and approaches $G$, as shown in Fig. 4b. This
behavior is expected, as the correct asymptotic behavior has been
incorportated explicitly into the approximant and convergence should be
uniform for large $\eta$.

(a)  (b)
---
Figure 4: (a) Comparison between the $N$-term series solution (10b) ($\\---$)
labeled as S$N$, the corresponding approximant (15) ($\\--$) labeled as A$N$,
and the numerical solution [35] ($\bullet$) to the Sakiadis problem (9b). (b)
Derivatives of approximant (15) with $N$=7 and $G_{\text{eff}}$ given by (14)
compared with the numerical solution [35] ($\bullet$).

For the Sakiadis problem, a relatively simple approximant that meets the
infinite condition in (9b), such as approximant (13), is adequate for
describing $f$ and its derivatives. However, in comparing the $G_{\text{eff}}$
curves of Figs. 3b and 4b, it is clear that an approximant such as (15) which
includes asymptotic corrections, is needed to predict $G$. As we shall see in
the next section, simple approximants can sometimes also be used to predict
higher-order asymptotic properties, provided that the approximant converges at
large enough $\eta$ for the asymptotic behavior to be reached.

Lastly, it is interesting to note that approximant (15) contains no
singularities, and thus cannot be used to predict the closest singularity
$\eta_{s}$ (and $S$) directly, as was done in Section 4.1.1 for approximant
(13). However, if we use either the new high precision values of $\kappa$ and
$C$ from approximant (15) or the numerical values listed in the table as
inputs to approximatnt (13), the location of the singularity converges to the
within the indicated digits to $\eta_{s}=-1.211393+3.88781$i and thus a radius
of convergence of $S\equiv|\eta_{s}|=4.07217$.

Although the implementation is different, the approximant (15) is similar in
style to that of the approximant of the Kidder boundary value problem given
recently in [42], where an exponential decay of the derivative at infinity is
incorporated by changing the independent variable of a Padé to reflect this
behavior. If convergence of this approximant for increasing order were
considered, this approach would fall within the framework of the asymptotic
approximants presented here.

The errors in the approximants of this section are given in Table 2 to enable
comparison with other approximate treatments of the Sakiadis problem; see for
example [38].

Table 2: The infinity norm of the error in the Sakiadis approximants and their 2nd derivatives (to within 1 significant digit), defined respectively as $E$=$||f_{A}-f_{\mathrm{numerical}}||_{\infty}$ and $E_{2}$=$||f_{A}^{\prime\prime}-f_{\mathrm{numerical}}^{\prime\prime}||_{\infty}$ on the interval $0\leq\eta\leq 20$. $N$ |  | $E$
---
from (13)
| $E_{2}$
---
from (13)
| $E$
---
from (15)
| $E_{2}$
---
from (15)
5 | 5$\times 10^{-1}$ | 6$\times 10^{-2}$ | 4$\times 10^{-2}$ | 2$\times 10^{-2}$
7 | 1$\times 10^{-1}$ | 1$\times 10^{-2}$ | 6$\times 10^{-3}$ | 2$\times 10^{-3}$
9 | 3$\times 10^{-2}$ | 1$\times 10^{-3}$ | 8$\times 10^{-4}$ | 2$\times 10^{-4}$
11 | 1$\times 10^{-2}$ | 7$\times 10^{-4}$ | 9$\times 10^{-5}$ | 2$\times 10^{-5}$
13 | $\\--$ | $\\--$ | 9$\times 10^{-6}$ | 2$\times 10^{-6}$
15 | $\\--$ | $\\--$ | 8$\times 10^{-7}$ | 2$\times 10^{-7}$
20 | $\\--$ | $\\--$ | 2$\times 10^{-9}$ | 3$\times 10^{-10}$
25 | $\\--$ | $\\--$ | 1$\times 10^{-11}$ | 5$\times 10^{-12}$
30 | $\\--$ | $\\--$ | 1$\times 10^{-11}$ | 4$\times 10^{-13}$

###  4.2 The Blasius problem

The Blasius problem is the archetypal boundary-layer problem found in most
undergraduate fluid mechanics books, describing the boundary layer due to a
moving fluid over a stationary flat plate [34]. The differential equation in
$f(\eta)$ is the same as the Sakiadis problem (and arises after a similarity
transform is applied to the governing equations),

$2f^{\prime\prime\prime}+ff^{\prime\prime}=0,$ (16a) while the boundary
conditions are now $f(0)=0,~{}f^{\prime}(0)=0,~{}f^{\prime}(\infty)=1.$ (16b)

The solution to (16b) may be constructed as a power series (10b) as in the
Sakiadis problem. However, one difference is that the Blasius series solution
skips every two terms, starting with the first two, which are zero as imposed
by the boundary conditions at $\eta$=0 in (16b). Again, the wall shear
parameter is defined as

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

which, along with the asymptotic properties [43]

$\lim_{\eta\to\infty}(f-\eta)\equiv B$ (17a)
$\lim_{\eta\to\infty}\left[\text{exp}\left(\frac{\eta^{2}}{4}+\frac{B\eta}{2}\right)f^{\prime\prime}\right]\equiv
Q,$ (17b)

will here be predicted by the asymptotic approximant, in a similar fashion as
done by [19] (for $\kappa$ using Padés). Recent numerical and analytical
predictions of these quantities, including a review of previous results may be
found in [44, 38, 45, 46].

The asymptotic behavior of the Blasius problem as $\eta\to\infty$ is given as
(see Appendix B.2)

$f\sim\eta+B+4Q\frac{\text{exp}[-\eta^{2}/4-B\eta/2]}{(\eta+B)^{2}}[1+O(\frac{1}{(\eta+B)^{2}})],~{}\textrm{as}\
\eta\rightarrow\infty.$ (18)

For comparison, the exponential correction analogous to (18) for the Sakiadis
problem is (36) in Appendix B.1. While we were able to incorporate (12)
directly in the Sakiadis approximant (15), this is not an option for the
Blasius problem since the asymptotic form (18) contains a singularity at
$\eta=-B$ and it is known that $B<0$ [43]. To clarify this issue, note that
(18) is fully valid as $\eta\to\infty$. However, our methodology for
constructing an approximant relies on the unification of asymptotic limits
such that the final form is capable of describing all $\eta>0$, and here is
where (18) poses an issue in its ability to inform such an approach.

While we cannot directly use the asymptotic expansion for the Blasius problem,
we can use a simple approximant that agrees asymptotically to
1${}^{\text{st}}$ order (i.e. with (16b) and (17a)) and then verify that it
approaches the exponential correction (17b) as $\eta\to\infty$. We use an
approximant of the form

$f_{A}=\eta+B-B\left(1+\sum_{n=1}^{N}A_{n}~{}\eta^{n}\right)^{-1},$ (19a)
which automatically satisfies both (16b) and (17a). As similarly discussed in
the context of the Sakiadis Approximant (13a), the above form is not a Padé
approximant. Firstly, $B$ is an unknown parameter in the approximant, and so
(19a) cannot be considered a Padé for $f-\eta-B$. Secondly, if one combines
the terms of (19a) through a common denominator, the coefficients of the
numerator will have an explicit dependence on those in the denominator,
whereas this is not the case for standard Padés.

The coefficients $A_{0}\dots A_{N}$, $B$, and $\kappa$ are now calculated such
that the $N$-term Taylor expansion of (19a) about $\eta=0$ is exactly equal to
the $N$-term truncation of (10b). Following the same procedure as in Section
4.1.1, we arrive at the following recursion for the coefficients, which are
now functions of $\kappa$ and $B$:

$A_{n>0}=\frac{1}{B}\sum_{j=1}^{n}\tilde{a}_{j}~{}A_{n-j},~{}~{}A_{0}=1,$
(19b)

where $\tilde{a}_{1}=-1$ and $\tilde{a}_{j>1}=a_{j}$. The coefficients $A_{N}$
and $A_{N-1}$ are now sacrificed to simultaneously predict $\kappa$ and $B$.
Setting

$A_{N}(\kappa,B)=0,~{}A_{N-1}(\kappa,B)=0$ (19c)

in (19b) leads to two nonlinear equations, whose most rapidly converging
$\kappa$ and $B$ roots are given in Table 3, along with the magnitude of the
singularity closest to $\eta=0$ in the approximant, denoted by $S$. The table
reports values up to $N=50$, beyond which convergence cannot be established in
the 8th digit. Our predictions converge to within 7 digits of the benchmark
value for $\kappa$, 6 digits for $B$, and 2 digits for the radius of
convergence $S$. Again, we determine $S$ as the magnitude of the closest
singularity, $\eta_{s}$, from the origin in the complex $\eta$ plane. Although
not shown in the table, the approximant (19) is consistent with the literature
[43] in that, as $N$ is increased, it predicts $\eta_{s}$ to lie on the
negative real axis.

Table 3: Predictions from the Blasius approximant (19). The infinity norm of the error in the approximant and its 2nd derivative are recorded below (to within 1 significant digit), defined respectively as $E$=$||f_{A}-f_{\mathrm{numerical}}||_{\infty}$ and $E_{2}$=$||f_{A}^{\prime\prime}-f_{\mathrm{numerical}}^{\prime\prime}||_{\infty}$ on the interval $0\leq\eta\leq 8.8$. $N$ | $\kappa$ | $B$ | $S$ | $E$ | $E_{2}$
---|---|---|---|---|---
5 | 0.12545065 | -6.08950239 | 9.85302 |  3$\times 10^{0}$ |  2$\times 10^{-1}$
10 | 0.30018461 | -2.00003122 | 4.47053 |  3$\times 10^{-1}$ |  3$\times 10^{-2}$
15 | 0.32626056 | -1.77168204 | 4.59756 |  5$\times 10^{-2}$ |  5$\times 10^{-3}$
20 | 0.33090243 | -1.73173129 | 4.92273 |  1$\times 10^{-2}$ |  1$\times 10^{-3}$
25 | 0.33181978 | -1.72323142 | 5.18569 |  2$\times 10^{-3}$ |  2$\times 10^{-4}$
30 | 0.33200793 | -1.72133641 | 5.43362 |  5$\times 10^{-4}$ |  5$\times 10^{-5}$
35 | 0.33204717 | -1.72090879 | 5.64460 |  1$\times 10^{-4}$ |  1$\times 10^{-5}$
40 | 0.33205518 | -1.72081494 | 5.69196 |  3$\times 10^{-5}$ |  2$\times 10^{-6}$
45 | 0.33205693 | -1.72079309 | 5.68929 |  5$\times 10^{-6}$ |  4$\times 10^{-7}$
50 | 0.33205731 | -1.72078801 | 5.68933 |  3$\times 10^{-7}$ |  3$\times 10^{-8}$
numerical [45, 47] | 0.332057336215196 | -1.7207876575205 | 5.6900380545 | 0 | 0

Once $B$ and $\kappa$ are substituted into (19b) for all non-zero coefficients
($A_{1}$ through $A_{N-2}$), the approximant series in (19) may be
constructed. The approximant (19) is compared with the series and numerical
solution of (16b) in Fig. 5. Like the simple Sakiadis approximant given in
Section (4.1.1), this basic Blasius approximant converges surprisingly well,
as can be seen in Fig. 5a. In Fig. 5b, we confirm that $f^{\prime}$ and
$f^{\prime\prime}$ are accurately obtained from the approximant, once $f$ has
reasonably converged (here, for $N=30$). Unlike the simple Sakiadis
approximant, the Blasius approximant is capable of picking up a higher-order
asymptotic quantity - namely $Q$ defined in (17b). This is seen in Fig. 5c,
where an effective $Q$, defined as

$Q_{\text{eff}}=\text{exp}\left(\frac{\eta^{2}}{4}+\frac{B\eta}{2}\right)f^{\prime\prime},~{}Q_{\text{eff}}\to
Q~{}\text{as}~{}\eta\to\infty$ (20)

is plotted versus $\eta$ for various series and approximant order. While the
approximant converges directly to a constant $Q$ value at large $\eta$, it is
striking that the series solution also reaches the $Q$ plateau prior to
diverging. This may be a result of the Blasius series having a region of
convergence that overlaps with asymptotic effects felt at relatively small
$\eta$. As seen in Fig. 5c, both the approximant and the numerical solution
converge to $Q\approx 0.1115$, which is consistent with the value
$Q$=0.111483755 obtained numerically by [45] and semi-analytically (i.e.
dependent on numerical $\kappa$) in [46]111$Q$ in [46] is defined as
$e^{-B^{2}/4}$ multiplied by the $Q$ defined here.; in both references, the
value is reported well-beyond these digits, but both round to the value above.
Since the behavior (20) is not explicitly incorporated into the approximant
(19), the convergence of $Q_{\text{eff}}$ is not uniform at large $\eta$. For
example, in Fig. 5c, the behavior of the $N=20,30,40$ approximants (labeled
A20, A30, A40) demonstrates the non-uniform convergence that will occur for
any order $N$ at some value of $\eta$.

(a)  (b)
---
(c)
Figure 5: (a) Comparison between the $N$-term series solution (10b) ($\\---$)
labeled as S$N$ and the corresponding approximant (19) ($\\--$) labeled as
A$N$ to the Blasius problem (16b). (b) Derivatives of approximant (19) with
$N$=30. (c) $Q_{\text{eff}}$ given by (20) as predicted by the series (10b)
and approximant (19). For comparison, the numerical solution from [45] is
shown in all above figures as $\bullet$s.

The errors in the approximant are given in the last two columns of Table 3 to
enable comparison with other approximate treatments of the Blasius problem
[19, 43, 44, 38, 46].

The virtue of approximant (19) is that it is globally accurate from small to
large $\eta$, is able to predict $\kappa$ and $B$ (i.e. does not rely on these
values as inputs), and has coefficients that are generated from a simple
recursion. A recent alternative globally accurate approximation to the Blasius
solution is provided in [46]. An attractive feature of this alternative
approximation is that it has a simpler form than the asymptotic approximant
(19). However, the form in [46] requires that $\kappa$ is known beforehand.
Nevertheless, one may directly compare the error in $f_{A}^{\prime\prime}$
given in Table 3 with that given in figure 2 of [46].

###  4.3 The Flierl-Petviashvili problem

The Flierl-Petviashvili (FP) equation in similarity variables $u(r)$ is used
to describe vortex solitons in the ocean, atmosphere [48], and Jupiter’s red
spot [49], and is given by

$u^{\prime\prime}+\frac{1}{r}u^{\prime}-u-u^{2}=0$ (21a)
$u^{\prime}(0)=0,~{}~{}u(\infty)=0.$ (21b)

The coefficients of the power series solution

$u=\sum\limits_{n=0}^{\infty}a_{n}r^{n}$ (22a) to (21b) are given by
$a_{n+2}=\frac{a_{n}+\sum\limits_{k=0}^{n}a_{k}a_{n-k}}{(n+2)^{2}},$ (22b)

which requires the specification of the first two coefficients to generate the
remaining even coefficients; all remaining odd terms are zero. The coefficient
$a_{1}$=0 is known from the boundary conditions (21b) and $a_{0}$ is an
unknown parameter,

$a_{0}\equiv z$

which, along with the asymptotic property

$\lim_{r\to\infty}(ue^{r}\sqrt{r})\equiv D$ (23)

will be predicted by the approximant, in a similar fashion as done by [19, 50]
(for $z$ using Padés). The property (23) above follows from an asymptotic
expansion of the FP problem as $r\to\infty$ (derived in Appendix B.3):

$u\sim D\frac{e^{-r}}{\sqrt{r}}[1+O(\frac{1}{r})],~{}\textrm{ as
}r\rightarrow\infty.$ (24)

As was the case for the Blasius problem, a singularity in the asymptotic form
(here at $r$=0 in the above) prevents direct incorporation of (24) to
construct an approximant. However, like the Blasius approximant of Section
4.2, we shall again use a simple approximant that agrees asymptotically to
zeroth order (i.e. with the 2nd condition of (21b)) and then verify that it
approaches the correction (24) as $r\to\infty$.

We consider a simple approximant of the form

$u_{A}=\frac{z}{1+\sum\limits_{n=1}^{N}A_{n}~{}r^{n}},$ (25a) which
automatically satisfies both boundary conditions of (21b), including the
condition $u(\infty)=0$ not explicitly captured by (22). Note that (25a) is a
Padé. The simple structure of this Padé, however, allows for the unknowns to
be calculated without inverting a matrix. The coefficients $A_{0}\dots A_{N}$
and $z$ in (25a) are now calculated such that the $N$-term Taylor expansion of
(25a) about $r=0$ is exactly equal to the $N$-term truncation of (22).
Following the same procedure as in Section 4.1.1, we arrive at the following
recursion for the coefficients, which are now functions of $z$:
$A_{n>0}=-\frac{1}{z}\sum_{j=1}^{n}a_{j}~{}A_{n-j},~{}~{}A_{0}=1.$ (25b) The
coefficients above mimic the series (22) in that it skips odd coefficients.
The unknown $z$ in (25b) is determined by sacrificing the degree of freedom
normally used to compute $A_{N}$. This is equivalent to setting $A_{N}=0$ in
(25b), leading to $\sum_{j=1}^{N}a_{j}A_{N-j}=0.$ (25c)

Note that since the $a_{j}$ coefficients are polynomials in $z$, (25c) is also
a polynomial in $z$. The most rapidly converging $z$-roots of (25c) (for
increasing $N$) are given in Table 4. Once $z$ is substituted into (25b) for
all non-zero coefficients ($A_{1}$ through $A_{N-1}$), the approximant series
in (25a) may be constructed. The magnitude of the singularity closest to $r=0$
in the approximant (25) is also given in Table 4, denoted by $S$. Although not
shown in the table, the approximant (25a) predicts that, within the reported
digits, the singularity appears to lie on the positive imaginary axis.

Table 4: Predictions from the Flierl-Petviashvili approximant (25a). The numerical method for computing the values in the bottom row is described in Appendix C. The infinity norm of the error in the approximant is also recorded below (to within 1 significant digit), defined as $E$=$||u_{A}-u_{\mathrm{numerical}}||_{\infty}$ on the interval $0\leq r\leq 10$.  $N$ | $z$ | $S$ | $E$
---|---|---|---
4 | -1.5 | 2.82843 | 9$\times 10^{-1}$
6 | -2.14039 | 2.57284 | 3$\times 10^{-1}$
8 | -2.34792 | 2.47227 | 4$\times 10^{-2}$
10 | -2.38564 | 2.60878 | 6$\times 10^{-3}$
12 | -2.39117 | 2.59927 | 7$\times 10^{-4}$
14 | -2.39185 | 2.61160 | 2$\times 10^{-4}$
16 | -2.39196 | 2.61154 | 1$\times 10^{-5}$
numerical | -2.3919564032 | 2.611541077 | 0

The approximant (25) is compared with the series and numerical solution of
(21b) in Fig. 6. As for the simple Sakiadis approximant given in Section 4.1.1
and the Blasius approximant of Section 4.2, the FP approximant also converges
uniformly, as can be seen in Fig. 6a. In Fig. 6b, we confirm that $u^{\prime}$
is accurately obtained from the approximant once $u$ has reasonably converged
(here, for $N=12$). Like the Blasius approximant, the FP approximant is
capable of picking up a higher-order asymptotic quantity - here, $D$ defined
in (23). This is seen in Fig. 6b, where an effective $D$, defined as

$D_{\text{eff}}=ue^{r}\sqrt{r},~{}D_{\text{eff}}\to
D~{}\text{as}~{}r\to\infty$ (26)

is plotted versus $r$ for various approximant order. As $N$ is increased, the
approximant $D_{\text{eff}}$ appears to be converging to the numerical
solution, which is in turn converging to a value of

$D\approx-10.7$

which may serve as a useful metric for future analytical and numerical
solutions of the FP problem. Again, Fig. 6 shows that convergence in
$D_{\text{eff}}$ is not uniform at large $r$, as expected by the assumed form
of the approximant.

(a)  (b)
---
Figure 6: (a) Comparison between the $N$-term series solution (22) ($\\---$)
labeled as S$N$ and the corresponding approximant (25) ($\\--$) labeled as
A$N$ to the Flierl-Petviashvilli problem (21b). (b) Derivatives of approximant
(25) with $N$=12 and $D_{\text{eff}}$ given by (26). For comparison, the
numerical solution (described in Appendix C) is shown in all above figures as
$\bullet$s.

The errors in the approximant are given in the rightmost column of Table 4 to
enable comparison with other approximate treatments of the FP equation; see
for example [19].

##  5 Summary

In this work, we formalize a new approach to sum series. We provide both a
methodology and examples that demonstrate how asymptotic approximants may be
constructed. The key feature is that such approximants are designed so that
the asymptotic behaviors in two regions of the domain may be joined. This
approach has seen recent success in the analytic continuation of “virial”
series that describe the pressure-density-temperature dependence of various
model fluids, where no differential equation is available and usually only the
first 3 to 12 terms of these (often divergent) power series are known [16, 17,
18].

Here, we provide an additional application for asymptotically consistent
approximants - namely, a remedy for divergent power series that arise as
solutions to nonlinear ordinary differential equations. The analytic forms
provided here enable symbolic differentiation and thus allow analytical
evaluation of the flow field (shear stress, vorticity, etc.) to any desired
resolution at low computational expense. Additionally, the approximant
coefficients are relatively straightforward to obtain, as they often only
require the implementation of a few series identities, enabling them to be
generated typically via a simple recurrence relation.

Note that improvements to the forms given here are easily made, depending on
how many parameters are to be embedded in (or predicted by) the approximant.
Here, we use the approximant to not only find a solution, but also provide the
as-of-yet most accurate and precise values of the wall shear, asymptotic
constants, and singularity of smallest magnitude for the Sakiadis boundary
layer. If better estimates become available from numerical (or other)
techniques, these may be used as inputs to the approximant to improve global
accuracy. For the Blasius problem and Flierl-Petviashvilli problem, more
accurate values for some of these properties are, in fact, available and could
have been used as inputs to the approximants; we chose instead to predict
these quantities to demonstrate that the method is self contained and does not
rely on external data.

The results presented here motivate the development and application of
asymptotic approximants to other problems where divergent, truncated, and
underspecified series arise. Considering the capabilities afforded and the
ease with which asymptotic approximants are generated, they deserve
consideration as an alternative to the more traditional Padé approximants in
problems of mathematical physics.

## References

  * [1] E. C. Titchmarsh.  The Theory of Functions.  Oxford University Press, 2 edition, 1968.
  * [2] A. J. Masters.  Virial expansions.  J. Phys.: Condens. Matter, 20(283102), 2008.
  * [3] C. M. Bender and S. A. Orszag.  Advanced Mathematical Methods for Scientists and Engineers I: Asymptotic Methods and Perturbation Theory.  McGraw-Hill, 1978.
  * [4] G. A Baker Jr.  Quantative Theory of Critical Phenomena.  Academic Press, London, 1990.
  * [5] J. S. R. Chisholm.  Generalisations of Padé approximants.  Circuits Syst. Signal Process, 1:279, 1982.
  * [6] A. J. Guttmann.  Asymptotic analysis of power series expansions.  In Phase Transitions and Critical Phenomenon Vol. 13, pages 1–234. Academic Press, NY, 1989.
  * [7] A. J. Guttmann and I. Jensen.  Series analysis.  In Polygons, Polyominoes, and Polycubes, pages 181–202. Springer, 2009.
  * [8] G. A. Baker Jr. and P. Graves-Morris.  Padé approximants.  Cambridge University Press, 2 edition, 1996.
  * [9] I. V. Andrianov and J. Awrejcewicz.  New trends in asymptotic approaches: summation and interpolation methods.  Appl Mech. Rev., 54(1), 2001.
  * [10] P. A. Frost and E. Y. Harper.  An extended padé procedure for constructing global approximations from asymptotic expansions: an explication with examples.  SIAM Rev., 18(1):62–91, 1976.
  * [11] G. A Baker Jr. and J. L. Gammel.  The Padé approximant.  J. Math. Anal. Appl., 2:21–30, 1961.
  * [12] M. Van Dyke.  Analysis of improvement of perturbation series.  Q. J. Mech. Appl. Math., 27(4):423–450, 1974.
  * [13] N. Clisby and B. M. McCoy.  Ninth and tenth order virial coefficients for hard spheres in $d$ dimensions.  J. Stat. Phys., 122(1):15–55, 2006.
  * [14] A. Guerrero and A. Bassi.  On Padé approximants to virial series.  J. Chem. Phys., 129(044509), 2008.
  * [15] T. B. Tan, A. J. Schultz, and D. A. Kofke.  Virial coefficients, equation of state, and solid–fluid coexistence for the soft sphere model.  Mol. Phys., 109(1):123–132, 2011.
  * [16] N. S. Barlow, A. J. Schultz, S. J. Weinstein, and D. A. Kofke.  An asymptotically consistent approximant method with application to soft- and hard-sphere fluids.  J. Chem. Phys., 137:204102, 2012.
  * [17] N. S. Barlow, A. J. Schultz, S. J. Weinstein, and D. A. Kofke.  Critical isotherms from virial series using asymptotically consistent approximants.  AIChE J., 60(9):3336–3349, 2014.
  * [18] N. S. Barlow, A. J. Schultz, S. J. Weinstein, and D. A. Kofke.  Communication: Analytic continutation of the virial series through the critical point using parametric approximants.  J. Chem. Phys., 143:071103:1–5, 2015.
  * [19] J. P. Boyd.  Padé approximant algorithm for solving nonlinear ordinary differential equation boundary value problems on an unbounded domain.  Comp. Phys., 11:299–303, 1997.
  * [20] E. A. Mason and T. H. Spurling.  The Virial Equation of State.  Pergammon Press, Oxford, 1969.
  * [21] D. A. McQuarrie.  Statistical Mechanics.  University Science Books, 2000.
  * [22] J. K. Singh and D. A. Kofke.  Mayer sampling: Calculation of cluster integrals using free-energy perturbation methods.  Phys. Rev. Lett., 92:220601, 2004.
  * [23] Richard J Wheatley.  Calculation of High-Order Virial Coefficients with Applications to Hard and Soft Spheres.  Physical Review Letters, 110(20):200601, May 2013.
  * [24] N. F. Carnahan and K. E. Starling.  Equation of state for nonattracting rigid spheres.  J. Chem. Phys., 51(2):635–636, 1969.
  * [25] J.-P. Hansen and I. R. McDonald.  Theory of Simple Liquids.  Academic Press, London, 3rd edition, 2006.
  * [26] P. Henrici.  Automatic computations with power series.  JACM, 3:10–15, 1956.
  * [27] F. J. Rogers and D. A. Young.  New, thermodynamically consistent, integral equation for simple fluids.  Phys. Rev. A, 30(2):999–1007, 1984.
  * [28] G. A Baker Jr.  Application of the Padé approximant method to the investigation of some magnetic properties of the Ising model.  Phys. Rev., 124(3):768–774, 1961.
  * [29] C. J. Thompson.  Mathematical Statistical Mechanics, chapter 6.5: Numerical Analysis of the Three-dimensional Ising Model.  Princeton University Press, New Jersey, 1972.
  * [30] M. E. Fisher.  Notes, definitions, and formulas for critical point singularities.  In M. S. Green and J. V. Sengers, editors, Critical phenomena, Proceedings of a conference held in Washington, pages 21–26, April 1965.
  * [31] H Behnejad, J. V. Sengers, and M. A. Anisimov.  Thermodynamic behaviour of fluids near critical points.  In Applied Thermodynamics of Fluids, pages 321–367. Royal Society of Chemistry, 2010.
  * [32] J. Pérez-Pellitero, P. Ungerer, G. Orkoulas, and A. Mackie.  Critical point estimation of the Lennard-Jones pure fluid and binary mixtures.  J. Chem. Phys., 125:054515, 2006.
  * [33] B. C. Sakiadis.  Boundary-layer behavior on continuous solid surfaces: II the boundary layer on a continuous flat surface.  AlChE J., 7:221–225, 1961.
  * [34] H. Blasius.  Grenzschichten in flussigkeiten mit kleiner reibung.  Zeitschrift fur Mathematik und Physik, 56:1–37, 1908.
  * [35] R. Cortell.  Numerical comparisons of Blasius and Sakiadis flows.  MATEMATIKA, 26(2):187–196, 2010.
  * [36] S. A. Eftekhari and A. A. Jafari.  Numerical solution of general boundary layer problems by the method of differential quadrature.  Sci. Iran. B, 20(4):1278–1301, 2013.
  * [37] R. Fazio.  The iterative transformation method for the Sakiadis problem.  Comp. Fluids, 106:196–200, 2015.
  * [38] D Xu and X. Guo.  Application of fixed point method to obtain semi-analytical solution to Blasius flow and its variation.  Applied Mathematics and Computation, 224:791–802, 2013.
  * [39] T. A. Driscoll, N. Hale, and L. N. Trefethen, editors.  Chebfun Guide.  Pafnuty Publications, Oxford, 2014.
  * [40] H. I. Andersson and J. B. Aarseth.  Sakiadis flow with variable fluid properties revisited.  I. J. Eng. Sci., 45(2-8):554–561, 2007.
  * [41] N. Macon and A. Spitzbart.  Inverses of Vandermonde matrices.  Am. Math. Monthly, 65(2):95–100, 1958.
  * [42] R. Iacono and J. P. Boyd.  The Kidder equation: $u_{xx}+2xu_{x}/\sqrt{1-\alpha u}=0$.  Stud. Appl. Math., 135:63–85, 2014.
  * [43] J. P. Boyd.  The Blasius function in the complex plane.  Exper. Math., 8(4):381–394, 1999.
  * [44] B. Yun.  Intuitive approach to the approximate analytical solution for the Blasius problem.  Applied Mathematics and Computation, 215:3489–3494, 2010.
  * [45] S. Anil Lal and Neeraj Paul M.  An accurate taylors series solution with high radius of convergence for the Blasius function and parameters of asymptotic variation.  JAFM, 7(4):557–564, 2014.
  * [46] R. Iacono and J. P. Boyd.  Simple analytic approximations for the blasius problem.  Physica D, 310:72–78, 2015.
  * [47] J. P. Boyd.  The Blasius function: Computations before computers, the value of tricks, undergraduate projects, and open research problems.  SIAM Rev., 50(4):791–804, 2008.
  * [48] G. R. Flierl.  Baroclinic solitary waves with radial symmetry.  Dyn. Atmos. Oceans, 3:15–38, 1979.
  * [49] V. I. Petviashvili.  Red spot of jupiter and the drift soliton in a plasma.  JETP Lett., 32:619–622, 1981.
  * [50] A. M. Wazwaz.  The Volterra integro-differential forms of the singular Flierl-Petviashvili and the Lane-Emden equations with boundary conditions.  Rom. J. Phys., 58:685–693, 2013.
  * [51] R. V. Churchill.  Complex Variables.  McGraw-Hill, 1948.
  * [52] M. Abramowitz and I Stegun.  Handbook of Mathematical Functions, page 298.  Dover, 1972.
  * [53] P. R. Turner.  Guide to scientific computing.  CRC press, 2 edition, 2001.
  * [54] M. Maleki, I. Hashim, and S. Abbasbandy.  Analysis of IVPs and BVPs on semi-infinite domains via collocation methods.  Journal of Applied Mathematics, 2012, 2012.
  * [55] C. Domb and M. F. Sykes.  On the susceptibility of a ferromagnetic above the curie point.  Proc. Roy. Soc. London A, 240(1221):214–228, 1957.

##  Appendix A Useful Series Formulae

The following relations may be used to develop a recursion for the
coefficients of an asymptotic approximant, allowing one to avoid solving an
algebraic system. The first relation is the well-known Cauchy product of two
series [51]:

$\sum_{n=0}^{N}a_{n}x^{n}\sum_{n=0}^{N}A_{n}x^{n}=\sum_{n=0}^{N}\left(\sum_{j=0}^{n}a_{j}A_{n-j}\right)x^{n}.$
(27)

Setting both sides of (27) equal to one and rearranging, the recursion leads
to a representation for the $N$-term expansion of the reciprocal of a series:

$\left(\sum_{n=0}^{N}a_{n}x^{n}\right)^{-1}=\sum_{n=0}^{N}A_{n}x^{n},$ (28a)
where
$A_{n>0}=-\frac{1}{a_{0}}\sum_{j=1}^{n}a_{j}A_{n-j},~{}A_{0}=\frac{1}{a_{0}}.$
(28b)

The generalization of (28) for the $N$-term expansion of a series raised to
any real power $s$ is given by J. C. P. Miller’s formula [26]:

$\left(\sum_{n=0}^{N}a_{n}x^{n}\right)^{s}=\sum_{n=0}^{N}A_{n}x^{n},$ (29a)
where
$A_{n>0}=\frac{1}{n~{}a_{0}}\sum_{j=1}^{n}(js-n+j)a_{j}A_{n-j},~{}A_{0}=(a_{0})^{s}.$
(29b)

##  Appendix B  Asymptotic Expansions

###  B.1 Sakiadis Problem

We wish to determine the $\eta\rightarrow\infty$ asymptotic behavior of the
solution $f(\eta)$ to the Sakiadis equation (9b). For purposes of the
analysis, the infinite constraint (11) alone suffices. We write:

$f=C+h(\eta),\ \textrm{with}\ h\rightarrow 0\ \textrm{as}\
\eta\rightarrow\infty$ (30)

where $C$ is a constant and $h(\eta)$ is a function to be determined. The form
(30) is substituted into equation (9a) to obtain

$2h^{\prime\prime\prime}+(C+h)h^{\prime\prime}=0.$ (31)

The equation (31) may be simplified by noting that terms quadratic in $h$ are
small relative to other terms as $\eta\to\infty$, i.e.:

$2h^{\prime\prime\prime}+Ch^{\prime\prime}=0,~{}\ \textrm{as}\
\eta\rightarrow\infty.$ (32)

The solution of this equation is

$h\sim Ge^{-C\eta/2},~{}\ \textrm{as}\ \eta\rightarrow\infty,$ (33)

where $G$ is a constant. Higher order corrections can be obtained by assuming
the form

$h\sim Ge^{-C\eta/2}+D(\eta),~{}\ \textrm{as}\ \eta\rightarrow\infty,$ (34)

and applying the method of dominant balance [3] when substituted into (31),
which leads to

$D(\eta)\sim\frac{G^{2}}{4C}e^{-C\eta},~{}\ \textrm{as}\
\eta\rightarrow\infty.$ (35)

The asymptotic behavior is written concisely by combining equation (30), (34),
and (35) to yield:

$f\sim C+Ge^{-C\eta/2}+\frac{G^{2}}{4C}e^{-C\eta}+O(e^{-3C\eta/2}),\
\textrm{as}\ \eta\rightarrow\infty.$ (36)

The above process for obtaining corrections may be repeated and leads to a
series of exponentials with arguments $-nC\eta/2$ for $n=1,2,3,4\dots$.

###  B.2 Blasius Problem

We wish to determine the $\eta\rightarrow\infty$ asymptotic behavior of the
solution $f(\eta)$ to the Blasius problem (16b). For purposes of the analysis,
the infinite constraint (17a) alone suffices. We write:

$f=\eta+B+g(\eta),\ \textrm{with}\ g\rightarrow 0\ \textrm{as}\
\eta\rightarrow\infty,$ (37)

where B is an unknown constant and $g(\eta)$ is a function to be determined.
The form (37) is substituted into (16b) to obtain

$2g^{\prime\prime\prime}+(\eta+B+g)g^{\prime\prime}=0.$ (38)

The equation (38) may be simplified by noting that terms quadratic in $g$ are
small relative to other terms as $\eta\to\infty$, i.e.:

$2g^{\prime\prime\prime}+(\eta+B)g^{\prime\prime}=0,\ \textrm{as}\
\eta\rightarrow\infty.$ (39)

The solution of this equation is

$g^{\prime\prime}\sim Qe^{[-(\eta^{2}+2B\eta)/4]},~{}\textrm{as}\
\eta\rightarrow\infty,$ (40)

where Q is an unknown constant. We can then write:

$g^{\prime}\sim
Qe^{(B^{2}/4)}\int_{\infty}^{\eta}e^{[-(\xi+B^{2})/4]}d\xi,~{}\textrm{as}\
\eta\rightarrow\infty$ (41)

and after application of integration by parts, we obtain:

$g^{\prime}\sim-2Qe^{(B^{2}/4)}\frac{e^{[-(\eta+B^{2})/4]}}{\eta+B}[1+O(\frac{1}{(\eta+B)^{2}})],~{}\textrm{as}\
\eta\rightarrow\infty.$ (42)

The process can then be repeated to find an expression for $g$, i.e.:

$g\sim
2Qe^{(B^{2}/4)}\int_{\infty}^{\eta}\frac{e^{[-(\xi+B^{2})/4]}}{\xi+B}d\xi,~{}\textrm{as}\
\eta\rightarrow\infty,$ (43)

from which integration by parts yields

$g\sim
4Qe^{(B^{2}/4)}\frac{e^{[-(\eta+B^{2})/4]}}{(\eta+B)^{2}}[1+O(\frac{1}{(\eta+B)^{2}})],~{}\textrm{as}\
\eta\rightarrow\infty,$ (44)

which can be combined with the original assumed form (37) to obtain the
asymptotic behavior:

$f\sim\eta+B+4Q\frac{\exp[-\eta^{2}/4-B\eta/2]}{(\eta+B)^{2}}[1+O(\frac{1}{(\eta+B)^{2}})],~{}\textrm{as}\
\eta\rightarrow\infty.$ (45)

###  B.3 Flierl-Petviashvili Problem

We wish to determine the $r\rightarrow\infty$ asymptotic behavior of the
solution $u(r)$ to the Flierl-Petviashvili problem (21b). For purposes of the
analysis, the infinite constraint in (21b) alone suffices. Since $u\rightarrow
0$ as $r\rightarrow\infty$, it is seen by inspection of (21b) that $u^{2}<<u$
as $r\rightarrow\infty$, leading to the linear equation

$r^{2}u^{\prime\prime}+ru-r^{2}u=0,~{}\textrm{ as }r\rightarrow\infty.$ (46)

The solution of (46) is a modified Bessel function of zeroth order [52](p. 374
equation 9.6.1). The leading-order behavior of this function is given by
equation 9.7.2 (p. 378) of [52] as

$u\sim D\frac{e^{-r}}{\sqrt{r}}[1+O(\frac{1}{r})],~{}\textrm{ as
}r\rightarrow\infty,$ (47)

where $D$ is an unknown constant.

##  Appendix C Numerical Solution of the Flierl-Petviashvili Problem

The numerical solution to the Flierl-Petviashvili equation used in Table 4 is
obtained by first applying the transformation $\alpha=1-e^{-r}$ to (21b),
leading to

$\displaystyle(1-\alpha)^{2}\frac{d^{2}u}{d\alpha^{2}}-[(1-\alpha)+\frac{1+\alpha}{\ln(1-\alpha)}]\frac{du}{d\alpha}-u-u^{2}=0$
$\displaystyle\frac{du}{d\alpha}(0)=0,~{}u(1)=0.$ (48)

Note that, like the original FP equation (21b), the solution to (48) is not
unique. There is a trivial solution $u=0$ and a non-trivial solution - we seek
the latter. The “Shooting Method” [53] (i.e. iterating on $u(0)=z$ guesses) is
used to effectively convert (48) into an initial value problem, which is
solved via 4th order Runge-Kutta. The boundary condition $u(1)=0$, which
informs the shooting iteration, is replaced with $u(1-\varepsilon)=0$, where
$\varepsilon<<1$. Table 5 lists numerical predictions of $z$ for decreasing
$\varepsilon$. For each value of $\varepsilon$, the numerical $\alpha$ step
size used in the Runge-Kutta implementation was successively decreased until
convergence was established to within the reported digits of the table.

Table 5: Predictions of $u(0)\equiv z$ of (48), obtained numerically by shooting to the boundary condition $u(1-\varepsilon)=0$. The reported digits are within the convergence tolerance of the shooting iteration. For comparison, the previous benchmark values reported by [19] and [54] are $z=-2.3919564$ and $-2.391956403$, respectively.  $\varepsilon$ | $z$
---|---
$10^{-1}$ | -2.6972451986
$10^{-2}$ | -2.4011647700
$10^{-3}$ | -2.3921055641
$10^{-4}$ | -2.3919584270
$10^{-5}$ | -2.3919564287
$10^{-6}$ | -2.3919564035
$10^{-7}$ | -2.3919564032
$10^{-8}$ | -2.3919564032

Using the high precision $z$ value of -2.3919564032 in Table 5, the magnitude
$S$ of the closest singularity to the origin in the complex $r$-plane is
computed by applying the ratio-test to (22b) and making a Domb-Sykes plot of
the ratios versus inverse coefficient order [55], shown in Fig. 7 and taken up
to $n=2000$. A linear fit is then made near $1/(n+1)$=0 with an $n=\infty$
intercept of $S=2.611541077$, where the digits have been conservatively
truncated since the residual of the linear fit is of order $10^{-12}$.

Figure 7: Domb-Sykes ratio plot of the coefficients in (22b) using
$a_{0}\equiv z=-2.3919564032$ from Table 5. The intercept at $1/(n+1)$=0 is
2.611541077.