The functions \( J_0(x) \), \( J_1(x) \), and \( J_n(x) \) are the solutions of Bessel's differential equation
\[ x^2 \frac{d^2 y}{dx^2} + x \frac{dy}{dx} + (x^2 - \nu^2)\, y = 0 \]
that remain finite at the origin. For non-negative integer order \( n \) they admit the series
\[ J_n(x) = \sum_{m=0}^{\infty} \frac{(-1)^m}{m!\,(m+n)!} \left(\frac{x}{2}\right)^{2m+n}. \]
\( J_\nu(x) \) is oscillatory and decays like \( \sqrt{2/(\pi x)}\,\cos\!\left(x - \nu\pi/2 - \pi/4\right) \) as \( x \to \infty \).
Function shape
Domain and special values
| \( x \) | \( J_0(x) \) | \( J_1(x) \) |
|---|---|---|
| \( 0 \) | \( 1 \) | \( 0 \) |
| \( -x \) | \( J_0(x) \) (even) | \( -J_1(x) \) (odd) |
| \( +\infty \) | \( 0 \) | \( 0 \) |
The J family is defined for all real \( x \). besseljn(n, x) follows the parity of n: even n gives an even function, odd n an odd function.
Usage
Accuracy and performance
besselj0/besselj1 use minimax polynomials near the origin, root- and extremum-centred polynomials on the middle range, and a Hankel asymptotic expansion for large \( x \) (the SQ2OPI/PIO4 phase form). On the benchmark machine they run roughly 2× faster than the gfortran intrinsics bessel_j0/bessel_j1 while matching them to within a few ULP. besseljn composes the fixed-order kernels through a stable recurrence.
The figure below shows the relative error of this library against an mpmath arbitrary-precision reference, in ULPs of real64:
