(15C)Bessel Functions, Arbitrary Order for HP-15C

[Namir]

Bessel Function J, Arbitrary Order for HP-15C

Based on the HP-25 code by Peter Henrici in his book "Computational Analysis with the HP-25 Pocket Calculator", page 253. The pdf file for the book is available for download in Eric rechlin's librarry. The books shows the equations used and the flowchart for the program.

The program uses the following equations:

Jv(x) = a * sum((b(i)),for n= 0 to infinity

where,

a = (x/2)^v / gamma(v+1)
b(0) = 1
b(n) = -(x/2)^2/n/(v+n)*b(n-1) for n = 1, 2, 3, ...

The range for x values is 0<= x <= 10, according to Henrici.


Memory Map
==========

Code:

R0 = v
R1 = gamma(v+1)
R2 = 1e-12
R3 = a
R4 = sum
R5 = b
R6 = n
R7 = x/2

Program Listing
===============

Code:


01 LBL A
02 STO 0
03 X!
04 STO 1
05 X<>Y
06 2
07 /
08 STO 7
09 EEX
10 CHS
11 1
12 2
13 STO 2
14 RCL 7
15 RCL 0
16 Y^X
17 RCL 1
18 /
19 STO 3
20 1
21 STO 4
22 STO 5
23 STO 6
24 LBL 0
25 RCL 7
26 X^2
27 RCL 6
28 RCL 0
29 +
30 RCL 6
31 *
32 /
33 CHS
34 STO* 5
35 RCL 5
36 STO+ 4
37 ABS
38 RCL 2
39 X>=Y? (Test 9)
40 GTO 1
41 1
42 STO+ 6
43 GTO 0
44 LBL 1
45 RCL 3
46 RCL 4
47 *
48 RTN

Example
=======
To calculate J<1.5>(5), perform the following:

1) Set calculate to [FIX] [8].
2) Enter 5 (for x) and press [ENTER].
3) Enter the order of 1.5 and the press [f][A].
4) The program performs the calculations and stops displaying -0.16965131.

[Eddie W. Shore]

Question, Namir: Should line 38 be RCL 1 or RCL 2 instead?

[Namir]

(11-25-2023 12:32 AM)Eddie W. Shore Wrote:  Question, Namir: Should line 38 be RCL 1 or RCL 2 instead?

It should be RCL 2 (which I changed to). Thanks for catching the source code typo my friend. I also corrected other typos. Sorry!

Namir

[Eddie W. Shore]

(11-25-2023 02:09 AM)Namir Wrote:  

(11-25-2023 12:32 AM)Eddie W. Shore Wrote:  Question, Namir: Should line 38 be RCL 1 or RCL 2 instead?

It should be RCL 2 (which I changed to). Thanks for catching the source code typo my friend. I also corrected other typos. Sorry!

Namir

The program works! Thank you Namir!

[Namir]

No, thank YOU Eddie!!!

[C.Ret]

Thanks to you two.

I wanted to check the results that Namir obtains with those that I obtain with a code that I used before.



I then realized that my program for calculating the values of \( B_n(x) \) only works for integer values of \(n\).
I may therefore have to modify it by adding an INT instruction before storing \(n\) in the R0 register.

My code uses the integration function, which is more explicit than the nifty algorithm found by Namir, but takes much more time. It is based on the following relationship:



Obviously, digital integration takes a very long time on my HP-15C from 1985. But it's no longer very serious, today only the old HP-15Cs are struggling. The time factor is no longer a concern on the HP-15C LE or CE and other DM-15s.

As the calculation time no longer matters, I can try to use a second integration for the corrective term when calculating \(B_\alpha(x) \) with non-integer \(\alpha\).



As a short circuit, my code tests the fractional part of the argument x and jumps to label 2 if it does not exist.
The second integral grow up to infinity ! But don't be crazy, note that \(sinh(9)\simeq4\;052\) is already almost too large since already \(e^{-4\;052}\ll10^{-10}\) and 'larger' infinity limit may lead to erroneous integration by the \(\int_y^x\) operator.

The calculations are very long, fortunately, I took the precaution of saving the result in register R2 because most of the time, my HP-15C turn off before I wake up to read the result. As for the first version, the argument \(x\) is stored in register R1 and the order \(n\) or \(\alpha\) of the Bessel function in register R0.


The code completed by the calculation of the two integrations is not much shorter nor more implicit than the code proposed by Namir. I give it here just to illustrate my points, I don't think it will be of use to anyone. Especially since it takes an infinitely long time to display the slightest result.



Furthermore, integration towards infinity is not, by far, the fastest operation that can be performed on an HP-15C. Using an iterative algorithm like the one Namir propose is therefore an excellent and fully justified idea.


Thank all of you for sharing this code which is very effective and efficient, especially for long-term users of an HP-15c from 1985 as me.