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

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(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. Math. Software，1986。

Answer 2

您使用的前向递归关系不稳定。要了解原因，请考虑 BesselJ(n,x) 的值变得越来越小，大约减小了 1/2n 倍。您可以通过查看 J 的泰勒级数的第一项来看到这一点。

因此，您所做的就是从一个稍小的数字的倍数中减去一个大的数字，以获得一个更小的数字。从数字上看，这不会很好地发挥作用。

这么看吧。我们知道结果的数量级为10^-25。您从 1 数量级的数字开始。因此，为了从中获得哪怕一位准确的数字，我们必须知道前两个数字，精度至少为 25 位。我们显然不这么认为，而且重复率实际上是有分歧的。

使用相同的递推关系，从高阶到低阶，是稳定的。当您从 J(20,1) 和 J(19,1) 的正确值开始时，您也可以完全准确地计算低至 0 的所有阶数。为什么这有效？因为现在每一步的数字都在变大。您从一个较大数字的精确倍数中减去一个非常小的数字，以获得一个更大的数字。

Answer 3

您只需修改下面用于球面贝塞尔函数的代码即可。它经过充分测试，适用于所有参数和顺序范围。抱歉，这是 C# 语言

     public static Complex bessel(int n, Complex z)
    {
        if (n == 0) return sin(z) / z;
        if (n == 1) return sin(z) / (z * z) - cos(z) / z;

        if (n <= System.Math.Abs(z.real))
        {
            Complex h0 = bessel(0, z);
            Complex h1 = bessel(1, z);
            Complex ret = 0;
            for (int i = 2; i <= n; i++)
            {
                ret = (2 * i - 1) / z * h1 - h0;
                h0 = h1;
                h1 = ret;
                if (double.IsInfinity(ret.real) || double.IsInfinity(ret.imag)) return double.PositiveInfinity;
            }
            return ret;
        }
        else
        {
            double u = 2.0 * abs(z.real) / (2 * n + 1);

            double a = 0.1;
            double b = 0.175;
            int v = n - (int)System.Math.Ceiling((System.Math.Log(0.5e-16 * (a + b * u * (2 - System.Math.Pow(u, 2)) / (1 - System.Math.Pow(u, 2))), 2)));

            Complex ret = 0;
            while (v > n - 1)
            {
                ret = z / (2 * v + 1.0 - z * ret);
                v = v - 1;
            }


            Complex jnM1 = ret;
            while (v > 0)
            {
                ret = z / (2 * v + 1.0 - z * ret);
                jnM1 = jnM1 * ret;
                v = v - 1;
            }
            return jnM1 * sin(z) / z;
        }
    }