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 leadingorder solution $u_0$ contains two branch points, namely
with the boundary conditions
Differential equation (1.1) possesses a closedform 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
By substituting this expression into (1.1) and matching powers of $\epsilon $ in the limit that $\epsilon \to 0$ , we find that
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 leadingorder solution (1.4) of the differential equation. We will find that the exponentially small oscillations appear as special curves in the complex xplane, 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 leadingorder 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 leadingorder solution $u_0$ analytically [Reference Chapman, Trinh and Witelski9, Reference Deng and Lustri14–Reference Deng, Lustri and Porter16]. In this case, the leadingorder 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 leadingorder behaviour.
We will sample the true leadingorder 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 leadingorder 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 leadingorder 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 illposed, 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,
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 highlevel 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
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 leastsquares 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 ${ {fdn} }_{X^{(m)}}$ subject to ${{w}}_m=1$ , where ${ {\cdot } }_{X^{(m)}}$ is the discrete 2norm over $X^{(m)}$ and ${ {\cdot } } _m$ is the discrete $L^{2}$ on mvectors. If the 2norm 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
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
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 leadingorder 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 leadingorder behaviour is given in (1.4).
Matching at $\mathcal {O}(\epsilon ^n)$ as $\epsilon \to 0$ gives
This recurrence relation can be used to calculate $u_n$ exactly, giving
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.
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
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
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
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
Using (2.2) and (2.1) to simplify gives
We will solve this differential equation using variation of parameters. The first step is to determine the homogeneous solutions to (2.7), given by
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
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
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,
Applying this transformation to (2.7) and rearranging gives
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
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.
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
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
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
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,
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 finiteamplitude oscillations present as $x \to \infty $ . This condition indicates that $u_{\mathrm {exp}} = 0$ on the lefthand 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 righthand side of the Stokes line, the exponentially small solution contribution $u_{\mathrm {exp}}$ must be
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 leadingorder 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 leadingorder 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
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.
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 leadingorder solution $u_0$ , and the approximated leadingorder 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}$ .
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.
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
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
Solving this recurrence relation gives
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.
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 upperhalf 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 lefthand side of the Stokes line, we have $\hat {R}_N = 0$ . On the righthand side of the Stokes line, we find that
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
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 righthand side of the Stokes line that intersects the real line at $x = 0.0015$ .
Note that the exponential contributions in (3.6) associated with upperhalf 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 leadingorder 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 leadingorder 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.
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.
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 leadingorder 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
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 leadingorder 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 leadingorder 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:
with the boundary conditions in (1.2). The constant $81/1024$ was chosen to make the exponential asymptotic analysis particularly convenient. The true leadingorder behaviour is given by
We can sample this leadingorder 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 leadingorder 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
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
This expression can then be used as the leadingorder 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.
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.
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 higherorder 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 leadingorder 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 leadingorder behaviour. Instead of using the exact leadingorder solution $u_0$ , we sampled the leadingorder behaviour on a discrete set of points and used the AAA algorithm to produce a rational approximation for the leadingorder 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 leadingorder 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 Boreltransformed 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 leadingorder 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.