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

Question

I tried to implement bessel function using that 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 MATLAB's bessel function to compare it with this one, I get too high different values. For example if I type Bessel(20) it gives me 3.1689e+005 as result, if instead I type bessel(20,1) it gives me 3.8735e-025 , a totally different result.

Answer 1

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

Consider the following comparison:

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')

So you can see how the computed values start to differ significantly after 9.

According to MATLAB:

BESSELJ uses a MEX interface to a Fortran library by D. E. Amos.

and gives the following as references for their implementation:

D. E. Amos, "A subroutine package for Bessel functions of a complex argument and nonnegative order", Sandia National Laboratory Report, SAND85-1018, May, 1985.

D. E. Amos, "A portable package for Bessel functions of a complex argument and nonnegative order", Trans. Math. Software, 1986.

Answer 2

The forward recurrence relation you are using is not stable. To see why, consider that the values of BesselJ(n,x) become smaller and smaller by about a factor 1/2n. You can see this by looking at the first term of the Taylor series for J.

So, what you're doing is subtracting a large number from a multiple of a somewhat smaller number to get an even smaller number. Numerically, that's not going to work well.

Look at it this way. We know the result is of the order of 10^-25. You start out with numbers that are of the order of 1. So in order to get even one accurate digit out of this, we have to know the first two numbers with at least 25 digits precision. We clearly don't, and the recurrence actually diverges.

Using the same recurrence relation to go backwards, from high orders to low orders, is stable. When you start with correct values for J(20,1) and J(19,1), you can calculate all orders down to 0 with full accuracy as well. Why does this work? Because now the numbers are getting larger in each step. You're subtracting a very small number from an exact multiple of a larger number to get an even larger number.

Answer 3

You can just modify the code below which is for the Spherical bessel function. It is well tested and works for all arguments and order range. I am sorry 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;
        }
    }