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

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

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，这是完全不同的结果。

最佳答案

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

考虑以下比较:

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之后的计算值开始出现显着差异。

根据MATLAB:



并提供以下内容作为其实现的引用:

关于math - 使用Jm + 1 = 2mj(m)-j(m-1)公式在MATLAB中计算贝塞尔函数，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/7810086/