# Exponential Asymptotics using Numerical Rational Approximation in Linear
Differential Equations

Christopher J. Lustri111Electronic address<EMAIL_ADDRESS>School of Mathematics and Statistics, The University of Sydney, New South
Wales, 2006, Australia  Theoretical Sciences Visiting Program, Okinawa
Institute of Science and Technology Graduate University, Onna, 904-0495,
Japan.  Samuel Crew  Faculty of Computer Science, Ruhr Universität Bochum,
Universitätsstraße, Bochum, 44799, Germany  Theoretical Sciences Visiting
Program, Okinawa Institute of Science and Technology Graduate University,
Onna, 904-0495, Japan.  S. Jonathan Chapman  Mathematical Institute,
University of Oxford, Oxford, OX2 6GG, United Kingdom

###### Abstract

Singularly-perturbed ordinary differential equations often exhibit Stokes’
phenomenon, which describes the appearance and disappearance of oscillating
exponentially small terms across curves in the complex plane known as Stokes
curves. These curves originate at singular points in the leading-order
solution to the differential equation. In many important problems, it is
impossible to obtain a closed-form expression for these leading-order
solutions, and it is therefore challenging to locate these singular points. We
present evidence that the analytic leading-order solution of a linear
differential equation can be replaced with a rational approximation based on a
numerical leading-order solution using the adaptive Antoulas-Anderson (AAA)
method. We show that the subsequent exponential asymptotic analysis accurately
predicts the exponentially small behaviour present in the solution. We explore
the limitations of this approach, and show that for sufficiently small values
of the asymptotic parameter, this approach breaks down; however, the range of
validity may be extended by increasing the number of poles in the rational
approximation. We finish by presenting a related nonlinear problem and
discussing the challenges that arise when attempting to apply this method to
nonlinear problems.

##  1 Introduction

The behaviour known as “Stokes’ phenomenon” was first observed by George
Gabriel Stokes in his analysis of solutions to the Airy equation [25, 26]. The
Airy function displays oscillatory behaviour for large negative arguments, and
decays exponentially for large positive arguments. Stokes wished to understand
the transition between these two regimes by studying the large-argument
asymptotic behaviour in the complex plane. In doing so, Stokes discovered that
exponentially-small correction terms in the asymptotic expansion appear or
disappear as special curves in the complex plane are crossed. These curves are
known as ‘Stokes curves’.

New asymptotic techniques were developed in order to study the switch in
behaviour across Stokes curves. Berry [2, 3, 4], and Berry & Howls [5, 6]
showed that truncating a a divergent asymptotic series optimally results in
asymptotic error of exponentially small order. Hence, rescaling the problem to
study this truncation error directly allows for the calculation of these
exponentially small terms. This idea was used to develop new asymptotic
methods for studying exponentially small behaviour, known as exponential
asymptotics, asymptotics-beyond-all-orders or hyperasymptotics [6, 21]. We
will denote exponentially small contributions to the asymptotic approximation
of some solution $u$ as $u_{\mathrm{exp}}$.

Exponential asymptotic techniques typically require the explicit calculation
of the leading-order solution of the differential equation, which we will
denote as $u_{0}$. Stokes curves originate at singularities of $u_{0}$ [17].
Additionally, the asymptotic size of $u_{\mathrm{exp}}$ as $\epsilon\to 0$
typically depends on the location and strength of the singularity. For many
problems, knowing the analytical form of the leading-order solution near
singularities provides enough information to calculate the Stokes curves and
approximate the exponential contributions [7].

For many practical problems, it is impossible to calculate the leading-order
solution analytically [14, 15, 16, 9]. In this case, the leading-order
solution must instead be approximated numerically. There is growing evidence
that numerical methods for analytic continuation, including numerically
stepping into the complex plane [9] and rational approximation methods [16]
can be used in conjunction with exponential asymptotics.

In general, the leading-order solution to an ordinary differential equation
can be computed on a real domain using any one of a variety of standard
numerical methods [19]. Using the result of numerical computations to describe
the complex plane behaviour of these solutions is challenging. It is well-
known that analytic continuation is ill-posed, and that small changes to the
function being continued, such as numerical error associated with
computational methods, can lead to significant changes in the analytically-
continued result.

Despite this, significant progress has been made on performing analytic
continuation using numerical methods; this is typically performed by requiring
the output function to take a particular algebraic form. The most commonly-
used methods for numerical complex analytic continuation are methods that give
a rational function as an output. These rational approximation methods include
Padé approximation [1] and the adaptive Antoulas–Anderson (AAA) algorithm
[20].

We wish to consider settings in which the solution has been numerically
approximated on some discrete set of points along the real axis which we then
seek to analytically continue into the complex plane. This setup makes the AAA
algorithm particularly well-suited for our purpose. This algorithm is
described in Section 1.1, and produces a rational function, denoted here as
$\hat{u}_{0}$ which approximates the leading-order behaviour. Note that all
singularities in $\hat{u}_{0}$ must be simple poles; if $u_{0}$ contains other
singularities, such as branch points, the branch will be typically be
approximated as an accumulation of simple poles instead [18, 27, 29]. Recall
that the strength of the singularity in $u_{0}$ is generically different to
the simple poles present in $\hat{u}_{0}$, and the asymptotic behaviour of
$u_{\mathrm{exp}}$ depends on the strength of the singularity. This argument
indicates that any exponential term obtained using this leading-order
approximation must display different asymptotic behaviour to
${u}_{\mathrm{exp}}$ in the limit that $\epsilon\to 0$.

In this study, we focus our attention on a singularly-perturbed linear
ordinary differential equation in the limit that $\epsilon\to 0$, whose
leading-order solution $u_{0}$ contains two branch points, namely

$\epsilon^{2}u^{\prime\prime}(x)+u(x)=\frac{1}{\sqrt{x+\mathrm{i}}}+\frac{1}{\sqrt{x-\mathrm{i}}},$
(1)

with the boundary conditions

$\lim_{x\to-\infty}u(x)=0,\qquad\lim_{x\to-\infty}u^{\prime}(x)=0.$ (2)

We note that differential equation (1) possesses a closed-form solution in
terms of the exponential integral function,

$\displaystyle u(x)=\,$ $\displaystyle
A(\epsilon)\cos\left(\frac{x}{\epsilon}\right)+B(\epsilon)\sin\left(\frac{x}{\epsilon}\right)$
$\displaystyle-\frac{\mathrm{i}\sqrt{x+\mathrm{i}}}{2\epsilon}\mathrm{e}^{(1-\mathrm{i}x)/\epsilon}\int_{1}^{\infty}\frac{\mathrm{e}^{-(1-\mathrm{i}x)t/\epsilon}}{\sqrt{t}}\,\mathrm{d}t+\frac{\mathrm{i}\sqrt{x+\mathrm{i}}}{2\epsilon}\mathrm{e}^{-(1-\mathrm{i}x)/\epsilon}\int_{1}^{\infty}\frac{\mathrm{e}^{(1-\mathrm{i}x)t/\epsilon}}{\sqrt{t}}\,\mathrm{d}t+\mathrm{c.c.},$
(3)

where $A$ and $B$ may be determined from the boundary conditions. It is
possible to determine $A$ and $B$ by performing a careful analysis on the
integral expressions in the limit that $x\to-\infty$, but we will instead take
a different approach, and determine an asymptotic expansion for the solution
in the limit that $\epsilon\to 0$.

We will perform a typical exponential asymptotic analysis on (1)–(2), using
the method proposed in [22], to determine ${u}_{\mathrm{exp}}$ for this
problem. We will then sample $u_{0}$ on a discrete set of points and use this
as the input for a rational approximation $\hat{u}_{0}$. We will apply the
same exponential asymptotic method on this leading-order expression to obtain
$\hat{u}_{\mathrm{exp}}$. By comparing the analytic form of
${u}_{\mathrm{exp}}$ and $\hat{u}_{\mathrm{exp}}$, we will see that the two
expressions have different asymptotic behaviour in the limit $\epsilon\to 0$.

We will then show that, despite this difference, $\hat{u}_{\mathrm{exp}}$ is
able to accurately approximate $u_{\mathrm{exp}}$ for a range of values of
$\epsilon$, and that this range can be extended by increasing the accuracy of
the AAA approximation (ie. reducing the $L^{2}$ error threshold on the set of
sample points). Finally, we will provide evidence that the threshold value of
$\epsilon$ below which $\hat{u}_{\mathrm{exp}}$ is inaccurate is proportional
to the difference between the true branch point location in $u_{0}$ and the
pole nearest to this point in $\hat{u}_{0}$. This may be thought of as the
approximation error of the true branch point location from the rational
approximation algorithm.

###  1.1 Numerical Analytic Continuation

Rational approximation methods typically approximate some function $f(x)$ as a
rational expression,

$\displaystyle f(x)\approx\frac{n(x)}{d(x)},$ (4)

where $n(x)$ and $d(x)$ are polynomials.

Historically, the most widely-used method for obtaining rational
approximations is Padé approximation, which is constructed using the Taylor
coefficients of $f(x)$ at a particular point. We want to study data that takes
a different form; where values of $f(x)$ are calculated numerically on some
discrete support set of points. While it is possible to approximate Taylor
coefficients using values of the function on a discrete set of points, there
are other methods more well-suited to input data with this form.

On such is the AAA approximation, which works by iteratively solving a
minimisation problem to produce an optimal rational approximation. The
explanation that follows is a high-level description of the AAA algorithm; for
a detailed explanation, see [20].

The algorithm is based on expressing the rational approximation in the form

$\displaystyle
f(x)\approx\frac{n(x)}{d(x)}={\sum_{j=1}^{m}\frac{w_{j}f_{j}}{x-x_{j}}}\,\Bigg{/}\,{\sum_{j=1}^{m}\frac{w_{j}}{x-x_{j}}},$
(5)

The points $x_{j}$ for $j=1,2,\ldots,m$, are known as support points. They are
drawn from the sample set $X$, which consists of the points on which $f(x)$ is
known; we define $f_{j}=f(x_{j})$. The algorithm solves an optimization
problem in order to determine the weights $w_{j}$, and then generates another
support point $x_{j+1}$. This process then repeats iteratively.

The weights $w_{j}$ are generated by solving a linear least-squares problem
over a set of points in the restricted domain
$X^{(m)}=X\setminus\\{x_{1},\ldots,x_{m}\\}$; these points are labelled
$X_{i}^{(m)}$. The weight vector $w=(w_{1},w_{2},\dots,w_{m})^{T}$ is chosen
to minimize $||fd-n||_{X^{(m)}}$ subject to $||w||_{m}=1$, where
$||\cdot||_{X^{(m)}}$ is the discrete 2-norm over $X^{(m)}$ and
$||\cdot||_{m}$. is the discrete $L^{2}$ on $m$-vectors. If the 2-norm of the
approximation residuals on $X^{(m)}$ are beneath some specified tolerance
(relative to the maximum value taken by $|f(x_{j})|$), the algorithm
terminates.

If the algorithm does not terminate, the value of $x_{m+1}$ is determined for
the next iteration. It is found by choosing the value of $x_{m+1}\in X^{(m)}$
that maximizes the quantity

$\displaystyle\sum_{j=1}^{m}\frac{w_{j}f(x_{m+1})}{x_{m+1}-x_{j}}-\sum_{j=1}^{m}\frac{w_{j}f_{j}}{x_{m+1}-x_{j}}.$
(6)

The process then repeats until the termination criteria is reached.

The key information for our purposes is that the method takes in the function
values on a set of points, and returns a rational approximation that minimizes
the $L^{2}$-norm of the approximation error on the data support set, and that
it is an iterative process that increases the number of poles in the solution
until some predetermined tolerance is met. We often write the approximation in
the form

$f(x)\approx\sum_{j=0}^{m}\frac{a_{j}}{x-p_{j}},$ (7)

where $a_{j}$ are the residues of the function at each pole location $p_{j}$.

As previously alluded to, one potential obstacle in applying rational
approximation methods is that these methods produce meromorphic functions; all
singularities in the approximation are isolated simple poles. Previous studies
of AAA approximation [18, 27, 29] showed that the rational function
approximation of a target function with a branch point contains an exponential
clustering of poles (and zeroes) in the approximation approaching the branch
point. This configuration of poles accurately approximates the effects of the
branch cut. The distribution of poles when approximating branch cuts in
rational approximations has been studied rigourously in the related problem of
Padé approximation [24].

In the following we show that if we approximate the leading-order solution of
an ordinary differential equation with a rational function of the form (7), it
is possible to use exponential asymptotics on this rational function in order
to approximate the exponentially small terms that appear in the solution due
to Stokes’ phenomenon.

##  2 Exponential Asymptotic Analysis

###  2.1 Asymptotic Power Series

We propose an asymptotic series solution to (1)–(2) of the form

$u(x;\epsilon)\sim\sum_{n=0}^{\infty}\epsilon^{2n}u_{n}(x).$ (8)

By substituting this into the ordinary differential equation (1) and matching
powers of $\epsilon$ in the limit that $\epsilon\to 0$, we find that the
leading-order behaviour is given by

$\displaystyle
u_{0}=\frac{1}{\sqrt{x+\mathrm{i}}}+\frac{1}{\sqrt{x-\mathrm{i}}}.$ (9)

Matching at $\mathcal{O}(\epsilon^{n})$ as $\epsilon\to 0$ gives

$\displaystyle u_{n}=-u_{n-1}^{\prime\prime},\qquad n\geq 1.$ (10)

This recurrence relation can be used calculate $u_{n}$ exactly, giving

$u_{n}=\frac{(4n)!(-1)^{n}}{2^{4n}(2n)!(x-\mathrm{i})^{2n+{1/2}}}+\mathrm{c.c.},$
(11)

The series (8) with terms given by (11) is divergent for any fixed choice of
$x$. Figure 1 depicts $\epsilon^{n}u_{n}$ at $x=0$ with $\epsilon=0.1$. The
magnitude of the terms decreases until $n=5$, after which the size of
successive terms increases, causing the series to diverge.

$-5$$3$$-4$$-3$$-2$$-1$$0$$n$$1$$2$$\log_{10}|\epsilon^{2n}u_{n}|$ at
$x=0$$0$$5$$10$$15$Minimum value Figure 1:  Magnitude of series terms from (8)
evaluated at $x=0$ for $\epsilon=0.1$. As $n$ increases, the terms become
smaller until a minimum value is reached at $n=5$, after which the terms
increase in size due to the factorial contribution to the numerator of $u_{n}$
in (11) The series (8) must therefore be divergent, and the optimal truncation
point occurs at the minimum value.

Typically, truncating the series at the value of $n$ that minimises
$\epsilon^{2n}|u_{n}|$ produces the most accurate approximation that can be
achieved by the series. This value of $n$ is often known as the “optimal
truncation point”, which we denote as $N_{\mathrm{opt}}$.

Divergence of an asymptotic series, such as that seen in Figure 1, indicates
that the solution contains exponentially small terms that cannot be described
by the algebraic series. Behaviour that occurs on this scale is smaller than
any algebraic term in the limit that $\epsilon\to 0$, and hence cannot be
captured by the series expression.

Exponential asymptotic methods are built on the observation that truncating a
series optimally at $n=N_{\mathrm{opt}}$ produces a remainder that is
exponentially small in the asymptotic limit [3]. By rescaling to study the
exact truncation remainder, we can directly calculate exponentially small
components of the asymptotic expansion.

###  2.2 Finding the Stokes Curves

The analysis presented here is a standard application of the exponential
asymptotic method developed in [22] in order to study Stokes’ phenomenon.
Nonetheless, this section will provide a relatively complete summary of the
procedure in order to make the analysis accessible.

The purpose of this analysis is to obtain the leading exponentially small
correction term to the series (8), and to therefore reveal the effects of
Stokes switching in the solution. We will concentrate on the contribution
arising from the term explicitly shown in (11), and note that a complex
conjugate contribution is also present in the solution, which we will include
in the final expression.

We now truncate the power series (8) after $N$ terms, giving

$u(x;\epsilon)=\sum_{n=0}^{N-1}\epsilon^{2n}u_{n}(x)+R_{N}(x),$ (12)

where $R_{N}$ is the remainder. The optimal truncation point
$N_{\mathrm{opt}}$ corresponds to minimising $\epsilon^{2N}|u_{N}|$. The most
straightforward way to find $N_{\mathrm{opt}}$ is to determine the value at
which consecutive terms have equal magnitude, or

$\left|\frac{(4n)!(-1)^{n}\epsilon^{2n}}{2^{4n}(2n)!(x-\mathrm{i})^{2n+{1/2}}}\right|=\left|\frac{(4n+4)!(-1)^{n+1}\epsilon^{2n+2}}{2^{4n+4}(2n+2)!(x-\mathrm{i})^{2n+{5/2}}}\right|,$
(13)

Optimal truncation occurs after a large number of terms in the limit that
$\epsilon\to 0$. We can make use of this to find that the optimal truncation
point must satisfy $N_{\mathrm{opt}}\sim|x-\mathrm{i}|/2\epsilon$ as
$\epsilon\to 0$. We therefore set

$N_{\mathrm{opt}}=\frac{|x-\mathrm{i}|}{2\epsilon}+\omega,$ (14)

where $\omega\in[0,1)$ is chosen such that $N_{\mathrm{opt}}$ is an integer.
This expression for $N_{\mathrm{opt}}$ is consistent with Figure 1 which has
$|x-\mathrm{i}|=1$ and $\epsilon=0.1$, corresponding to $N_{\mathrm{opt}}=5$.
We will eventually use this value for $N$ in (12), but for algebraic
simplicity we will not make the substitution immediately.

We substitute the truncated series (12) into the governing equation (1) to
obtain

$\sum_{n=0}^{N-1}\epsilon^{2n+2}u_{n}^{\prime\prime}+\sum_{n=0}^{N-1}\epsilon^{2n}u_{n}+\epsilon^{2}R_{N}^{\prime\prime}+R_{N}=\frac{1}{\sqrt{x+\mathrm{i}}}+\frac{1}{\sqrt{x-\mathrm{i}}}.$
(15)

Using (11) and (10) to simplify gives

$\epsilon^{2}R_{N}^{\prime\prime}+R_{N}=\epsilon^{2N}u_{N}=\frac{(4N)!(-1)^{N}\epsilon^{2N}}{2^{4N}(2N)!(x-\mathrm{i})^{2N+{1/2}}},$
(16)

We will solve this differential equation using variation of parameters. The
first step is to determine the homogeneous solutions to (16), given by

$R_{N}(x)=K_{1}\mathrm{e}^{\mathrm{i}x/\epsilon}\qquad\mathrm{and}\qquad
R_{N}(x)=K_{2}\mathrm{e}^{-\mathrm{i}x/\epsilon},$ (17)

where $K_{1}$ and $K_{2}$ are arbitrary constants. The next step is to permit
the arbitrary constants to vary in $x$, and then consider the full
differential equation in (16). We will find that for (16) the second solution
is the relevant one, so we set

$R_{N}(x)=K_{2}(x)\mathrm{e}^{-\mathrm{i}x/\epsilon}=\mathcal{S}(x)\mathrm{e}^{-\mathrm{i}(x-\mathrm{i})/\epsilon},$
(18)

where $K_{2}(x)=\mathcal{S}(x)\mathrm{e}^{-1/\epsilon}$, with the scaling
chosen so that the exponent in (18) is equal to 0 when $x=\mathrm{i}$. When
written in this form, $\mathcal{S}(x)$ is known as the Stokes multiplier. It
is constant except within a narrow region surrounding a Stokes curve. In this
region, $\mathcal{S}$ varies rapidly from zero on one side of the curve to a
nonzero value on the other. This rapid variation generates Stokes switching in
the solution.

Substituting (18) into (16) and simplifying gives

$-2\mathrm{i}\epsilon\mathcal{S}^{\prime}\mathrm{e}^{-\mathrm{i}(x-\mathrm{i})/\epsilon}=\frac{(4N)!(-1)^{N}\epsilon^{2N}}{2^{4N}(2N)!(x-\mathrm{i})^{2N+{1/2}}}.$
(19)

Recall that the optimal truncation given in (14) depends on $|x-\mathrm{i}|$.
This suggests that the transformation
$\mathrm{i}(x-\mathrm{i})=r\mathrm{e}^{\mathrm{i}\theta}$ will simplify this
expression in a useful way, giving $N_{\mathrm{opt}}=r/2\epsilon+\omega$. As
in [22], we then fix $r$ and consider only the faster variation that occurs in
the angular direction $\theta$. From the chain rule,

$\frac{\mathrm{d}}{\mathrm{d}x}=\frac{\mathrm{e}^{-\mathrm{i}\theta}}{r}\frac{\mathrm{d}}{\mathrm{d}\theta}.$
(20)

Applying this transformation to (16) and rearranging gives

$\frac{\mathrm{d}\mathcal{S}}{\mathrm{d}\theta}=-\frac{\mathrm{i}\sqrt{\mathrm{i}}r\mathrm{e}^{\mathrm{i}\theta}(4N)!\epsilon^{2N-1}}{2^{4N+1}(2N)!(r\mathrm{e}^{\mathrm{i}\theta})^{2N+1/2}}\exp\left(\frac{r\mathrm{e}^{\mathrm{i}\theta}}{\epsilon}\right)\qquad\mathrm{as}\qquad\epsilon\to
0.$ (21)

We now finally substitute in the optimal truncation value
$N_{\mathrm{opt}}=r/2\epsilon+\omega$. Making use of Stirling’s formula in the
limit that $\epsilon\to 0$ gives

$\frac{\mathrm{d}\mathcal{S}}{\mathrm{d}\theta}\sim-\frac{\mathrm{i}\sqrt{2\mathrm{i}r}}{\epsilon}\exp\left(\frac{r}{\epsilon}\left(\mathrm{e}^{\mathrm{i}\theta}-1-\mathrm{i}\theta\right)+\mathrm{i}\theta\left(\frac{1}{2}+\omega\right)\right)\qquad\mathrm{as}\qquad\epsilon\to
0.$ (22)

As promised, this expression is exponentially small, except in the
neighbourhood of $\theta=0$, where the exponential term is algebraic. Hence,
$\mathcal{S}$ varies rapidly in this region, indicating that $\theta=0$ is a
Stokes curve.

Recall that $\mathrm{i}(x-\mathrm{i})=r\mathrm{e}^{\mathrm{i}\theta}$, so
$\theta=0$ corresponds to $x-\mathrm{i}$ taking negative imaginary values, or
a vertical line extending downwards from the branch point $x=\mathrm{i}$. If
we perform a corresponding analysis for the complex conjugate term in (11), we
identify that there is a second Stokes curve extending vertically upwards from
the branch point $x=-\mathrm{i}$. This Stokes structure is illustrated in
Figure 2.

$\mathrm{Re}(x)$$\mathrm{Im}(x)$i$-$i$u_{\mathrm{exp}}=0$$u_{\mathrm{exp}}\sim\sqrt{\frac{4\pi}{\epsilon}}\mathrm{e}^{-1/\epsilon}\cos\left(\frac{x}{\epsilon}+\frac{\pi}{4}\right)$Stokes
curvesBranch cutBranch pointNo exp.termsExp. smallcorrections Figure 2:
Stokes’ phenomenon in the exact solution to (1) with boundary conditions (2).
The leading-order solution contains branch points at $x=\pm\mathrm{i}$, with
the branch cuts extending vertically away from the real axis. The branch
points generate Stokes curves, which connect the two points. On the left-hand
side of the Stokes curves, there are no exponential contributions. On the
right-hand side of the Stokes curves, the solution contains exponentially
small oscillations of the form given in (27).

If we cross the Stokes line along $\mathrm{Re}(x)=0$, we expect that there
will be an exponentially small jump in the asymptotic solution.

###  2.3 Exponentially Small Contribution

To determine the behaviour that appears as the Stokes curve is crossed, we
define a new inner variable $\theta=\epsilon^{1/2}\vartheta$, which gives

$\frac{\mathrm{d}\mathcal{S}}{\mathrm{d}\vartheta}\sim-\mathrm{i}\sqrt{\frac{2\mathrm{i}r}{\epsilon}}\mathrm{e}^{-r\vartheta^{2}/2}\qquad\mathrm{as}\qquad\epsilon\to
0.$ (23)

By comparing $\vartheta$ with the original coordinates, we see that $x<0$
corresponds to $\vartheta\to-\infty$, while $x>0$ corresponds to
$\vartheta\to\infty$. Hence, the jump across the Stokes curve is given by

$[\mathcal{S}]_{-}^{+}=\lim_{\vartheta\to+\infty}\mathcal{S}-\lim_{\vartheta\to-\infty}\mathcal{S}\sim-\mathrm{i}\sqrt{\frac{2\mathrm{i}r}{\epsilon}}\int_{-\infty}^{\infty}\mathrm{e}^{-r\vartheta^{2}/2}\,\mathrm{d}\vartheta=-2\mathrm{i}\sqrt{\frac{\pi\mathrm{i}}{\epsilon}}.$
(24)

If the value of $\mathcal{S}$ changes by (24) as the Stokes curve is crossed,
the exponentially small contribution to the solution for $x>0$ is given by

$[R_{N}]_{-}^{+}\sim-2\mathrm{i}\sqrt{\frac{\pi\mathrm{i}}{\epsilon}}\mathrm{e}^{-\mathrm{i}(x-\mathrm{i})/\epsilon}.$
(25)

If we add in the complex conjugate term, we obtain the jump in the
exponentially small terms $u_{\mathrm{exp}}$ as the Stokes curve is crossed
from left to right,

$[u_{\mathrm{exp}}]_{-}^{+}\sim-2\mathrm{i}\sqrt{\frac{\pi\mathrm{i}}{\epsilon}}\,\mathrm{e}^{-\mathrm{i}(x-\mathrm{i})/\epsilon}+\mathrm{c.c.}=4\sqrt{\frac{\pi}{\epsilon}}\,\mathrm{e}^{-1/\epsilon}\cos\left(\frac{x}{\epsilon}+\frac{\pi}{4}\right),\qquad\mathrm{as}\qquad\epsilon\to
0.$ (26)

Note that this expression has constant amplitude. By comparing this with the
boundary conditions in (2), we see that there cannot be any finite-amplitude
oscillations present as $x\to-\infty$. This indicates that
$u_{\mathrm{exp}}=0$ on the left-hand side of the Stokes curve ($x<0$), and
that the exponentially small behaviour appears as the Stokes curve is crossed
from left to right. On the right-hand side of the Stokes curve, the
exponentially small solution contribution $u_{\mathrm{exp}}$ must be

$u_{\mathrm{exp}}\sim
4\sqrt{\frac{\pi}{\epsilon}}\,\mathrm{e}^{-1/\epsilon}\cos\left(\frac{x}{\epsilon}+\frac{\pi}{4}\right),\qquad
x>0.$ (27)

This behaviour is illustrated in Figure 2.

##  3 AAA Approximation

The analysis in Section 2 was a relatively straightforward application of the
method of [22]. It centered on the observation that singularities of the
leading-order solution in the complex plane, such as the branch points at
$x=\pm\mathrm{i}$ in (9), generate Stokes curves. We can then optimally
truncate and rescale the equation to study the exponentially small
contributions in the neighbourhood of these terms.

###  3.1 AAA Approximation for $u_{0}$

In order to simulate a problem where we only possess a numerical description
of the leading-order solution $u_{0}$ from (9), we first define a set of
points $x_{j}$ separated by some $\Delta x$. We then evaluate (9) at each
point to obtain $u_{0}(x_{j})$. We will use this set of points $x_{j}$ and
function values $u_{0}(x_{j})$ as the basis for a AAA rational approximation.

We apply the AAA algorithm to obtain a rational approximation
$\hat{u}_{0}(x)$, given by

$u_{0}\approx\hat{u}_{0}=\sum_{r=0}^{m}\frac{a_{r}}{x-p_{r}}.$ (28)

For example, if we choose $x_{j}$ to be evenly distributed in the interval
$[-4,4]$ with $\Delta x=0.1$, applying the AAA algorithm with an error
tolerance of $10^{-12}$ gives a rational function with 15 poles. The poles and
residues are are shown in Table 1. We have given each pair of poles a
designation, so that we may refer to them later.

Pole Index: $r$ | Pole Pair Name | Pole locations: $p_{r}$ | Pole residues: $a_{r}$
---|---|---|---
1, 2 | Pair 1 | $-0.0015\pm 1.0178\mathrm{i}$ | $0.1256\mp 0.1155\mathrm{i}$
3, 4 | Pair 2 | $-0.0142\pm 1.1647\mathrm{i}$ | $0.1319\mp 0.1204\mathrm{i}$
5, 6 | Pair 3 | $-0.0444\pm 1.4872\mathrm{i}$ | $0.1458\mp 0.1312\mathrm{i}$
7, 8 | Pair 4 | $-0.1052\pm 2.0575\mathrm{i}$ | $0.1725\mp 0.1513\mathrm{i}$
9, 10 | Pair 5 | $-0.2369\pm 3.0523\mathrm{i}$ | $0.2283\mp 0.1898\mathrm{i}$
11, 12 | Pair 6 | $-0.6072\pm 5.0113\mathrm{i}$ | $0.3780\mp 0.2791\mathrm{i}$
13, 14 | Pair 7 | $-2.5372\pm 10.7492\mathrm{i}$ | $1.1038\mp 0.6053\mathrm{i}$
15 | Unpaired | $-6.6066+0.0000\mathrm{i}$ | $-0.0003+0.0000\mathrm{i}$
Table 1: Poles and residues in the AAA approximation of $u_{0}$, which
contains branch points at $x=\pm\mathrm{i}$. The approximation $\hat{u}_{0}$
was produced by sampling $u_{0}$ on the domain $x\in[-4,4]$ at intervals of
$\Delta x=0.1$ with an error tolerance of $10^{-12}$. The poles in
$\hat{u}_{0}$ accumulate at $x=\pm\mathrm{i}$, or the true branch points of
$u_{0}$. An unpaired pole appears on the real axis outside the sampling
interval, which we discard as a numerical artefact. The remaining poles occur
in complex conjugate pairs. The pole $p_{1}$ is the nearest pole to the true
branch point at $x=\mathrm{i}$, and this pole will play an important role in
subsequent analysis.

Note that there is a pole located on the real axis at $x\approx-6.066$. The
AAA algorithm sometimes generates poles on the real axis that lie outside of
the approximation interval. These poles have no practical effect on the
solution, and we will omit their contributions from the sum (28). The
remaining poles appear as complex conjugate pairs, as the solution is real
when $x$ is real.

The pole pairs are not restricted to the imaginary axis, which is where we
defined the branch cuts to be in the true solution $u_{0}$ (9). In fact, the
AAA algorithm places the poles along a curve, representing the branch cut,
which we are not free to prescribe. This is a common feature of rational
approximation methods, and was studied in depth for Padé approximation in
[24], using a quantity known as the “condenser capacity” in the approximation.
More recently, an explanation of the branch curve selection for AAA that built
on these ideas was presented in [28].

It is possible to adjust the rational approximation algorithm, for example by
using the minimax algorithm from the chebfun package, to obtain an
approximation which forces the poles to align along the imaginary axis.

In Figure 3, we present the real and imaginary parts of the true leading-order
solution $u_{0}$, and the approximated leading-order solution $\hat{u}_{0}$.
Note that the true solution possesses vertical branch cuts originating at
$\pm\mathrm{i}$, while the approximated solution possesses simple poles that
simulate the effect of a branch cut. The solutions appear visually identical
except in a region surrounding the branch cut/poles. In this, and all
subsequent analysis, we denote $p_{1}$ as the pole nearest to $x=\mathrm{i}$.

(a) Re($u_{0}$)

(b) Im($u_{0}$)

(c) Re($\hat{u}_{0}$)

(d) Im($\hat{u}_{0}$)

Figure 3: The real and imaginary parts of the true leading-order solution
$u_{0}$ and the approximated leading-order solution $\hat{u}_{0}$ described in
Table 1. The two expressions are visually indistinguishable except on the
imaginary axis. The function $u_{0}$ contains vertical branch cuts. The
function $\hat{u}_{0}$ is a rational function which can only contain simple
poles; these poles are arranged in such a way that they approximate the effect
of a branch cut in the solution.

In Figure 4, we present the error between the two on a logarithmic plot. This
confirms that the approximation error is small except in a region surrounding
the branch cut/poles, where the error becomes significant. Given that the
approximation is highly accurate everywhere on the sample domain except near
the branch cut, we hope that the exponential asymptotic analysis will produce
equivalently accurate results.

(a) $\log_{10}|u_{0}-\hat{u}_{0}|$ Figure 4: This figure shows the error
$|u_{0}-\hat{u}_{0}|$, using $u_{0}$ and $\hat{u}_{0}$ from Figure 3. Note
that the error is extremely small except in a region near the imaginary axis,
where the true branch cut lies.

###  3.2 Asymptotic Power Series

In effect, we are using the AAA approximation in place of the inhomogeneous
term in (1). Hence, we are determining the asymptotic behaviour of

$\epsilon^{2}\hat{u}^{\prime\prime}(x)+\hat{u}(x)=\sum_{r=0}^{m}\frac{a_{r}}{x-p_{r}}.$
(29)

By expanding the solution as a series (8) and matching powers of $\epsilon$,
we can again obtain a recurrence relation for the series terms $u_{n}$. The
recurrence relation gives

$\displaystyle\hat{u}_{0}(x)=\sum_{r=0}^{m}\frac{a_{r}}{x-p_{r}},\qquad\hat{u}_{n}(x)=-\hat{u}_{n-1}^{\prime\prime}(x).$
(30)

Solving this recurrence relation gives

$\hat{u}_{n}=\sum_{r=0}^{m}\frac{a_{r}(2n)!}{(x-p_{r})^{2n+1}}.$ (31)

Each of the poles $p_{r}$ will generate a Stokes curve, and lead to an
exponentially small contribution appearing in the solution.

###  3.3 Stokes Curves and Exponentially Small Terms

To determine the location of the Stokes curves, and the quantity that appears
as they are crossed, we proceed with an almost identical analysis to that of
Sections 2.2 and 2.3, with the only significant difference being the form of
the series terms in (31).

This analysis reveals that there are exponentially small contributions in the
solution associated with each pole. These contributions appear across Stokes
curves that extends vertically from the corresponding pole and intersect the
real axis. This behaviour is shown for the example problem in Figure 5. We
will later find that the most significant exponential contributions in this
example appear across the Stokes curves generated by pole pairs 1, 2, and 3.

$\mathrm{Re}(x)$$\mathrm{Im}(x)$Pair 7Pair 6Pair
5$2.5$Inset$-0.2$$0.5$$0.1$See Inset$1$$0$Pair 1Pair 2Pair 3Pair
4$-1$$-2$$-3$$5$$10$$-5$$-10$Stokes curves Figure 5: Stokes structure in the
solution to (29), which uses $\hat{u}_{0}$ as the leading-order solution. The
solution contains 7 pairs of poles, with locations given in Table 1. Each pair
of poles generates Stokes curves that extend vertically from the poles,
intersecting the real axis. As each Stokes curve is crossed from left to
right, an exponentially small asymptotic contribution appears in the solution.
Note that the poles accumulate near the true branch points of $u_{0}$ at
$x=\pm\mathrm{i}$. We will later find that the largest exponential
contributions arise from the poles that are nearest to $x=\pm\mathrm{i}$.

Using an essentially identical set of steps to Section 2.3, we can determine
the form of the exponentially small remainder. We will determine the
contribution that appears across the Stokes curve generated by a pole in the
upper-half plane located at at $x=p_{r}$. The Stokes curve extends vertically
downwards from the pole along $\mathrm{Re}(x)=\mathrm{Re}(p_{r})$ and
intersects the real axis at $x=\mathrm{Re}(p_{r})$.

We denote the optimally-truncated remainder as $\hat{R}_{N}$. On the left-hand
side of the Stokes curve, we have $\hat{R}_{N}=0$. On the right-hand side of
the Stokes curve, we find that

$\hat{R}_{N}\sim\frac{2\pi
a_{r}}{\epsilon}\,\mathrm{e}^{-\mathrm{i}(x-p_{r})/\epsilon}\quad\mathrm{as}\quad\epsilon\to
0.$ (32)

Once all of the Stokes curves have been crossed from left to right, we find
that the combined exponentially small contribution, which we denote as
$\hat{u}_{\mathrm{exp}}$, is therefore given by

$\hat{u}_{\mathrm{exp}}\sim\sum_{r=0}^{m}\frac{2\pi
a_{r}}{\epsilon}\,\mathrm{e}^{-\mathrm{i}(x-p_{r})/\epsilon}\quad\mathrm{as}\quad\epsilon\to
0.$ (33)

For example, for the set of poles presented in Table 1, we find that all of
the exponentially small contributions are present in the solution on the
right-hand side of the Stokes curve that intersects the real line at
$x=-0.0015$.

Note that the exponential contributions in (33) associated with upper-half
plane poles at $x=p_{r}$ decay exponentially as $\mathrm{Im}(p_{r})$ increases
(with the complex conjugate contributions exhibiting corresponding exponential
decay). This suggests that the poles nearest to the real axis will produce the
largest exponential contributions.

##  4 Comparison

We have now calculated the exponentially small contributions that appear in
the asymptotic solution to the differential equation (1) with boundary
condition (2), given by $u_{\mathrm{exp}}$ in (27). We have also calculated
the exponentially small contributions that appear in the asymptotic solution
to the approximate problem (29) with the same boundary condition, given by
$\hat{u}_{\mathrm{exp}}$ in (33). The remaining question is whether the
approximate exponential contributions are able to accurately approximate the
true exponential contributions.

$-0.10$$-0.05$$0.10$$0.05$$0$$0$$5$$1$$2$$3$$4$$x$$u_{\mathrm{exp}}$Pairs
4–7Pair 3Pair 2Pair 1$\hat{u}_{\mathrm{exp}}$$u_{\mathrm{exp}}$ Figure 6:
Exponentially small oscillations in the asymptotic solution to (1) and (29)
for $\epsilon=0.2$. The true exponential contribution $u_{\mathrm{exp}}$ is
shown as a black dashed curve. The approximate exponential contribution
$\hat{u}_{\mathrm{exp}}$ is shown as a red curve. This contribution was
generated using the poles and residues from Table 2. The two curves are
visually indistinguishable. The contribution $\hat{u}_{\mathrm{exp}}$ is the
sum of contributions from each of the seven pole pairs. These contributions
are shown individually as blue curves; it is apparent that the largest
contributions arise from the pole pairs that are nearest to $x=\pm\mathrm{i}$
(ie. pole pairs 1, 2, and 3), with the amplitude of the contributions decaying
as the distance of the pair from $x=\pm\mathrm{i}$ increases.

Note that (27) and (33) cannot be identical as $\epsilon\to 0$, due to the
following differences:

  * •

The algebraic prefactor of $u_{\mathrm{exp}}$ is proportional to
$\epsilon^{-1/2}$, while the algebraic prefactor of the approximate
exponential behaviour $\hat{u}_{\mathrm{exp}}$ is proportional to
$\epsilon^{-1}$. This difference in exponent will be generic: the AAA rational
function must have simple poles that generate algebraic prefactors that are
proportional to $\epsilon^{-1}$, while more general singular points lead to
prefactors whose algebraic power depends on the strength of the singularity.

  * •

The exponential scaling of each expression is determined by the location of
singular points in the leading-order solutions $u_{0}$ (9) and $\hat{u}_{0}$
(30). The branch points in $u_{0}$ that produce (27) are located at
$x=\pm\mathrm{i}$, while the poles in $\hat{u}_{0}$ that produce (33) are
located at $x=p_{r}$, where $|\mathrm{Im}(p_{r})|>1$. Hence the two
expressions have different exponential decay as $\epsilon\to 0$, which can be
seen by directly comparing the two expressions.

These differences appear to suggest that the two expressions
$u_{\mathrm{exp}}$ and $\hat{u}_{\mathrm{exp}}$ cannot possibly describe the
same behaviour, due to the difference in the strength and position of the
singular points between the exact and approximate leading-order behaviour.
Indeed, the two expressions must be different in the limit that $\epsilon\to
0$, as all of the exponential terms in (33) are smaller than the exponential
terms in (27) in the asymptotic limit (even though the algebraic prefactor is
larger).

Despite this, we find that $\hat{u}_{\mathrm{exp}}$ is able to accurately
approximate $u_{\mathrm{exp}}$ for some range of $\epsilon$. We compare the
two contributions for our sample problem with $\epsilon=0.2$ over a segment of
the positive real axis in Figure 6. The contributions $u_{\mathrm{exp}}$ and
$\hat{u}_{\mathrm{exp}}$ appear indistinguishable in this figure.

In addition to the exponentially small contributions, Figure 6 also shows the
contribution to $\hat{u}_{\mathrm{exp}}$ from each pole pair, and reveals that
the largest contributions come from pole pairs 1–3. This indicates that pairs
which are nearest to the branch points at $x=\pm\mathrm{i}$ also have the most
significant effect on the exponentially-small behaviour of the solution, which
is consistent with the exponential decay of the solution contributions as
$|\mathrm{Im}(p_{r})|$ increases.

$1$$0.5$$0$$0.05$$0.10$$0$$0.15$$\mathrm{Tol.}=10^{-14}$, 10
pairs$\mathrm{Tol.}=10^{-13}$, 9 pairs$\mathrm{Tol.}=10^{-12}$, 9
pairs$\mathrm{Tol.}=10^{-11}$, 8 pairs$\mathrm{Tol.}=10^{-10}$, 7
pairs$\mathrm{Tol.}=10^{-9}$, 7 pairs$\mathrm{Tol.}=10^{-8}$, 6
pairs$\mathrm{Tol.}=10^{-7}$, 5 pairs$\mathrm{Tol.}=10^{-8}$, 4
pairs$\epsilon$Amplitude of $\hat{u}_{\mathrm{exp}}$/Amplitude of
$u_{\mathrm{exp}}$$\epsilon=|p_{\mathrm{n}}-\mathrm{i}|$ Figure 7: Ratio of
the amplitudes of the approximate exponential contribution
$\hat{u}_{\mathrm{exp}}$ and the true exponential contribution
$u_{\mathrm{exp}}$ as $\epsilon$ is varied. If $\hat{u}_{\mathrm{exp}}$ is an
accurate approximation of $u_{\mathrm{exp}}$, this ratio is close to 1. The
two expressions have different behaviour as $\epsilon\to 0$, so it is
impossible for the approximation to remain accurate indefinitely as $\epsilon$
is decreased. This is apparent in the figure as each approximation is accurate
until some lower threshold value of $\epsilon$ is reached, beyond which the
approximation becomes inaccurate. Increasing the tolerance, and hence the
number of poles in the AAA approximation, reduces this threshold value of
$\epsilon$. We also identify the value of $\epsilon$ equal to
$|p_{1}-\mathrm{i}|$ on each curve. This value of $\epsilon$ corresponds to a
roughly constant value of the error in each curve, suggesting that the
approximation error of the exponential terms is related the quantity
$|p_{1}-\mathrm{i}|$.

In Figure 7, we compare the amplitude of $u_{\mathrm{exp}}$ and
$\hat{u}_{\mathrm{exp}}$ over a range of $\epsilon$ values. We depict several
different curves, each of which shows the results for a different AAA
algorithm tolerance.

For each curve, the ratio is approximately 1 (indicating that the
approximation is accurate) until some threshold value of $\epsilon$ is reached
from above; below this threshold value the ratio tends to 0. This is
consistent with the algebraic expressions for the two terms, as the
exponential contributions in $\hat{u}_{\mathrm{exp}}$ are smaller than the
contributions in $u_{\mathrm{exp}}$ as $\epsilon\to 0$.

Reducing the error tolerance of the AAA algorithm has the effect of decreasing
the value of $\epsilon$ for which the approximation of the exponential terms
is accurate. Furthermore, this reduction does not happen in a uniform fashion.
Changing the tolerance from $10^{-9}$ to $10^{-10}$ or from $10^{-12}$ to
$10^{-13}$ does not change the number of poles in the approximate leading-
order behaviour, and we see that this change has little effect on the
threshold value of $\epsilon$.

It is apparent from Figure 7 that higher AAA tolerances produce approximations
with more poles, that are accurate for smaller values of $\epsilon$. In Table
2, we show the effect of the error in branch cut prediction, or the difference
in location between the true branch point (at $x=\mathrm{i}$) and the nearest
pole in the approximation (at $x=p_{1}$), or $|p_{1}-\mathrm{i}|$. We compare
this with the relative error in the amplitude, measured by

$\textrm{Relative Error}=1-\frac{\textrm{Amplitude of
}\hat{u}_{\mathrm{exp}}}{\textrm{Amplitude of }{u}_{\mathrm{exp}}}.$ (34)

AAA Tolerance | No. of Pairs | $\quad|p_{1}-\mathrm{i}|\quad$ | Rel. Err. at $\epsilon=|p_{1}-\mathrm{i}|$
---|---|---|---
$10^{-6}$ | 4 | $0.0450$ | $0.1607$
$10^{-7}$ | 5 | $0.0327$ | $0.1642$
$10^{-8}$ | 6 | $0.0264$ | $0.1619$
$10^{-9}$ | 7 | $0.0193$ | $0.1656$
$10^{-10}$ | 7 | $0.0178$ | $0.1643$
$10^{-11}$ | 8 | $0.0151$ | $0.1670$
$10^{-12}$ | 9 | $0.0121$ | $0.1670$
$10^{-13}$ | 9 | $0.0121$ | $0.1670$
$10^{-14}$ | 10 | $0.0096$ | $0.1671$
Table 2: This table shows the distance between the true branch point in
$u_{0}$ at $x=\mathrm{i}$ and the nearest pole in $\hat{u}_{0}$, denoted as
$p_{1}$, for the different AAA approximations presented in Figure 7. As the
tolerance of the approximation increases, the distance $|p_{1}-\mathrm{i}|$
decreases. The relative approximation error is evaluated for
$\epsilon=|p_{1}-\mathrm{i}|$, and the values are roughly constant. This
suggests that the approximation error depends on the quantity
$|p_{1}-\mathrm{i}|$, or how accurately the AAA approximation predicts the
location of the true branch point.

Table 2 suggests that, no matter what AAA tolerance is chosen, the relative
error remains relatively similar at $\epsilon=|p_{1}-\mathrm{i}|$. The points
$\epsilon=|p_{1}-\mathrm{i}|$ are also marked in Figure 7. This figure shows
that for $\epsilon\ll|p_{1}-\mathrm{i}|$, the approximated amplitudes are
inaccurate while for $\epsilon\gg|p_{1}-\mathrm{i}|$, the approximate
amplitudes are accurate. We conjecture that if there is a singularity in the
true leading-order behaviour at $x=x_{\mathrm{s}}$ and the nearest pole to
this point in the rational approximation is at $x=p_{1}$, the exponential
asymptotic predictions will be accurate as long as
$|p_{1}-x_{\mathrm{s}}|\ll\epsilon$.

Obviously, this method is most valuable for problems in which $x_{\mathrm{s}}$
is not known beforehand, so $|p_{1}-x_{\mathrm{s}}|$ is not known. It is
possible to estimate $x_{\mathrm{s}}$ by carefully examining the accumulation
rate of poles in the complex plane [29]. We intend to study this in more
detail in future work, and to explain the dependence of the threshold value on
$|p_{1}-x_{\mathrm{s}}|$.

##  5 A Nonlinear Differential Equation

From this analysis, it appears that rational approximation methods can be used
to apply exponential asymptotic methods to linear ordinary differential
equations, as long as care is taken to ensure that $\epsilon$ is not too
small. The next natural question is whether we can extend these methods to
nonlinear ordinary differential equations, using the exponential asymptotic
method of [7].

This is a challenging problem, because we can no longer add the pole
contributions in $\hat{u}_{\mathrm{exp}}$ independently. From [30], we know
that singularities in nonlinear problems that are near to each other (ie.
within a neighbourhood whose width is proportional to a particular power of
$\epsilon$) can interact, and that the resultant asymptotic behaviour changes
due to these interactions. This is likely to occur for leading-order solutions
generated using rational approximation methods due to the accumulation of
poles near the true singular point. In future work, we will determine the
asymptotic corrections to $u_{\mathrm{exp}}$ due to pole interactions in
nonlinear problems.

In many problems it is possible to determine the strength of the singular
points in $u_{0}$ using asymptotic balancing arguments, even if it is not
possible to easily determine their location. In this case, another possible
method for studying these problems is to apply a map to the sampled data so
that the AAA algorithm is being fitted to a function with simple poles, then
to invert the map so that it contains singularities with the correct order.

Consider the following simple example:

$\epsilon^{2}u^{\prime\prime}(x)+u(x)^{2}=\frac{81}{1024}\left(\frac{1}{\sqrt{x-\mathrm{i}}}+\frac{1}{\sqrt{x+\mathrm{i}}}\right)^{2},$
(35)

with the boundary conditions in (2). The constant $81/1024$ was chosen to make
the exponential asymptotic analysis particularly convenient. The true leading-
order behaviour is given by

$u_{0}=\frac{9}{32}\left(\frac{1}{\sqrt{x-\mathrm{i}}}+\frac{1}{\sqrt{x+\mathrm{i}}}\right).$
(36)

We can sample this leading-order solution on a set of points $x_{j}$ to obtain
data points $u_{0}(x_{j})$, as in Section 3.1, and use this as the basis for
an approximation $\hat{u}_{0}$. Using this, we can compute the exponential
contributions that appear across each Stokes curve in the solution
independently and add them together. If we do so, however, we obtain
predictions of the exponentially small behaviour that are inaccurate, because
the analysis does not take into account nonlinear interactions between pole
contributions.

We instead sample the leading-order solution on $x_{j}$ and then square the
output, to obtain $u_{0}(x_{j})^{2}$, which is equivalent to taking samples of

$u_{0}^{2}=\frac{81}{1024}\left(\frac{1}{x-\mathrm{i}}+\frac{1}{x+\mathrm{i}}+\frac{2}{\sqrt{x^{2}+1}}\right).$
(37)

We can obtain a AAA rational approximation for $u_{0}^{2}$ based on this data,
which we denote $\widehat{u_{0}^{2}}$. Finally, we can take the square root of
this rational approximation to obtain

$u_{0}\approx\left(\widehat{u^{2}_{0}}\right)^{1/2}=\sqrt{\sum_{r=0}^{m}\frac{a_{r}}{x-p_{r}}}.$
(38)

This expression can then be used as the leading-order solution for an
exponential asymptotic analysis using the methods from [7].

$0$$x$$u_{\mathrm{exp}}$$5\times 10^{-4}$$-5\times 10^{-4}$$2.5\times
10^{-4}$$-2.5\times 10^{-4}$$7.5\times 10^{-4}$$-7.5\times
10^{-4}$$-20$$20$$-10$$0$$10$$u_{\mathrm{exp}}$$\hat{u}_{\mathrm{exp}}$ Figure
8: Comparison of the true exponential terms and the AAA-approximated
exponential terms for the nonlinear differential equation (35) with
$\epsilon=0.1$. The AAA-approximated exponential terms were obtained using the
leading-order from (38). In order to show that this method captures nonlinear
effects, we show the expression on the entire real axis, but note that this
contribution is only actually present in the asymptotic solution for $x>0$.
The two contributions are visually similar but not completely identical, due
to nonlinear effects caused by the subdominant branch points at
$x=\pm\mathrm{i}$ in (37).

In Figure 8, we present the true exponentially small terms $u_{\mathrm{exp}}$
and the approximated exponentially small terms $\hat{u}_{\mathrm{exp}}$
obtained using this method for $\epsilon=0.1$. The rational approximation used
was generated using sample points on $x\in[-10,10]$ with $\Delta x=0.1$. We
present the exponential contribution on $x\in[-20,20]$ to demonstrate that
this method accurately captures nonlinear wave effects, but we actually expect
this contribution to appear as the Stokes curve intersecting $x=0$ is crossed
from left to right, and therefore only be present in the asymptotic solution
for $x>0$. The qualitative match seen in this figure shows that the
exponentially-small terms in (35) can be accurately approximated using this
method.

The accuracy of this approximation compared to naively approximating the
exponential terms based on $\hat{u}_{0}$ can be understood by comparing the
poles and residues of the AAA approximation of $u_{0}$ and $u_{0}^{2}$, using
the parameters from Figure 8. These are presented in Table 3.

In $u_{0}$, the only true singularities are the branch points at
$x=\pm\mathrm{i}$. In Table 3, we see that the residue of each pole of
$\hat{u}_{0}$ is similar in magnitude, as these poles approximate the effect
of the branch cuts. In $u_{0}^{2}$, there are simple poles at
$x=\pm\mathrm{i}$, as well as branches that also originate at these points.
The dominant behaviour for this expression is the simple poles. In Table 3, we
see that the poles located closest to $\pm\mathrm{i}$ have a larger residue,
with the other poles having a smaller residue. This is because the poles
nearest to $\pm\mathrm{i}$ in the AAA approximation are reproducing the simple
pole behaviour in the true expression for $u_{0}^{2}$, and the remaining poles
are approximating the subdominant branch cut behaviour.

| $\hat{u}_{0}$ | $\left(\widehat{u^{2}_{0}}\right)^{1/2}$
---|---|---
| $p_{r}$ |  $|a_{r}|$ | $p_{r}$ | $|a_{r}|$
Pair 1 | $0.0020\pm 1.0172\mathrm{i}$ |  $0.2173$ | $0.0018\pm 1.0004\mathrm{i}$ | $0.2845$
Pair 2 | $0.0180\pm 1.1582\mathrm{i}$ |  $0.2210$ | $-0.0067\pm 1.1204\mathrm{i}$ | $0.1342$
Pair 3 | $0.0507\pm 1.4604\mathrm{i}$ |  $0.2287$ | $-0.0207\pm 1.3724\mathrm{i}$ | $0.1329$
Pair 4 | $0.1014\pm 1.9686\mathrm{i}$ |  $0.2404$ | $-0.0381\pm 1.7959\mathrm{i}$ | $0.1335$
Pair 5 | $0.1728\pm 2.7620\mathrm{i}$ |  $0.2566$ | $-0.0550\pm 2.4573\mathrm{i}$ | $0.1330$
Pair 6 | $0.2726\pm 3.9795\mathrm{i}$ |  $0.2793$ | $-0.0665\pm 3.4701\mathrm{i}$ | $0.1374$
Table 3: Comparison of the six pole pairs nearest to $x=\pm\mathrm{i}$
contained in a rational approximation for $u_{0}$ from (36) with the six
nearest pole pairs in the approximation for $u_{0}^{2}$ from (37). The poles
in $\hat{u}_{0}$ have residues with similar magnitude, as the singularities in
$u_{0}$ are branch points at $x=\pm\mathrm{i}$. In
$(\widehat{u}_{0}^{2})^{1/2}$, the poles nearest to $x=\pm\mathrm{i}$ have a
larger residue than the remaining poles, as it approximating the contribution
from the simple poles in (37). The remaining poles mimic the behaviour of the
subdominant branch cut.

Because the strongest singularities in $u_{0}^{2}$ from (37) are simple poles,
the AAA approximation can be used as the basis for an accurate exponential
asymptotic analysis. However, this expression does still contain branch cuts
originating at $x=\pm\mathrm{i}$. This produces a string of additional poles
in the AAA approximation. Hence, any higher-order corrections to the
exponential terms, such as those studied in [23, 8], will be inaccurate unless
pole interaction effects are taken into account. A more thorough analysis of
nonlinear differential equations using rational approximation methods is
beyond the scope of this article and will be the subject of future work.

##  6 Conclusions and Discussion

In this study, we applied exponential asymptotic methods from [22] to study
Stokes’ phenomenon in a linear ordinary differential equation in the
small-$\epsilon$ limit. The leading-order solution $u_{0}$ (9) contains two
branch points, located at $x=\pm\mathrm{i}$. These branch points produce
Stokes curves, and oscillating exponentially small terms appear in the
solution as these Stokes curves are crossed at $x=0$ on the real axis. We
calculated the form of these exponentially small oscillations,
$u_{\mathrm{exp}}$ (27).

We then repeated this process using a rational approximation for the leading-
order behaviour. Instead of using the exact leading-order solution, we sampled
the leading-order behaviour on a discrete set of points and used the AAA
algorithm to produce a rational approximation for the leading-order solution
based on the sampled data, $\hat{u}_{0}$ (28). We then performed an
exponential asymptotic analysis and found that each pair of poles in the
solution produced Stokes curves, each of which generated exponentially small
oscillations. Taking the sum of these oscillations produced the total
exponentially small behaviour $\hat{u}_{\mathrm{exp}}$ (33).

Finally, we compared $u_{\mathrm{exp}}$ and $\hat{u}_{\mathrm{exp}}$, and
found that $\hat{u}_{\mathrm{exp}}$ is able to accurately approximate
$u_{\mathrm{exp}}$ for nonzero values of $\epsilon$, despite having different
asymptotic behaviour in the limit that $\epsilon\to 0$. If $\epsilon$ is too
small, however, the approximation is inaccurate. By reducing the tolerance of
the AAA algorithm, and therefore increasing the accuracy of the rational
approximation, we can reduce the threshold value of $\epsilon$ beyond which
the approximation loses accuracy.

Empirical tests suggest that this threshold value of $\epsilon$ is
proportional to the distance between the branch point in $u_{0}$ and the
nearest pole in $\hat{u}_{0}$. This is a measure of how accurately the
rational approximation predicts the true singularity location. While the true
branch point (or other singularity) in $u_{0}$ cannot typically be calculated
in practice, it is possible to estimate the true location by studying the rate
at which the poles in $\hat{u}_{0}$ accumulate [29].

There are two promising avenues available for extending this method to study
nonlinear ordinary differential equations. The first is to determine the
asymptotic corrections to the approximated exponentials that are caused by
nonlinear interactions between nearby simple poles. The second is to use
asymptotic arguments to determine the strength of singularities in the
leading-order behaviour, and apply the AAA algorithm to a mapped version of
the data that contains simple poles instead of other singularities. We briefly
showed in Figure 8 that this method can be used successfully for a nonlinear
ordinary differential equation, and plan to explore this and related problems
in future work.

This approach has connections to recent work [10, 11] where Padé approximants
are generated from numerical truncated power series data (as is typically
available in Borel plane analyses [13]). These approximations generate
accumulations of simple poles near branch points in the original Borel-
transformed function. The authors use conformal mappings to convert these into
poles, which may be more easily studied. It is likely that a variant of the
the conformal map techniques developed in [11] could be applied to the AAA
approximation to $u_{0}(z)$ developed in the present work, as well as more
complicated leading-order behaviour.

Finally, it is important to determine how robust this method is to noisy input
data. The effect of noise on numerical rational approximation has been studied
in [12], the authors study the effects of noise on conformal maps generated
using Padé approximation, and in [31], in which the authors study the effects
of noise on rational approximations generated using the AAA algorithm. It
would be valuable to use similar methods to study the effect of noise on
exponential asymptotic analyses.

##  7 Acknowledgements

CJL has been supported by Australian Research Council Discovery Project
DP190101190. The authors would like to thank the Isaac Newton Institute (INI)
for Mathematical Sciences for support and hospitality during the programme
“Applicable resurgent asymptotics: towards a universal theory”, where part of
the work on this paper was undertaken. The INI programme was supported by the
EPSRC grant no EP/R014604/1. The authors would like to thank the Okinawa
Institute of Science and Technology (OIST) for support and hospitality. Part
of this research was conducted while visiting OIST through the Theoretical
Sciences Visiting Program (TSVP). CJL would like to thank Nick L. Trefethen
and John R. King for valuable discussions.

## References

  * [1] G. A. Baker Jr, A. George, and J. L. Gammel.  The Padé approximant.  J. Math. Anal. and Appl., 2(1):21–30, 1961.
  * [2] M. V. Berry.  Stokes’ phenomenon; smoothing a Victorian discontinuity.  Pub. Math. de L’IHÉS, 68:211–221, 1988.
  * [3] M. V. Berry.  Uniform asymptotic smoothing of Stokes’s discontinuties.  Proc. Roy. Soc. Lond. A, 422(1862):7–21, 1989.
  * [4] M. V. Berry.  Asymptotics, superasymptotics, hyperasymptotics.  In H. Segur, S. Tanveer, and H. Levine, editors, Asymptotics Beyond All Orders, pages 1–14. Plenum Publishing Corporation, Amsterdam, The Netherlands, 1991.
  * [5] M. V. Berry and C. J. Howls.  Hyperasymptotics.  Proc. Roy. Soc. Lond. A, 430(1880):653–668, 1990.
  * [6] M. V. Berry and C. J. Howls.  Hyperasymptotics for integrals with saddles.  Proc. Roy. Soc. Lond. A, 434(1892):657–675, 1991.
  * [7] S. J. Chapman, J. R. King, and K. L. Adams.  Exponential asymptotics and Stokes lines in nonlinear ordinary differential equations.  Proc. Roy. Soc. Lond. A, 454(1978):2733–2755, 1998.
  * [8] S. J. Chapman and D. B. Mortimer.  Exponential asymptotics and stokes lines in a partial differential equation.  Proc. Roy. Soc. Lond. A, 461(2060):2385–2421, 2005.
  * [9] S. J. Chapman, P. H. Trinh, and T. P. Witelski.  Exponential asymptotics for thin film rupture.  SIAM J. Appl. Math., 73:232–253, 2013.
  * [10] O. Costin and G. V. Dunne.  Physical resurgent extrapolation.  Physics Letters B, 808:135627, sep 2020.
  * [11] O. Costin and G. V. Dunne.  Uniformization and constructive analytic continuation of taylor series.  Communications in Mathematical Physics, 392(3):863–906, 2022.
  * [12] O. Costin, G. V. Dunne, and M. Meynig.  Noise effects on padé approximants and conformal maps.  Journal of Physics A: Mathematical and Theoretical, 55(46):464007, 2022.
  * [13] S. Crew and P. H. Trinh.  Resurgent aspects of applied exponential asymptotics.  arXiv preprint arXiv:2208.07290, 2022.
  * [14] G. Deng, C. J. Lustri, and M. A. Porter.  Nanoptera in weakly nonlinear woodpile chains and diatomic granular chains.  SIAM J. Appl. Dyn. Syst., 20(4):2412–2449, 2021.
  * [15] G. Deng and C.J. Lustri.  Nanoptera in nonlinear woodpile chains with zero precompression.  Physica D, 429:133053, 2022.
  * [16] G. Deng and C.J. Lustri.  Exponential asymptotics of woodpile chain nanoptera using numerical analytic continuation.  Stud. Appl. Math., 150(2):520–557, 2023.
  * [17] R. B. Dingle.  Asymptotic Expansions: Their Derivation and Interpretation.  Academic Press, New York, NY, USA, 1973.
  * [18] A. Gopal and L. N. Trefethen.  Solving laplace problems with corner singularities via rational functions.  SIAM J. on Numer. Anal., 57(5):2074–2094, 2019.
  * [19] F. B. Hildebrand.  Introduction to numerical analysis.  Courier Corporation, 1987.
  * [20] Y. Nakatsukasa, O. Sète, and L. N. Trefethen.  The AAA algorithm for rational approximation.  SIAM J. Sci. Comput., 40:1494–1522, 2018.
  * [21] A. B. Olde Daalhuis.  Hyperasymptotics and the Stokes’ phenomenon.  Proceedings of the Royal Society of Edinburgh Section A: Mathematics, 123(4):731–743, 1993.
  * [22] A. B. Olde Daalhuis, S. J. Chapman, J. R. King, J. R. Ockendon, and R. H. Tew.  Stokes phenomenon and matched asymptotic expansions.  SIAM J. Appl. Math., 55(6):1469–1483, 1995.
  * [23] J. Shelton, S. Crew, and P. H. Trinh.  Exponential asymptotics and higher-order Stokes phenomenon in singularly perturbed ODEs.  arXiv preprint arXiv:2303.07866, 2023.
  * [24] H. Stahl.  The convergence of Padé approximants to functions with branch points.  J. Approx. Theory, 91(2):139–204, 1997.
  * [25] G. G. Stokes.  On the discontinuity of arbitrary constants which appear in divergent developments.  Trans. Cam. Phil. Soc., 10:106–128, 1864.
  * [26] G. G. Stokes.  On the discontinuity of arbitrary constants that appear as multipliers of semi-convergent series.  Acta Mathematica, 26(1):393–397, 1902.
  * [27] L. N. Trefethen.  Numerical conformal mapping with rational functions.  Comput. Meth. Funct. Th., 20:369–387, 2020.
  * [28] L. N. Trefethen.  Numerical analytic continuation.  Japan J. Ind. Appl. Math., pages 1–50, 2023.
  * [29] L. N. Trefethen, Y. Nakatsukasa, and J. A. Weideman.  Exponential node clustering at singularities for rational approximation, quadrature, and PDEs.  Numer. Math., pages 1–28, 2021.
  * [30] P. H. Trinh and S. J. Chapman.  Exponential asymptotics with coalescing singularities.  Nonlinearity, 28(5):1229, 2015.
  * [31] H. Wilber, A. Damle, and A. Townsend.  Data-driven algorithms for signal processing with trigonometric rational functions.  SIAM Journal on Scientific Computing, 44(3):C185–C209, 2022.