1. Introduction
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
$$ \begin{align} \epsilon^2 u"(x) + u(x) = \frac{1}{\sqrt{x+{i}}} + \frac{1}{\sqrt{x-{i}}}, \end{align} $$
with the boundary conditions
$$ \begin{align} \lim_{x\to-\infty} u(x) = 0, \quad \lim_{x\to-\infty} u'(x) = 0. \end{align} $$
Differential equation (1.1) possesses a closed-form solution in terms of the exponential integral function, but we will instead take an alternative approach and determine an asymptotic expansion for the solution in the limit that
$\epsilon \to 0$
. To determine this expansion, we write a series solution to (1.1)–(1.2) in the form
$$ \begin{align} u(x; \epsilon) \sim \sum_{n=0}^{\infty} \epsilon^{2n}u_n(x). \end{align} $$
By substituting this expression into (1.1) and matching powers of
$\epsilon $
in the limit that
$\epsilon \to 0$
, we find that
$$ \begin{align} u_0(x) = \frac{1}{\sqrt{x+{i}}} + \frac{1}{\sqrt{x-{i}}}. \end{align} $$
The solution
$u(x)$
can be shown to contain nondecaying oscillations in the region
$x> 0$
. These oscillations cannot be described using the power series (1.3), as the amplitude of the oscillations is exponentially small in the limit
$\epsilon \to 0$
. Exponential asymptotic techniques are asymptotic methods that may be used to study behaviour that occurs on this scale [Reference Berry2–Reference Berry and Howls6, Reference Olde Daalhuis, Chapman, King, Ockendon and Tew21].
We will perform an exponential asymptotic analysis on (1.1)–(1.2), using the method proposed in [Reference Olde Daalhuis, Chapman, King, Ockendon and Tew21], to calculate the form of the exponentially small oscillations, denoted
${u}_{\mathrm {exp}}$
. As is typical of exponential asymptotic techniques, this method requires the explicit calculation of the leading-order solution (1.4) of the differential equation. We will find that the exponentially small oscillations appear as special curves in the complex x-plane, which originate at singularities of
$u_0$
in the complex plane [Reference Dingle17] and are known as “Stokes lines”, are crossed. It is important to note that the asymptotic size of
$u_{\mathrm {exp}}$
as
$\epsilon \to 0$
typically depends on the location and strength of the singularity. Knowing the analytical form of the leading-order solution
$u_0$
near singularities often provides enough information to calculate the Stokes lines and approximate the exponential contributions [Reference Chapman, King and Adams7].
For many practical problems, it is impossible to calculate the leading-order solution
$u_0$
analytically [Reference Chapman, Trinh and Witelski9, Reference Deng and Lustri14–Reference Deng, Lustri and Porter16]. 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 [Reference Chapman, Trinh and Witelski9] and rational approximation methods [Reference Deng and Lustri15], can be used in conjunction with exponential asymptotics.
The purpose of this paper is to determine whether numerical rational approximation can be combined with exponential asymptotics to correctly approximate exponentially small asymptotic effects. 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 adaptive Antoulas–Anderson (AAA) algorithm [Reference Nakatsukasa, Sète and Trefethen20] 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.
We will sample the true leading-order solution
$u_0$
on a discrete set of points and use these data as the input for an AAA 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 significantly 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 (that is, reducing the
$L^2$
error threshold on the set of sample points).
The difference in asymptotic behaviour between
$u_{\mathrm {exp}}$
and
$\hat {u}_{\mathrm {exp}}$
is a consequence of the fact that all singularities in
$\hat {u}_0$
must be simple poles; if
$u_0$
contains other singularities, such as branch points, the branch will typically be approximated as an accumulation of simple poles instead [Reference Gopal and Trefethen18, Reference Trefethen24, Reference Trefethen, Nakatsukasa and Weideman26]. Hence, the strength of the singularity in
$u_0$
is generically different to the simple poles present in
$\hat {u}_0$
. Recall that the asymptotic behaviour of
$u_{\mathrm {exp}}$
depends on the strength of the singularity. It is therefore remarkable that, despite such a numerical approximation having completely different singularity types and locations, the asymptotic behaviour
$u_{\mathrm {exp}}$
can be well approximated by
$\hat {u}_{\mathrm {exp}}$
.
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 difference may be thought of as the approximation error of the true branch point location from the rational approximation algorithm.
1.1. Numerical analytic continuation
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 [Reference Hildebrand19]. 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, can lead to significant changes in the analytically continued result.
Despite this challenge, significant progress has been made on performing analytic continuation using numerical methods, which typically require 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 [Reference Baker and Gammel1] and the AAA algorithm [Reference Nakatsukasa, Sète and Trefethen20].
Rational approximation methods typically approximate some function
$f(x)$
as a rational expression,
$$ \begin{align} f(x) \approx\frac{n(x)}{d(x)}, \end{align} $$
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 take 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.
One such method 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 [Reference Nakatsukasa, Sète and Trefethen20].
The algorithm is based on expressing the rational approximation in the form
$$ \begin{align} f(x) \approx\frac{n(x)}{d(x)}={\sum_{j=1}^m\frac{w_jf_j}{x-x_j}} \, \bigg/ \, {\sum_{j=1}^m\frac{w_j}{x-x_j}}. \end{align} $$
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 optimisation problem 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,\ldots ,w_m)^T$
is chosen to minimise
${|| {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 maximises the quantity
$$ \begin{align} \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}. \end{align} $$
The process then repeats until the termination criteria are 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 minimises 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
$$ \begin{align} f(x) \approx \sum_{j=0}^{m} \frac{a_j}{x - p_j}, \end{align} $$
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 [Reference Gopal and Trefethen18, Reference Trefethen24, Reference Trefethen, Nakatsukasa and Weideman26] 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 rigorously in the related problem of Padé approximation [Reference Stahl23].
In the following study, we show that if we approximate the leading-order solution of an ordinary differential equation with a rational function of the form (1.8), it is possible to use exponential asymptotics on this rational function to approximate the exponentially small terms that appear in the solution due to Stokes’ phenomenon.
1.2. Paper outline
In Section 2, we perform an exponential asymptotic analysis of the system (1.1)–(1.2) to determine the Stokes switching behaviour in the solution. In Section 3, we use an AAA rational approximation to solve (1.1)–(1.2) with
$\epsilon = 0$
, and use this as the basis for an exponential asymptotic analysis. In Section 4, we compare the results from the two preceding sections, and show that the AAA rational approximation can correctly describe the Stokes switching behaviour. Finally, in Section 5, we discuss how this method can be applied to study a nonlinear ordinary differential equation, and describe the additional challenges to this method caused by nonlinearity.
2. Exponential asymptotic analysis
2.1. Asymptotic power series
We propose an asymptotic series solution to (1.1)–(1.2) of the form presented in (1.3). By substituting this expression into the ordinary differential equation (1.1) and matching powers of
$\epsilon $
in the limit that
$\epsilon \to 0$
, we find that the leading-order behaviour is given in (1.4).
Matching at
$\mathcal {O}(\epsilon ^n)$
as
$\epsilon \to 0$
gives
$$ \begin{align} u_n = -u_{n-1}", \quad n \geq 1. \end{align} $$
This recurrence relation can be used to calculate
$u_n$
exactly, giving
$$ \begin{align} u_n = \frac{(4n)!(-1)^n}{2^{4n}(2n)!(x-{i})^{2n+{1/2}}} + \mathrm{c.c.}, \end{align} $$
where c.c. denotes the complex conjugate contribution. The series (1.3) with terms given by (2.2) 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.
Figure 1 Magnitude of series terms from (1.3) 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 (2.2). The series (1.3) 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 [Reference Berry3]. By rescaling to study the exact truncation remainder, we can directly calculate exponentially small components of the asymptotic expansion.
2.2. Finding the Stokes lines
The analysis presented here is a standard application of the exponential asymptotic method developed in [Reference Olde Daalhuis, Chapman, King, Ockendon and Tew21] to study Stokes’ phenomenon. Nonetheless, this section will provide a relatively complete summary of the procedure to make the analysis accessible.
The purpose of this analysis is to obtain the leading exponentially small correction term to the series (1.3), 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 (2.2), 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 (1.3) after N terms, giving
$$ \begin{align} u(x; \epsilon) = \sum_{n=0}^{N-1} \epsilon^{2n}u_n(x) + R_N(x), \end{align} $$
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
$$ \begin{align} \bigg|\frac{(4n)!(-1)^n \epsilon^{2n}}{2^{4n}(2n)!(x-{i})^{2n+{1/2}}} \bigg| = \bigg|\frac{(4n+4)!(-1)^{n+1}\epsilon^{2n+2}}{2^{4n+4}(2n+2)!(x-{i})^{2n+{5/2}}}\bigg|. \end{align} $$
Optimal truncation occurs after a large number of terms in the limit that
$\epsilon \to 0$
. We can make use of this fact to find that the optimal truncation point must satisfy
$N_{\mathrm {opt}} \sim |x-{i}|/2\epsilon $
as
$\epsilon \to 0$
. We therefore set
$$ \begin{align} N_{\mathrm{opt}} = \frac{|x-{i}|}{2\epsilon} + \omega, \end{align} $$
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 - {i}| = 1$
and
$\epsilon = 0.1$
, corresponding to
$N_{\mathrm {opt}} = 5$
. We will eventually use this value for N in (2.3), but for algebraic simplicity, we will not make the substitution immediately.
We substitute the truncated series (2.3) into the governing equation (1.1) to obtain
$$ \begin{align} \sum_{n=0}^{N-1} \epsilon^{2n+2} u_n" + \sum_{n=0}^{N-1}\epsilon^{2n} u_n + \epsilon^2 R_N" + R_N = \frac{1}{\sqrt{x+{i}}} + \frac{1}{\sqrt{x-{i}}}. \end{align} $$
Using (2.2) and (2.1) to simplify gives
$$ \begin{align} \epsilon^2 R_N" + R_N = \epsilon^{2N}u_N = \frac{ (4N)!(-1)^N\epsilon^{2N}}{2^{4N}(2N)!(x-{i})^{2N+{1/2}}}. \end{align} $$
We will solve this differential equation using variation of parameters. The first step is to determine the homogeneous solutions to (2.7), given by
$$ \begin{align} R_N(x) = K_1 {e}^{{i} x/\epsilon}\quad\mathrm{and}\quad R_N(x) = K_2 {e}^{-{i} x/\epsilon}, \end{align} $$
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 (2.7). We will find that for (2.7), the second solution is the relevant one, so we set
$$ \begin{align} R_N(x) = K_2(x){e}^{-{i} x/\epsilon} = \mathcal{S}(x){e}^{-{i} (x-{i})/\epsilon}, \end{align} $$
where
$K_2(x) = \mathcal {S}(x){e}^{-1/\epsilon }$
, with the scaling chosen so that the exponent in (2.9) is equal to 0 when
$x = {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 line. 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 (2.9) into (2.7) and simplifying gives
$$ \begin{align} -2 {i} \epsilon \mathcal{S}' {e}^{-{i}(x-{i})/\epsilon} = \frac{ (4N)!(-1)^N\epsilon^{2N}}{2^{4N}(2N)!(x-{i})^{2N+{1/2}}}. \end{align} $$
Recall that the optimal truncation given in (2.5) depends on
$|x-{i}|$
. This observation suggests that the transformation
${i}(x - {i}) = r{e}^{{i}\theta }$
will simplify this expression in a useful way, giving
$N_{\mathrm {opt}} = r/2\epsilon + \omega $
. As in [Reference Olde Daalhuis, Chapman, King, Ockendon and Tew21], we then fix r and consider only the faster variation that occurs in the angular direction
$\theta $
. From the chain rule,
$$ \begin{align} \frac{d}{dx} = \frac{{e}^{-{i}\theta}}{r}\frac{d}{d\theta}. \end{align} $$
Applying this transformation to (2.7) and rearranging gives
$$ \begin{align} \frac{d\mathcal{S}}{d\theta} = -\frac{{i}\sqrt{{i}}r{e}^{{i}\theta}(4N)!\epsilon^{2N-1}}{2^{4N+1} (2N)!(r{e}^{{i}\theta})^{2N+1/2}}\exp \bigg(\frac{r{e}^{{i}\theta}}{\epsilon}\bigg) \quad \text{as } \epsilon \to 0. \end{align} $$
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
$$ \begin{align} \frac{d\mathcal{S}}{d\theta} \sim -\frac{{i}\sqrt{2{i} r}}{\epsilon}\exp\bigg(\frac{r}{\epsilon}({e}^{{i}\theta} - 1 - {i}\theta) + {i}\theta\bigg(\frac{1}{2} + \omega\bigg)\bigg)\quad\text{as } \epsilon \to 0. \end{align} $$
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 line.
Recall that
${i}(x-{i}) = r{e}^{{i}\theta }$
, so
$\theta = 0$
corresponds to
$x - {i}$
taking negative imaginary values, or a vertical line extending downwards from the branch point
$x = {i}$
. If we perform a corresponding analysis for the complex conjugate term in (2.2), we identify that there is a second Stokes line extending vertically upwards from the branch point
$x = -{i}$
. This Stokes structure is illustrated in Figure 2.
Figure 2 Stokes’ phenomenon in the exact solution to (1.1) with boundary conditions (1.2). The leading-order solution
$u_0$
contains branch points at
$x = \pm {i}$
, with the branch cuts extending vertically away from the real axis. The branch points generate Stokes lines, which connect the two points. On the left-hand side of the Stokes lines, there are no exponential contributions. On the right-hand side of the Stokes lines, the solution contains exponentially small oscillations of the form given in (2.18).
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 line is crossed, we define a new inner variable
$\theta = \epsilon ^{1/2}\vartheta $
, which gives
$$ \begin{align} \frac{d\mathcal{S}}{d\vartheta} \sim -{i}\sqrt{\frac{2{i} r}{\epsilon}}{e}^{-r\vartheta^2/2} \quad \text{as } \epsilon \to 0. \end{align} $$
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 line is given by
$$ \begin{align} [\mathcal{S}]_-^+ = \lim_{\vartheta\to +\infty} \mathcal{S} -\lim_{\vartheta\to -\infty} \mathcal{S} \sim -{i}\sqrt{\frac{2{i} r}{\epsilon}}\int_{-\infty}^{\infty}{e}^{-r\vartheta^2/2}{\,{d}} \vartheta = -2 {i}\sqrt{\frac{\pi{i}}{\epsilon}}. \end{align} $$
If the value of
$\mathcal {S}$
changes by (2.15) as the Stokes line is crossed, the exponentially small contribution to the solution for
$x> 0$
is given by
$$ \begin{align} [R_N]_-^+ \sim -2 {i}\sqrt{\frac{\pi{i}}{\epsilon}}{e}^{-{i}(x-{i})/\epsilon}. \end{align} $$
If we add in the complex conjugate term, we obtain the jump in the exponentially small terms
$u_{\mathrm {exp}}$
as the Stokes line is crossed from left to right,
$$ \begin{align} [u_{\mathrm{exp}}]_-^+ \sim -2 {i}\sqrt{\frac{\pi{i}}{\epsilon}}\,{e}^{-{i}(x-{i})/\epsilon} + \mathrm{c.c.} = 4\sqrt{\frac{\pi}{\epsilon}}\,{e}^{-1/\epsilon}\cos\bigg(\frac{x}{\epsilon}+\frac{\pi}{4}\bigg) \quad \text{as } \epsilon \to 0. \end{align} $$
Note that this expression has constant amplitude. By comparing this observation with the boundary conditions in (1.2), we see that there cannot be any finite-amplitude oscillations present as
$x \to -\infty $
. This condition indicates that
$u_{\mathrm {exp}} = 0$
on the left-hand side of the Stokes line (
$x < 0$
), and that the exponentially small behaviour appears as the Stokes line is crossed from left to right. On the right-hand side of the Stokes line, the exponentially small solution contribution
$u_{\mathrm {exp}}$
must be
$$ \begin{align} u_{\mathrm{exp}} \sim 4\sqrt{\frac{\pi}{\epsilon}}\,{e}^{-1/\epsilon}\cos\bigg(\frac{x}{\epsilon}+\frac{\pi}{4}\bigg),\quad x> 0. \end{align} $$
This behaviour is illustrated in Figure 2.
3. AAA approximation
The analysis in Section 2 was a relatively straightforward application of the method of [Reference Olde Daalhuis, Chapman, King, Ockendon and Tew21]. It centred on the observation that singularities of the leading-order solution
$u_0$
in the complex plane, such as the branch points at
$x = \pm {i}$
in (1.4), generate Stokes lines. 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$
To simulate a problem where we only possess a numerical description of the leading-order solution
$u_0$
from (1.4), we first define a set of points
$x_j$
separated by some
$\Delta x$
. We then evaluate (1.4) 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 an AAA rational approximation.
We apply the AAA algorithm to obtain a rational approximation
$\hat {u}_0(x)$
, given by
$$ \begin{align} u_0 \approx \hat{u}_0 = \sum_{r=0}^{m} \frac{a_r}{x - p_r}. \end{align} $$
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 shown in Table 1. We have given each pair of poles a designation, so that we may refer to them later.
Table 1 Poles and residues in the AAA approximation of
$u_0$
, which contains branch points at
$x = \pm {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 {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 = {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 in (3.1). 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$
in (1.4). In fact, the AAA algorithm places the poles along a curve, representing the branch cut, which we are not free to prescribe. This limitation is a common feature of rational approximation methods, and was studied in depth for Padé approximation in [Reference Stahl23], 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 [Reference Trefethen25].
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 {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 = {i}$
.
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. (Colour available online.)
In Figure 4, we present the error between the two on a logarithmic plot. The data in this figure confirm 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.
Figure 4 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. (Colour available online.)
3.2. Asymptotic power series
In effect, we are using the AAA approximation in place of the inhomogeneous term in (1.1). Hence, we are determining the asymptotic behaviour of
$$ \begin{align} \epsilon^2 \hat{u}"(x) + \hat{u}(x) = \sum_{r=0}^{m} \frac{a_r}{x - p_r}. \end{align} $$
By expanding the solution as a series (1.3) and matching powers of
$\epsilon $
, we can again obtain a recurrence relation for the series terms
$u_n$
. The recurrence relation gives
$$ \begin{align} \hat{u}_0(x) = \sum_{r=0}^{m} \frac{a_r}{x - p_r}, \quad \hat{u}_n(x) = -\hat{u}_{n-1}"(x). \end{align} $$
Solving this recurrence relation gives
$$ \begin{align} \hat{u}_n = \sum_{r=0}^{m} \frac{a_r (2n)!}{(x - p_r)^{2n+1}}. \end{align} $$
Each of the poles
$p_r$
will generate a Stokes line, and lead to an exponentially small contribution appearing in the solution.
3.3. Stokes lines and exponentially small terms
To determine the location of the Stokes lines, 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 (3.4).
This analysis reveals that there are exponentially small contributions in the solution associated with each pole. These contributions appear across Stokes lines that extend 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 lines generated by pole pairs 1, 2 and 3.
Figure 5 Stokes structure in the solution to (3.2), which uses
$\hat {u}_0$
as the leading-order solution. The solution contains seven pairs of poles, with locations given in Table 1. Each pair of poles generates Stokes lines that extend vertically from the poles, intersecting the real axis. As each Stokes line 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 {i}$
. We will later find that the largest exponential contributions arise from the poles that are nearest to
$x = \pm {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 line generated by a pole in the upper-half plane located at at
$x = p_r$
. The Stokes line 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 line, we have
$\hat {R}_N = 0$
. On the right-hand side of the Stokes line, we find that
$$ \begin{align} \hat{R}_N \sim \frac{2 \pi a_r}{\epsilon}\,{e}^{-{i}(x-p_r)/\epsilon} \quad \text{as } \epsilon \to 0. \end{align} $$
Once all of the Stokes lines 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
$$ \begin{align} \hat{u}_{\mathrm{exp}} \sim \sum_{r=0}^{m}\frac{2 \pi a_r}{\epsilon}\,{e}^{-{i}(x-p_r)/\epsilon} \quad \text{as } \epsilon \to 0. \end{align} $$
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 line that intersects the real line at
$x = -0.0015$
.
Note that the exponential contributions in (3.6) 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 observation 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.1) with boundary condition (1.2), given by
$u_{\mathrm {exp}}$
in (2.18). We have also calculated the exponentially small contributions that appear in the asymptotic solution to the approximate problem (3.2) with the same boundary condition, given by
$\hat {u}_{\mathrm {exp}}$
in (3.6). The remaining question is whether the approximate exponential contributions are able to accurately approximate the true exponential contributions.
Note that (2.18) and (3.6) 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$
in (1.4) and
$\hat {u}_0$
in (3.3). The branch points in
$u_0$
that produce (2.18) are located at
$x = \pm {i}$
, while the poles in
$\hat {u}_0$
that produce (3.6) 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 (3.6) are smaller than the exponential terms in (2.18) in the asymptotic limit (even though the algebraic prefactor is larger).
Despite this difference in asymptotic behaviour, 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.
Figure 6 Exponentially small oscillations in the asymptotic solution to (1.1) and (3.2) 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 {i}$
(that is, pole pairs 1, 2 and 3), with the amplitude of the contributions decaying as the distance of the pair from
$x = \pm {i}$
increases. (Colour available online.)
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 observation indicates that pairs which are nearest to the branch points at
$x = \pm {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.
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.
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 behaviour 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 - {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 to the quantity
$|p_1 - {i}|$
. (Colour available online.)
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 behaviour 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, which 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 = {i}$
) and the nearest pole in the approximation (at
$x = p_1$
), or
$|p_1 - {i}|$
. We compare this value with the relative error in the amplitude, measured by
$$ \begin{align} \textrm{Relative Error} = 1 - \frac{\textrm{Amplitude of } \hat{u}_{\mathrm{exp}}}{\textrm{Amplitude of } {u}_{\mathrm{exp}}}. \end{align} $$
Table 2 The distance between the true branch point in
$u_0$
at
$x = {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 - {i}|$
decreases. The relative approximation error is evaluated for
$\epsilon = |p_1 - {i}|$
, and the values are roughly constant. This observation suggests that the approximation error depends on the quantity
$|p_1 - {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 - {i}|$
. The points
$\epsilon = |p_1 - {i}|$
are also marked in Figure 7. This figure shows that for
$\epsilon \ll |p_1 - {i}|$
, the approximated amplitudes are inaccurate, while for
$\epsilon \gg |p_1 - {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 [Reference Trefethen, Nakatsukasa and Weideman26]. We intend to study this accumulation 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 [Reference Chapman, King and Adams7].
Resolving this question is a challenging problem, because we can no longer add the pole contributions in
$\hat {u}_{\mathrm {exp}}$
independently. From [Reference Trinh and Chapman27], we know that singularities in nonlinear problems that are near to each other (that is, 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 behaviour is likely to occur for leading-order solutions
$\hat {u}_0$
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:
$$ \begin{align} \epsilon^2 u"(x) + u(x)^2 = \frac{81}{1024}\bigg(\frac{1}{\sqrt{x-{i}}} + \frac{1}{\sqrt{x+{i}}}\bigg)^2, \end{align} $$
with the boundary conditions in (1.2). The constant
$81/1024$
was chosen to make the exponential asymptotic analysis particularly convenient. The true leading-order behaviour is given by
$$ \begin{align} u_0 = \frac{9}{32}\bigg(\frac{1}{\sqrt{x-{i}}} + \frac{1}{\sqrt{x+{i}}}\bigg). \end{align} $$
We can sample this leading-order solution
$u_0$
on a set of points
$x_j$
to obtain data points
$u_0(x_j)$
, as in Section 3.1, and use these data as the basis for an approximation
$\hat {u}_0$
. Using this approximation, we can compute the exponential contributions that appear across each Stokes line 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
$u_0$
on
$x_j$
and then square the output, to obtain
$u_0(x_j)^2$
, which is equivalent to taking samples of
$$ \begin{align} u_0^2 = \frac{81}{1024}\bigg(\frac{1}{x-{i}} + \frac{1}{x+{i}} + \frac{2}{\sqrt{x^2+1}}\bigg). \end{align} $$
We can obtain an AAA rational approximation for
$u_0^2$
based on these data, which we denote
$\widehat {u_0^2}$
. Finally, we can take the square root of this rational approximation to obtain
$$ \begin{align} u_0 \approx (\widehat{u^2_0})^{1/2} = \sqrt{\sum_{r=0}^m \frac{a_r}{x - p_r}}. \end{align} $$
This expression can then be used as the leading-order solution for an exponential asymptotic analysis using the methods from [Reference Chapman, King and Adams7]. 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 line 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 (5.1) can be accurately approximated using this method.
Figure 8 Comparison of the true exponential terms and the AAA-approximated exponential terms for the nonlinear differential equation (5.1) with
$\epsilon = 0.1$
. The AAA-approximated exponential terms were obtained using the leading-order from (5.4). 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 {i}$
in (5.3).
The accuracy of this approximation compared with 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 {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 {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 {i}$
have a larger residue, with the other poles having a smaller residue. This behaviour occurs because the poles nearest to
$\pm {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.
Table 3 Comparison of the six pole pairs nearest to
$x = \pm {i}$
contained in a rational approximation for
$u_0$
from (5.2) with the six nearest pole pairs in the approximation for
$u_0^2$
from (5.3). The poles in
$\hat {u}_0$
have residues with similar magnitude, as the singularities in
$u_0$
are branch points at
$x = \pm {i}$
. In
$(\widehat {u}_0^2)^{1/2}$
, the poles nearest to
$x = \pm {i}$
have a larger residue than the remaining poles, as it approximates the contribution from the simple poles in (5.3). The remaining poles mimic the behaviour of the subdominant branch cut.
Because the strongest singularities in
$u_0^2$
from (5.3) 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 {i}$
. These branch cuts produce strings of additional poles in the AAA approximation. Hence, any higher-order corrections to the exponential terms, such as those studied in [Reference Chapman and Mortimer8, Reference Shelton, Crew and Trinh22], 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 [Reference Olde Daalhuis, Chapman, King, Ockendon and Tew21] to study Stokes’ phenomenon in a linear ordinary differential equation in the small-
$\epsilon $
limit. The leading-order solution
$u_0$
in (1.4) contains two branch points, located at
$x = \pm {i}$
. These branch points produce Stokes lines, and oscillating exponentially small terms appear in the solution as these Stokes lines are crossed at
$x = 0$
on the real axis. We calculated the form of these exponentially small oscillations,
$u_{\mathrm {exp}}$
in (2.18).
We then repeated this process using a rational approximation for the leading-order behaviour. Instead of using the exact leading-order solution
$u_0$
, 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$
in (3.1). We then performed an exponential asymptotic analysis and found that each pair of poles in the solution produced Stokes lines, each of which generated exponentially small oscillations. Taking the sum of these oscillations produced the total exponentially small behaviour
$\hat {u}_{\mathrm {exp}}$
in (3.6).
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 result provides 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 [Reference Trefethen, Nakatsukasa and Weideman26].
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 [Reference Costin and Dunne10, Reference Costin and Dunne11] where Padé approximants are generated from numerical truncated power series data (as is typically available in Borel plane analyses [Reference Crew and Trinh13]). 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 [Reference Costin and Dunne11] could be applied to the AAA approximation of
$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 [Reference Costin, Dunne and Meynig12], where the authors study the effects of noise on conformal maps generated using Padé approximation, and in [Reference Wilber, Damle and Townsend28], where 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.
Acknowledgements
C. J. Lustri 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 Engineering and Physical Sciences Research Council (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). C. J. Lustri would like to thank L. Nick Trefethen and John R. King for valuable discussions.
Open access
