我尝试使用该公式实现贝塞尔函数,这是代码:
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 ,这是一个完全不同的结果。
最佳答案
这种递归关系在数学上很好,但在使用 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 之后开始出现显着差异。
根据 MATLAB:
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.
关于math - 在MATLAB中使用Jm+1=2mj(m) -j(m-1)公式计算bessel函数,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/7810086/