math - 在MATLAB中使用Jm+1=2mj(m) -j(m-1)公式计算bessel函数

Question

我尝试使用该公式实现贝塞尔函数，这是代码:

function result=Bessel(num);


if num==0
    result=bessel(0,1);
elseif num==1
    result=bessel(1,1);
else
    result=2*(num-1)*Bessel(num-1)-Bessel(num-2);
end;

但是如果我使用 MATLAB 的 bessel 函数将它与这个函数进行比较，我会得到太大的不同值。例如，如果我输入 Bessel(20) 它会给我 3.1689e+005 作为结果，如果我输入 bessel(20,1) 它会给我 3.8735e-025 ，这是一个完全不同的结果。

Answer 1

这种递归关系在数学上很好，但在使用 float 的有限精度表示实现算法时在数值上不稳定。

考虑以下比较:

x = 0:20;
y1 = arrayfun(@(n)besselj(n,1), x);   %# builtin function
y2 = arrayfun(@Bessel, x);            %# your function
semilogy(x,y1, x,y2), grid on
legend('besselj','Bessel')
title('J_\nu(z)'), xlabel('\nu'), ylabel('log scale')

因此您可以看到计算值在 9 之后开始出现显着差异。

根据 MA​​TLAB:

BESSELJ uses a MEX interface to a Fortran library by D. E. Amos.

并给出以下作为其实现的引用:

D. E. Amos, "A subroutine package for Bessel functions of a complex argument and nonnegative order", Sandia National Laboratory Report, SAND85-1018, May, 1985.

D. E. Amos, "A portable package for Bessel functions of a complex argument and nonnegative order", Trans. Math. Software, 1986.