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

TL;DR

This paper introduces a method called asymptotic approximants to accurately sum divergent or incomplete power series solutions by leveraging known asymptotic behavior, demonstrated on nonlinear physics problems.

Contribution

The paper formalizes and provides an algorithm for constructing asymptotic approximants that bridge asymptotic regions, enabling accurate solutions from limited series data.

Findings

Successfully constructed approximants for three nonlinear physics problems.

Predicted unknown coefficients and properties with high accuracy.

Provided new benchmark values for the Sakiadis boundary layer.

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.

Tables

Table 1. Table 1: Predictions from approximants for the Sakiadis problem. For comparison, the numerical values of κ 𝜅 \kappa and C 𝐶 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 η ∞ subscript 𝜂 \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 𝑆 S were obtained by inserting the numerical κ 𝜅 \kappa and C 𝐶 C values into ( 13 ), locating the singularity of smallest magnitude, and systematically increasing the approximant order (up to N = 200 𝑁 200 N=200 ), where the value of S 𝑆 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

Table 2. Table 2: The infinity norm of the error in the Sakiadis approximants and their 2 nd derivatives (to within 1 significant digit), defined respectively as E 𝐸 E = ‖ f A − f numerical ‖ ∞ subscript norm subscript 𝑓 𝐴 subscript 𝑓 numerical ||f_{A}-f_{\mathrm{numerical}}||_{\infty} and E 2 subscript 𝐸 2 E_{2} = ‖ f A ′′ − f numerical ′′ ‖ ∞ subscript norm superscript subscript 𝑓 𝐴 ′′ superscript subscript 𝑓 numerical ′′ ||f_{A}^{\prime\prime}-f_{\mathrm{numerical}}^{\prime\prime}||_{\infty} on the interval 0 ≤ η ≤ 20 0 𝜂 20 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}$

Table 3. Table 3: Predictions from the Blasius approximant ( 19 ). The infinity norm of the error in the approximant and its 2 nd derivative are recorded below (to within 1 significant digit), defined respectively as E 𝐸 E = ‖ f A − f numerical ‖ ∞ subscript norm subscript 𝑓 𝐴 subscript 𝑓 numerical ||f_{A}-f_{\mathrm{numerical}}||_{\infty} and E 2 subscript 𝐸 2 E_{2} = ‖ f A ′′ − f numerical ′′ ‖ ∞ subscript norm superscript subscript 𝑓 𝐴 ′′ superscript subscript 𝑓 numerical ′′ ||f_{A}^{\prime\prime}-f_{\mathrm{numerical}}^{\prime\prime}||_{\infty} on the interval 0 ≤ η ≤ 8.8 0 𝜂 8.8 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

Table 4. 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 𝐸 E = ‖ u A − u numerical ‖ ∞ subscript norm subscript 𝑢 𝐴 subscript 𝑢 numerical ||u_{A}-u_{\mathrm{numerical}}||_{\infty} on the interval 0 ≤ r ≤ 10 0 𝑟 10 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

Table 5. Table 5: Predictions of u ​ ( 0 ) ≡ z 𝑢 0 𝑧 u(0)\equiv z of ( 48 ), obtained numerically by shooting to the boundary condition u ​ ( 1 − ε ) = 0 𝑢 1 𝜀 0 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 𝑧 2.3919564 z=-2.3919564 and − 2.391956403 2.391956403 -2.391956403 , respectively.

$\epsilon$	$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

Equations

$$f=\sum _{n=0}^{\infty }{a}_n(x-x_0{)}^n,$$

$$f\sim C{f}_a(x)\text{as}x\to x_a,$$

$$P=kT\sum _{n=1}^N{B}_n(T){\rho }^n,B_1=1$$

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

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

$$Z_A={\left(1+\sum _{n=1}^{N-1}{A}_n(h){\tilde{\rho }}^n\right)}^{\frac{h/3}{N-1}}$$

$${\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.$$

$$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.$$

$$\left(\frac{P}{P_c}-1\right)\sim \pm D{\left|\frac{\rho }{{\rho }_c}-1\right|}^{\delta },T=T_c,\rho \to {\rho }_c,$$

$$P_A=P_c-\sum _{n=0}^N{A}_n{\rho }^n{\left(1-\frac{\rho }{{\rho }_c}\right)}^{\delta }.$$

$$\left(P_c-k{T}_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$$

$$A_{n>0}=\frac{P_c\Gamma (\delta +n)}{n!{\rho }_c^n\Gamma (\delta )}-\frac{k{T}_c}{\Gamma (\delta )}\sum _{j=0}^{n-1}\frac{B_{n-j}\Gamma (\delta +j)}{{\rho }_c^{j}j!},A_0=P_c.$$

$$P_c-\frac{k{T}_{c}N!}{\Gamma (\delta +N)}\sum _{j=1}^N\frac{\Gamma (\delta +N-j)}{(N-j)!}{B}_j{\rho }_c^j=0$$

$$f_S=\sum _{n=0}^N{a}_n(x-x_0{)}^n.$$

$$f_A=f_A\left(\sum _{n=0}^N{A}_n[g(x){]}^n\right).$$

$$2f^{\prime \prime \prime }+f{f}^{\prime \prime }=0$$

$$f(0)=0,f^{\prime }(0)=1,f^{\prime }(\infty )=0.$$

$$f=\sum _{n=0}^{\infty }{a}_n{\eta }^n$$

$$a_{n+3}=\frac{-\sum _{j=0}^n(j+1)(j+2)a_{j+2}{a}_{n-j}}{2(n+1)(n+2)(n+3)}.$$

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

$$\underset{\eta \to \infty }{\lim }f\equiv C.$$

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

$$f_A=C-C{\left(1+\sum _{n=1}^N{A}_n{\eta }^n\right)}^{-1},$$

$${\left[1-\frac{1}{C}\sum _{n=0}^N{a}_n{\eta }^n\right]}^{-1}=1+\sum _{n=1}^N{A}_n{\eta }^n.$$

$$A_{n>0}=\frac{1}{C}\sum _{j=1}^n{a}_j{A}_{n-j},A_0=1.$$

$$\sum _{j=1}^N{a}_j{A}_{N-j}=0,\sum _{j=1}^{N-1}{a}_j{A}_{N-1-j}=0.$$

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

$$f_A=C+\sum _{n=1}^N{A}_n{e}^{-nC\eta /2},$$

$$\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],$$

$$A_N(\kappa ,C)=0,A_{N-1}(\kappa ,C)=0$$

$$\kappa =-0.443748313369$$

$$C=1.61612544681$$

$$G=-2.1313459241$$

$$2f^{\prime \prime \prime }+f{f}^{\prime \prime }=0,$$

Full text

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)$:

[TABLE]

and an asymptotic behavior

[TABLE]

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:

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

[TABLE]

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

[TABLE]

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:

[TABLE]

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

[TABLE]

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.

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]

[TABLE]

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

[TABLE]

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.

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:

[TABLE]

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:

[TABLE]

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)$:

[TABLE]

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

[TABLE]

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

[TABLE]

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$:

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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$:

[TABLE]

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

[TABLE]

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).

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$

[TABLE]

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.

4.1.2 Asymptotic approximant with exponential corrections

We now consider an approximant,

[TABLE]

[TABLE]

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

[TABLE]

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

[TABLE]

and

[TABLE]

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

[TABLE]

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$.

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].

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),

[TABLE]

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

[TABLE]

which, along with the asymptotic properties [43]

[TABLE]

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)

[TABLE]

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

[TABLE]

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$:

[TABLE]

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

[TABLE]

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.

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

[TABLE]

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$.

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

[TABLE]

The coefficients of the power series solution

[TABLE]

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,

[TABLE]

which, along with the asymptotic property

[TABLE]

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):

[TABLE]

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

[TABLE]

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.

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

[TABLE]

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

[TABLE]

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.

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.

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]:

[TABLE]

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:

[TABLE]

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]:

[TABLE]

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:

[TABLE]

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

[TABLE]

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.:

[TABLE]

The solution of this equation is

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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:

[TABLE]

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

[TABLE]

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.:

[TABLE]

The solution of this equation is

[TABLE]

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

[TABLE]

and after application of integration by parts, we obtain:

[TABLE]

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

[TABLE]

from which integration by parts yields

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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.

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}$.