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

Question

我尝试使用该公式实现bessel函数,这是代码:

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函数将它与这个函数进行比较,我会得到太高的不同值.例如,如果我键入贝塞尔(20)它给我3.1689e + 005作为结果,如果我输入bessel(20,1)它给我3.8735e-025,一个完全不同的结果.

Answer 1

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

考虑以下比较：

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 使用 DE Amos 的 Fortran 库的 MEX 接口。

并给出以下实施参考：

DE Amos，“复杂参数和非负阶贝塞尔函数的子程序包”，桑迪亚国家实验室报告，SAND85-1018，1985 年 5 月。

DE Amos，“复杂参数和非负阶贝塞尔函数的可移植包”，Trans。数学。软件，1986。