Search code examples
mathmatlabrecurrencebessel-functions

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

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.

Solution

  • recurrence_bessel

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

    comparison_plot

    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.