§9.17 Methods of Computation

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 $\mathrm{Ai}\left(z\right)$ and $\mathrm{Bi}\left(z\right)$ this accuracy can be increased considerably by use of the exponentially-improved forms of expansion supplied in §9.7(v).

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 $\mathrm{Ai}\left(z\right)$, for example, this means that in the sectors we may integrate along outward rays from the origin with initial values obtained from §9.2(ii). But when the integration has to be towards the origin, with starting values of $\mathrm{Ai}\left(z\right)$ and ${\mathrm{Ai}}^{\prime }\left(z\right)$ computed from their asymptotic expansions. On the remaining rays, given by $\mathrm{ph}z=\pm \frac{1}{3}\pi$ and $\pi$, 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 $\mathrm{Gi}\left(x\right)$ on the positive real axis. In these cases boundary-value methods need to be used instead; see §3.7(iii).

Among the integral representations of the Airy functions the Stieltjes transform (9.10.18) furnishes a way of computing $\mathrm{Ai}\left(z\right)$ 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 $\mathrm{Ai}\left(z\right)$ and ${\mathrm{Ai}}^{\prime }\left(z\right)$ for $z\in ℂ$. 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 ${\mathrm{Ai}}^{\prime }\left(z\right)$ 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).

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

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).