Computation of Asymptotic Expansions of Turning Point Problems via Cauchy’s Integral Formula: Bessel Functions

Linear second-order differential equations having a large real parameter and turning point in the complex plane are considered. Classical asymptotic expansions for solutions involve the Airy function and its derivative, along with two infinite series, the coefficients of which are usually difficult to compute. By considering the series as asymptotic expansions for two explicitly defined analytic functions, Cauchy’s integral formula is employed to compute the coefficient functions to a high order of accuracy. The method employs a certain exponential form of Liouville–Green expansions for solutions of the differential equation, as well as for the Airy function. We illustrate the use of the method with the high accuracy computation of Airy-type expansions of Bessel functions of complex argument.

the use of the method with the high accuracy computation of Airy-type expansions of

Introduction
In this paper, we study linear second-order differential equations having a simple turning point. Specifically, we consider the differential equation where u is positive and large, f (z) has a simple zero (turning point) at z = z 0 (say), and f (z) and g(z) are analytic in an unbounded domain containing the turning point. This is a classical problem, with applications to numerous special functions. To obtain asymptotic solutions, the Liouville transformation is applied, where either sign in front of the integral can be chosen. As a result, we transform (1.1) to the form where (1.4) The lower integration limit in (1.2) ensures that the turning point z = z 0 of (1.1) is mapped to the turning point ζ = 0 of (1.3). Throughout this paper, we shall assume that this turning point is bounded away from any other turning points or singularities of (1.1); equivalently, ψ(ζ ) is analytic for 0 ≤ |ζ | < R for some positive R which is independent of u. When the turning point z 0 is real and f (z) is real on a real interval around z 0 , the sign in (1.2) is usually chosen in such a way that the new variable ζ is real when z is real and in a neighborhood of the turning point.
In general, these coefficients are difficult to compute, primarily due to the requirement of repeated integrations. They also show cancellations near the turning point. For complex ζ close to 0, one can compute these coefficients in a numerically stable way by considering power series expansions for the coefficients, as done in [1,2], where they are expanded in powers of ω = √ 1 − z 2 . For other computational approaches to compute the coefficients, and in particular for real values of ζ , see [12].
The purpose of this paper is to provide a more simple means of computing a large number of these terms. We shall employ Cauchy's integral formula to do so, and our results will be valid for real and complex ζ lying in a bounded (but not necessarily small) domain containing the turning point ζ = 0. Our approach can potentially be extended to other situations, including the cases of simple poles [9,Chap. 12], and coalescing turning points.
We illustrate the use of the method with the high accuracy computation of Airy-type expansions of Bessel functions of complex argument.

General Method
We first present Liouville-Green expansions for solutions of (1.1), a certain form of which will be required for our method. These only involve elementary (exponential) functions but are not valid at the turning point. The appropriate Liouville-Green transformation is given by [9,Chap. 10,§2]; namely, ξ = 2 3 ζ 3/2 , V = f 1/4 (z)w. (2.1) With these, equation (1.1) is transformed to equation where 3) The branch for the first of (2.1) will be dependent on the solutions under consideration, as described below. Note that as ζ completes one circuit about the turning point ζ = 0, the variable ξ correspondingly crosses more than one Riemann sheet.
It turns out that solutions where asymptotic expansions appear inside exponentials are more convenient for our purposes. Specifically, from [9, Chap. 10, Ex. 2.1], we have solutions In these, the coefficients are given by and Primes are derivatives with respect to ξ . The integration constants in (2.5) will be discussed below, and we find that for the odd coefficients E 2 j+1 (ξ ) ( j = 0, 1, 2, . . .), they must be suitably chosen. Explicit error bounds for ε ± n (u, ξ), which verify the asymptotic validity of the expansions (2.4), are given in [5]. In particular, for arbitrary δ > 0, under certain conditions on ψ(ξ ), we have that ε ± n (u, ξ) = e ±uξ O u −n as ξ → ∞ in certain domains Ξ ± as described in [9,Chap. 10,§3]. These are the same as those for the corresponding asymptotic solutions of the more common form It is the relation (2.7) that is the reason why the expansions (2.4) are numerically advantageous: the coefficients F s (ξ ) can all be determined explicitly without resorting to integration. Furthermore, from (2.5) we observe that only one integration is required (numerical or explicit) to evaluate each E s (ξ ), as opposed to repeated integrals for computing the coefficients A s (ξ ) in (2.8).
Remarkably, it turns out that integration is not required to evaluate the even terms E 2 j (ξ ) ( j = 1, 2, 3, . . .). To see this, consider the Wronskian of the solutions V ± n (u, ξ) given by (2.4). Since this is a constant (by Abel's theorem), we infer that which, on taking logarithms, yields We then asymptotically expand the RHS of this relation in inverse powers of u 2 , and equate the coefficients of both sides. As a result, we find that (to within an arbitrary additive constant in each instance) (2.10) and so on. In particular, the even coefficients E 2 j (ξ ) ( j = 1, 2, 3, . . .) are explicitly given in terms of F 2k+1 (ξ ) (k = 0, 1, 2, . . . , j − 1) (which in turn are given by (2.6) and (2.7)). At this stage, we consider Liouville-Green solutions of (1.3). Comparing (1.2), (2.1), and (2.4), we obtain three asymptotic solutions W j (u, ζ ) ( j = 0, ±1) of (1.3) which are recessive (respectively) for ζ ∈ S j , possessing the asymptotic expansions and as u → ∞, uniformly for ζ ∈ Ω j (say). For j = 0, ±1, the branches are such that Re e −i jπ ζ 3/2 ≥ 0 for ζ ∈ S j , and Re e −i jπ ζ 3/2 ≤ 0 for ζ / ∈ S j . As is typically the case in practice, for each j we assume that Ω j ∩ S j is unbounded.
We remark that, on account of the analyticity of ψ(ζ ) in the disk 0 ≤ |ζ | < R, the expansions (2.11) and (2.12) certainly hold in the bounded sector δ ≤ |ζ | < R, |arg(ζ e −2πi j/3 )| ≤ π − δ, where here and elsewhere δ denotes an arbitrary small positive constant. We also note that the recessive property at ζ = ∞ in Ω j ∩S j uniquely defines W j (u, ζ ) up to a multiplicative constant. Indeed, we have that W j (u, ζ ) = c 2n+1, j W 2n+1, j (u, ζ ) for some constants c 2n+1, j , although we shall not use these relations. Now, since no two from these three solutions are linearly dependent, we can assume they satisfy a connection formula of the form for certain constants λ −1 and λ 1 (which may of course depend on u). The factor i is for convenience. We note that each W j (u, ζ ) ( j = 0, ±1), being a solution of (1.1), is analytic in a neighborhood of the turning point. Based on (1.5), and following [4,12], we thus can define functions A(u, z) and B(u, z), analytic at z = z 0 (ζ = 0), implicitly by the pair of equations 1 2π 1/2 u 1/6 W 0 (u, ζ ) = Ai 0 u 2/3 ζ A(u, z) + Ai 0 u 2/3 ζ B(u, z) (2.14) and e πi/6 λ 1 2π 1/2 u 1/6 W 1 (u, ζ ) = Ai 1 u 2/3 ζ A(u, z) + Ai 1 u 2/3 ζ B(u, z), (2.15) where (as shown below) the multiplicative constants on the LHS of both equations have been chosen to yield the appropriate behavior of A(u, z) and B(u, z) as u → ∞. We remark that for computational purposes, it is more convenient to consider A(u, z) and B(u, z) as functions of z, although in deriving their asymptotic expansions, we shall regard them as functions of ζ as necessary.
Next, from the connection formula (2.13) and the corresponding well-known connection formula for the Airy functions, we derive from (2.14) and (2.15) the following Airy function representation for the solution which is recessive in S −1 : We shall show, using (2.14), (2.15), and (2.16), that A(u, z) and B(u, z) are slowly varying relative to the fast variation of the Airy functions in a full neighborhood of the turning point. Specifically, on referring to (1.5), A(u, z) and B(u, z) will admit the following asymptotic expansions as u → ∞: uniformly with respect to ζ lying in a certain unbounded domain that contains the disk 0 ≤ |ζ | < R.
We now use (A.3) and (A.4) (see Appendix), along with the corresponding expansions (which are valid for arg ζ e 2πi/3 ≤ π − δ): and B(u, z) ∼ 1 These asymptotic expansions certainly hold for (at least) −π + δ ≤ arg(ζ ) ≤ 1 3 π − δ, 0 < |ζ | < R. Similar slowly-varying expansions can be obtained for all other values of arg(ζ ) near the turning point by solving for A(u, z) and B(u, z) from a suitably chosen pair from the three equations (2.14), (2.15), and (2.16). Each of these expansions can be expressed in the same form as above, except that the coefficients E s (ξ ) may differ by the integration constants in (2.5).
In fact, if these integration constants are arbitrarily chosen, we find that the RHS of (2.20) and (2.21) generally each has a branch point at ζ = 0 (z = z 0 ). Hence, being multi-valued (in fact unbounded) for small ζ , the expansions are generally only valid for restricted arg(ζ ) near the turning point. It is only for specific integration constants in (2.5), at least for the odd coefficients E 2 j+1 (ξ ) ( j = 0, 1, 2 . . .), that the expansions are single-valued, and hence by a continuity argument, (2.20) and (2.21) hold for 0 < |ζ | < R, with arg(ζ ) unrestricted. This essentially corresponds to the choice of the lower integration constant in (1.6), which (as remarked above) ensures that each coefficient B s (ζ ) is analytic at ζ = 0.
Recalling that ξ = 2 3 ζ 3/2 , it is straightforward to show that the expansions (2.20) and (2.21 are meromorphic: here and throughout, meromorphic means with respect to ζ , and at ζ = 0. If we assume for the moment that this is true for ζ 1/2 E 2 j+1 (ξ ) and E 2 j (ξ ), then on re-expanding (2.20) and (2.21) into the forms (2.17), we deduce that the coefficients A s (ζ ) and B s (ζ ) in the latter expansions must be single-valued. By a uniqueness argument, these coefficients must satisfy (1.6) and (1.7), and moreover, single-valuedness means the lower integration limit in the former must indeed be 0. But in this case, we know that A s (ζ ) and B s (ζ ) are actually analytic at ζ = 0. We deduce that if ζ 1/2 E 2 j+1 (ξ ) and E 2 j (ξ ) are meromorphic, then (2.20) and (2.21) have a removable singularity when they are re-expanded into the forms (2.17). With an appropriate choice of integration constant for each E 2 j+1 (ξ ) ( j = 0, 1, 2, . . .), we now show that this is indeed so.
To establish the desired meromorphicity of ζ 1/2 E 2 j+1 (ξ ) and E 2 j (ξ ), let us consider these separately. First, from (2.5), we see that ζ 1/2 E 2 j+1 (ξ ) is meromorphic if and only if the Laurent-type expansion satisfies α 2 j+1 = 0. Clearly, the integration constants in (2.5) can be selected in order for this to be so, and we assume this from now on.
Next consider the even terms. Again from (2.5), we observe that must either be meromorphic or have a logarithmic singularity at ζ = 0. However, the latter possibility can immediately be discarded by referring to (2.9) along with the meromorphicity of F 2 j+1 (ξ ). We note that the integration constants in (2.24) do not affect the meromorphicity of E 2 j (ξ ), and hence they can be arbitrarily chosen.
In summary, we have established that α 2 j+1 = 0 in (2.23) is a necessary and sufficient condition for the expansions (2.20) and (2.21) to be valid in a neighborhood of the turning point with arg(ζ ) unrestricted.
In practice, there are several ways of ensuring the correct choice of E 2 j+1 (ξ ). If each of these functions can be explicitly determined from (2.5), then a symbolic algebra system could be used to determine the expansion (2.23), and if α 2 j+1 = 0, one can simply subtract this constant from the function given by the anti-derivative derived from (2.5).
If explicit integration of (2.5) is not possible, and quadrature is required, then from (2.23), we observe that Since the coefficients will be computed numerically on a loop surrounding ζ = 0, the values on the RHS of (2.25) can also be computed, provided ξ and ξ * correspond to the initial and terminal points of the loop. Hence, as in the case of the explicitly known coefficients, these numerically computed values can be subtracted from the initial values of E 2 j+1 (ξ ) evaluated from (2.5).
In the common situation where ψ(ζ ) is real in some real interval (−a, a), where a > 0, this numerical method can be simplified somewhat. Specifically, if we choose the lower limit of integration in (2.5) to be at a real point on our path of integration, ξ = ξ(ζ 0 ) say, where 0 < ζ 0 < a, then we similarly find from (2.23) that and we can then proceed as above.
We summarize the principal result as follows.
Theorem 2.1 For the differential equation

)-(2.7), where the integration constants for the odd coefficients in (2.5) must be selected so that
can alternatively be evaluated by expanding the RHS of (2.9) in inverse powers of u 2 , and equating the coefficients of both sides. This avoids integration for evaluating these terms. We also note that the expansions (2.20) and (2.21) are valid in the same domains as those for the corresponding asymptotic solutions of [9, Chap. 11, Theorem 9.1] and can be unbounded provided f (z) and g(z) have the appropriate behavior at ∞ (see [9, Chap. 11, Sect.

9.3]).
Our focus is to utilize the expansions (2.20) and (2.21) to efficiently compute A(u, z) and u 4/3 B(u, z) to O u −2m for some prescribed m. As we shall see in the application to Bessel functions below, in general it may be more convenient to consider scaled functions where φ(u, z) is some suitably chosen function that is analytic in a domain Ω (say) in which the expansion (2.20) is valid. Our aim is to employ the Cauchy integrals where L is a positively orientated closed loop lying in Ω and surrounding t = z and t = z 0 , and therefore (2.20) and (2.21) can be inserted in the respective integrands. As we showed above, the asymptotic expansions (2.20) and (2.21) are in theory valid close to the turning point ξ = ζ = 0 (z = z 0 ) as u → ∞. However, for fixed u, the terms in the series become unbounded as ξ → 0, and therefore these series are numerically unsatisfactory near the turning point. It is for this reason that we use the Cauchy integral formulas (2.31) for their numerical approximation. At this point, it is important to realize that we will need to compute the coefficients E s (ξ ) in (2.20) and (2.21) at specific fixed points on L. Since this computation is done once and for all, the computational efficiency in the computation of these coefficients is not so important. And, as noted earlier, if the integrals in (2.5) cannot be explicitly evaluated, only one numerical evaluation of an integral is required for each odd coefficient.
For our purposes, it will be sufficient to consider a circular path of integration with center at z c , not necessarily with z c = z 0 , but which in any case must contain the turning point z 0 . For computing the integral, we therefore use the parametrization t (θ ) = z c + Re iθ , θ ∈ [0, 2π ], and we have and likewise for B(u, z). The function F(θ ) is periodic, and we are integrating over one period. Since the function is infinitely differentiable as a function of θ , we can expect that the trapezoidal rule will give good convergence [8,Thm. 5.6]. We write and therefore Numerical experiments show that N does not need to be large and that, for instance in our application to Bessel's equation (see § 5.2), N = 500 is enough for more than 15-digits accuracy in a wide region around the turning point (namely for z c = 2, R = 1.8); the number can be reduced for smaller regions, and for z c = 1, R = 0.5, 15 digits accuracy is reached with N = 150. This is consistent with the expected good performance of the trapezoidal rule.
For the derivatives of solutions w j (u, z) (say) of (1.1), suppose scaled coefficient functions A(u, z) and B(u, z) are defined by for some appropriate connection coefficients λ ±1 . We can then proceed similarly as before, by defining corresponding coefficients C(u, z) and D(u, z) by Then, solving for these coefficients, we obtain the same expressions as for the previous coefficients, but with w j (u, z) substituted by their derivatives. Then, to obtain similar expansions to (2.20) and (2.21), we would also need the Liouville-Green expansions for the derivatives of these functions, which can be obtained by differentiation of the corresponding expansions for the functions themselves. A preferable way to compute the new coefficients is to differentiate (2.34). After differentiating with respect to z, using the Airy differential equation to eliminate the second derivative of the Airy functions, then solving and comparing to (2.35), we obtain the relations To compute the derivatives of the coefficients, we can also use Cauchy's integral formula to write (2.38) Then, from (2.30), (2.31), and (2.36)-(2.38), we obtain and We then use the computed values of A(u, t) and B(u, t) on L and proceed as above. Again, the trapezoidal rule is a good choice; in fact, it is in some sense optimal [3]. We remark that it is desirable that (2.36) and (2.37) be numerically stable, in the sense that, for large u, there is no cancellation in leading order terms in the two terms in either representation. The choice of scaling function φ(u, z) in (2.30) should be such that this is indeed the case.

Bessel's Equation: Preliminary Transformations
We illustrate the new technique using the Airy function asymptotic expansions for Bessel functions. The first step in doing so is to apply the Liouville transformations described in § §1 and 2 to Bessel's equation. To this end, we first note that functions Here ν plays the role of our large positive parameter u, and z is complex. From (1.2), and taking the negative sign, let and The transformed variable ζ is real for real z ∈ (0, 1) (ζ ∈ (0, +∞)), and ζ(z) can be defined by analytic continuation in the whole complex plane cut along the negative real axis. This transformation, and its correspondence in the z-plane, is depicted in Fig. 2.
Observe that if we are assuming that the principal values are taken, the expression (3.1) should not be used in all the complex plane, and, for instance, for real z > 1, we would have complex values of ζ , when ζ should be real for real z > 0. This problem is solved by taking the following: which in fact is an alternative when z > 1 (and in fact a larger region). We can also write this in terms of logarithms as We get (1.3), with u replaced by ν, and From [9, Chap. 11, §10], we identify the asymptotic solutions by (3.5) and
In the general case, the coefficientsÊ s (z) can be computed numerically by quadrature via (3.7). But for Bessel's equation, it turns out that the coefficients can be explicitly computed, and in particular, they have the expression where P s (x) are polynomials of degree s in x.
Before establishing this, we note for the odd terms that where the term in the square brackets is meromorphic at z = 1. Hence, on expanding around the turning point ζ = 0 (z = 1), we deduce that (2.23) holds with α 2 j+1 = 0, as desired.
In order to establish (3.10), we first show that the coefficientsF s (z) have the form with Q s (t) a polynomial of degree not larger than s in t.
We first observe thatF 1 (z) has this form. Next, we substitute (3.11) in (3.9), and then, for (3.11) and (3.9) to hold, we have that the polynomials must satisfy (3.12) which, starting from Q 1 (t) = (1/2 + t/8), shows that Q n is of degree not larger than n and that (3.11) holds with polynomials given by (3.12). Now, the coefficientsÊ s can be computed from where the starting point a is chosen depending on the solution to be matched. These integrals can be explicitly computed. We have and because Q s is a polynomial of degree not larger than s, we are left with integrals of the form with k < a −1 and R k a polynomial of degree not larger than k which can be computed by differentiating (3.13). Because k < a − 1, we have H (∞) = 0. The polynomials P s in (3.10) have the properties: where C 2s+1 are the coefficients in the Stirling asymptotic series We found that the property P 2s+1 (0) = C 2s+1 holds by comparing (3.4) with a similar expansion for J ν (νz) but matching the solution at z = 0, using that J ν (νz) ∼ (νz/2) ν / (ν + 1) as z → 0. The first few polynomials are x(x + 4), x 32 + 288x + 232x 2 + 13x 3 .
Next, using as z → ∞, we match the solutions that are recessive in the upper half z-plane, yielding (3.14) as ν → ∞, uniformly for z lying in a domain that contains the first quadrant (excluding a neighborhood of the turning point z = 1).

(3.15)
We remark that the expansions (3.14) and (3.15) are a reformulation of the Debye expansions [11, §10.19(ii)] for z ∈ (0, 1), which also holds for certain complex values. Debye expansions valid for 1 < z < ∞ are also given in this reference. Comparing these two representations with (3.4) and (3.5), we perceive that in (2.30), the scaling function here is given by

(4.4)
A similar relation for −2Y ν (νz) is obtained by replacing Ai and its derivative by Bi and its derivative in (4.4).
As described earlier, the idea for computing the slowly-varying coefficients A(ν, z) and B(ν, z) in a region containing the turning point z = 1 1 is to invoke Cauchy's integral formula (2.31). In this, the integration is taken along a positively oriented closed loop containing z in its interior, but not the origin; additionally, it should not cross the negative real axis in the z-plane; in this way, we can guarantee that in the ζ plane, we stay away from the shaded region in Fig. 2 (right) and therefore that all the functions appearing in the integration are analytic in a domain containing the path of integration and its interior (and therefore Cauchy's integral formula holds for this integration path). For our purposes, it will be sufficient to consider a circular path of integration in (2.32), which can be centered or not at z = 1, but which in any case must contain z = 1, but not z = 0.
From (2.21), (2.30), and (4.3), we have the following expansion, which we use to compute B(ν, z) on the path of integration: where, formally, where a 1 = a 2 = 5 72 , and subsequent terms given by (A.2). We observe that the cancellation of the coefficient as ν → ∞ is located in the sinh term. In order to avoid loss of accuracy, it is better to write where the function sinhc(x) = sinh(x)/x can be computed for small x with the Maclaurin series Several properties that could be expected for the coefficient B(ν, z) can be seen to hold explicitly from (4.5): 1. It is real for real values of z.
For the function A(ν, z), we have the expected properties from (4.6): 1. It is real for real values of z.
For the derivatives H (1) ν (νz) and H (2) ν (νz), we proceed as before, and we express them in the form and Solving (4.8) and (4.9) for C(ν, z) and D(ν, z), and considering the Airy-type expansions for the derivatives [11, 10.20.9], we find that as ν → ∞. From the behavior noted above of A(ν, z) and B(ν, z) for large ν (which also hold for their derivatives), we see that in the dominant coefficient D(ν, z), the first term in (2.37) is the largest, which is indeed O(ν 1/3 ), while for the coefficient C(ν, z), both terms in (2.36) are O(ν −1/3 ), which is the correct order of this coefficient. Therefore, once the cancellation for the coefficient B(ν, z) is avoided, no cancellations occur as ν → ∞, and we expect both (2.36) and (2.37) to be numerically stable.
There are two possible approaches in numerically evaluating the coefficients C(ν, z) and D(ν, z). With A(ν, z) and B(ν, z) computed on the path of integration as described above, the first method is to numerically evaluate the integrals (2.39) and (2.40).
Alternatively, one can use (2.36) and (2.37), with A(ν, z) and B(ν, z), and their derivatives, computed via Cauchy integral formulas. Then, ζ and ζ can be evaluated directly from (3.1), (3.2), and (3.3). There is, however, a possible loss of accuracy as z → 1 in the computation of ζ (z) but it can be easily eliminated. We have and both the numerator and the denominator tend to 0 as z → 1, which implies loss of accuracy. In order to avoid this, we put δ = √ 1 − z 2 , and when |δ| is small, we consider the Maclaurin series for Returning to A(ν, z) and B(ν, z), we end this section describing a more direct but less stable method for their computation (again, on the path of integration of (2.31)). From (4.1) and (4.2), we have the exact representations Now, because we are assuming that a method to compute the Airy functions is available (we are precisely considering expansions in terms of Airy functions), we could compute these coefficients on the Cauchy contour without the need to substitute the Airy functions by their Liouville-Green expansions, as done before. For examples of Fortran implementations of complex Airy functions, see [1,6,7]. As before, the integration path is chosen in such a way that the functions appearing in (4.11) and (4.12) can be computed without recourse to the Airy-type expansion. With respect to the cylinder functions involved in the computation of the coefficients, we use (3.14) and (3.15), along with H (2) In order to compute the coefficients in a numerically stable way along the contour of integration, we need to verify that the two terms are not suffering cancellations, which would happen if the two terms are exponentially large and cancel each other. If both Hankel functions are large, one should be replaced by the J Bessel function. For example, we see that inside the eye-shaped curve (Fig 2), the expression for B(ν, z) in (4.12) is unstable because two exponentially large quantities are subtracting: all terms in this expression are dominant inside the eye-shaped curve. Outside the eye, this expression is in principle stable because in each term a dominant function multiplies a recessive function.
Inside the eye, we can write a satisfactory expression using H This expression not only can be used inside the eye-shaped region, but also for the rest of the half-plane (z) ≥ 0; however, close to the real axis when z > 1, we get higher accuracy from Liouville-Green expansions using (4.12). In our numerical algorithms, we use (4.12) when (z) > 1 and (4.13) in the rest of this half-plane. But, as described before, if we also expand the Airy functions instead of computing them separately, we obtain asymptotic expansions for the coefficient away from the turning point, and switching from one expression to the other for the coefficients is not needed in this case. See Sect. 4.
Similarly for A(ν, z), we have We could do analogous substitutions to get a satisfactory formula when z ≤ 0. However this is not really necessary because the coefficients A(ν, z) and B(ν, z) are real on the real line and, as commented before, analytic in a domain containing the integration path. Therefore, by the Schwarz reflection principle, in this domain, A(ν,z) = A(ν, z) and B(ν,z) = B(ν, z). From a computational point of view, this reduces by one half the complexity of evaluating the Cauchy integral, provided we take a symmetric contour with respect to the real axis, as we will do. This more direct approach has some disadvantages with respect to the expansions in terms of stability. First, we notice that, even when the solutions are chosen adequately, there still remains some cancellation in the B(ν, z) coefficient for large orders.
Additionally, there is also some accuracy degradation in the computation of Airy functions for large arguments due to unavoidable loss of accuracy in the computation of exponentials of large argument, as we next describe.

Numerical Results
We now give several numerical illustrations for the performance of both the Liouville-Green expansions and the approximation around the turning point. For this purpose, we have coded our algorithms in Fortran 90 and compared with the values given by Amos' algorithm (which has typically 13-14 digits accuracy). Amos' algorithm uses a variety of methods for computing Bessel functions, depending on the values of ν and z, and in particular Airy-type expansions close to the turning point. We have also tested our methods using Maple T M . The comparison with Amos' algorithm is used as an exhaustive testbench of the numerical stability of our approximations in fixed precision arithmetic; both Amos' program and our implementation are fast enough to provide many thousands of function values in just a second. The tests in variable precision with Maple T M are much slower and much less exhaustive, but they will allow us to explore higher accuracies.

Testing of the New Liouville-Green Approximations
Here we numerically evaluate the new expansions (3.14) and (3.15). We notice that the analytic continuation formulas of [11, §10.11] can be used to compute cylinder functions for z < 0 from the values for z > 0. For instance, we have J ν ze ±πi = e ±νπi J ν (z).
Therefore, testing for z > 0 is enough. However, we are also considering the case z < 0 for the case of J ν (z) in order to show the full validity region of the expansion (3.15).
A test of the performance of the truncated expansion J ν (νz) ≈ 1 2πν can be seen in Fig. 3. The figure shows the comparison of the function values obtained with n = 14 in (5.1) against those obtained with Amos' algorithm [1] for computing J ν (νz) with ν = 100. Fig. 3 shows the comparison for random values of the variable z generated in the domain −2 < z < 2, −2 < z < 2, respectively. The points where the relative error is greater than 10 −12 are plotted in the figures. As the figure shows, and as expected, the Liouville-Green expansion (3.15) loses accuracy close to the turning point and for real values of z with |z| > 1. It is worth noting than in the neighborhood of the turning point, we can consider our expansions with coefficients computed via Cauchy integrals that we are discussing next, while for z > 1, we can compute Fig. 3 Comparison of the function values obtained with the expansion (3.15) against those obtained with Amos' algorithm [1] for computing J ν (νz) with ν = 100 and −2 < z < 2, −2 < z < 2. The points where the relative error is greater than 10 −12 are plotted Fig. 4 Comparison of the function values obtained with the expansion (3.14) against those obtained with Amos' algorithm [1] for computing H (1) ν (νz) with ν = 100 and 0 < z < 2, 0 < z < 2. The points where the relative error is greater than 10 −12 are plotted J ν (νz) using its relation with Hankel functions and the Liouville-Green expansions for these functions or, alternatively, we can use (4.4) with coefficients computed from its asymptotic approximation. It appears then that for ν ≥ 10 it is possible to compute J ν (νz) in the whole complex z-plane with around 15 digits accuracy by only resorting to asymptotic approximations.
Finally, the performance of the Liouville-Green approximation for z > 0, with again n = 14, is illustrated in Fig. 4. As expected, the expansion fails in the proximity of z = 1.

Airy-Type Expansions via Cauchy's Integral Formula
We now test the accuracy in the computation of the Airy-type expansion (4.1) for H (1) ν (νz). We employ two different approaches in the approximation of the coefficient functions A(ν, z) and B(ν, z), both of which use the Cauchy integral formulas (2.31). In the first approach, on the Cauchy contour, we use the Liouville-Green expansions only for the cylinder functions, and we assume an algorithm for computing complex Airy functions is available; we have used both Amos' algorithm [1] and our own algorithm [7], with similar results. Thus, in (4.1), A(ν, z) and B(ν, z) (expressed by the appropriate Airy and cylinder functions) are approximated by using the Liouville-Green expansions (5.1) and (5.2) in the integrands of (2.31), but with the Airy functions computed from the algorithms cited.
The second approach uses the Liouville-Green approximations for both the Airy functions and the cylinder functions in approximating A(ν, z) and B(ν, z) in (4.1) (although again we do use Airy algorithms to compute Ai −1 ν 2/3 ζ and Ai −1 ν 2/3 ζ when computing H (1) ν (νz) with (5.3)). Hence in this case, from (4.1), (4.5) and (4.6), we have for z inside L, Two different tests of the accuracy of the Airy-type asymptotic expansion for H (1) ν (νz) (with coefficients evaluated by Cauchy's integral formula), for z fixed or ν fixed, are considered. The results are shown in Fig. 5, for which the Airy function has been evaluated both using [7] and [1], with similar results. In all cases, we take m = 7, which means we are considering terms up to O(ν −n ), n = 2m = 14. Of  (1) ν (ν(1 + 0.1i)) against those obtained with Amos' algorithm [1]. In the evaluation of the coefficients, the Airy functions are also computed with [1] course, similar tests can be made for other solutions of the Bessel equation (for instance J ν (νz) or Y ν (νz)), and the accuracy results are very similar, as can be expected since all solutions are computed with the same coefficients.
In Fig. 5, the relative error in the comparison of the function values obtained with the Airy-type expansion for H (1) ν (ν(1 + 0.1i)), 2 < ν < 400, against those obtained with Amos' algorithm is shown. Several different sources of error are apparent in the figure: the error for ν small due to the Liouville-Green approximations for cylinder functions used in the Airy-type expansions, the small increase in the error due to rounding for ν large caused by the B(ν, z) coefficient, and by the unavoidable loss of accuracy in the computation of Airy functions for large arguments.
In Fig. 6, we show the same results as in Fig. 5, but using the asymptotic expansions for the coefficients given in Sect. 4. We observe a clear improvement over Fig. 5, particularly for low ν. As ν becomes larger, there is a slow increase in accuracy loss (smaller than in the previous case). This is due to the increasing inaccuracies in the computation of Airy functions as the argument becomes larger, which give rise to errors in the computation of the Hankel function when using (4.1). These inaccuracies are unavoidable in finite precision arithmetic and, as described in [7], can only be removed by considering scaled functions (with the dominant exponential factor exactly scaled out). Figure 7 shows the same tests, but focusing in smaller values of ν and for two selections in the number of terms for the expansion. We observe two tendencies, particularly for n = 2m = 18: a decrease of the error as ν increases due to the fact that the expansion becomes more accurate, and a slow increase of the error due to the finite precision computation of the Airy functions in the expression for H (1) ν (νz); this error increase is not attributable to our approach. These tendencies give an optimal accuracy for ν close to 5. As commented, the accuracy for larger ν in finite precision arithmetic can be improved by considering scaled functions.
In Fig. 8, we plot the points where the relative error is greater than 10 −13 in the comparison of the function values obtained with the Airy-type expansion for H (1) 10 (10z) against those obtained with Amos' algorithm. The random z points in the test have been generated inside the semi-circle limited by the integration contour used to compute  (1) ν (ν(1 + 0.1i)) against those obtained with Amos' algorithm [1]. The asymptotic expansion of the coefficients is considered over the Cauchy contour. We take n = 14 = 2m in Eqs.   10 (10z) against those obtained with Amos' algorithm [1]. The points where the relative error is greater than 10 −13 are plotted. We take n = 2m = 14 the coefficient functions A(ν, z) and B(ν, z) of the Airy-type asymptotic expansion (a circle of center z c = 2 and radius R = 1.8). As can be seen in the figure, the accuracy obtained with the expansion is better than 10 −13 in a large portion of the domain, although, as expected, it worsens when approaching the integration contour.
This figure illustrates the accuracy in the discretization of the Cauchy integral in a large region, but with z not too close to the contour of integration.
Preliminary tests show that the time spent for the computation of the Liouville-Green expansions in §3 and the Airy-type expansions in §4 is similar to the time for the implementations of Debye and Airy-type expansions in Amos' algorithm [1]. However, as expected, when Cauchy's integral formula is used, our approach becomes slower for computing a single function value, because we need to compute the coefficients a number of times over the contour.
The advantage of the Cauchy approach is that the most costly computation is the evaluation of the coefficients E j (ξ ) and the numerators in the sums of Eqs. (2.20) and (2.21), but this computation is done once and for all over the contour. Our approach permits the selection of a convenient Cauchy contour depending on the application. If many function values are needed in a certain z-region, a good approach is to precompute the coefficients in a circuit containing this region (if it is possible and the shaded regions of Fig. 2 can be avoided), and then applying the discretized version of the Cauchy integral formula (Eq. (2.33)) as many times as needed (and the recomputation of the sums in Eqs. (2.20) and (2.21) if ν is also varying).
At this point, it is important to stress that the main interest of the Cauchy technique lies in the computability more than in the efficiency (although it is efficient). It provides a new and direct method of computation of the coefficients of Airy-type expansions with potential applications to more complicated cases.
In addition to the tests in fixed precision arithmetic, we have performed additional test in variable precision using Maple T M . First, we have checked that the coefficients in the Airy-type expansions are computed in a numerically stable way with our scheme and that the remaining error degradation in Fig. 6 is due to loss of accuracy in the computation of the Airy functions when using (5.3). For this purpose, we have used Maple T M for computing the coefficients with a fixed number of digits (we take 16 digits) and then computed the Hankel function using (5.3) with a sufficiently high number of digits. When we test the value of the Hankel function with, say, 50 digits, we observe that the error is always close to 16 digits accuracy, which shows that the coefficients have been consistently computed with that accuracy and without accuracy loss for high ν. This is illustrated in Fig. 9.     In all the computations, we have used 500 points in the contour of integration. The error associated with the discretization of the Cauchy integral is so small that it has no impact on the previous results. For observing this error, a higher number of digits should be considered. Fig. 10 shows the results when 500 points over the Cauchy contour are considered and the computations are done with 50 digits.
As can be expected, because the integrand is periodic, the trapezoidal rule has exponential convergence. Indeed, we have checked that smallest reachable error roughly depends on the number of points over the Cauchy contour as 10 −N /20 .
We have checked numerically that when the effect of discretization of the Cauchy integral is negligible, the relative error when n terms in the asymptotic expansion of the coefficients are considered, the relative error in the computation of H (1) ν (νz) close to the turning point varies as C n+1 /ν n+1 , where C n+1 does not depend on ν and increases with n. This is the expected behavior for an asymptotic expansion for large ν. Some estimations of these constants close to z = 1 are shown in Table 1.
It would be important to be able to establish strict error bounds and to compare them with these experimental errors. This should be achievable by bounding the errors in the Liouville-Green expansions used for the coefficients, similar to the error bounds given in [5], and again appealing to Cauchy's integral formula. This will be considered in a subsequent paper.