5

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

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

4

3 回答 3

3

复发_贝塞尔

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

考虑以下比较：

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 年。

于 2011-10-27T05:46:47.147 回答
2

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

所以，你正在做的是从一个较小的数字的倍数中减去一个大的数字，以获得一个更小的数字。从数字上看，这不会很好地工作。

这样看。我们知道结果是10^-25. 您从 1 级的数字开始。因此，为了从中获得一个准确的数字，我们必须知道前两个数字的精度至少为 25 位。我们显然没有，而且复发实际上是发散的。

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

于 2011-10-27T07:38:31.217 回答
0

您可以只修改下面用于球形贝塞尔函数的代码。它经过了很好的测试，适用于所有参数和订单范围。我很抱歉它在 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;
        }
    }
于 2013-03-12T13:20:49.590 回答