1

我有一个涉及 0 阶球面贝塞尔函数的问题。我编写了自己的球面贝塞尔函数：

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

这似乎与我所有测试用例的 Mathematicas 内置一致。问题在nu或x =0。Mathematica 正确返回 1，但我的 MATLAB 脚本返回 NaN。我怎样才能更正我的代码，这样如果我输入一个say数组，x = 0:1:5我的第一个值就会得到1，而不是

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

我的自定义函数是解决这个问题的最佳方法吗？

谢谢

4

1 回答 1

3

x实际上，接近的浮点值0也返回Nan。实际上，您可能要担心三种情况。x = ±∞ 的结果NaN也是如此。这是一个处理这些的矢量化函数：

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

开发此类功能时的一个有用资源是http://functions.wolfram.com。第一类球面贝塞尔函数的页面有很多有用的关系。

于 2013-12-11T17:09:20.463 回答