A Simple Approximation for the Airy Function

Besides the Normal Distribution Function, I occasionally need the Airy Function Ai(x): it arises in perturbation theory and some other contexts. This function is most definitely not part of most standard numerics libraries! While high-quality implementations of it are part of most “serious” numerics libraries (such as SciPy or GSL), these libraries are not always available or convenient.

Here is a really simple numerical approximation: it is numerically “good enough” for casual work, and simple enough to be implemented on the fly as needed.

As usual with such strictly numerical approximation schemes, the formulas are quite opaque and the implementation contains its share of magic constants.

from math import sqrt, sin, cos, exp

def _ap(x):
    p0 = -0.002800908
    p1 =  0.326662423
    P0 = -0.007232251
    P1 = -0.044567423

    z = sqrt(0.0425 + x**3.0)
    z1 = z**(7.0/6.0)
    z2 = z**(2.0/3.0)
    return ( (p0 + p1*z) + (x/z2)*(P0 + P1*z) )*exp(-2.0*z/3.0)/z1

def _am(x):
    q0 = -0.043883564
    q1 =  0.3989422
    Q0 = -0.013883003
    Q1 = -0.3989422

    x0 = (-x)**3.0
    z = sqrt(0.37 + x0)
    x0 = sqrt(x0)

    z1 = z**(7.0/6.0)
    z2 = z**(5.0/6.0)
    return (q0+q1*z)*cos(2.0*x0/3.0)/z1 + x*(Q0+Q1*z)*sin(2.0*x0/3.0)/(z2*x0)

def airy(x):
    return _am(x) if x < 0 else _ap(x)

As I said in the introduction, the approximation is “good enough”, but not great: the largest absolute error is about $\pm 0.003$ for values of $x$ near zero, equivalent to $-1 \ldots +2$ percent relative error in that region. But for casual work, the simplicity of the implementation may well make up for the lack of precision.

Resources:

This approximation is presented in the paper:

Two-point Quasi-Fractional Approximation to the Airy Function Ai(x) by Pablo Martín, Ricardo Pérez, Antonio L. Guerrero; Journal of Computational Physics, Volume 99 (1992) 337

Beware of a misprint in one of the formulas; this misprint has been corrected in the Python implementation above.