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

by Namir Shammas

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

Algorithm

 The algorithm is:

 Jn(x) = (x/2)^n ∑ [(-x^2/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 = -X^2/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 hot 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

BACK

Copyright (c) Namir Shammas. All rights reserved.