About the Project
9 Airy and Related FunctionsComputation

§9.17 Methods of Computation

  1. §9.17(i) Maclaurin Expansions
  2. §9.17(ii) Differential Equations
  3. §9.17(iii) Integral Representations
  4. §9.17(iv) Via Bessel Functions
  5. §9.17(v) Zeros

§9.17(i) Maclaurin Expansions

Although the Maclaurin-series expansions of §§9.4 and 9.12(vi) converge for all finite values of z, they are cumbersome to use when |z| is large owing to slowness of convergence and cancellation. For large |z| the asymptotic expansions of §§9.7 and 9.12(viii) should be used instead. Since these expansions diverge, the accuracy they yield is limited by the magnitude of |z|. However, in the case of Ai(z) and Bi(z) this accuracy can be increased considerably by use of the exponentially-improved forms of expansion supplied in §9.7(v).

§9.17(ii) Differential Equations

A comprehensive and powerful approach is to integrate the defining differential equation (9.2.1) by direct numerical methods. As described in §3.7(ii), to ensure stability the integration path must be chosen in such a way that as we proceed along it the wanted solution grows at least as fast as all other solutions of the differential equation. In the case of Ai(z), for example, this means that in the sectors 13π<|phz|<π we may integrate along outward rays from the origin with initial values obtained from §9.2(ii). But when |phz|<13π the integration has to be towards the origin, with starting values of Ai(z) and Ai(z) computed from their asymptotic expansions. On the remaining rays, given by phz=±13π and π, integration can proceed in either direction.

For further information see Lozier and Olver (1993) and Fabijonas et al. (2004). The former reference includes a parallelized version of the method.

In the case of the Scorer functions, integration of the differential equation (9.12.1) is more difficult than (9.2.1), because in some regions stable directions of integration do not exist. An example is provided by Gi(x) on the positive real axis. In these cases boundary-value methods need to be used instead; see §3.7(iii).

§9.17(iii) Integral Representations

Among the integral representations of the Airy functions the Stieltjes transform (9.10.18) furnishes a way of computing Ai(z) in the complex plane, once values of this function can be generated on the positive real axis. For details, including the application of a generalized form of Gaussian quadrature, see Gordon (1969, Appendix A) and Schulten et al. (1979).

Gil et al. (2002c) describes two methods for the computation of Ai(z) and Ai(z) for z. In the first method the integration path for the contour integral (9.5.4) is deformed to coincide with paths of steepest descent (§2.4(iv)). The trapezoidal rule (§3.5(i)) is then applied. The second method is to apply generalized Gauss–Laguerre quadrature (§3.5(v)) to the integral (9.5.8). For the second method see also Gautschi (2002a). The methods for Ai(z) are similar.

For quadrature methods for Scorer functions see Gil et al. (2001), Lee (1980), and Gordon (1970, Appendix A); but see also Gautschi (1983).

§9.17(iv) Via Bessel Functions

In consequence of §9.6(i), algorithms for generating Bessel functions, Hankel functions, and modified Bessel functions (§10.74) can also be applied to Ai(z), Bi(z), and their derivatives.

§9.17(v) Zeros

Zeros of the Airy functions, and their derivatives, can be computed to high precision via Newton’s rule (§3.8(ii)) or Halley’s rule (§3.8(v)), using values supplied by the asymptotic expansions of §9.9(iv) as initial approximations. This method was used in the computation of the tables in §9.9(v). See also Fabijonas et al. (2004).

For the computation of the zeros of the Scorer functions and their derivatives see Gil et al. (2003c).