The Museum of HP Calculators

HP-71B Program to Calculate Bessel Functions Jn(x)

This program is by Namir Shammas and is used here by permission.

This program is supplied without representation or warranty of any kind. Namir Shammas and The Museum of HP Calculators therefore assume no responsibility and shall have no liability, consequential or otherwise, of any kind arising from the use of this program material or any part thereof.

The following program calculates the Bessel function Jn(x).

Algorithm

 The algorithm is:

 Jn(x) = (x/2)n ∑ [(-x2/4)k /(k! Γ(n + k + 1))]  for k = 0, 1, 2, …, ∞

The P-Code

Give arguments X and N, and tolerance T for terms used in the summation:

  1. Sum = 0
  2. k = 0
  3. p = -X2/4
  4. f1 = 1, f2 = 1, f3 = n!
  5. Term = f1 /(f2 * f3)
  6. Sum = Sum + Term
  7. Increment k
  8. f1 = f1 * p
  9. f2 = f2 * k
  10. f3 = f3 * (n + k)
  11. if |Term| > Tolerance then resume at step 5
  12. Return Jn(x) = (X/2)N * Sum

Test Data

The following table shows data for the Bessel function. Use values in this table to test the program.

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

 Here is an example that shows how to calculate J2(1):

> [RUN]
X? 1[END LINE]
N 2[END LINE]
J = 0.440051  

Here is the BASIC listing:

100 INPUT "X? "; X
110 INPUT "N? "; N
120 S = 0
130 P = -(X ^ 2 / 4)
140 F = 1
150 G = 1
160 REM calculate factorial of N
170 H = 1
180 FOR K = 2 TO N
190 H = K * H
200 NEXT K
210 K = 0
220 REM START LOOP
230 T = F / (G * H)
240 S = S + T
250 IF ABS(T) < 0.000001 THEN 310
260 K = K + 1
270 F = F * P
280 G = G * K
290 H = H * (N + K)
300 GOTO 230
310 DISP "J ="; (X / 2) ^ N * S
320 END

Go back to the software library
Go back to the main exhibit hall