EXPONENTIAL ASYMPTOTICS USING NUMERICAL RATIONAL APPROXIMATION IN LINEAR DIFFERENTIAL EQUATIONS

Abstract Singularly perturbed ordinary differential equations often exhibit Stokes’ phenomenon, which describes the appearance and disappearance of oscillating exponentially small terms across curves in the complex plane known as Stokes lines. These curves originate at singular points in the leading-order solution to the differential equation. In many important problems, it is impossible to obtain a closed-form expression for these leading-order solutions, and it is therefore challenging to locate these singular points. We present evidence that the analytic leading-order solution of a linear differential equation can be replaced with a numerical rational approximation using the adaptive Antoulas–Anderson (AAA) method. Despite such an approximation having completely different singularity types and locations, we show that the subsequent exponential asymptotic analysis accurately predicts the exponentially small behaviour present in the solution. For sufficiently small values of the asymptotic parameter, this approach breaks down; however, the range of validity may be extended by increasing the number of poles in the rational approximation. We present a related nonlinear problem and discuss the challenges that arise due to nonlinear effects. Overall, our approach allows for the study of exponentially small asymptotic effects without requiring an exact analytic form for the leading-order solution; this permits exponential asymptotic methods to be used in a much wider range of applications.


Introduction
The behaviour known as "Stokes' phenomenon" was first observed by George Gabriel Stokes in his analysis of solutions to the Airy equation [25,26].The Airy function displays oscillatory behaviour for large negative arguments, and decays exponentially for large positive arguments.Stokes wished to understand the transition between these two regimes by studying the large-argument asymptotic behaviour in the complex plane.In doing so, Stokes discovered that exponentially-small correction terms in the asymptotic expansion appear or disappear as special curves in the complex plane are crossed.These curves are known as 'Stokes curves'.
New asymptotic techniques were developed in order to study the switch in behaviour across Stokes curves.Berry [2][3][4], and Berry & Howls [5,6] showed that truncating a a divergent asymptotic series optimally results in asymptotic error of exponentially small order.Hence, rescaling the problem to study this truncation error directly allows for the calculation of these exponentially small terms.This idea was used to develop new asymptotic methods for studying exponentially small behaviour, known as exponential asymptotics, asymptotics-beyond-all-orders or hyperasymptotics [6,21].We will denote exponentially small contributions to the asymptotic approximation of some solution u as uexp.
Exponential asymptotic techniques typically require the explicit calculation of the leading-order solution of the differential equation, which we will denote as u0.Stokes curves originate at singularities of u0 [17].Additionally, the asymptotic size of uexp as ϵ → 0 typically depends on the location and strength of the singularity.For many problems, knowing the analytical form of the leading-order solution near singularities provides enough information to calculate the Stokes curves and approximate the exponential contributions [7].
For many practical problems, it is impossible to calculate the leading-order solution analytically [9,[14][15][16].In this case, the leading-order solution must instead be approximated numerically.There is growing evidence that numerical methods for analytic continuation, including numerically stepping into the complex plane [9] and rational approximation methods [16] can be used in conjunction with exponential asymptotics.
In general, the leading-order solution to an ordinary differential equation can be computed on a real domain using any one of a variety of standard numerical methods [19].Using the result of numerical computations to describe the complex plane behaviour of these solutions is challenging.It is wellknown that analytic continuation is ill-posed, and that small changes to the function being continued, such as numerical error associated with computational methods, can lead to significant changes in the analytically-continued result.
Despite this, significant progress has been made on performing analytic continuation using numerical methods; this is typically performed by requiring the output function to take a particular algebraic form.
The most commonly-used methods for numerical complex analytic continuation are methods that give a rational function as an output.These rational approximation methods include Padé approximation [1] and the adaptive Antoulas-Anderson (AAA) algorithm [20].
We wish to consider settings in which the solution has been numerically approximated on some discrete set of points along the real axis which we then seek to analytically continue into the complex plane.This setup makes the AAA algorithm particularly well-suited for our purpose.This algorithm is described in Section 1.1, and produces a rational function, denoted here as û0 which approximates the leading-order behaviour.Note that all singularities in û0 must be simple poles; if u0 contains other singularities, such as branch points, the branch will be typically be approximated as an accumulation of simple poles instead [18,27,29].Recall that the strength of the singularity in u0 is generically different to the simple poles present in û0, and the asymptotic behaviour of uexp depends on the strength of the singularity.This argument indicates that any exponential term obtained using this leading-order approximation must display different asymptotic behaviour to uexp in the limit that ϵ → 0.
In this study, we focus our attention on a singularly-perturbed linear ordinary differential equation in the limit that ϵ → 0, whose leading-order solution u0 contains two branch points, namely with the boundary conditions We note that differential equation ( 1) possesses a closed-form solution in terms of the exponential integral function, where A and B may be determined from the boundary conditions.It is possible to determine A and B by performing a careful analysis on the integral expressions in the limit that x → −∞, but we will instead take a different approach, and determine an asymptotic expansion for the solution in the limit that ϵ → 0.
We will perform a typical exponential asymptotic analysis on (1)-( 2), using the method proposed in [22], to determine uexp for this problem.We will then sample u0 on a discrete set of points and use this as the input for a rational approximation û0.We will apply the same exponential asymptotic method on this leading-order expression to obtain ûexp.By comparing the analytic form of uexp and ûexp, we will see that the two expressions have different asymptotic behaviour in the limit ϵ → 0.
We will then show that, despite this difference, ûexp is able to accurately approximate uexp for a range of values of ϵ, and that this range can be extended by increasing the accuracy of the AAA approximation (ie.reducing the L 2 error threshold on the set of sample points).Finally, we will provide evidence that the threshold value of ϵ below which ûexp is inaccurate is proportional to the difference between the true branch point location in u0 and the pole nearest to this point in û0.This may be thought of as the approximation error of the true branch point location from the rational approximation algorithm.

Numerical Analytic Continuation
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 takes a different form; where values of f (x) are calculated numerically on some discrete support set of points.While it is possible to approximate Taylor coefficients using values of the function on a discrete set of points, there are other methods more well-suited to input data with this form.
On such is the AAA approximation, which works by iteratively solving a minimisation problem to produce an optimal rational approximation.The explanation that follows is a high-level description of the AAA algorithm; for a detailed explanation, see [20].
The algorithm is based on expressing the rational approximation in the form The points xj for j = 1, 2, . . ., 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 fj = f (xj).The algorithm solves an optimization problem in order to determine the weights wj, and then generates another support point xj+1.This process then repeats iteratively.
The weights wj are generated by solving a linear least-squares problem over a set of points in the restricted domain X (m) = X \ {x1, . . ., xm}; these points are labelled X The process then repeats until the termination criteria is reached.
The key information for our purposes is that the method takes in the function values on a set of points, and returns a rational approximation that minimizes the L 2 -norm of the approximation error on the data support set, and that it is an iterative process that increases the number of poles in the solution until some predetermined tolerance is met.We often write the approximation in the form where aj are the residues of the function at each pole location pj.
As previously alluded to, one potential obstacle in applying rational approximation methods is that these methods produce meromorphic functions; all singularities in the approximation are isolated simple poles.Previous studies of AAA approximation [18,27,29] showed that the rational function approximation of a target function with a branch point contains an exponential clustering of poles (and zeroes) in the approximation approaching the branch point.This configuration of poles accurately approximates the effects of the branch cut.The distribution of poles when approximating branch cuts in rational approximations has been studied rigourously in the related problem of Padé approximation [24].
In the following we show that if we approximate the leading-order solution of an ordinary differential equation with a rational function of the form (7), it is possible to use exponential asymptotics on this rational function in order to approximate the exponentially small terms that appear in the solution due to Stokes' phenomenon.
2 Exponential Asymptotic Analysis

Asymptotic Power Series
We propose an asymptotic series solution to (1)-( 2) of the form By substituting this into the ordinary differential equation (1) and matching powers of ϵ in the limit that ϵ → 0, we find that the leading-order behaviour is given by the factorial contribution to the numerator of u n in (11) The series (8) must therefore be divergent, and the optimal truncation point occurs at the minimum value.
Matching at O(ϵ n ) as ϵ → 0 gives This recurrence relation can be used calculate un exactly, giving The series (8) with terms given by ( 11) is divergent for any fixed choice of x.Typically, truncating the series at the value of n that minimises ϵ 2n |un| 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 Nopt.
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 ϵ → 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 = Nopt produces a remainder that is exponentially small in the asymptotic limit [3].By rescaling to study the exact truncation remainder, we can directly calculate exponentially small components of the asymptotic expansion.

Finding the Stokes Curves
The analysis presented here is a standard application of the exponential asymptotic method developed in [22] in order to study Stokes' phenomenon.Nonetheless, this section will provide a relatively complete summary of the procedure in order to make the analysis accessible.
The purpose of this analysis is to obtain the leading exponentially small correction term to the series (8), and to therefore reveal the effects of Stokes switching in the solution.We will concentrate on the contribution arising from the term explicitly shown in (11), and note that a complex conjugate contribution is also present in the solution, which we will include in the final expression.
We now truncate the power series (8) after N terms, giving where RN is the remainder.The optimal truncation point Nopt corresponds to minimising ϵ 2N |uN |.The most straightforward way to find Nopt is to determine the value at which consecutive terms have equal magnitude, or (4n Optimal truncation occurs after a large number of terms in the limit that ϵ → 0. We can make use of this to find that the optimal truncation point must satisfy Nopt ∼ |x − i|/2ϵ as ϵ → 0. We therefore set where ω ∈ [0, 1) is chosen such that Nopt is an integer.This expression for Nopt is consistent with Figure 1 which has |x − i| = 1 and ϵ = 0.1, corresponding to Nopt = 5.We will eventually use this value for N in ( 12), but for algebraic simplicity we will not make the substitution immediately.
We substitute the truncated series ( 12) into the governing equation ( 1) to obtain Using ( 11) and ( 10) to simplify gives We will solve this differential equation using variation of parameters.The first step is to determine the homogeneous solutions to (16), given by where K1 and K2 are arbitrary constants.The next step is to permit the arbitrary constants to vary in x, and then consider the full differential equation in (16).We will find that for ( 16) the second solution is the relevant one, so we set where K2(x) = S(x)e −1/ϵ , with the scaling chosen so that the exponent in ( 18) is equal to 0 when x = i.
When written in this form, S(x) is known as the Stokes multiplier.It is constant except within a narrow region surrounding a Stokes curve.In this region, S varies rapidly from zero on one side of the curve to a nonzero value on the other.This rapid variation generates Stokes switching in the solution.
Substituting ( 18) into ( 16) and simplifying gives Recall that the optimal truncation given in ( 14) depends on |x − i|.This suggests that the transformation i(x − i) = re iθ will simplify this expression in a useful way, giving Nopt = r/2ϵ + ω.As in [22], we then fix r and consider only the faster variation that occurs in the angular direction θ.From the chain rule,  The branch points generate Stokes curves, which connect the two points.On the left-hand side of the Stokes curves, there are no exponential contributions.On the right-hand side of the Stokes curves, the solution contains exponentially small oscillations of the form given in (27).
Applying this transformation to (16) and rearranging gives We now finally substitute in the optimal truncation value Nopt = r/2ϵ + ω.Making use of Stirling's formula in the limit that ϵ → 0 gives As promised, this expression is exponentially small, except in the neighbourhood of θ = 0, where the exponential term is algebraic.Hence, S varies rapidly in this region, indicating that θ = 0 is a Stokes curve.
Recall that i(x − i) = re iθ , so θ = 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 (11), we identify that there is a second Stokes curve extending vertically upwards from the branch point x = −i.This Stokes structure is illustrated in Figure 2.
If we cross the Stokes line along Re(x) = 0, we expect that there will be an exponentially small jump in the asymptotic solution.

Exponentially Small Contribution
To determine the behaviour that appears as the Stokes curve is crossed, we define a new inner variable By comparing ϑ with the original coordinates, we see that x < 0 corresponds to ϑ → −∞, while x > 0 corresponds to ϑ → ∞.Hence, the jump across the Stokes curve is given by If the value of S changes by (24) as the Stokes curve 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 uexp as the Stokes curve is crossed from left to right, Note that this expression has constant amplitude.By comparing this with the boundary conditions in (2), we see that there cannot be any finite-amplitude oscillations present as x → −∞.This indicates that uexp = 0 on the left-hand side of the Stokes curve (x < 0), and that the exponentially small behaviour appears as the Stokes curve is crossed from left to right.On the right-hand side of the Stokes curve, the exponentially small solution contribution uexp must be This behaviour is illustrated in Figure 2.

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

AAA Approximation for u 0
In order to simulate a problem where we only possess a numerical description of the leading-order solution u0 from (9), we first define a set of points xj separated by some ∆x.We then evaluate (9) at each point to obtain u0(xj).We will use this set of points xj and function values u0(xj) as the basis for a AAA rational approximation.
We apply the AAA algorithm to obtain a rational approximation û0(x), given by For example, if we choose xj to be evenly distributed in the interval [−4, 4] with ∆x = 0.1, applying the AAA algorithm with an error tolerance of 10 −12 gives a rational function with 15 poles.The poles and residues are are shown in Table 1.We have given each pair of poles a designation, so that we may refer to them later.
Note that there is a pole located on the real axis at x ≈ −6.066.The AAA algorithm sometimes generates poles on the real axis that lie outside of the approximation interval.These poles have no practical effect on the solution, and we will omit their contributions from the sum (28).The remaining poles appear as complex conjugate pairs, as the solution is real when x is real.
The pole pairs are not restricted to the imaginary axis, which is where we defined the branch cuts to be in the true solution u0 (9).In fact, the AAA algorithm places the poles along a curve, representing the branch cut, which we are not free to prescribe.This is a common feature of rational approximation error tolerance of 10 −12 .The poles in û0 accumulate at x = ±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.
methods, and was studied in depth for Padé approximation in [24], using a quantity known as the "condenser capacity" in the approximation.More recently, an explanation of the branch curve selection for AAA that built on these ideas was presented in [28].
It is possible to adjust the rational approximation algorithm, for example by using the minimax algorithm from the chebfun package, to obtain an approximation which forces the poles to align along the imaginary axis.
In Figure 3, we present the real and imaginary parts of the true leading-order solution u0, and the approximated leading-order solution û0.Note that the true solution possesses vertical branch cuts originating at ±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 p1 as the pole nearest to x = i.
In Figure 4, we present the error between the two on a logarithmic plot.This confirms that the approximation error is small except in a region surrounding the branch cut/poles, where the error becomes significant.Given that the approximation is highly accurate everywhere on the sample domain except near the branch cut, we hope that the exponential asymptotic analysis will produce equivalently accurate results.

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

Stokes Curves and Exponentially Small Terms
To determine the location of the Stokes curves, and the quantity that appears as they are crossed, we proceed with an almost identical analysis to that of Sections 2.2 and 2.3, with the only significant difference being the form of the series terms in (31).
This analysis reveals that there are exponentially small contributions in the solution associated with each pole.These contributions appear across Stokes curves that extends vertically from the corresponding pole and intersect the real axis.This behaviour is shown for the example problem in Figure 5.We will later find that the most significant exponential contributions in this example appear across the Stokes curves generated by pole pairs 1, 2, and 3.
Using an essentially identical set of steps to Section 2.3, we can determine the form of the exponentially small remainder.We will determine the contribution that appears across the Stokes curve generated by a pole in the upper-half plane located at at x = pr.The Stokes curve extends vertically downwards from the pole along Re(x) = Re(pr) and intersects the real axis at x = Re(pr).
We denote the optimally-truncated remainder as RN .On the left-hand side of the Stokes curve, we have RN = 0. On the right-hand side of the Stokes curve, we find that RN ∼ 2πar ϵ e −i(x−pr )/ϵ as ϵ → 0.
Once all of the Stokes curves have been crossed from left to right, we find that the combined exponentially small contribution, which we denote as ûexp, is therefore given by ûexp ∼ m r=0 2πar ϵ e −i(x−pr )/ϵ as ϵ → 0. ( For example, for the set of poles presented in Table 1, we find that all of the exponentially small  (29), which uses û0 as the leading-order solution.The solution contains 7 pairs of poles, with locations given in Table 1.Each pair of poles generates Stokes curves that extend vertically from the poles, intersecting the real axis.As each Stokes curve is crossed from left to right, an exponentially small asymptotic contribution appears in the solution.Note that the poles accumulate near the true branch points of u 0 at x = ±i.We will later find that the largest exponential contributions arise from the poles that are nearest to x = ±i.contributions are present in the solution on the right-hand side of the Stokes curve that intersects the real line at x = −0.0015.
Note that the exponential contributions in (33) associated with upper-half plane poles at x = pr decay exponentially as Im(pr) increases (with the complex conjugate contributions exhibiting corresponding exponential decay).This suggests that the poles nearest to the real axis will produce the largest exponential contributions.

Comparison
We have now calculated the exponentially small contributions that appear in the asymptotic solution to the differential equation ( 1) with boundary condition (2), given by uexp in (27).We have also calculated the exponentially small contributions that appear in the asymptotic solution to the approximate problem (29) with the same boundary condition, given by ûexp in (33).The remaining question is whether the approximate exponential contributions are able to accurately approximate the true exponential contributions.  1) and ( 29) for ϵ = 0.2.The true exponential contribution u exp is shown as a black dashed curve.The approximate exponential contribution û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 û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 = ±i (ie.pole pairs 1, 2, and 3), with the amplitude of the contributions decaying as the distance of the pair from x = ±i increases.
Note that ( 27) and ( 33) cannot be identical as ϵ → 0, due to the following differences: • The algebraic prefactor of uexp is proportional to ϵ −1/2 , while the algebraic prefactor of the approximate exponential behaviour ûexp is proportional to ϵ −1 .This difference in exponent will be generic: the AAA rational function must have simple poles that generate algebraic prefactors that are proportional to ϵ −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 u0 (9) and û0 (30).The branch points in u0 that produce (27) are located at x = ±i, while the poles in û0 that produce (33) are located at x = pr, where |Im(pr)| > 1.
Hence the two expressions have different exponential decay as ϵ → 0, which can be seen by directly comparing the two expressions.
These differences appear to suggest that the two expressions uexp and û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 ϵ → 0, as all of the exponential terms in (33) are smaller than the exponential terms in (27) in the asymptotic limit (even though the algebraic prefactor is larger).
Despite this, we find that ûexp is able to accurately approximate uexp for some range of ϵ.We compare the two contributions for our sample problem with ϵ = 0.2 over a segment of the positive real axis in Figure 6.The contributions uexp and ûexp appear indistinguishable in this figure.
In addition to the exponentially small contributions, Figure 6 also shows the contribution to ûexp from each pole pair, and reveals that the largest contributions come from pole pairs 1-3.This indicates that pairs which are nearest to the branch points at x = ±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 |Im(pr)| increases.expressions have different behaviour as ϵ → 0, so it is impossible for the approximation to remain accurate indefinitely as ϵ is decreased.This is apparent in the figure as each approximation is accurate until some lower threshold value of ϵ 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 ϵ.We also identify the value of ϵ equal to |p 1 − i| on each curve.This value of ϵ corresponds to a roughly constant value of the error in each curve, suggesting that the approximation error of the exponential terms is related the quantity |p 1 − i|.
In Figure 7, we compare the amplitude of uexp and ûexp over a range of ϵ 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 ϵ is reached from above; below this threshold value the ratio tends to 0. This is consistent with the algebraic expressions for the two terms, as the exponential contributions in ûexp are smaller than the contributions in uexp as ϵ → 0.
Reducing the error tolerance of the AAA algorithm has the effect of decreasing the value of ϵ 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 ϵ.
It is apparent from Figure 7 that higher AAA tolerances produce approximations with more poles, that are accurate for smaller values of ϵ.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 = p1), or |p1 − i|.We compare this with the relative error in the amplitude, measured by Relative Error = 1 − Amplitude of ûexp Amplitude of uexp .
Table 2 suggests that, no matter what AAA tolerance is chosen, the relative error remains relatively similar at ϵ = |p1 − i|.The points ϵ = |p1 − i| are also marked in Figure 7.This figure shows that for ϵ ≪ |p1 − i|, the approximated amplitudes are inaccurate while for ϵ ≫ |p1 − i|, the approximate amplitudes are accurate.We conjecture that if there is a singularity in the true leading-order behaviour at x = xs and the nearest pole to this point in the rational approximation is at x = p1, the exponential Obviously, this method is most valuable for problems in which xs is not known beforehand, so |p1 −xs| is not known.It is possible to estimate xs by carefully examining the accumulation rate of poles in the complex plane [29].We intend to study this in more detail in future work, and to explain the dependence of the threshold value on |p1 − xs|.

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 ϵ is not too small.The next natural question is whether we can extend these methods to nonlinear ordinary differential equations, using the exponential asymptotic method of [7].
This is a challenging problem, because we can no longer add the pole contributions in ûexp independently.From [30], we know that singularities in nonlinear problems that are near to each other (ie. within a neighbourhood whose width is proportional to a particular power of ϵ) can interact, and that the resultant asymptotic behaviour changes due to these interactions.This is likely to occur for leadingorder solutions generated using rational approximation methods due to the accumulation of poles near the true singular point.In future work, we will determine the asymptotic corrections to uexp due to pole interactions in nonlinear problems.
In many problems it is possible to determine the strength of the singular points in u0 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.The AAA-approximated exponential terms were obtained using the leading-order from (38).In order to show that this method captures nonlinear effects, we show the expression on the entire real axis, but note that this contribution is only actually present in the asymptotic solution for x > 0. The two contributions are visually similar but not completely identical, due to nonlinear effects caused by the subdominant branch points at x = ±i in (37).
Consider the following simple example: with the boundary conditions in (2).The constant 81/1024 was chosen to make the exponential asymptotic analysis particularly convenient.The true leading-order behaviour is given by We can sample this leading-order solution on a set of points xj to obtain data points u0(xj), as in Section 3.1, and use this as the basis for an approximation û0.Using this, we can compute the exponential contributions that appear across each Stokes curve in the solution independently and add them together.
If we do so, however, we obtain predictions of the exponentially small behaviour that are inaccurate, because the analysis does not take into account nonlinear interactions between pole contributions.
We instead sample the leading-order solution on xj and then square the output, to obtain u0(xj) 2 , which is equivalent to taking samples of We can obtain a AAA rational approximation for u 2 0 based on this data, which we denote u 2 0 .Finally, we can take the square root of this rational approximation to obtain This expression can then be used as the leading-order solution for an exponential asymptotic analysis using the methods from [7].In Figure 8, we present the true exponentially small terms uexp and the from (36) with the six nearest pole pairs in the approximation for u 2 0 from (37).The poles in û0 have residues with similar magnitude, as the singularities in u 0 are branch points at x = ±i.In ( u 2 0 ) 1/2 , the poles nearest to x = ±i have a larger residue than the remaining poles, as it approximating the contribution from the simple poles in (37).The remaining poles mimic the behaviour of the subdominant branch cut.approximated exponentially small terms ûexp obtained using this method for ϵ = 0.1.The rational approximation used was generated using sample points on x ∈ [−10, 10] with ∆x = 0.1.We present the exponential contribution on x ∈ [−20 , 20] to demonstrate that this method accurately captures nonlinear wave effects, but we actually expect this contribution to appear as the Stokes curve intersecting x = 0 is crossed from left to right, and therefore only be present in the asymptotic solution for x > 0. The qualitative match seen in this figure shows that the exponentially-small terms in (35) can be accurately approximated using this method.
The accuracy of this approximation compared to naively approximating the exponential terms based on û0 can be understood by comparing the poles and residues of the AAA approximation of u0 and u 2 0 , using the parameters from Figure 8.These are presented in Table 3.
In u0, the only true singularities are the branch points at x = ±i.In Table 3, we see that the residue of each pole of û0 is similar in magnitude, as these poles approximate the effect of the branch cuts.In u 2 0 , there are simple poles at x = ±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 ±i have a larger residue, with the other poles having a smaller residue.This is because the poles nearest to ±i in the AAA approximation are reproducing the simple pole behaviour in the true expression for u 2 0 , and the remaining poles are approximating the subdominant branch cut behaviour.
Because the strongest singularities in u 2 0 from (37) are simple poles, the AAA approximation can be used as the basis for an accurate exponential asymptotic analysis.However, this expression does still contain branch cuts originating at x = ±i.This produces a string of additional poles in the AAA approximation.Hence, any higher-order corrections to the exponential terms, such as those studied in [8,23], 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.

Conclusions and Discussion
In this study, we applied exponential asymptotic methods from [22] to study Stokes' phenomenon in a linear ordinary differential equation in the small-ϵ limit.The leading-order solution u0 (9) contains two branch points, located at x = ±i.These branch points produce Stokes curves, and oscillating exponentially small terms appear in the solution as these Stokes curves are crossed at x = 0 on the real axis.We calculated the form of these exponentially small oscillations, uexp (27).
We then repeated this process using a rational approximation for the leading-order behaviour.Instead of using the exact leading-order solution, we sampled the leading-order behaviour on a discrete set of points and used the AAA algorithm to produce a rational approximation for the leading-order solution based on the sampled data, û0 (28).We then performed an exponential asymptotic analysis and found that each pair of poles in the solution produced Stokes curves, each of which generated exponentially small oscillations.Taking the sum of these oscillations produced the total exponentially small behaviour ûexp (33).
Finally, we compared uexp and ûexp, and found that ûexp is able to accurately approximate uexp for nonzero values of ϵ, despite having different asymptotic behaviour in the limit that ϵ → 0. If ϵ 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 ϵ beyond which the approximation loses accuracy.
Empirical tests suggest that this threshold value of ϵ is proportional to the distance between the branch point in u0 and the nearest pole in û0.This is a measure of how accurately the rational approximation predicts the true singularity location.While the true branch point (or other singularity) in u0 cannot typically be calculated in practice, it is possible to estimate the true location by studying the rate at which the poles in û0 accumulate [29].
There are two promising avenues available for extending this method to study nonlinear ordinary differential equations.The first is to determine the asymptotic corrections to the approximated exponentials that are caused by nonlinear interactions between nearby simple poles.The second is to use asymptotic arguments to determine the strength of singularities in the leading-order behaviour, and apply the AAA algorithm to a mapped version of the data that contains simple poles instead of other singularities.We briefly showed in Figure 8 that this method can be used successfully for a nonlinear ordinary differential equation, and plan to explore this and related problems in future work.This approach has connections to recent work [10,11] where Padé approximants are generated from numerical truncated power series data (as is typically available in Borel plane analyses [13]).These approximations generate accumulations of simple poles near branch points in the original Borel-transformed function.The authors use conformal mappings to convert these into poles, which may be more easily studied.It is likely that a variant of the the conformal map techniques developed in [11] could be applied to the AAA approximation to u0(z) developed in the present work, as well as more complicated leading-order behaviour.
Finally, it is important to determine how robust this method is to noisy input data.The effect of noise on numerical rational approximation has been studied in [12], the authors study the effects of noise on conformal maps generated using Padé approximation, and in [31], in which the authors study the effects of noise on rational approximations generated using the AAA algorithm.It would be valuable to use similar methods to study the effect of noise on exponential asymptotic analyses.

.
The weight vector w = (w1, w2, . . ., wm) T is chosen to minimize ||f d − n|| X (m) subject to ||w||m = 1, where || • || X (m) is the discrete 2-norm over X (m) and ||•||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 (xj)|), the algorithm terminates.If the algorithm does not terminate, the value of xm+1 is determined for the next iteration.It is found by choosing the value of xm+1 ∈ X (m) that maximizes the quantity

Figure 1 :
Figure 1: Magnitude of series terms from (8) evaluated at x = 0 for ϵ = 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

Figure 1
depicts ϵ n un at x = 0 with ϵ = 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 2 :
Figure 2: Stokes' phenomenon in the exact solution to (1) with boundary conditions (2).The leading-order solution contains branch points at x = ±i, with the branch cuts extending vertically away from the real axis.

Figure 3 :
Figure3: The real and imaginary parts of the true leading-order solution u 0 and the approximated leadingorder solution û0 described in Table1.The two expressions are visually indistinguishable except on the imaginary axis.The function u 0 contains vertical branch cuts.The function û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.

Figure 4 :
Figure 4: This figure shows the error |u 0 − û0 |, using u 0 and û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.

Figure 5 :
Figure 5: Stokes structure in the solution to (29), which uses û0 as the leading-order solution.The solution

Figure 7 :
Figure 7: Ratio of the amplitudes of the approximate exponential contribution ûexp and the true exponential contribution u exp as ϵ is varied.If ûexp is an accurate approximation of u exp , this ratio is close to 1.The two

Table 1 :
Pole Index: r Pole Pair Name Pole locations: p r Poles and residues in the AAA approximation of u 0 , which contains branch points at x = ±i.The approximation û0 was produced by sampling u 0 on the domain x ∈ [−4, 4] at intervals of ∆x = 0.1 with an

Table 2 :
This table shows the distance between the true branch point in u 0 at x = i and the nearest pole in û0 , denoted as p 1 , for the different AAA approximations presented in Figure7.As the tolerance of the approximation increases, the distance |p 1 − i| decreases.The relative approximation error is evaluated for ϵ = |p 1 − i|, and the values are roughly constant.This 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 3 :
Comparison of the six pole pairs nearest to x = ±i contained in a rational approximation for u 0