Relating Airy and Bessel functions

The Endeavour 2013-09-01

The Airy functions Ai(x) and Bi(x) are independent solutions to the differential equation

For negative x they act something like sin(x) and cos(x). For positive x they act something like exp(x) and exp(-x). This isn’t surprising if you look at the differential equation. If you replace x with a negative constant, you sines and cosines, and if you replace it with a positive constant, you get positive and negative exponentials.

The Airy functions can be related to Bessel functions as follows:

$\mathrm{Ai}(x) = \left\{ \begin{array}{ll} \frac{1}{3}\sqrt{\phantom{-}x} \left(I_{-1/3}(\hat{x}) - I_{1/3}(\hat{x})\right) & \mbox{if } x > 0 \\ \\ \frac{1}{3}\sqrt{-x} \left(J_{-1/3}(\hat{x}) + J_{1/3}(\hat{x})\right) & \mbox{if } x < 0 \end{array} \right.$

and

$\mathrm{Bi}(x) = \left\{ \begin{array}{ll} \sqrt{\phantom{-}x/3} \left(I_{-1/3}(\hat{x}) + I_{1/3}(\hat{x})\right) & \mbox{if } x > 0 \\ \\ \sqrt{-x/3} \left(J_{-1/3}(\hat{x}) - J_{1/3}(\hat{x})\right) & \mbox{if } x < 0 \end{array} \right.$

Here J is a “Bessel function of the first kind” and I is a “modified Bessel function of the first kind.” Also

$\hat{x} = \frac{2}{3} \left(\sqrt{|x|}\right)^3$

To verify the equations above, and to show how to compute these functions in Python, here’s some code.

The SciPy function airy computes both functions, and their first derivatives, at once. I assume that’s because it doesn’t take much longer to compute all four functions than to compute one. The code for Ai2 and Bi2 below uses np.where instead of if ... else so that it can operate on NumPy vectors all at once. You can plot Ai and Ai2 and see that the two curves lie on top of each other. The same holds for Bi and Bi2.

from scipy.special import airy, jv, ivfrom numpy import sqrt, wheredef Ai(x):    (ai, ai_prime, bi, bi_prime) = airy(x)    return aidef Bi(x):   (ai, ai_prime, bi, bi_prime) = airy(x)    return bidef Ai2(x):    third = 1.0/3.0    hatx = 2*third*(abs(x))**1.5    return where(x > 0,        third*sqrt( x)*(iv(-third, hatx) - iv(third, hatx)),        third*sqrt(-x)*(jv(-third, hatx) + jv(third, hatx)))def Bi2(x):    third = 1.0/3.0    hatx = 2*third*(abs(x))**1.5    return where(x > 0,        sqrt( x/3.0)*(iv(-third, hatx) + iv(third, hatx)),        sqrt(-x/3.0)*(jv(-third, hatx) - jv(third, hatx)))

There is a problem with Ai2 and Bi2: they return nan at 0. A more careful implementation would avoid this problem, but that’s not necessary since these functions are only for illustration. In practice, you’d simply use airy and it does the right thing at 0.