Special Functions: Bessel, Airy & Orthogonal Polynomials

Overview

This page describes the numerical evaluation algorithms for three closely related families of functions: the Bessel-type functions (cylinder functions), the Airy functions, and the classical orthogonal polynomials. Each arises as a solution of a second-order linear ordinary differential equation, and all are evaluated with the same toolkit of series expansions, asymptotic expansions and three-term recurrences.

• Bessel-type functions: solutions of the Bessel equation $x^2 y'' + x y' + (x^2 - \nu^2) y = 0$. They appear in wave, diffusion and potential problems in cylindrical coordinates. The family includes the first/second kind $J_\nu, Y_\nu$, the modified functions $I_\nu, K_\nu$, the spherical Bessel functions and the Hankel functions.
• Airy functions: the solutions $\mathrm{Ai}, \mathrm{Bi}$ of $y'' = x y$. They appear in asymptotics near turning points (WKB), in the caustics of optics and in uniform approximations of quantum mechanics. They can be expressed through Bessel functions of order $\pm 1/3$.
• Orthogonal polynomials: sequences of polynomials orthogonal with respect to a given weight function. Legendre appears in the angular part of spherical harmonics, Hermite in the quantum harmonic oscillator, and Laguerre in the radial wavefunctions of the hydrogen atom.

Related API: SpecialFunctions. Standard real-argument Bessel and spherical Bessel functions delegate to the C++17 standard library (std::cyl_bessel_j and the like), while complex arguments and arbitrary precision (Complex<double>, Complex<Float>) are evaluated by sangi's own series.

Bessel Functions $J_\nu, Y_\nu, I_\nu, K_\nu$

Four kinds of cylinder functions are handled:

• First kind $J_\nu(x)$ and second kind (Neumann) $Y_\nu(x)$ — the two independent solutions of the ordinary Bessel equation
• Modified first kind $I_\nu(x)$ and modified second kind $K_\nu(x)$ — solutions of the modified Bessel equation $x^2 y'' + x y' - (x^2 + \nu^2) y = 0$

Small arguments: power series

Near the origin the Taylor series converges quickly. The first kind is

$$J_\nu(z) = \left(\frac{z}{2}\right)^{\nu} \sum_{m=0}^{\infty} \frac{(-1)^m}{m!\,\Gamma(m+\nu+1)} \left(\frac{z}{2}\right)^{2m}$$

The modified first kind has no sign alternation:

$$I_\nu(z) = \left(\frac{z}{2}\right)^{\nu} \sum_{m=0}^{\infty} \frac{1}{m!\,\Gamma(m+\nu+1)} \left(\frac{z}{2}\right)^{2m}$$

sangi's complex integer-order implementation evaluates this series recursively. Each term is updated from the previous one by

$$\text{term}_m = \text{term}_{m-1} \cdot \frac{\mp (z/2)^2}{m\,(m+n)}$$

so there is no need to recompute factorials or the gamma function at each step ($J$ carries $-(z/2)^2$ in the numerator, $I$ carries $+(z/2)^2$). The sum is truncated at the relative tolerance $|\text{term}_m| < \varepsilon\,|\text{sum}|$.

Integer order $Y_n, K_n$: series with a logarithmic term

The integer-order second-kind functions diverge logarithmically at the origin, so a series containing $\ln(z/2)$ is used (A&S 9.1.11 / 9.6.11). The structure of $Y_n$ is

$$Y_n(z) = \frac{2}{\pi} J_n(z)\left(\gamma + \ln\frac{z}{2}\right) - \frac{1}{\pi}\left(\frac{z}{2}\right)^{-n}\!\!\sum_{k=0}^{n-1}\frac{(n-k-1)!}{k!}\left(\frac{z^2}{4}\right)^{k} - \frac{1}{\pi}\left(\frac{z}{2}\right)^{n}\!\!\sum_{k=0}^{\infty}\frac{(-1)^k (H_k + H_{k+n})}{k!\,(k+n)!}\left(\frac{z^2}{4}\right)^{k}$$

Here $\gamma$ is the Euler-Mascheroni constant and $H_k = \sum_{j=1}^{k} 1/j$ is the harmonic number. $K_n$ takes a similar form containing $\ln(z/2)$ and $(H_k + H_{k+n} - 2\gamma)$ (with different signs and coefficients). sangi builds up the harmonic numbers and gamma terms recursively and reuses the $J_n$ / $I_n$ above as leading terms.

Large arguments: asymptotic expansion

In the region of large argument the oscillating / decaying asymptotic forms are used. The first kind is

$$J_\nu(x) \sim \sqrt{\frac{2}{\pi x}}\,\cos\!\left(x - \frac{\nu\pi}{2} - \frac{\pi}{4}\right), \qquad x \to \infty$$

This is a divergent series, so it is truncated just before the term magnitude reaches its minimum (optimal truncation). sangi's real-order Bessel functions delegate to the standard library, so this region splitting is handled inside that implementation.

Backward recurrence (Miller's algorithm)

To obtain $J_n$ for all orders up to an integer at once, the three-term recurrence

$$J_{\nu+1}(x) = \frac{2\nu}{x} J_\nu(x) - J_{\nu-1}(x)$$

becomes unstable when run in the direction of increasing order (forward). $J_n$ and $Y_n$ satisfy the same recurrence, and in the forward direction the $Y_n$ component (which grows exponentially with order) is excited by rounding error and swamps the solution.

The stable direction is backward. Miller's algorithm starts from a sufficiently high order $N \gg n$ with the assumptions $J_{N+1} = 0,\ J_N = \varepsilon$ and runs

$$J_{\nu-1}(x) = \frac{2\nu}{x} J_\nu(x) - J_{\nu+1}(x)$$

decreasing $\nu = N, N-1, \ldots$, and finally normalizes everything with the summation rule $J_0(x) + 2\sum_{k\geq 1} J_{2k}(x) = 1$. Because the error component decays in the backward direction, the ratio $J_n / J_0$ is recovered correctly regardless of the initial values.

Wronskian normalization and derivatives

The Wronskian of the first and second kind is

$$W\{J_\nu, Y_\nu\}(x) = J_\nu(x) Y_\nu'(x) - J_\nu'(x) Y_\nu(x) = \frac{2}{\pi x}$$

which can be used to set the scale of the independent solutions and to check their mutual consistency. Derivatives are obtained from neighbouring orders by recurrence (DLMF 10.6.1, 10.29.1, etc.):

$$J_\nu'(x) = \tfrac{1}{2}\big(J_{\nu-1}(x) - J_{\nu+1}(x)\big), \qquad I_\nu'(x) = \tfrac{1}{2}\big(I_{\nu-1}(x) + I_{\nu+1}(x)\big)$$

$Y_\nu'$ and $K_\nu'$ are analogous ($K$ satisfies $K_\nu'(x) = -\tfrac12(K_{\nu-1}+K_{\nu+1})$). sangi's besselJPrime, besselIPrime and so on implement these recurrences directly.

Spherical Bessel and Hankel functions

The spherical Bessel functions are half-integer-order Bessel functions that close in elementary functions:

$$j_n(x) = \sqrt{\frac{\pi}{2x}}\, J_{n+1/2}(x), \qquad y_n(x) = \sqrt{\frac{\pi}{2x}}\, Y_{n+1/2}(x)$$

Starting from $j_0(x) = \sin x / x$ and $y_0(x) = -\cos x / x$, the recurrence (DLMF 10.51.1)

$$j_{n+1}(x) = \frac{2n+1}{x} j_n(x) - j_{n-1}(x)$$

is applied in the forward direction (practically stable for small orders and moderate arguments). The Hankel functions (third kind) represent travelling waves as complex combinations built from the first and second kind:

$$H_n^{(1)}(z) = J_n(z) + i\,Y_n(z), \qquad H_n^{(2)}(z) = J_n(z) - i\,Y_n(z)$$

The spherical Hankel functions $h_n^{(1,2)} = j_n \pm i\,y_n$ are formed in the same way. sangi provides these as hankelH1, hankelH2, sphericalHankelH1 and so on.

Airy Functions $\mathrm{Ai}, \mathrm{Bi}$

The Airy functions are the two independent solutions of the equation $y'' = x y$. $\mathrm{Ai}(x)$ decays as $x \to +\infty$, while $\mathrm{Bi}(x)$ diverges.

Maclaurin series near the origin

Take the two power-series solutions of the Airy equation:

$$f(x) = 1 + \frac{x^3}{2\cdot 3} + \frac{x^6}{2\cdot 3\cdot 5\cdot 6} + \cdots, \qquad g(x) = x + \frac{x^4}{3\cdot 4} + \frac{x^7}{3\cdot 4\cdot 6\cdot 7} + \cdots$$

These are entire functions expressible through the hypergeometric function $_0F_1$ ($f = {}_0F_1(;\tfrac23;\tfrac{x^3}{9})$, etc.), and each term is updated by the recurrence

$$f_{\text{term}_{k+1}} = f_{\text{term}_{k}} \cdot \frac{x^3}{(3k+2)(3k+3)}, \qquad g_{\text{term}_{k+1}} = g_{\text{term}_{k}} \cdot \frac{x^3}{(3k+3)(3k+4)}$$

The Airy functions are the linear combinations of these two solutions:

$$\mathrm{Ai}(x) = \mathrm{Ai}(0)\, f(x) + \mathrm{Ai}'(0)\, g(x), \qquad \mathrm{Bi}(x) = \mathrm{Bi}(0)\, f(x) + \mathrm{Bi}'(0)\, g(x)$$

The combination coefficients are given by the gamma function:

$$\mathrm{Ai}(0) = \frac{1}{3^{2/3}\,\Gamma(2/3)}, \quad \mathrm{Ai}'(0) = -\frac{1}{3^{1/3}\,\Gamma(1/3)}, \quad \mathrm{Bi}(0) = \frac{1}{3^{1/6}\,\Gamma(2/3)}, \quad \mathrm{Bi}'(0) = \frac{3^{1/6}}{\Gamma(1/3)}$$

Since $f, g$ are entire, the series converges over the whole plane and applies directly to complex arguments. The derivatives $\mathrm{Ai}', \mathrm{Bi}'$ are evaluated as series for $f', g'$ in the same framework. The arbitrary-precision version computes $\Gamma(1/3), \Gamma(2/3)$ and $3^{1/3}, 3^{1/6}$ at the working precision before forming the combination coefficients.

Representation through Bessel functions

The Airy functions can be expressed through Bessel functions of order $\pm 1/3$, which makes their relationship explicit:

$$\mathrm{Ai}(x) = \frac{1}{\pi}\sqrt{\frac{x}{3}}\, K_{1/3}\!\left(\tfrac{2}{3}x^{3/2}\right) \quad (x>0)$$

$$\mathrm{Ai}(-x) = \frac{\sqrt{x}}{3}\left[J_{1/3}\!\left(\tfrac{2}{3}x^{3/2}\right) + J_{-1/3}\!\left(\tfrac{2}{3}x^{3/2}\right)\right] \quad (x>0)$$

This shows that the Airy functions belong to the same cylinder-function framework as the Bessel family, and at the same time provides a way to reuse the asymptotics of the modified and ordinary Bessel functions for large-argument evaluation.

Asymptotics for large $|x|$

With $\zeta = \tfrac{2}{3}|x|^{3/2}$, the asymptotic forms are used in the region away from the origin:

$$\mathrm{Ai}(x) \sim \frac{e^{-\zeta}}{2\sqrt{\pi}\, x^{1/4}} \quad (x\to +\infty), \qquad \mathrm{Ai}(-x) \sim \frac{1}{\sqrt{\pi}\, x^{1/4}}\sin\!\left(\zeta + \tfrac{\pi}{4}\right) \quad (x\to +\infty)$$

The behaviour characteristic of the Airy functions — exponential decay on the positive side and oscillation on the negative side (the transition across the turning point) — appears here.

Kelvin, Struve & Anger-Weber

This section treats the function families that derive from, or are adjacent to, the Bessel functions.

Kelvin functions ber, bei, ker, kei

The Kelvin functions are defined as the real and imaginary parts of Bessel functions of complex argument (order 0):

$$\mathrm{ber}(x) + i\,\mathrm{bei}(x) = J_0\!\left(x\,e^{3i\pi/4}\right), \qquad \mathrm{ker}(x) + i\,\mathrm{kei}(x) = K_0\!\left(x\,e^{i\pi/4}\right)$$

sangi's implementation calls $J_0$ / $K_0$ of complex argument directly and extracts the real and imaginary parts. $\mathrm{ber}, \mathrm{bei}$ are even functions defined for all real numbers, while $\mathrm{ker}, \mathrm{kei}$ are defined only for $x>0$ ($\mathrm{ker}$ diverges logarithmically as $x\to 0^+$, and $\mathrm{kei}(0) = -\pi/4$). They appear in the skin effect of electrical engineering and in the magnetic-field calculation of magnetic cores inside solenoid coils.

Struve functions $H_\nu, L_\nu$

The Struve functions are particular solutions of the inhomogeneous Bessel equation. For small arguments they use the series (DLMF 11.2.1):

$$H_\nu(z) = \left(\frac{z}{2}\right)^{\nu+1}\sum_{k=0}^{\infty} \frac{(-1)^k}{\Gamma(k+\tfrac32)\,\Gamma(\nu+k+\tfrac32)}\left(\frac{z}{2}\right)^{2k}$$

The modified Struve function $L_\nu$ has the same form without the sign alternation. For large arguments the asymptotic expansion (DLMF 11.6.1)

$$H_\nu(z) \sim Y_\nu(z) + \frac{1}{\pi}\sum_{k=0}^{N} \frac{\Gamma(k+\tfrac12)}{\Gamma(\nu+\tfrac12-k)}\left(\frac{z}{2}\right)^{\nu-2k-1}$$

is used ($L_\nu$ has $I_\nu$ as its leading term with a different sign), and the asymptotic series is truncated once the term magnitudes start to grow. sangi switches between the series and the asymptotic form at the boundary $|z| \lesssim 16$. At present only real arguments with $z \geq 0$ are supported.

Anger function $J_\nu$ and Weber function $E_\nu$

The Anger and Weber functions are defined by integrals (DLMF 11.10):

$$\mathbf{J}_\nu(z) = \frac{1}{\pi}\int_0^{\pi}\cos(\nu t - z\sin t)\,dt, \qquad \mathbf{E}_\nu(z) = \frac{1}{\pi}\int_0^{\pi}\sin(\nu t - z\sin t)\,dt$$

When the order $\nu$ is an integer $n$ the Anger function coincides with the Bessel function $J_n$, but for non-integer orders it forms an independent family. Since the integrand is a smooth analytic function, sangi integrates the defining integral directly by the composite Simpson rule (the series representation is not used, because handling the integer-order limit is cumbersome).

Orthogonal Polynomials

The classical orthogonal polynomials are sequences of polynomials orthogonal over an interval with respect to a given weight function $w(x)$. All of them can be evaluated stably with a three-term recurrence (a consequence of Favard's theorem).

Legendre polynomials $P_n(x)$

Weight $w(x) = 1$, interval $[-1, 1]$. Evaluated with Bonnet's recurrence:

$$(n+1) P_{n+1}(x) = (2n+1)\, x\, P_n(x) - n\, P_{n-1}(x), \qquad P_0 = 1,\ P_1 = x$$

Chebyshev polynomials $T_n(x), U_n(x)$

The first kind $T_n$ (weight $1/\sqrt{1-x^2}$) and the second kind $U_n$ (weight $\sqrt{1-x^2}$) share the same recurrence and differ only in the initial values:

$$T_{n+1}(x) = 2x\, T_n(x) - T_{n-1}(x), \quad T_0=1,\ T_1=x$$

$$U_{n+1}(x) = 2x\, U_n(x) - U_{n-1}(x), \quad U_0=1,\ U_1=2x$$

The trigonometric representations $T_n(\cos\theta) = \cos(n\theta)$ and $U_n(\cos\theta) = \sin((n+1)\theta)/\sin\theta$ also hold. The recurrence is stable for $|x| \leq 1$ and can be used as is for $|x| > 1$.

Hermite polynomials (physicists' / probabilists')

The recurrences differ between the physicists' $H_n$ (weight $e^{-x^2}$, eigenfunctions of the quantum harmonic oscillator) and the probabilists' $He_n$ (weight $e^{-x^2/2}$):

$$H_{n+1}(x) = 2x\, H_n(x) - 2n\, H_{n-1}(x), \quad H_0=1,\ H_1=2x$$

$$He_{n+1}(x) = x\, He_n(x) - n\, He_{n-1}(x), \quad He_0=1,\ He_1=x$$

The two are related by $H_n(x) = 2^{n/2} He_n(x\sqrt{2})$.

Laguerre polynomials $L_n(x), L_n^{\alpha}(x)$

Weight $e^{-x}$ (the generalized version has $x^\alpha e^{-x}$), interval $[0, \infty)$. They appear in the radial wavefunctions of the hydrogen atom.

$$(k+1) L_{k+1}^{\alpha}(x) = (2k+1+\alpha-x)\, L_k^{\alpha}(x) - (k+\alpha)\, L_{k-1}^{\alpha}(x), \quad L_0^{\alpha}=1,\ L_1^{\alpha}=1+\alpha-x$$

Setting $\alpha = 0$ reduces this to the ordinary Laguerre polynomials $L_n(x)$.

Jacobi polynomials $P_n^{(\alpha,\beta)}(x)$ and Gegenbauer $C_n^{\lambda}(x)$

The Jacobi polynomials are orthogonal with respect to the weight $(1-x)^{\alpha}(1+x)^{\beta}$ ($\alpha, \beta > -1$) and contain many of the above in a unified way:

$$P_n^{(0,0)} = P_n \ (\text{Legendre}), \quad P_n^{(-1/2,-1/2)} \propto T_n, \quad P_n^{(1/2,1/2)} \propto U_n, \quad P_n^{(\alpha,\alpha)} \propto C_n^{\alpha+1/2}$$

Evaluation runs the three-term recurrence of DLMF 18.9.2 (whose coefficients depend on $\alpha, \beta, n$). The Gegenbauer (ultraspherical) polynomials $C_n^{\lambda}$ are proportional to the symmetric Jacobi polynomials and have their own recurrence

$$n\, C_n^{\lambda}(x) = 2(n+\lambda-1)\, x\, C_{n-1}^{\lambda}(x) - (n+2\lambda-2)\, C_{n-2}^{\lambda}(x)$$

sangi provides Chebyshev, Hermite and Laguerre as dedicated implementations (lower overhead than going through Jacobi) and uses jacobiP for general $\alpha, \beta$.

Clenshaw algorithm

When evaluating a series of orthogonal polynomials $S(x) = \sum_{k=0}^{N} c_k\, p_k(x)$ as a whole, the Clenshaw algorithm is more efficient and stable than building each $p_k(x)$ and multiplying and adding. For a recurrence $p_{k+1} = (\alpha_k x + \beta_k) p_k - \gamma_k p_{k-1}$, compute the auxiliary sequence $b_k$ in the backward direction:

$$b_k = c_k + (\alpha_k x + \beta_k)\, b_{k+1} - \gamma_{k+1}\, b_{k+2}, \qquad b_{N+1} = b_{N+2} = 0$$

and $S(x)$ is then expressed with only $b_0, b_1$ and the low-order $p_0, p_1$. This is the standard technique for evaluating Chebyshev series and spherical-harmonic expansions.

Associated Legendre and spherical harmonics

The associated Legendre functions $P_n^m(x)$ appear in the angular part of spherical harmonics and are built up by recurrence starting from the diagonal terms:

$$P_m^m(x) = (2m-1)!!\,(1-x^2)^{m/2}, \quad P_{m+1}^m(x) = x(2m+1) P_m^m(x)$$

$$(n-m+1) P_{n+1}^m(x) = (2n+1)\, x\, P_n^m(x) - (n+m)\, P_{n-1}^m(x)$$

The normalized spherical harmonics (real part) are given by

$$Y_n^m(\theta) = \sqrt{\frac{2n+1}{4\pi}\cdot\frac{(n-m)!}{(n+m)!}}\; P_n^m(\cos\theta)$$

sangi's sphLegendre computes the factorial ratio $(n-m)!/(n+m)!$ in log space to avoid overflow before forming the normalization coefficient. This implementation does not include the Condon-Shortley phase (consistent with the C++17 standard library std::assoc_legendre).

Comparison Table

The common design principles are to update each series term recursively to avoid recomputing factorials and gamma values, to run recurrences in the stable direction (backward for Bessel, forward for orthogonal polynomials), and to switch methods by region (series for small arguments, asymptotics for large arguments).

References

• Abramowitz, M. & Stegun, I. A. (1964). Handbook of Mathematical Functions. National Bureau of Standards. (§9 Bessel functions, §10 Airy functions, §12 Struve and Anger-Weber, §22 orthogonal polynomials)
• Olver, F. W. J., Lozier, D. W., Boisvert, R. F. & Clark, C. W. (eds.) (2010). NIST Handbook of Mathematical Functions. Cambridge University Press. (DLMF, dlmf.nist.gov §10, §11, §18)
• Olver, F. W. J. (1974). Asymptotics and Special Functions. Academic Press.
• Watson, G. N. (1944). A Treatise on the Theory of Bessel Functions. 2nd ed. Cambridge University Press.
• Press, W. H., Teukolsky, S. A., Vetterling, W. T. & Flannery, B. P. (2007). Numerical Recipes: The Art of Scientific Computing. 3rd ed. Cambridge University Press. (§6 special functions, Miller's algorithm, Clenshaw recurrence)
• Gil, A., Segura, J. & Temme, N. M. (2007). Numerical Methods for Special Functions. SIAM.