Matlab custom spherical bessel function NaN at 0

Question

I have a problem involved spherical Bessel functions of order 0. I wrote my own spherical Bessel function:

function js = sphbesselj(nu,x)
js = sqrt(pi ./(2* x)) .* besselj(nu + 0.5, x);
end

which seems to agree with Mathematicas inbuilt one for all my test cases. The problem is at nu or x =0. Mathematica correctly returns 1, but my MATLAB scrip returns NaN. How could I correct my code, so that if I feed in an array of say x = 0:1:5 I get the 1 for my first value, instead of

>> sphbesselj(0,x)
ans =
   NaN    0.8415    0.4546    0.0470   -0.1892   -0.1918

Is my custom function the best way to go about this?

Thanks

Dividing by zero is always going to give you problems. Why not just catch the x=0 case and modify accordingly — mathematician1975, Dec 11 '13 at 09:56

horchler · Accepted Answer · 2013-12-12T16:36:10.777

3

In fact, floating point values of x close to 0 also return Nan. You actually have three cases to possibly worry about. x = ±∞ results in NaN as well. Here's a vectorized function that handles those:

function js = sphbesselj(nu,x)
    if isscalar(nu) && isscalar(x)
        js = 0;
    elseif isscalar(nu)
        js = zeros(size(x));
        nu = js+nu;
    else
        js = zeros(size(nu));
        x = js+x;
    end
    x0 = (abs(x) < realmin);
    x0nult0 = (x0 & nu < 0);
    x0nueq0 = (x0 & nu == 0);
    js(x0nult0) = Inf;          % Re(Nu) < 0, X == 0
    js(x0nueq0) = 1;            % Re(Nu) == 0, X == 0
    i = ~x0nult0 & ~x0nueq0 & ~(x0 & nu > 0) & (abs(x) < realmax);
    js(i) = sign(x(i)).*sqrt(pi./(2*x(i))).*besselj(nu(i)+0.5,x(i));
end

A useful resource when developing such functions is http://functions.wolfram.com. The page on the spherical Bessel function of the first kind has many useful relationships.

edited Dec 12 '13 at 16:36

answered Dec 11 '13 at 17:09

horchler

18,384
4
37
73

Hi, Thanks, this function is very clever and I learnt some new tricks from it. I have run into a new problem now though. The function returns 1 when x is any negative number. I can't quite work out how to edit it to re allow negative values again. Thanks! – Steve Hatcher Dec 12 '13 at 08:38
@SteveHatcher: You're correct. I've updated the code. That issue was fixed by `x0 = (abs(x) < realm in);`. I also made two other small changes to match Mathematica's `SphericalBesselJ`. In particular, the sign of the output seems to depend on the sign of `x`. This doesn't seem to be indicated by the definition of the function, but it makes sense as this is a spherical function. – horchler Dec 12 '13 at 16:40
Awesome thanks so much - that works perfectly! Yes Mathematica does seem to presume more about situations to make things simpler. I am still surprised MATLAB does not have spherical bessel built in functions. – Steve Hatcher Dec 13 '13 at 01:44

Matlab custom spherical bessel function NaN at 0

1 Answers1