Bessel Function (4 versions)

[Namir]

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