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

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?

• A good place to start is to investigate how people implement the evaluation of the Bessel function itself. If you understand how that is done, then you will also be able to implement $f(x)$ in a similar way. In your case, the approach will likely involve some sort of Taylor expansion or, more likely a Taylor expansion of $f(x) e^{x^2/\zeta}$ which you then multiply by $e^{-x^2/\zeta}$. Commented Mar 19, 2022 at 3:31

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:
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}}$$