Bessel

Bessel Functions for the HP-41

Overview

1°) Bessel & modified Bessel Functions of the first kind J_n(x) & I_n(x)
2°) Bessel Functions J_n(x) of integer order ( n >= 0 )
3°) Bessel & modified Bessel Functions of the second kind Y_n(x) & K_n(x) ( non-integer order )
4°) Bessel & modified Bessel Functions of the second kind Y_n(x) & K_n(x) ( integer order )
5°) Continued Fractions for J_n(x) & Y_n(x)
6°) An Asymptotic Expansion for J_n(x) and Y_n(x)
7°) An Asymptotic Expansion for K_n(x)
8°) Integrals of J_n(x) & I_n(x)
9°) Integral of J_n(x) ( integer order )
10°) Spherical Bessel Functions

a) Spherical Bessel Functions of the first kind
b) Spherical Bessel Functions of the second ( and first ) kind

11°) Zeros of Bessel Functions, Mac Mahon's Asymptotic Expansions

a) Zeros of J_n(x) & Y_n(x)
b) Zeros of the first derivatives J'_n(x) & Y'_n(x)

12°) Bessel Functions - Complex order & argument

   a) Bessel Functions of the 1st kind
   b) Bessel Functions of the 2nd kind
   c) Bessel Functions of the 2nd kind ( real integer order )

1°) Bessel & modified Bessel Functions of the first kind J_n(x) & I_n(x)

Formulae:

J_n(x) = (x/2)ⁿ [ 1/Gam(n+1) + (-x²/4)¹/ (1! Gam(n+2) ) + .... + (-x²/4)^k/ (k! Gam(n+k+1) ) + .... ] n # -1 ; -2 ; -3 ; ....

I_n(x) = (x/2)ⁿ [ 1/Gam(n+1) + (x²/4)¹/ (1! Gam(n+2) ) + .... + (x²/4)^k/ (k! Gam(n+k+1) ) + .... ] n # -1 ; -2 ; -3 ; ....

Data Registers: /
Flag: F01
Subroutine: "GAM" or "GAM+" ( cf "Gamma Function for the HP-41" )

-Synthetic registers N & O may be replaced by any unused data registers.

01 LBL "JNX"
02 2
03 /
04 STO N
05 X<>Y
06 STO O
07 1
08 ENTER^
09 STO M

10 ENTER^
11 LBL 01
12 X<> Z
13 FC? 01
14 CHS
15 RCL N
16 X^2
17 *
18 RCL M

19 ST/ Y
20 RCL O
21 +
22 /
23 ISG M
24 CLX
25 ST+ Y
26 X<> Z
27 X#Y?

28 GTO 01
29 RCL N
30 RCL O
31 Y^X
32 *
33 X<> O
34 1
35 +
36 XEQ "GAM"

37 RCL O
38 X<>Y
39 /
40 CLA
41 END

( 70 bytes / SIZE 000 )

STACK	INPUTS	OUTPUTS
Y	n	f(x)
X	x	f(x)

where f(x) = J_n(x) if flag F01 is clear and f(x) = I_n(x) if flag F01 is set.

Example: Calculate J_0.7(1.9) and I_0.7(1.9)

   CF 01
   0.7   ENTER^
   1.9   XEQ "JNX"   yields   J_0.7(1.9) = 0.584978102   ( in 12 seconds )

   SF 01
   0.7   ENTER^
   1.9   XEQ "JNX"   yields   I_0.7(1.9) = 1.727630598   ( in 12 seconds )

Notes:

-If n is a positive integer, J_-n(x) = (-1)ⁿ J_n(x) and I_-n(x) = I_n(x)
J₀(0) = 1 but this program doesn't work if x = n = 0 ( DATA ERROR line 31 )
-Add X=Y? X#0? GTO 00 SIGN X<>Y LBL 00 after line 30 and it will work in this case too.

-If you have the M-Code routine "HGF+" in your HP-41 ( cf "Hypergeometric Functions" ), the routine may be simplified as follows:
-Line 21 is - if possible - an M-Code routine that gives 0^0 = 1, see for instance "A few M-Code Routines for the HP-41"

01 LBL "JNX"
02 X<>Y
03 STO 02
04 1
05 +
06 STO 01
07 CLX
08 2
09 /
10 STO 00
11 CLST
12 SIGN
13 CHS
14 RCL 00
15 X^2
16 FC? 01
17 CHS
18 HGF+
19 RCL 00
20 RCL 02
21 Y^X
22 *
23 END

( 33 or 34 bytes / SIZE 003 )

STACK	INPUTS	OUTPUTS
Y	n	/
X	x	f(x)

where f(x) = J_n(x) if flag F01 is clear and f(x) = I_n(x) if flag F01 is set.

Example:

   CF 01
   0.7   ENTER^
   1.9   XEQ "JNX"   yields   J_0.7(1.9) = 0.584978103   ( in 4 seconds )

   SF 01
   0.7   ENTER^
   1.9   XEQ "JNX"   yields   I_0.7(1.9) = 1.727630603   ( in 4 seconds )

-However, with or without "HGF+", J_n(x) is not accurately computed for large arguments, "JNX1" is better:

2°) Bessel Functions J_n(x) of integer order ( n >= 0 )

-The following program is only a modification of a program given by Keith Jarett in "HP-41 Extended Functions made easy"
-It uses the relations: J₀(x) + 2 J₂(x) + 2 J₄(x) + 2 J₆(x) + ...... = 1 and J_n-1(x) + J_n+1(x) = (2n/x) J_n(x)

Data Registers: R00 and R02 to R05: temp R01 = x
Flags: /
Subroutines: /

01 LBL "JNX1"
02 X#0?
03 GTO 00
04 X#Y?
05 RTN
06 SIGN
07 RTN
08 LBL 00
09 STO 01
10 ABS

11 5
12 +
13 X<>Y
14 STO 00
15 X<Y?
16 X<>Y
17 INT
18 ST+ X
19 STO 03
20 ST- 00

21 CLST
22 STO 04
23 STO 05
24 SIGN
25 LBL 01
26 STO Z
27 RCL 03
28 ST+ X
29 *
30 RCL 01

31 /
32 X<>Y
33 -
34 ISG 00
35 STO 02
36 ST+ 04
37 RCL 05
38 X<> 04
39 STO 05
40 RDN

41 DSE 03
42 GTO 01
43 RCL 05
44 ST+ X
45 X<>Y
46 -
47 RCL 02
48 X<>Y
49 /
50 END

( 70 bytes / SIZE 006 )

STACK	INPUTS	OUTPUTS
Y	n	/
X	x	J_n(x)

Examples:

2 ENTER^
10 XEQ "JNX1" produces J₂(10) = 0.254630314 ( in 17 seconds )

3 ENTER^
100 R/S >>>> J₃(100) = 7.628420178 10^-2 ( in 2mn01s )

-Unlike "JNX" , "JNX1" yields accurate results for large x-values

-If you replace lines 45 thru 49 by RCL Y - ST/ 02 ST/ Z / RCL 02
you'll get J₁(x) in Z-register , J₀(x) in Y-register and J_n(x) in X-register.

3°) Bessel & modified Bessel Functions of the second kind Y_n(x) & K_n(x) ( non-integer order )

Formulae: Y_n(x) = ( J_n(x) cos(n(pi)) - J_-n(x) ) / sin(n(pi)) ; K_n(x) = (pi/2) ( I_-n(x) - I_n(x) ) / sin(n(pi)) n # .... -3 ; -2 ; -1 ; 0 ; 1 ; 2 ; 3 ....

Data Registers: R00 = x , R01 = n , R02 = temp

Flag: F01 If F01 is clear, "YNX" calculates Y_n(x)
If F01 is set, "YNX" --------- K_n(x)
Subroutine: "JNX"

01 LBL "YNX"
02 DEG
03 STO 00
04 X<>Y
05 STO 01
06 CHS
07 X<>Y

08 XEQ "JNX"
09 STO 02
10 RCL 01
11 RCL 00
12 XEQ "JNX"
13 PI
14 R-D

15 RCL 01
16 *
17 STO Z
18 FS? 01
19 CLX
20 COS
21 *

22 RCL 02
23 -
24 X<>Y
25 SIN
26 /
27 PI
28 2

29 /
30 CHS
31 FC? 01
32 ST/ X
33 *
34 END

( 54 bytes / SIZE 003 )

STACK	INPUTS	OUTPUTS
Y	n	/
X	x	f(x)

with f(x) = Y_n(x) if CF 01 ; f(x) = K_n(x) if SF 01

Example: Evaluate Y_1.4(3) and K_1.4(3)

CF 01 1.4 ENTER^ 3 XEQ "YNX" >>>> Y_1.4(3) = 0.137821836 ( in 31 seconds )
SF 01 1.4 ENTER^ 3 R/S >>>> K_1.4(3) = 0.046088036 --------------

4°) Bessel & modified Bessel Functions of the second kind Y_n(x) & K_n(x) ( integer order )

Formulae: with psi(x) = the digamma function, we have ( positive n ):

Y_n(x) = -(1/pi) (x/2)^-n SUM_{k=0,1,....,n-1} (n-k-1)!/(k!) (x²/4)^k + (2/pi) ln(x/2) J_n(x) - (1/pi) (x/2)ⁿ SUM_k=0,1,..... ( psi(k+1) + psi(n+k+1) ) (-x²/4)^k / (k!(n+k)!)

K_n(x) = (1/2) (x/2)^-n SUM_{k=0,1,....,n-1} (n-k-1)!/(k!) (-x²/4)^k - (-1)ⁿ ln(x/2) I_n(x) + (1/2) (-1)ⁿ (x/2)ⁿ SUM_k=0,1,..... ( psi(k+1) + psi(n+k+1) ) (x²/4)^k / (k!(n+k)!)

Data Registers: R00 = x/2 ; R01 = n ; R02 to R05: temp

Flag: F01 If F01 is clear, "YNX1" calculates Y_n(x)
If F01 is set, "YNX1" --------- K_n(x)

Subroutine: "JNX"

01 LBL "YNX1"
02 STO 00
03 2
04 ST/ 00
05 X<> Z
06 STO 01
07 X<>Y
08 XEQ "JNX"
09 RCL 00
10 X^2
11 LN
12 *
13 1
14 FS? 01
15 CHS
16 RCL 01
17 Y^X
18 *
19 STO 02
20 RCL 01
21 STO 03
22 CLX

23 LBL 01
24 RCL 01
25 RCL 03
26 -
27 FACT
28 RCL 03
29 1
30 -
31 X<0?
32 CLST
33 FACT
34 ST/ Y
35 CLX
36 RCL 00
37 ST* X
38 FS? 01
39 CHS
40 LASTX
41 Y^X
42 *
43 +
44 DSE 03

45 GTO 01
46 RCL 00
47 RCL 01
48 Y^X
49 /
50 ST- 02
51 1.15443133
52 STO 05
53 RCL 01
54 LBL 02
55 X#0?
56 1/X
57 ST- 05
58 X<> L
59 DSE X
60 GTO 02
61 1
62 STO 03
63 STO 04
64 RCL 05
65 LBL 03
66 RCL 04

67 RCL 00
68 X^2
69 FC? 01
70 CHS
71 *
72 RCL 03
73 ST/ Y
74 RCL 01
75 +
76 /
77 STO 04
78 RCL 05
79 LASTX
80 1/X
81 -
82 RCL 03
83 1/X
84 -
85 STO 05
86 *
87 +
88 ISG 03

89 CLX
90 X#Y?
91 GTO 03
92 RCL 01
93 FACT
94 /
95 RCL 00
96 FS? 01
97 CHS
98 RCL 01
99 Y^X
100 *
101 RCL 02
102 +
103 FC? 01
104 PI
105 FS? 01
106 -2
107 /
108 END

( 151 bytes / SIZE 006 )

STACK	INPUTS	OUTPUTS
Y	n	/
X	x	f(x)

with f(x) = Y_n(x) if CF 01 ; f(x) = K_n(x) if SF 01 ; n = 0 , 1 , 2 , ........ , 69

Example: Evaluate Y₂(3) and K₂(3)

CF 01 2 ENTER^ 3 XEQ "YNX1" >>>> Y₂(3) = -0.160400393 ( in 23 seconds )
SF 01 2 ENTER^ 3 R/S >>>> K₂(3) = 0.061510458 --------------

Note:

-If n is a positive integer, Y_-n(x) = (-1)ⁿ Y_n(x) and K_-n(x) = K_n(x)

5°) Continued Fractions for J_n(x) & Y_n(x)

-Unlike "JNX" & "YNX" ( with CF 01 ) and "YNX1" , the following program produces accurate results for large arguments.

Formulae:

With p + i.q = -1/(2x) + i + (i/x) [ ( 0.5² - n² )/( 2x + 2i + ( 1.5² - n² )/( 2x + 4i + ..... ) ) ]
and g_n = -1/(((2n + 2)/x) - 1/(((2n + 4)/x) - ..... ))

J_n(x) = sign(D) [ ( 2q/(x.Pi) ) / ( q² + ( p - g_n - n/x )² ) ]^1/2 where D = the denominator of the second continued fraction
Y_n(x) = [ ( p - g_n - n/x )/q ] J_n(x)

Data Registers: R00 thru R15 are used by "CFRZ" R01 = x , R16 = n
Flags: /
Subroutines: "CFR" & "CFRZ" ( cf "Continued Fractions for the HP-41" )

01 LBL "JYNX"
02 STO 01
03 X<>Y
04 STO 16
05 "T"
06 ASTO 00
07 CLST
08 RCL 01
09 XEQ "CFRZ"
10 RCL 01
11 ST/ Z
12 /
13 1
14 +
15 X<>Y
16 CHS
17 RCL 01

18 ST+ X
19 1/X
20 -
21 STO 09
22 X<>Y
23 STO 10
24 "U"
25 ASTO 00
26 CLX
27 CLA
28 RCL 01
29 XEQ "CFR"
30 CHS
31 STO 08
32 RCL 09
33 +
34 RCL 16

35 RCL 01
36 /
37 -
38 STO 11
39 RCL 10
40 R-P
41 LAST X
42 ST+ X
43 PI
44 RCL 01
45 *
46 /
47 SQRT
48 X<>Y
49 /
50 RCL 04
51 SIGN

52 *
53 STO 12
54 RCL 11
55 *
56 RCL 10
57 /
58 RCL 12
59 RTN
60 LBL "T"
61 RCL 03
62 ST+ X
63 RCL 01
64 ST+ X
65 RCL 03
66 .5
67 -
68 X^2

69 RCL 16
70 X^2
71 -
72 0
73 X<>Y
74 RTN
75 LBL "U"
76 RCL 02
77 RCL 16
78 +
79 ST+ X
80 RCL 01
81 /
82 1
83 CHS
84 END

( 125 bytes / SIZE 017 )

STACK	INPUTS	OUTPUTS
Y	n >= 0	Y_n(x)
X	x > 0	J_n(x)

Examples:

10 ENTER^ XEQ "JYNX" >>>> J₁₀(10) = 0.207486107 X<>Y Y₁₀(10) = -0.359814151 ( in 2mn27s )

3.14 ENTER^
100 XEQ "JYNX" >>>> J_3.14(100) = 0.079535723 X<>Y Y_3.14(100) = 0.006582327 ( in 4mn14s )

-The methoddoesn't work if n is a negative integer.
-If n < 0 we can use the relations J_n = J_-n cos n.Pi + Y_-n sin n.Pi and Y_n = -J_-n sin n.Pi + Y_-n cos n.Pi
-If x < 0 the results are generally complex.

6°) An Asymptotic Expansion for J_n(x) and Y_n(x)

Formulae:

J_n(x) = (2/(pi*x))^1/2 ( P.cos(x-(2n+1)pi/4) - Q.sin(x-(2n+1)pi/4) ) and Y_n(x) = (2/(pi*x))^1/2 ( P.sin(x-(2n+1)pi/4) + Q.cos(x-(2n+1)pi/4) )

where P ~ 1 - (4n²-1)(4n²-9)/(2!(8x)²) + (4n²-1)(4n²-9)(4n²-25)(4n²-49)/(4!(8x)⁴) - ......

and Q ~ (4n²-1)/(8x) - (4n²-1)(4n²-9)(4n²-25)/(3!(8x)³) + ......

Data Registers: R00 = x ; R01 = n ; R02 to R06: temp
Flags: /
Subroutines: /

01 LBL "AEJY"
02 DEG
03 STO 00
04 X<>Y
05 STO 01
06 ST+ X
07 X^2
08 STO 04
09 1
10 STO 02
11 -
12 RCL 00
13 8
14 *
15 STO 05
16 ST* 05
17 /
18 STO 03

19 XEQ 01
20 STO 06
21 CLX
22 STO 02
23 SIGN
24 STO 03
25 XEQ 01
26 GTO 02
27 LBL 01
28 RCL 03
29 RCL 02
30 2
31 +
32 STO 02
33 ST+ X
34 3
35 -
36 X^2

37 RCL 04
38 -
39 *
40 RCL 04
41 RCL 02
42 ST+ X
43 DSE X
44 X^2
45 -
46 *
47 RCL 02
48 ST/ Y
49 DSE X
50 /
51 RCL 05
52 /
53 STO 03
54 +

55 X#Y?
56 GTO 01
57 RTN
58 LBL 02
59 RCL 00
60 RCL 01
61 ST+ X
62 1
63 +
64 4
65 /
66 PI
67 *
68 -
69 R-D
70 STO 04
71 X<>Y
72 P-R

73 RCL 04
74 RCL 06
75 P-R
76 ST+ T
77 RDN
78 -
79 2
80 PI
81 /
82 RCL 00
83 /
84 SQRT
85 ST* Z
86 *
87 END

( 112 bytes / SIZE 007 )

STACK	INPUTS	OUTPUTS
Y	n	Y_n(x)
X	x > 0	J_n(x)

Example: 4 ENTER^ 100 XEQ "AEJY" >>>> J₄(100) = 0.026105809
X<>Y Y₄(100) = -0.075430120 ( in 10 seconds )

-The program will not work if x is too "small" or if n is too "large"
-If you have the M-Code routine HGF+ in your HP-41, you can use the following variant:

01 LBL "AEJY"
02 DEG
03 STO 00
04 X<>Y
05 ST+ X
06 STO 06
07 1
08 +
09 STO 01
10 LASTX
11 RCL 06
12 -
13 STO 02
14 LASTX
15 3
16 +

17 STO 03
18 LASTX
19 RCL 06
20 -
21 STO 04
22 .5
23 STO 05
24 X^2
25 ST* 01
26 ST* 02
27 ST* 03
28 ST* 04
29 1/X
30 1
31 RCL 00
32 X^2

33 1/X
34 CHS
35 STO 06
36 HGF+
37 STO 07
38 RCL 01
39 STO 08
40 RCL 02
41 STO 09
42 1
43 ST+ 01
44 ST+ 02
45 ST+ 05
46 4
47 X<>Y
48 RCL 06

49 HGF+
50 RCL 08
51 RCL 09
52 *
53 ST+ X
54 *
55 RCL 00
56 ST/ Y
57 PI
58 RCL 08
59 *
60 -
61 R-D
62 STO 06
63 X<>Y
64 P-R

65 RCL 06
66 RCL 07
67 P-R
68 X<> Z
69 -
70 X<> Z
71 +
72 PI
73 RCL 00
74 *
75 2
76 /
77 SQRT
78 ST/ Z
79 /
80 END

( 105 bytes / SIZE 010 )

STACK	INPUTS	OUTPUTS
Y	n	Y_n(x)
X	x > 0	J_n(x)

Examples:

    4     ENTER^
100   XEQ "AEJY" >>>>      J₄(100) = 0.026105809
                                 X<>Y     Y₄(100) = -0.075430120    ( in 6 seconds )

-Likewise, we find

J_PI(11.6) = 0.238578119
Y_PI(11.6) = 0.002890137 ( in 13 seconds )

-But for n = PI and x = 11.5 , the series diverge too soon and we have an infinite loop ( press any key to stop the loop - or ENTER^ ON )

7°) An Asymptotic Expansion for K_n(x)

Formula: K_n(x) ~ (pi/(2x))^1/2 e^-x ( 1 + (4n²-1)/(8x) + (4n²-1)(4n²-9)/(2!(8x)²) + (4n²-1)(4n²-9)(4n²-25)/(3!(8x)³) + ...... )

Data Registers: R00 = x ; R01 = 4n² ; R02 , R03: temp
Flags: /
Subroutines: /

01 LBL "AEK"
02 STO 00
03 X<>Y
04 ST+ X
05 X^2
06 STO 01
07 SIGN
08 STO 02

09 STO 03
10 LBL 01
11 RCL 03
12 RCL 01
13 RCL 02
14 ST+ X
15 DSE X
16 X^2

17 -
18 *
19 RCL 00
20 RCL 02
21 *
22 8
23 *
24 /

25 STO 03
26 +
27 ISG 02
28 CLX
29 X#Y?
30 GTO 01
31 PI
32 RCL 00

33 ST+ X
34 /
35 SQRT
36 *
37 RCL 00
38 E^X
39 /
40 END

( 54 bytes / SIZE 004 )

STACK	INPUTS	OUTPUTS
Y	n	/
X	x	K_n(x)

Examples:

2 ENTER^ 10 XEQ "AEK" >>>> K₂(10) = 0.000,021,509,817,05 ( in 15 seconds )
1.4 ENTER^ 19 R/S >>>> K_1.4(19) = 1.683198846 10^-9 ( in 6 seconds )

-The program will not work if x is too "small" or if n is too "large".
-For instance, n = 2 , x = 7 dont work.
-However, if we stop the iterations before the terms begin to increase, we find K₂(7) = 0.0005545622
-Divergence is more typical than convergence for asymptotic series...

Remark: I_n(x) is accurately computed by "JNX" with flag F01 set, but if you need an asymptotic expansion for this function:

replace line 39 with *
line 36 with /
line 34 with * and add CHS after line 24 ( for instance, I_1.4(19) = 15,597,340 ( in 8 seconds ) )

-The M-Code routine HGF+ may also be used to simplify this program ( cf "Hypergeometric Functions for the HP-41" )

01 LBL "AEK"
02 STO 00
03 X<>Y
04 STO 01
05 CHS
06 .5
07 ST+ 01
08 +
09 STO 02
10 LASTX
11 1/X
12 0
13 RCL 00
14 ST+ X
15 STO 03
16 1/X
17 CHS
18 HGF+
19 PI
20 RCL 03
21 /
22 SQRT
23 *
24 RCL 00
25 CHS
26 E^X
27 *
28 END

( 40 bytes )

STACK	INPUTS	OUTPUTS
Y	n	/
X	x	K_n(x)

Example:

PI 10.1 XEQ "AEK" >>>> K_Pi(10.1) = 0.00002545492107 ---Execution time = 6s---

-If x is not large enough, the divergence occurs too soon and we fall in an infinite loop
-Press any key ( or ENTER^ ON ) to stop the loop.

-For instance, with n = PI and x = 10 , this program doesn't work.

8°) Integrals of J_n(x) & I_n(x)

-Here, we integrate the series expansions used by "JNX" ( n # -1 ; -2 ; -3 ; .... )

Data Registers: /

Flag: F01 CF 01 for §₀^x J_n(t).dt
SF 01 for §₀^x I_n(t).dt

Subroutine: "1/G+" ( Reciprocal of the Gamma Function )

-Synthetic registers M , N & O may be replaced by any unused data registers.
-Lines 46-47-48 may be replaced by XEQ "GAM+" RCL O X<>Y /

01 LBL "ITIJ"
02 STO N
03 X<>Y
04 STO O
05 2
06 ST/ N
07 SIGN
08 STO M
09 +
10 /

11 ENTER^
12 ENTER^
13 LBL 01
14 X<> Z
15 FC? 01
16 CHS
17 RCL N
18 X^2
19 *
20 RCL M

21 ST/ Y
22 RCL O
23 +
24 ST/ Y
25 RCL M
26 +
27 1
28 -
29 ST* Y
30 2

31 +
32 /
33 ISG M
34 CLX
35 ST+ Y
36 X<> Z
37 X#Y?
38 GTO 01
39 RCL N
40 RCL O

41 Y^X
42 *
43 X<> O
44 1
45 +
46 XEQ "1/G+"
47 RCL O
48 *
49 CLA
50 END

( 85 bytes / SIZE 000 )

STACK	INPUTS	OUTPUTS
Y	n	/
X	x	§₀^x f_n(t).dt

where f_n(t) = J_n(t) if flag F01 is clear
and f_n(t) = I_n(t) if flag F01 is set

Examples:

          CF 01
1.4   ENTER^
   3     XEQ"ITIJ" >>>>   §₀³ J_1.4(x).dx = 1.049262785 ( in 14s )

          SF 01
1.4   ENTER^
   3     XEQ"ITIJ" >>>>   §₀³ I_1.4(x).dx = 2.918753200 ( in 13s )

-If flag F01 is clear, this program cannot be used for large arguments ( like "JNX" )

9°) Integral of J_n(x) ( integer order )

Formula: §₀^x J_n(t).dt = 2 ( J_n+1(x) + J_n+3(x) + J_n+5(x) + ........ )

Data Registers: R00 thru R05 ( or R06 ) are used by "JNX1" ( R01 = x ) R07: temp R08: partial sums and eventually, (1/2) §₀^x J_n(t).dt
Flag: F07
Subroutine: "JNX1"

01 LBL "ITJ"
02 STO 01
03 X<>Y
04 STO 07
05 CLX
06 STO 08

07 DSE 07
08 LBL 01
09 SF 07
10 LBL 02
11 RCL 07
12 2

13 +
14 STO 07
15 RCL 01
16 XEQ "JNX1"
17 RCL 08
18 +

19 STO 08
20 LASTX
21 X#Y?
22 GTO 01
23 FS?C 07
24 GTO 02

25 ST+ X
26 END

( 45 bytes / SIZE 009 )

STACK	INPUTS	OUTPUTS
Y	n	/
X	x	§₀^x J_n(t).dt

Examples:

1 ENTER^
3 XEQ "ITJ" >>>> §₀³ J₁(x).dx = 1.260051955 ( in 2mn01s ) - "ITIJ" yields the same result in 14 seconds -

0 ENTER^
10 R/S >>>> §₀¹⁰ J₀(x).dx = 1.067011304 ( in 5mn54s ) - "ITIJ" yields 1.067011425 ( error ~ 10^-7 ) in 23 seconds -

50 ENTER^
30 R/S >>>> §₀³⁰ J₅₀(x).dx = 1.478729947 10^-8 ( in 12mn29s )

-For small arguments, "ITJ" is much slower than "ITIJ"
-There is also a risk of OUT OF RANGE for large x-values,
but if you modify "JNX1" as explained on the right of its listing ( cf paragraph 2°) )
you'll get for instance §₀¹⁰⁰ J₅₀(x).dx = 1.088806747 ( in about 10 seconds with a V41 emulator in turbo mode )

10°) Spherical Bessel Functions

a) Spherical Bessel Functions of the first kind

-Spherical Bessel functions of the first kind are defined by j_n(x) = (Pi/(2x))^1/2 J_n+1/2(x) where n is an integer
-So, we can compute these functions by one of theprograms above.
-However, another approach is used hereafter:

Formulae: j_n-1(x) = j_n(x) (2n+1)/x - j_n+1(x) this recurrence relation is used , starting with j_m(x) = 1 , j_m+1(x) = 0 for some large m
-Then, the results are normalized by the sum 1 j₀²(x) + 3 j₁²(x) + 5 j₂²(x) + ....... = 1
-We could also normalize by calculating j₀(x) = (sin x)/x but this would be inaccurate if sin x ~ 0

Data Registers: R00 , R02 , R03: temp R01 = x , R04 = Sum
Flags: /
Subroutines: /

-The sum in register R04 may easily exceed 9.999999999 10⁹⁹
-So, in order to avoid an OUT OF RANGE,
add    ENTER^   ABS   RCL 05   X>Y?   GTO 01   ST/ 02   ST/ Z   ST/ T   X^2   ST/ 04   LBL 01   R^ R^   after line 36
and add   E40   STO 05 after line 23
-Thus you'll get, for instance,   j₁₀₀(50) = 1.019012264 10^-22 ( in 3mn22s )

01 LBL "SB1"
02 X#0?
03 GTO 00
04 X#Y?
05 RTN
06 SIGN
07 RTN
08 LBL 00
09 STO 01
10 ABS
11 5

12 +
13 X<>Y
14 STO 00
15 X<Y?
16 X<>Y
17 INT
18 ST+ X
19 STO 03
20 ST- 00
21 ST+ X
22 STO 04

23 CLST
24 SIGN
25 ST+ 04
26 LBL 01
27 RCL 03
28 ST+ X
29 1
30 +
31 RCL Y
32 *
33 RCL 01

34 /
35 R^
36 -
37 ISG 00
38 STO 02
39 ENTER^
40 X^2
41 RCL 03
42 ST+ X
43 DSE X
44 *

45 ST+ 04
46 RDN
47 DSE 03
48 GTO 01
49 RCL 02
50 RCL 04
51 SQRT
52 /
53 END

( 74 bytes / SIZE 005 )

STACK	INPUTS	OUTPUTS
Y	n >= 0	/
X	x	j_n(x)

Examples:

2 PI XEQ "SB1" >>>> j₂(Pi) = 0.303963551 ( execution time = 14 seconds )

10 ENTER^
2 R/S >>>> j₁₀(2) = 6.825300865 10^-8 ( in 17 seconds )

100 ENTER^ R/S >>>> j₁₀₀(100) = 0.01088047702 ( in 3mn )

b) Spherical Bessel Functions of the second ( and first ) kind

-Spherical Bessel functions of the second kind are defined by y_n(x) = (Pi/(2x))^1/2 Y_n+1/2(x) where n is an integer
-But the following program uses another method:

Formulae:

y_n+1(x) = y_n(x) (2n+1)/x - y_n-1(x) , j_n+1(x) = j_n(x) (2n+1)/x - j_n-1(x)
y₀(x) = -(cos x)/x , y₁(x) = -(cos x)/x² - (sin x)/x , j₀(x) = (sin x)/x , j₁(x) = (sin x)/x² - (cos x)/x

Data Registers:   R00 = k.nnn ; R01 = x
Flag: F01    CF 01 to calculate y_n(x)
                    SF 01 to calculate   j_n(x)    Set flag F01 only if x is not significantly smaller than n - the reccurence relation is unstable otherwise.
Subroutines: /

-Line 21 may be replaced by RTN and line 33 may be deleted.

01 LBL "SB2"
02 RAD
03 STO 00
04 X<>Y
05 E3
06 /
07 STO 01

08 SIGN
09 P-R
10 RCL 00
11 ST/ Z
12 /
13 FC? 01
14 CHS

15 DEG
16 FS? 01
17 X<>Y
18 LBL 01
19 ISG 01
20 FS? 30
21 GTO 02

22 STO Z
23 RCL 01
24 INT
25 ST+ X
26 DSE X
27 *
28 RCL 00

29 /
30 X<>Y
31 -
32 GTO 01
33 LBL 02
34 END

( 53 bytes / SIZE 002 )

STACK	INPUTS	OUTPUTS
Y	0 <= n < 1000	/
X	x # 0	f_n(x)

where f_n(x) = y_n(x) if CF 01
and f_n(x) = j_n(x) if SF 01 ( only if x is not significantly smaller than n )

Examples:

a)         CF 01
     2     ENTER^
   3.14 XEQ "SB2" >>>>   y₂(3.14) = -0.222053752 ( in 2 seconds )

           CF 01
    30   ENTER^
     5   XEQ "SB2" >>>>   y₅(30) = -7.760717577 10¹⁸ ( in 15 seconds )

b)         SF 01
      4    ENTER^
    100 XEQ "SB2" >>>>   j₄(100) = -4.179461836 10^-3 ( in 3 seconds )

SF 01
100 ENTER^ XEQ "SB2" >>>> j₁₀₀(100) = 1.088047702 10^-2 ( in 47 seconds )

-If n < 0 , we can use the relation y_n(x) = (-1)ⁿ⁺¹ j_-n-1(x) which actually holds for all integers n ( n < 0 , n = 0 , n > 0 )

11°) Zeros of Bessel Functions: Mac Mahon's Asymptotic Expansions

a) Zeros of J_n(x) & Y_n(x)

-The s-th positive zero of J_n(x) is denoted j_n,s and the s-th positive zero of Y_n(x) is denoted y_n,s

-This routine employs the first 5 terms of Mac Mahon's asymptotic expansions ( s >> n ):

j_n,s ~ b - (µ-1) / (8.b) - (4/3) (µ-1) (7µ-31) / (8.b)³ - (32/15) (µ-1) (83µ²-982µ+3779) / (8.b)⁵
- (64/105) (µ-1) (6949µ³-153855µ²+1585743µ-6277237) / (8.b)⁷ - ................

where b = ( s + n/2 - 1/4 ) PI and µ = 4 n²

y_n,s ~ the same formula with b' = ( s + n/2 - 3/4 ) PI instead of b

Data Registers: R00 temp

R01 = b R03 = 8.b R05 = (8.b)² R07 = last term of the series for J
R02 = b' R04 = 8.b' R06 = (8.b')² R08 = last term of the series for Y
Flags: /
Subroutines: /

01 LBL "JYZ"
02 X<>Y
03 STO 02
04 2
05 /
06 +
07 4
08 1/X
09 -
10 STO 01
11 .5
12 -
13 PI
14 ST* 01
15 *
16 X<> 02
17 ST+ X
18 X^2
19 ENTER^
20 ENTER^
21 STO 00
22 6949
23 *

24 153855
25 -
26 *
27 1585743
28 +
29 *
30 6277237
31 -
32 64
33 *
34 105
35 /
36 STO 07
37 STO 08
38 RCL 01
39 RCL 02
40 8
41 ST* Z
42 *
43 STO 04
44 X^2
45 STO 06
46 ST/ T

47 RDN
48 STO 03
49 X^2
50 STO 05
51 /
52 RCL 00
53 83
54 *
55 982
56 -
57 RCL 00
58 *
59 3779
60 +
61 32
62 *
63 15
64 /
65 ST+ Z
66 +
67 RCL 05
68 /
69 X<>Y

70 RCL 06
71 /
72 RCL 00
73 7
74 *
75 31
76 -
77 .75
78 /
79 ST+ Z
80 +
81 RCL 06
82 /
83 X<>Y
84 RCL 05
85 /
86 1
87 ST+ Y
88 ST+ Z
89 RCL 00
90 -
91 ST* 07
92 ST* 08

93 ST* Z
94 *
95 RCL 03
96 /
97 X<>Y
98 RCL 04
99 /
100 RCL 03
101 7
102 Y^X
103 ST/ 07
104 CLX
105 RCL 04
106 LASTX
107 Y^X
108 ST/ 08
109 CLX
110 RCL 02
111 +
112 X<>Y
113 RCL 01
114 +
115 END

( 172 bytes / SIZE 009 )

STACK	INPUTS	OUTPUTS
Y	n >= 0	y_n,s
X	s	j_n,s

Examples:

PI   ENTER^
12    XEQ "JYZ"   >>>>   j_Pi,12 = 41.73324496            R07 = -12 E-9
                              RDN   y_Pi,12 = 40.15792512            R08 = -16 E-9

0    ENTER^
4       R/S            >>>>    j_0,4 = 11.79153443            R07 = -58 E-9
                              RDN   y_0,4 = 10.22234503            R08 = -158 E-9

0    ENTER^
3       R/S            >>>>    j_0,3 = 8.653728            R07 = -508 E-9
                              RDN   y_0,3 = 7.086051            R08 = -2069 E-9

Notes:

-Registers R07 & R08 contain the 5th terms of the expansions: this may be used to give an idea of the accuracy.
-In the first 2 examples, almost all the decimals are correct. In the 3rd example, the results are exact to 6D
-Otherwise, you can delete lines 100 to 109 and lines 92-91-37-36

b) Zeros of the first derivatives J'_n(x) & Y'_n(x)

-The s-th positive zero of J'_n(x) is denoted j'_n,s and the s-th positive zero of Y'_n(x) is denoted y'_n,s

-This routine also uses the first 5 terms of Mac Mahon's asymptotic expansions ( s >> n ):

j'_n,s ~ b - (µ+3) / (8.b) - (4/3) (7µ²+82µ-9) / (8.b)³ - (32/15) (83µ³+2075µ²-3039µ+3537) / (8.b)⁵
- (64/105) (6949µ⁴+296492µ³-1248002µ²+7414380µ-5853627) / (8.b)⁷ - ................

where b = ( s + n/2 - 3/4 ) PI and µ = 4 n²

y'_n,s ~ the same formula with b' = ( s + n/2 - 1/4 ) PI instead of b

Data Registers: R00 temp

R01 = b R03 = 8.b R05 = (8.b)² R07 = last term of the series for J'
R02 = b' R04 = 8.b' R06 = (8.b')² R08 = last term of the series for Y'
Flags: /
Subroutines: /

01 LBL "dJYZ"
02 X<>Y
03 STO 01
04 2
05 /
06 +
07 4
08 1/X
09 -
10 STO 02
11 .5
12 -
13 PI
14 ST* 02
15 *
16 X<> 01
17 ST+ X
18 X^2
19 ENTER^
20 ENTER^
21 STO 00
22 6949
23 *
24 296492
25 +

26 *
27 1248002
28 -
29 *
30 7414380
31 +
32 *
33 5853627
34 -
35 64
36 *
37 105
38 /
39 CHS
40 STO 07
41 STO 08
42 RCL 01
43 RCL 02
44 8
45 ST* Z
46 *
47 STO 04
48 X^2
49 STO 06
50 ST/ T

51 RDN
52 STO 03
53 X^2
54 STO 05
55 /
56 RCL 00
57 83
58 *
59 2075
60 +
61 RCL 00
62 *
63 3039
64 -
65 RCL 00
66 *
67 3537
68 +
69 32
70 *
71 15
72 /
73 ST- Z
74 -
75 RCL 05

76 /
77 X<>Y
78 RCL 06
79 /
80 RCL 00
81 7
82 *
83 82
84 +
85 RCL 00
86 *
87 9
88 -
89 .75
90 /
91 ST- Z
92 -
93 RCL 06
94 /
95 X<>Y
96 RCL 05
97 /
98 RCL 00
99 3
100 +

101 ST- Z
102 -
103 RCL 03
104 /
105 X<>Y
106 RCL 04
107 /
108 RCL 03
109 7
110 Y^X
111 ST/ 07
112 CLX
113 RCL 04
114 LASTX
115 Y^X
116 ST/ 08
117 CLX
118 RCL 02
119 +
120 X<>Y
121 RCL 01
122 +
123 END

( 187 bytes / SIZE 009 )

STACK	INPUTS	OUTPUTS
Y	n >= 0	y'_n,s
X	s	j'_n,s

Examples:

PI   ENTER^
12    XEQ "dJYZ"   >>>>   j'_Pi,12 = 40.14532124            R07 = -57 E-9
                                RDN   y'_Pi,12 = 41.72112783            R08 = -43 E-9

0    ENTER^
4       R/S            >>>>    j'_0,4 = 10.17346816            R07 = 147 E-9
                              RDN   y'_0,4 = 11.74915483           R08 =   54 E-9

0    ENTER^
3       R/S            >>>>    j'_0,3 = 7.015587            R07 = 1930 E-9
                              RDN   y'_0,3 = 8.596006            R08 =   474 E-9

Notes:

12°) Bessel Functions - Complex order & arguments

a) Bessel Functions of the 1st kind

Data Registers: R00 thru R10: temp

Flag: F01 CF 01 to compute J_n(z)
SF 01 to compute I_n(z)

Subroutines: "GAMZ+" ( cf "Gamma Function for the HP-41" § 1°)-h) )
"Z^2" "Z*Z" "Z^Z" "Z/Z" ( cf "Complex Functions for the HP-41" )

"Z^2" "Z*Z" ... and so on ... may replaced by M-Code routines ( cf "A few M-Code Routines for the HP-41" )

01 LBL "JNZ"
02 STO 01
03 RDN
04 STO 02
05 RDN
06 STO 05
07 X<>Y
08 STO 06
09 2
10 ST/ 01
11 ST/ 02
12 RCL 02
13 RCL 01
14 XEQ "Z^2"
15 STO 09
16 X<>Y
17 STO 10

18 CLX
19 STO 04
20 STO 08
21 SIGN
22 STO 00
23 STO 03
24 STO 07
25 LBL 01
26 RCL 10
27 RCL 09
28 RCL 04
29 RCL 03
30 XEQ "Z*Z"
31 RCL 00
32 FC? 01
33 CHS
34 ST/ Z

35 /
36 RCL 05
37 RCL 00
38 +
39 RCL 06
40 X<>Y
41 XEQ "Z/Z"
42 STO 03
43 RCL 07
44 +
45 STO 07
46 LATSX
47 -
48 X<>Y
49 STO 04
50 RCL 08
51 +

52 STO 08
53 LASTX
54 -
55 R-P
56 ISG 00
57 CLX
58 X#0?
59 GTO 01
60 RCL 02
61 RCL 01
62 RCL 06
63 RCL 05
64 XEQ "Z^Z"
65 RCL 08
66 RCL 07
67 XEQ "Z*Z"
68 STO 07

69 X<>Y
70 STO 08
71 RCL 06
72 RCL 05
73 1
74 +
75 XEQ "GAMZ+"
76 RCL 08
77 RCL 07
78 R^
79 R^
80 XEQ "Z/Z"
81 END

( 125 bytes / SIZE 011 )

STACK	INPUTS	OUTPUTS
T	b	/
Z	a	/
Y	y	Im [ f(z) ]
X	x	Re [ f(z) ]

n = a + i b f(z) = J_n(z) if flag F01 is clear
z = x + i y f(z) = I_n(z) if flag F01 is set.

Examples: n = 1 + 4 i , z = 2 + 3 i

• CF 01

   4   ENTER^
   1   ENTER^
   3   ENTER^
   2   XEQ "JNZ"   >>>>   J_1+4i ( 2+3.i ) = 0.342267397 - 0.424349408 i        ---Execution time = 33s---

• SF 01

   4   ENTER^
   1   ENTER^
   3   ENTER^
   2   XEQ "JNZ"   >>>>   I_1+4i ( 2+3.i ) = 1.391608538 + 0.315041667 i        ---Execution time = 33s---

b) Bessel Functions of the 2nd kind

Data Registers: R00 thru R14: temp

Flag: F01 CF 01 to compute Y_n(z)
SF 01 to compute K_n(z)

Subroutines: "JNZ" ( cf §12-a) )
"Z*Z" "Z/Z" "COSZ" "SINZ" ( cf "Complex Functions for the HP-41" )

01 LBL "YNZ"
02 STO 11
03 X<>Y
04 STO 12
05 X<>Y
06 XEQ "JNZ"
07 FS? 01
08 GTO 00
09 STO 13
10 X<>Y
11 STO 14

12 RCL 06
13 RCL 05
14 PI
15 ST* Z
16 *
17 XEQ "COSZ"
18 RCL 14
19 RCL 13
20 XEQ "Z*Z"
21 LBL 00
22 STO 13

23 X<>Y
24 STO 14
25 RCL 06
26 CHS
27 RCL 05
28 CHS
29 RCL 12
30 RCL 11
31 XEQ "JNZ"
32 ST- 13
33 X<>Y

34 ST- 14
35 RCL 06
36 RCL 05
37 PI
38 ST* Z
39 *
40 XEQ "SINZ"
41 RCL 14
42 RCL 13
43 R^
44 R^

45 XEQ "Z/Z"
46 PI
47 2
48 /
49 FC? 01
50 SIGN
51 FC? 01
52 CHS
53 ST* Z
54 *
55 END

( 98 bytes / SIZE 015 )

STACK	INPUTS	OUTPUTS
T	b	/
Z	a	/
Y	y	Im [ f(z) ]
X	x	Re [ f(z) ]

n = a + i b f(z) = J_n(z) if flag F01 is clear
z = x + i y f(z) = I_n(z) if flag F01 is set.

Examples: n = 1 + 4 i , z = 2 + 3 i

• CF 01

   4   ENTER^
   1   ENTER^
   3   ENTER^
   2   XEQ "YNZ"   >>>>   Y_1+4i ( 2+3.i ) = -0.552969731 - 0.212458380 i        ---Execution time = 73s---

• SF 01

   4   ENTER^
   1   ENTER^
   3   ENTER^
   2   XEQ "YNZ"   >>>>   K_1+4i ( 2+3.i ) = 0.011395936 - 0.065079638 i        ---Execution time = 71s---

Note:

-This program does not work if n is a real integer.

c) Bessel Functions of the 2nd kind ( real integer order )

Data Registers: R00 thru R12: temp R01-R02 = the result , R05 = n

Flag: F01 CF 01 to compute Y_n(z)
SF 01 to compute K_n(z)

Subroutines: "JNZ" ( cf §12-a) )
"Z*Z" "Z^2" "Z^X" "LNZ" ( cf "Complex Functions for the HP-41" )

01 LBL "YNZ1"
02 X<>Y
03 STO 12
04 X<>Y
05 STO 11
06 0
07 RDN
08 XEQ "JNZ"
09 STO 01
10 X<>Y
11 STO 02
12 2
13 ST/ 11
14 ST/ 12
15 RCL 12
16 RCL 11
17 XEQ "LNZ"
18 1
19 FS? 01
20 CHS
21 STO 00
22 RCL 05
23 STO 10
24 Y^X
25 ST+ X
26 ST* Z
27 *
28 RCL 02
29 RCL 01
30 XEQ "Z*Z"
31 STO 01
32 X<>Y
33 STO 02
34 RCL 12
35 RCL 11

36 XEQ "Z^2"
37 STO 03
38 X<>Y
39 STO 04
40 CLX
41 STO 07
42 STO 08
43 LBL 01
44 RCL 05
45 RCL 10
46 -
47 FACT
48 RCL 10
49 1
50 -
51 X<0?
52 CLST
53 STO 09
54 FACT
55 /
56 RCL 00
57 RCL 09
58 Y^X
59 *
60 RCL 04
61 RCL 03
62 RCL 09
63 XEQ "Z^X"
64 R^
65 ST* Z
66 *
67 ST+ 07
68 X<>Y
69 ST+ 08
70 DSE 10

71 GTO 01
72 RCL 12
73 RCL 11
74 RCL 05
75 CHS
76 XEQ "Z^X"
77 RCL 08
78 RCL 07
79 XEQ "Z*Z"
80 ST- 01
81 X<>Y
82 ST- 02
83 1.15443133
84 STO 09
85 RCL 05
86 LBL 02
87 X#0?
88 1/X
89 ST- 09
90 X<> L
91 DSE X
92 GTO 02
93 RCL 00
94 ST* 11
95 ST* 12
96 CHS
97 ST* 03
98 ST* 04
99 CLX
100 STO 08
101 STO 06
102 SIGN
103 STO 10
104 RCL 09
105 RCL 05

106 FACT
107 1/X
108 STO 07
109 *
110 STO 00
111 LBL 03
112 RCL 04
113 RCL 03
114 RCL 08
115 RCL 07
116 XEQ "Z*Z"
117 STO 07
118 X<>Y
119 STO 08
120 RCL 10
121 RCL 10
122 RCL 05
123 +
124 *
125 ST/ 07
126 ST/ 08
127 RCL 09
128 LASTX
129 1/X
130 -
131 RCL 10
132 1/X
133 -
134 STO 09
135 ISG 10
136 CLX
137 RCL 07
138 *
139 RCL 00
140 +

141 STO 00
142 LASTX
143 -
144 ABS
145 RCL 08
146 RCL 09
147 *
148 RCL 06
149 +
150 STO 06
151 LASTX
152 -
153 ABS
154 +
155 X#0?
156 GTO 03
157 RCL 12
158 RCL 11
159 RCL 05
160 XEQ "Z^X"
161 RCL 06
162 RCL 00
163 XEQ "Z*Z"
164 ST+ 01
165 X<>Y
166 ST+ 02
167 PI
168 FS? 01
169 -2
170 ST/ 01
171 ST/ 02
172 RCL 02
173 RCL 01
174 END

( 262 bytes / SIZE 013 )

STACK	INPUTS	OUTPUTS
Z	n	/
Y	y	Im [ f(z) ]
X	x	Re [ f(z) ]

where z = x + i y , n is a non-negative integer

and f(z) = Y_n(z) if flag F01 is clear ; f(z) = K_n(z) if flag F01 is set.

Examples: n = 3 , z = 1 + 2 i

• CF 01

   3   ENTER^
   2   ENTER^
   1   XEQ "YNZ1"   >>>>   Y₃ ( 1+2.i ) = 0.290153294 - 0.212118771 i        ---Execution time = 54s---

• SF 01

   3   ENTER^
   2   ENTER^
   1   XEQ "YNZ1"   >>>>   K₃ ( 1+2.i ) = -0.681436426 + 0.625154655 i        ---Execution time = 54s---

-Likewise, Y₀ ( 1+2.i ) = 1.367418719 + 1.521506580 i & K₀ ( 1+2.i ) = -0.242345103 - 0.176267190 i

Note:

-If n is a negative integer, one can use: Y_n(x) = (-1)ⁿ Y_-n(x) and K_n(x) = K_-n(x)

References:

[1] Abramowitz and Stegun , "Handbook of Mathematical Functions" - Dover Publications - ISBN 0-486-61272-4
[2] Press , Vetterling , Teukolsky , Flannery , "Numerical Recipes in C++" - Cambridge University Press - ISBN 0-521-75033-4

hp41programs

Bessel Functions for the HP-41