Question about how HP Math pacs calculate Bessel Functions!

[Namir]

Hello All,

If you look in your HP-65 and HP-67 Math Pac manuals and the HP-41 High Level Math solution book, you may come across programs that calculate the Jn(x) and In(x) Bessel functions. These programs use recursive algorithms that calculate series of a specified number of terms m. The documentation calculates m using an interesting formula:

m = 2 * INT[(6 + max(n,z) + 9z/(z + 2))/2]

Where n is the order of the Bessel function and z = 3x/2.

My question, can anyone find the reference for HP's algorithms to calculate the Jn(x) and In(x) Bessel functions using the above equation?? The HP documentation mentions the "Handbook of Mathematical Functions" by Abramowitz et al. as a reference. I looked in that book and and the closest thing I found was at the bottom of page 385 (under the heading Numerical Methods). The handbook discusses informally the algorithm HP used but does not mention m and how it is calculated. Did I miss something? Or is the equation for m something HP came up with?

I guess you can say that this is a math challenge, with a new twist!

Namir

Edited: 3 Nov 2009, 4:20 p.m.

[Frank Rottgardt]

I don´t know which edition of Abramowitz you have excess to. This source cites the 1972 edition which would fit the HP-67 timewise.

MathWorld

Abramowitz, M. and Stegun, I. A. (Eds.). "Spherical Bessel Functions." §10.1 in Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, 9th printing. New York: Dover, pp. 437-442, 1972.

[Namir]

I also have the blue-soft cover version. On page 385 (Bessel functions of integer orders) the book discusses in general terms basically the same algorithm that HP uses, except there is no mention for of m and the equation to calculate it. The discussion does mention the use of recursive iteration to calculate non-normalized Bessel function values, and then briefly hints at how to normalize them.

The HP documentation is much clearer. I was able to implement the algorithm in Excel VBA in two versions: using arrays and without arrays. Both yield accurate results.

Namir

[Ángel Martin]

There's programs in the SandMath-5 module to calculate Jn, In, Kn, and Yn for both integer and non-integer indexes, and any real number. I used the usual approach given the relatively fast convergenge of the series.

Cheers,
ÁM

Edited: 4 Nov 2009, 3:35 a.m.

[Namir]

Angel,

Your SandMath ROM looks VERY nice!!!

You calculate the Bessel functions Jn(x), In(x), Kn(x), and Yn(x) in your SandMath ROM. Do you base your calculations for ALL FOUR Bessel functions on an algorithm that is a similar one that HP uses? If yes, how do you calculate Kn(x) and Yn(x)? Sorry for being nosy.

Cheers,

Namir

Edited: 4 Nov 2009, 7:03 a.m.

[Ángel Martin]

Hi Namir,

No complex stuff: it's basically just about summing terms until there's no significant contribution to the sum value. The process is sped-up by using a couple of MCODE functions involved in the function term calculation (like obviously GAMMA for non-integer orders, HARMN for the harmonic number, X=YR? for rounded comparisons, CHSYX, etc.), and a couple of re-arranging tricks of the series syntax to optimize the efficiency. It still takes its time, but this is a non-issue if you use V41's TURBO mode of course - it feels instantaneous response.

I considered an ALL-MCODE implementation but so far had no time to implemented, maybe some day...

I'm on the road for a couple of weeks but can send you more details when I'm back if you give me an email address to send it to.

Best,
ÁM.

Edited: 5 Nov 2009, 12:30 a.m.