Recurrence relations for Bessel functions show directional numerical instability

Recurrence relations for Bessel functions are a common tool in computational mathematics, but their stability depends on both the function type and the direction of computation. For Bessel functions of the first kind (Jn), computing higher-order values from lower-order ones (forward direction) amplifies rounding errors, while the reverse direction remains stable. The opposite holds for Bessel functions of the second kind (Yn), where forward computation is stable but reverse is not. This asymmetry arises because one solution decays while the other grows with increasing order, and rounding errors can inadvertently introduce the growing component into the decaying solution.

The instability stems from the behavior of the two independent solutions to second-order linear difference equations. When one solution decays and the other grows, numerical errors can dominate the decaying solution over time. For example, Jn(x) decays to zero as n increases, while Yn(x) diverges to negative infinity. If a small component of the growing solution is introduced due to rounding, it can overwhelm the desired decaying solution, leading to incorrect results. This issue is not unique to Bessel functions but applies to any recurrence where the ratio of solutions tends to zero or infinity.

Engineers must carefully select computation direction or employ stabilization techniques to mitigate these errors. Miller’s algorithm is one such method, designed to compute minimal (decaying) solutions accurately by working backward from a high-order estimate. Without such measures, applications relying on these recurrences, such as signal processing, physics simulations, or numerical analysis, risk accumulating errors that render results unusable. The choice of method depends on the specific function and the range of n being computed, requiring awareness of the underlying stability characteristics.