Speeding up the Bessel function

The Bessel function is used in a few scattering functions, for example in that of a cylinder. The implementation of the Bessel function in Matlab is quite slow, however, which annoyed me considerably as it slowed down my fitting functions. I therefore had a look at whether I could speed it up.

The paper by Frank B. Gross entitled “New approximations to J0 and J1 Bessel functions” (IEEE transactions on antennas and propagation vol. 43, no. 8, 1995, pp.904) seemed like a good place to start. It is mentioned there that there are two fast approximations already available, i.e. the “small argument approximation”, valid for x<0.5 and the “large argument approximation”, valid for x>5, which are very quick functions indeed. However, for 0.5<x<5, the only approximation available was a Taylor series.

dr. Gross then comes up with a new approximation expressed as a polynomial with integer coefficients. The only drawback of this implementation, is that for the first order Bessel function J1, it requires the calculation of the derivative at x, and would therefore require either more points, or would put stringent requirements on the order and spacing of x.

In a later paper by Millane and Eads (IEEE transactions on antennas and propagation, 2003, vol. 51, no. 6, pp. 1398) the approximation derived by Gross is further generalised to encompass all Bessel functions of integer order. More interestingly, it does away with the derivation. However, I could not get their approximation to result in the correct values.

The last paper that “solved” the issue is a paper from 2006 by Li (and Gross, Applied mathematics and computation 2006, 183, pp. 1220), which generalises the generalisation by Millane and Eads to obtain an approximation for the Bessel functions for all fractional order. This contains an equation which provides a good and quick approximation to the Bessel function.

When implemented, however, it could not beat the internal function in Matlab. The speed reached by the approximation was virtually identical to that of the internal Matlab function. Now, while this was not a great success, it can be considered a partial success, as my interpreted code was as fast as the optimized, internal code. Of course, it remains an approximation…

Finally, I tried simply using a linear interpolation between several results of the Bessel function between 0<x<5, and using the large argument approximation beyond that. This gives moderately good results, as it is two to three times faster as the internal Matlab function. While not much faster, it has certainly taught me a lot about Bessel function implementations. If there is anyone among you out there with knowledge on how to compile Matlab functions into optimized code, please let me know.

Cheers,

Brian.