使用JM+1 = 2MJ（M）-J（M -1）公式计算MATLAB中的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

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

考虑以下比较：

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

并将以下作为其实施的参考文献：

德阿莫斯（De Amos），“复杂论证和非负秩序贝塞尔（Bessel）功能的子例程包”，桑迪亚国家实验室报告，Sand85-1018，1985年5月。

de amos，“用于复杂参数和非负顺序的贝塞尔函数的便携式软件包”，trans。数学。软件，1986年。

Answer 2

您使用的正向复发关系不稳定。要查看为什么，请考虑Besselj（N，X）的值越来越小于一个因素 1/2n. 。您可以通过查看J. Taylor系列的第一学期来看到这一点。

因此，您要做的是从较小的数字的倍数中减去一个大数，以获取更小的数字。从数字上讲，这不会很好地工作。

这样看。我们知道结果是 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;
        }
    }