Implementing MATLAB's besselj function

Hi,

I've been trying to learn how to implement a Bessel function of the 1st kind. I understand that this function is provided in the MATLAB library of functions, "besselj(nu,Z)". The reason I am trying to implement it is to become more familiar with the formula, and how to program formula's in MATLAB.

I have come across the definition of the Bessel function of the 1st kind in various text, and also in the MATLAB documentation of the besselj function. When attempting to implement the formula myself, and then plot against MATLAB's implementation of the function I hit some errors. I'm wondering if anyone here can perhaps explain in a bit more detail how MATLAB have implemented this function?

When I view the graph of the in-built function 'besselj', regardless of the values of 'z', the graph continues fine, as it should. (You have no control of the value 'k' when calling MATLAB's besselj function). When I view the graph of my own implementation of the formula (defined in the MATLAB definition section of the documentation for besselj) I get strange discontinuities for large 'z' values. Whether I use a large or small value for 'k' it never converges to be the same as the output from the besselj function alone. If I have a small upper value of 'z' set to say 20, I can easily get a match between the output of besselj and my own implementation.

I suspect there may be a way of defining an optimum value for 'k' depending on the largest 'z' value?! Can someone perhaps clarify this for me?

Thank you!

Answer 1

[Steven Lord]

Computing special functions like BESSELJ can be tricky, as Cleve describes in one of his Cleve's Corner articles. It's particularly involved for larger values of the inputs, again as Cleve discussed. While the documentation includes a definition of the Bessel function as an infinite sum, that's not necessarily the way we compute it. [I can't really be more specific, as I can't show you the code.] If you edit the besselj.m file, you can see a reference to an article about computing Bessel functions (and some of the issues involved) that may be of interest to you.

[Mark Thompson]

Hi Steven,

Thanks for your response. It would appear then that MATLAB's implementation of this function is a little more complex than just the equation depicted in the besselj documentation. This is a little misleading. I guess as users we just have to "assume" that MATLAB have implemented it correctly and that the results are trustworthy for our work? I'm sure it's fine...

[Steven Lord]

No, you don't have to assume. Trust but verify.

 besselj(0, (0:0.1:5).')

Compare that result to the first column of Table 9.1 in the standard reference manual Abramowitz & Stegun's Handbook of Mathematical Functions. That Wikipedia page includes links where you may legally download it, though it is also available in "dead tree" form from Dover. The results from MATLAB and the results in A&S should agree. Repeat for the rest of the columns in that table and the rest of the tables in that chapter. If A&S and MATLAB agree, they're both right or both wrong, and I claim that the chance that they're both right is much higher than the chance they're both wrong in the exact same way.

If they don't agree, see if A&S has had errata issued for that table and if that errata resolves the discrepancy. If it hasn't, check another source (symbolic computations using Symbolic Math Toolbox, another software package, one of the implementations cited in the companion to the NIST Handbook of Mathematical Functions such as W. J. Cody's FORTRAN code, another text, etc.) and/or contact Technical Support. [I would steer clear of Microsoft Excel as a check, though; a quick test found they only agree in the first six to eight decimal places with A&S and MATLAB.]

[Mark Thompson]

Thanks very much for your advice and suggestions Steven! I will continue exploring using the resources you have provided me.