+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Software Libraries (/forum-10.html)
+--- Forum: HP-41C Software Library (/forum-11.html)
+--- Thread: Bessel Function (4 versions) (/thread-223.html)
Bessel Function (4 versions) - Namir - 12-24-2013 12:38 PM
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 ifMemory 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
RTNUsage
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
|
\/ 0For 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
RTNUsage
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
|
\/ 0For 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
*
RTNUsage
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 * sNote: 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
*
RTNUsage
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.430171Code:
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.384612RE: Bessel Function (4 versions) - aurelio - 12-25-2013 02:31 PM
(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..
RE: Bessel Function (4 versions) - Namir - 12-28-2013 04:41 AM
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.