Robust ways of evaluating $j_n(x+iy)/e^y$

Question

For reasons to do with extending the range of applicability of this Mathematica package, I would like to have a good handle on the numerical evaluation of functions of the form $$ f(x) = \frac{j_n(\sqrt{x^2-\zeta^2})}{e^{\zeta}}, $$ where $j_n$ is a spherical Bessel function of the first kind, $x$ is real-valued and reasonably small (say, in the range 0-10 or so) and $\zeta$ is a large real number (larger than $x$, and say, around 50 to 100).

From a mathematical point of view, $f(x)$ isn't particularly problematic and it does not involve any exponential catastrophes, since for $\zeta \gg x$ we have $\sqrt{x^2-\zeta^2} \sim i(\zeta-x^2/2\zeta)$, and therefore the Bessel-function term, $$ j_n(\sqrt{x^2-\zeta^2}) \sim j_n(i(\zeta-x^2/\zeta)) %=i^n i_n^{(1)}(\zeta-x^2/2\zeta) \sim \frac{ i^n }{ 2\zeta-x^2/\zeta } e^{\zeta-x^2/\zeta}, $$ contains an asymptotic contribution of order $e^\zeta$ that cancels out with the denominator in $f(x)$.

.... however.

When I try to implement this in Mathematica, I run into its General::munfl error, which informs me that when trying to evaluate Exp[-ζ] the number is "is too small to represent as a normalized machine number; precision may be lost". Mathematica then sets $e^{-\zeta}$ to zero and kills the entire calculation.

What are robust ways to get around this problem, and to design a numerical approach that will handle this function well in this limit? Is there a well-established way to extract the exponential asymptotic from the Bessel function and leave behind a well-behaved remnant?

Answer 1

The evaluation of $f(x) = e^{-\zeta}j_{n}(\sqrt{x^{2}-\zeta^{2}})$ should not run into numerical issues when using a verbatim translation of the formula for $x$ in $[0,10]$ and $10 \lt \zeta \lt 700$ when evaluated in double precision (which is typically mapped to the IEEE-754 binary64 format on modern systems), as $\exp(-700)$ and $j_n(700i)$ are both representable as double-precision operands. Larger values of $\zeta$ can be accommodated as follows.

For the real numbers $x \lt \zeta$, we know that $\sqrt{x^{2}-\zeta^{2}}$ is an imaginary number $iy$, $y=\sqrt{\zeta^2-x^2}$. The spherical Bessel function $j_{0}$ can be used to compute all spherical Bessel functions of integer order $j_{n}$. In this case $j_0(iy)$ is real-valued, and for sufficiently large $y$, $$j_{0}(iy) = \frac{\sin iy}{iy} \rightarrow\frac{e^y}{2y}$$ It remains the computation of the ratio of spherical Bessel functions $j_n/j_0$ in order to compute $j_n(iy)$ from $j_0(iy)$. The following paper demonstrates that this can be accomplished accurately and reasonably efficiently with a concise continued fraction computation:

William J. Lentz, "Continued fraction calculation of spherical Bessel functions," Computers in Physics, Vol. 4, No. 4, Jul/Aug 1990, pp. 403-407 (online)

While Lentz's code guards against overflow in the continued fraction computation by using a separate scale factor, it does not guard against overflow in the computation of $j_0$. Since for large $y$ (some quick experiments suggest that $y \gt 20$ is sufficient for double-precision computation) we can express $j_0(iy)$ in the exponential form shown above, this is a trivial modification, by combining $e^{y}e^{-\zeta}$ into $e^{y-\zeta}$. A numerically advantageous way of computing $y-\zeta$ when $\zeta > 2x$ is $$y-\zeta = \frac{-x^2}{\zeta + \zeta \sqrt{1-\left(\frac{x}{\zeta}\right)^2}}$$