Calculate the Bessel function in MATLAB using the formula Jm + 1 = 2mj (m) -j (m-1)

Question

I tried to implement the bessel function using this formula, this is the code:

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;

But if I use the Bessel MATLAB function to compare with this, I get different values ​​that are too high. For example, if I find Bessel (20), it gives me the result 3.1689e + 005, if instead I type bessel (20,1), it gives me 3.8735e-025, a completely different result.

Answer 1

such recurrence relations are good in mathematics, but are numerically unstable when implementing algorithms using representations with limited precision of floating point numbers.

:

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 MEX Fortran . . .

:

. . , " ", , SAND85-1018, 1985 .

. . , " ", Trans. Math. Software, 1986.

Answer 2

, , . , , , J (n, x) 1/2n. , J.

, , , . .

. , 10^-25. , 1. , , 25 . , .

, , , . J (20,1) J (19,1), 0 . ? . , .

Answer 3

You can simply change the code below, which is intended for the spherical bessel function. It is well tested and works for all arguments and order. I'm sorry that it is in 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;
        }
    }