# HP Forums

Full Version: Bessel Function (4 versions)
You're currently viewing a stripped down version of our content. View the full version with proper formatting.
This document presents a set of programds that implement the Bessel Function Jn(x).

Version 1

The first version is based on an HP-25 program by Peter Henrici. The HP-41C implementation extends the original HP-25 program to calculate the Bessel functions to any non-negative integer order. The original HP-25 calculated J(0,x) through J(4,x). The HP-41C implementation also does that and then determines if the targeted order is in the range of 0 to 4. If, so, the program simply returns the result already available. If the targetd order is 5 and up, the program uses a recurive relation to calculate the higher order function starting with the values of J(3,x) and J(4,x) which are already available.

Algorithm

Given Bessel order N, argument x, and maximum number of iterations m:

Code:
```Let y0 = 1 Let y1, y2, y3, and y4 = 0 Let c = 0 Do   y4 = y3   y3 = y2   y2 = y1   y1 = y0      y0 = 2*N/x*y1 - y2   m = m - 1   If m = 0 then exit Do loop   if m is even then c = c + 2*y0 Loop c = c + y0 J0 = y0 / c J1 = y1 / c J2 = y2 / c J3 = y3 / c J4 = y4 / c ff N <= 4 then   Display J(N) else   i = 4    ' has valid range of 4 to N-1. When it is equal to N we stop iterating   Let J0 = J3   Let J1 = J4   Do     J2 = 2*i/x*J1 - J0     J0 = J1     J1 = J2     i = i + 1   Loop until i = N   Display J2 end if```

Memory Map

R00 = J0
R01 = J1
R02 = J2
R03 = J3
R04 = J4
R05 =
R06 = x
R07 = m
R08 = n
R09 = i

Listing

Code:
``` LBL "BESSL" LBL a "M?" PROMPT                # prompt for max iterations STO 07 LBL A "N^X?" PROMPT                # prompt for Bessel order N and argument x STO 06 X<>Y STO 08 RCL 07 X=0? 50 STO 07 2 STO* 07 0 STO 01 STO 05 1 STO 00 LBL 00 RCL 03 STO 04 RCL 02 STO 03 RCL 01 STO 02 CHS RCL 00 STO 01 RCL 07 RCL 06 / * + STO 00 RCL 07 2 - STO 07 X=0? GTO 01 4 / FRC X#0? GTO 00 RCL 00 2 * STO+ 05 GTO 00 LBL 01 RCL 00 STO+ 05 RCL 05 1/X 4E-3 Rv LBL 02 STO* IND T ISG T GTO 02 4 RCL 08 X>Y?                    # Bessel order > 4 GTO 03                # Calculate Bessel function for order greater than 4 RCL IND X            # Recall result already available in registers 0 to 4 "J=" ARCL X PROMPT                # show result for N <= 4 RTN LBL 03                # Calculate Bessel function for order greater than 4                 4 STO 09                # set i to 4 RCL 03 STO 00                # set J0 = J3 RCL 04 STO 01                # set J1 = J4 LBL 04                # loop to calculate J(i+1,x) RCL 09 2 * RCL 06 / RCL 01 * RCL 00 - STO 02                # calculate J(i+1,x) = 2*i/x*j(i,x) - J(i-1,x) RCL 01                # Shift registers for the next iteration STO 00 RCL 02 STO 01 RCL 09                 1 + STO 09 RCL 08 X>Y?                # Is N > i? GTO 04 RCL 02 "J=" ARCL X PROMPT RTN```

Usage

1. To store the maximum number of iterations, press [f][A]. The program displays the "M?" prompt to ask you for the maximum number of iterations. Enter that value and press [R/S]. If enter 0 or if you skip this value, the program assigns 50 to the maximum number of iterations.
2. To enter the order of the Bessel function and its argument, press [A]. The program displays the "N^X?" prompt to request these values. Enter the order of the Bessel function and its argument, and then press [R/S].
3. The program iterates and then displays the result as "J=<value>".
4. To enter a new function order and argument resume at step 2.

Version 2

The second version uses the definition of Bessel function that is based on an integral. The program uses the INTEG command found in the Advantage Module to perform required numerical integration. Using the INTEG command shortens the code significantly, while gaining speed, and without compromise of the accuracy of the results.

Algorithm

Code:
``` . . . . . . . . . /\ pi                   | J(n,x) = (1/pi) * | cos(n*t - x*sin(t)) dt                    |                 \/ 0```

For n being a non-negative integer.

Memory Map

R00 = x
R01 = n

Listing

Code:
```LBL "BESSL2" LBL A RAD "N^X?" PROMPT STO 00                    # store x X<>Y STO 1                        # store n "BSFX"                    # Set the name of the function to be integrated 0 PI                            # Set the integration limits INTEG                        # Integrate PI /                                # divide by pi to get the result "J=" ARCL X PROMPT RTN LBL "BSFX"            # The function used by the integral ENTER ENTER RCL 01 * X<>Y SIN RCL 00 * - COS RTN```

Usage

1. To enter the order of the Bessel function and its argument, press [A]. The program displays the "N^X?" prompt to request these values. Enter the order of the Bessel function and its argument, and then press [R/S].
2. The program calculates the Bessel function and then displays its value as "J=<value>".

Version 3

The third version uses the definition of Bessel function that is also based on an integral and also works with non-integer Bessel function orders. The program uses the INTEG command found in the Advantage Module to perform required numerical integration. Using the INTEG command shortens the code significantly, while gaining speed, and without compromise of the accuracy of the results.

Algorithm

Code:
``` . . . . . . . . . . . . . . . . . . . . . . . . .. /\ x                                                    | J(n,x) = 1/[2^(n-1) * Gamma(n+0.5)*sqrt(pi)*x^n] * | ((x^2 - t^2)^(n-0.5))*cos(t) dt                                                     |                                                  \/ 0```

For n being the order of the Bessel function.

Memory Map

R00 = x
R01 = n
R02 = argument for calculating the gamma function and for calculating the integral
R03 = used to calculate the result

Listing

Code:
```LBL "BESSL3" LBL A RAD "N^X?" PROMPT STO 00                    # store x X<>Y STO 1                        # store n "BSFX2"                    # Set the name of the function to be integrated 0 RCL 00                    # Set the integration limits INTEG                        # Integrate STO 03                    # Store integral 2 RCL 01 1 - Y^X STO/ 03 RCL 01 0.5 + XEQ 00                    # Calculate gamma(n+0.5) STO/ 03 PI SQRT STO/ 03 RCL 00 RCL 01 Y^X STO/ 03 RCL 03 "J=" ARCL X PROMPT RTN LBL 00                    # Calculate the Gamma function STO 02 0.5 - RCL 02 LN * RCL 02 - 12 RCL 02 * 1/X + RCL 02 3 Y^X 360 * 1/X - RCL 02 5 Y^X 1260 * 1/X + RCL 02 7 Y^X 1680 * 1/X - EXP PI 2 * SQRT * RTN LBL "BSFX2"            # The function used by the integral STO 02 X^2 CHS RCL 00 X^2 + RCL 01 0.5 - Y^X RCL 02 COS * RTN```

Usage

1. To enter the order of the Bessel function and its argument, press [A]. The program displays the "N^X?" prompt to request these values. Enter the order of the Bessel function and its argument, and then press [R/S].
2. The program calculates the Bessel function and then displays its value as "J=<value>".

Version 4

The fourth version uses the definition of Bessel function that uses a summation and one value for the Gamma function. The implementation calculates the Bessel function for any non-negative integer and non-integer order.

Algorithm

J(n,x) = (x/2^n / Gamma(n+1) * sum((-1)^i * (x/2)^(2*i) / (n+1) / (i!), for i = 0 to infinity)

For n being the order of the Bessel function.

Given tolerance, x, and order n:

Code:
```Let y = x/2 a = y / Gamma(n+1) b = 1 s = 1 m = 1 Do   b = b * y^2 / (m * (n + m))   s = s + b   m = m + 1 Loop Until |b| < tolerance J(n,x) = a * s```

Note: In Henrici's book the flow chart shows J being calculated as a * b. This is an error. The book's listing for the HP-25 multiplies the correct operands.

Memory Map

R00 = n
R01 = Gamma(n+1)
R02 = tolerance
R03 = a
R04 = sum
R05 = b
R06 = N
R07 = x/2
R08 = used to calculate the gamma function in LBL 00

Listing

Code:
```LBL "BESSL4" LBL a "TOLER?" PROMPT STO 02 LBL A "N^X?" PROMPT 2 / STO 07                    # store x/2 X<>Y STO 00                    # store n 1 + XEQ 00                    # Calculate gamma(n+1) STO 01                    # store gamma(n+1) RCL 07 RCL 00 Y^X RCL 01 / STO 03 1 STO 04 STO 05 STO 06 LBL 01 RCL 07 X^2 RCL 06 RCL 00 + RCL 06 * / CHS STO* 05 RCL 05 STO+ 04 ABS RCL 02 X>Y? GTO 02 1 STO+ 06 GTO 01 LBL 02 RCL 03 RCL 04 * "J=" ARCL X PROMPT RTN LBL 00                    # Calculate the Gamma function STO 08 0.5 - RCL 08 LN * RCL 08 - 12 RCL 08 * 1/X + RCL 08 3 Y^X 360 * 1/X - RCL 08 5 Y^X 1260 * 1/X + RCL 08 7 Y^X 1680 * 1/X - EXP PI 2 * SQRT * RTN```

Usage

1. To enter the tolerance value press [f][A]. The program uses the "TOLER?" prompt you to ask for the toelrance value. Enter this value and press [R/S]. The program resumes to step 2 without the need to press [A]/
2. To enter the order of the Bessel function and its argument, press [A]. The program displays the "N^X?" prompt to request these values. Enter the order of the Bessel function and its argument, and then press [R/S].
3. The program calculates the Bessel function and then displays its value as "J=<value>".

Test Data
Code:
``` X    J0(x)        J1(x)        J2(x)        J3(x) -------------------------------------------------------------------------- 0.0    1            0        0        0 0.1    0.997502    0.049938    0.001249    2.08E-05 0.2    0.990025    0.099501    0.004983    0.000166 0.3    0.977626    0.148319    0.011166    0.000559 0.4    0.960398    0.196027    0.019735    0.00132 0.5    0.93847    0.242268    0.030604    0.002564 0.6    0.912005    0.286701    0.043665    0.0044 0.7    0.881201    0.328996    0.058787    0.00693 0.8    0.846287    0.368842    0.075818    0.010247 0.9    0.807524    0.40595    0.094586    0.014434 1.0    0.765198    0.440051    0.114903    0.019563 1.1    0.719622    0.470902    0.136564    0.025695 1.2    0.671133    0.498289    0.159349    0.032874 1.3    0.620086    0.522023    0.183027    0.041136 1.4    0.566855    0.541948    0.207356    0.050498 1.5    0.511828    0.557937    0.232088    0.060964 1.6    0.455402    0.569896    0.256968    0.072523 1.7    0.397985    0.577765    0.281739    0.08515 1.8    0.339986    0.581517    0.306144    0.098802 1.9    0.281819    0.581157    0.329926    0.113423 2.0    0.223891    0.576725    0.352834    0.128943 2.1    0.166607    0.568292    0.374624    0.145277 2.2    0.110362    0.555963    0.395059    0.162325 2.3    0.05554    0.539873    0.413915    0.179979 2.4    0.002508    0.520185    0.43098    0.198115 2.5    -0.04838    0.497094    0.446059    0.2166 2.6    -0.0968    0.470818    0.458973    0.235294 2.7    -0.14245    0.441601    0.469562    0.254045 2.8    -0.18504    0.409709    0.477685    0.272699 2.9    -0.22431    0.375427    0.483227    0.291093 3.0    -0.26005    0.339059    0.486091    0.309063 3.1    -0.29206    0.300921    0.486207    0.326443 3.2    -0.32019    0.261343    0.483528    0.343066 3.3    -0.3443    0.220663    0.478032    0.358769 3.4    -0.3643    0.179226    0.469723    0.373389 3.5    -0.38013    0.137378    0.458629    0.38677 3.6    -0.39177    0.095466    0.444805    0.398763 3.7    -0.39923    0.053834    0.42833    0.409225 3.8    -0.40256    0.012821    0.409304    0.418026 3.9    -0.40183    -0.02724    0.387855    0.425044 4.0    -0.39715    -0.06604    0.364128    0.430171```

Code:
```n            x        J(n,x) -------------------------------------------- 0.5        1        0.671397 1.5        1        0.240298 2.5        1        0.049497 3.5        1        0.007186 4.5        1        0.000807 0.5        2        0.513016 1.5        2        0.491294 2.5        2        0.223925 3.5        2        0.068518 4.5        2        0.015887 0.5        3        0.065008 1.5        3        0.477718 2.5        3        0.412710 3.5        3        0.210132 4.5        3        0.077598 0.5        4           -0.301921 1.5        4        0.185286 2.5        4        0.440885 3.5        4        0.365820 4.5        4        0.199300 0.5        5           -0.342168 1.5        5           -0.169651 2.5        5        0.240377 3.5        5        0.410029 4.5        5        0.333663 0.5        6           -0.091015 1.5        6           -0.327930 2.5        6           -0.072950 3.5        6        0.267139 4.5        6        0.384612```
(12-24-2013 12:38 PM)Namir Wrote: [ -> ]This document presents a set of programds that implement the Bessel Function Jn(x).
The first version is based on an HP-25 program by Peter Henrici.

thank-you Namir, I have that book (computational analysis with the HP25 pocket calculator, isn't it?) and I was reading it in train just this morning.....
I started with programs designed for the numeric series...........

I would like to transport any of the programs suited for the HP25 and HP33 to the HP41C....what a big enterprise for a "old" beginner like me..
The older the programmable calculator is the easier it is to translate its commands to a newer one. Just today, I easily keyed in 2 HP-65 programs in an HP-41CX. I had no trouble entering and running the HP-65 program in the HP-41CX.
Reference URL's
• HP Forums: https://www.hpmuseum.org/forum/index.php
• :