Computing Bessel functions in complex interval arithmetic

Question

Does anyone know any resource of algorithms for evaluating the Bessel function of the first kind and second kind under interval-valued domains? In particular, I am interested in complex interval arithmetic. The only paper I have found is:

Zakovic, S.: Evaluation of Bessel-Ricatti Functions using Interval Arithmetic. Working Paper, Numerical Optimisation Centre, University of Hertfordshire, 1994.

which I cannot find for download online.

Thanks.

Answer 1

As I understand it, interval arithmetic maps intervals to intervals. So, if the domain is $[a,b]$, then $J_\alpha([a,b])$ is basically $$ \big[\min_{x\in[a,b]} J_\alpha(x),\max_{x\in[a,b]} J_\alpha(x)\big] $$

That is, you basically evaluate the function on the entire domain, and take the smallest and largest values of the function on that domain, and that is the range interval. Maybe this isn't rigorous, but is right for this case. The same concept applies to the other Bessel functions of interest, so long as you don't include the $Y_\alpha(0)$ case, which will require special handling in what I propose below.

So, if you have a numerical algorithm that can evaluate the Bessel function at points on an arbitrary interval, and you can numerically optimize over that interval (minimize and maximize), then you have an algorithm to do what you want. To give you some concrete ideas, here are some of the many approaches to evaluation for single values. Then, you can, for example, use an adaptive polynomial interpolation algorithm that represents the function on the domain as a polynomial accurate to machine precision. Once you have this representation, maxima and minima can be found by differentiation (a simple operation on the polynomial coefficients) and polynomial rootfinding (equivalent to numerically finding the eigenvalues of a matrix whose characteristic polynomial is the one of interest). For this last part, you'll have to form the something like a companion matrix, which generally is a matrix whose characteristic equation is the polynomial you want to find the roots of. Finally, then there are standard eigenvalue algorithms available.

To make things a little more concrete, if you have Matlab at your disposal, I think you could use the Chebfun package, which uses Chebyshev polynomial interpolation, to quite easily evaluate Bessel functions in interval arithmetic. I imagine it would be less than 10 lines of code, because Matlab has good Bessel function algorithms and Chebfun does the polynomial fitting, differentiation, and rootfinding using eigenvalue algorithms all automatically for you. If you don't have all that at your disposal, you could reimplement these ideas yourself, using the above code that's available as a reference.

Update

The following code in Matlab with Chebfun maps rectangular intervals to rectangular intervals that encompass the actual range.

%order of the bessel function we're working with
nu = 0;
%the complex interval
X = [-5,5];
Y = [-6,6];

f = chebfun2(@(z) besselj(nu,z),[X Y]);

rangeX = [min(min(real(f))) max(max(real(f)))];
rangeY = [min(min(imag(f))) max(max(imag(f)))];

fprintf('\\Re[J_\\nu(z)] \\in [%f,%f]\n',rangeX);
fprintf('\\Im[J_\\nu(z)] \\in [%f,%f]\n',rangeY);

It outputs the following

$$ \Re[J_\nu(z)] \in [-62.328459,67.234407] $$ $$ \Im[J_\nu(z)] \in [-65.803201,65.803201] $$

This is compact, but implements very powerful algorithms, as I described in the rest of this post. My original ideas about minima and maxima apply, but just to the real and imaginary parts of the range.