Threaded Mode | Linear Mode

Namir · (This post was last modified: 12-24-2013 12:46 PM by 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    1            0        0        0

1    0.997502    0.049938    0.001249    2.08E-05

2    0.990025    0.099501    0.004983    0.000166

3    0.977626    0.148319    0.011166    0.000559

4    0.960398    0.196027    0.019735    0.00132

5    0.93847    0.242268    0.030604    0.002564

6    0.912005    0.286701    0.043665    0.0044

7    0.881201    0.328996    0.058787    0.00693

8    0.846287    0.368842    0.075818    0.010247

9    0.807524    0.40595    0.094586    0.014434

0    0.765198    0.440051    0.114903    0.019563

1    0.719622    0.470902    0.136564    0.025695

2    0.671133    0.498289    0.159349    0.032874

3    0.620086    0.522023    0.183027    0.041136

4    0.566855    0.541948    0.207356    0.050498

5    0.511828    0.557937    0.232088    0.060964

6    0.455402    0.569896    0.256968    0.072523

7    0.397985    0.577765    0.281739    0.08515

8    0.339986    0.581517    0.306144    0.098802

9    0.281819    0.581157    0.329926    0.113423

0    0.223891    0.576725    0.352834    0.128943

1    0.166607    0.568292    0.374624    0.145277

2    0.110362    0.555963    0.395059    0.162325

3    0.05554    0.539873    0.413915    0.179979

4    0.002508    0.520185    0.43098    0.198115

5    -0.04838    0.497094    0.446059    0.2166

6    -0.0968    0.470818    0.458973    0.235294

7    -0.14245    0.441601    0.469562    0.254045

8    -0.18504    0.409709    0.477685    0.272699

9    -0.22431    0.375427    0.483227    0.291093

0    -0.26005    0.339059    0.486091    0.309063

1    -0.29206    0.300921    0.486207    0.326443

2    -0.32019    0.261343    0.483528    0.343066

3    -0.3443    0.220663    0.478032    0.358769

4    -0.3643    0.179226    0.469723    0.373389

5    -0.38013    0.137378    0.458629    0.38677

6    -0.39177    0.095466    0.444805    0.398763

7    -0.39923    0.053834    0.42833    0.409225

8    -0.40256    0.012821    0.409304    0.418026

9    -0.40183    -0.02724    0.387855    0.425044

0    -0.39715    -0.06604    0.364128    0.430171

Code:

n            x        J(n,x)

--------------------------------------------

5        1        0.671397

5        1        0.240298

5        1        0.049497

5        1        0.007186

5        1        0.000807

5        2        0.513016

5        2        0.491294

5        2        0.223925

5        2        0.068518

5        2        0.015887

5        3        0.065008

5        3        0.477718

5        3        0.412710

5        3        0.210132

5        3        0.077598

5        4           -0.301921

5        4        0.185286

5        4        0.440885

5        4        0.365820

5        4        0.199300

5        5           -0.342168

5        5           -0.169651

5        5        0.240377

5        5        0.410029

5        5        0.333663

5        6           -0.091015

5        6           -0.327930

5        6           -0.072950

5        6        0.267139

5        6        0.384612

aurelio · (This post was last modified: 12-25-2013 02:34 PM by aurelio.)

(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..

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.