The Museum of HP Calculators
The Museum of HP Calculators
This program is Copyright © 2006 by Jean-Marc Baillard and is used here by permission.
This program is supplied without representation or warranty of any kind. Jean-Marc Baillard 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.
Overview
1°) Bessel & modified Bessel Functions
of the first kind Jn(x) & In(x)
2°) Bessel Functions Jn(x) of
integer order ( n >= 0 ) ( improved )
3°) Bessel & modified Bessel Functions
of the second kind Yn(x) & Kn(x)
( non-integer order )
4°) Bessel & modified Bessel Functions
of the second kind Yn(x) & Kn(x)
( integer order )
5°) Continued Fractions for Jn(x)
& Yn(x) ( new )
6°) An Asymptotic Expansion for Jn(x)
and Yn(x)
7°) An Asymptotic Expansion for Kn(x)
8°) Integrals of Jn(x)
& In(x) ( new )
9°) Integral of Jn(x)
( integer order ) ( improved )
10°) Spherical Bessel Functions (
new )
a) Spherical Bessel Functions of the first
kind
b) Spherical Bessel Functions of the
second ( and first ) kind
1°) Bessel & modified Bessel Functions of the
first kind Jn(x) & In(x)
Formulae: Jn(x) = (x/2)n [ 1/Gam(n+1) + (-x2/4)1/ (1! Gam(n+2) ) + .... + (-x2/4)k/ (k! Gam(n+k+1) ) + .... ] n # -1 ; -2 ; -3 ; ....
In(x) = (x/2)n [ 1/Gam(n+1) + (x2/4)1/ (1! Gam(n+2) ) + .... + (x2/4)k/ (k! Gam(n+k+1) ) + .... ] n # -1 ; -2 ; -3 ; ....
Data Registers: /
Flag: F01
Subroutine: "GAM" ( Gamma Function
)
-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) = Jn(x) if flag F01 is clear and f(x) = In(x) if flag F01 is set.
Example: Calculate J0.7(1.9) and I0.7(1.9)
CF 01
0.7 ENTER^
1.9 XEQ "JNX" yields
J0.7(1.9) = 0.584978102 (
in 12 seconds )
SF 01
0.7 ENTER^
1.9 XEQ "JNX" yields
I0.7(1.9) = 1.727630598 (
in 12 seconds )
Notes:
-If n is a positive integer, J-n(x) = (-1)n
Jn(x) and I-n(x) =
In(x)
J0(0) = 1 but this program doesn't work if x = n =
0 ( DATA ERROR line 31 )
if you add X=Y? X#0? GTO 00 SIGN
X<>Y LBL 00 after line 30 , it will work in this case too.
Jn(x) is not accurately computed for large arguments,
"JNX1" is better:
2°) Bessel Functions Jn(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: J0(x) + 2 J2(x)
+ 2 J4(x) + 2 J6(x) + ...... = 1 and
Jn-1(x) + Jn+1(x) = (2n/x) Jn(x)
Data Registers: R00 and R02 to R05: temp
R01 = x
Flags: /
Subroutines: /
01 LBL "JNX1"
02 X#0?
Lines 02 to 08 are useful to yield J0(0) = 1 and
Jn(0) = 0 if n # 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
The sums in register R04-R05 may exceed 9.999999999 1099
33 -
So, in order to avoid an OUT OF RANGE,
34 ISG 00
add ENTER^ ABS RCL 06
X>Y? GTO 01 ST/ 02 ST/ Z
ST/ T ST/ 04 ST/ 05 LBL 01
R^ R^ after line 33
35 STO 02
and add E40 STO 06 after line 23
36 ST+ 04
37 RCL 05
Thus you'll get, for instance: J0(1000) =
0.02478668624 ( in 24mn26s )
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 | Jn(x) |
Examples:
2 ENTER^
10 XEQ "JNX1" produces
J2(10) = 0.254630314 ( in 17 seconds
)
3 ENTER^
100 R/S >>>>
J3(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 J1(x) in Z-register , J0(x)
in Y-register and Jn(x) in X-register.
3°) Bessel & modified Bessel Functions of the
second kind Yn(x) & Kn(x) (
non-integer
order )
Formulae: Yn(x)
= ( Jn(x) cos(n(pi)) - J-n(x) ) / sin(n(pi))
; Kn(x) = (pi/2) ( I-n(x)
- In(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 Yn(x)
If F01 is set, "YNX" ---------
Kn(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) = Yn(x) if CF 01 ; f(x) = Kn(x) if SF 01
Example: Evaluate Y1.4(3) and K1.4(3)
CF 01 1.4 ENTER^
3 XEQ "YNX" >>>> Y1.4(3)
= 0.137821836 ( in 31 seconds )
SF 01 1.4 ENTER^
3 R/S
>>>> K1.4(3) = 0.046088036
--------------
4°) Bessel & modified Bessel Functions of the
second kind Yn(x) & Kn(x) (
integer
order )
Formulae: with psi(x) = the digamma function, we have ( positive n ):
Yn(x) = -(1/pi) (x/2)-n SUMk=0,1,....,n-1 (n-k-1)!/(k!) (x2/4)k + (2/pi) ln(x/2) Jn(x) - (1/pi) (x/2)n SUMk=0,1,..... ( psi(k+1) + psi(n+k+1) ) (-x2/4)k / (k!(n+k)!)
Kn(x) = (1/2) (x/2)-n SUMk=0,1,....,n-1
(n-k-1)!/(k!) (-x2/4)k - (-1)n ln(x/2)
In(x) + (1/2) (-1)n (x/2)n SUMk=0,1,.....
( psi(k+1) + psi(n+k+1) ) (x2/4)k / (k!(n+k)!)
Data Registers: R00 = x/2 ; R01 = n
; R02 to R05: temp
Flag: F01 If
F01 is clear, "YNX1" calculates Yn(x)
If F01 is set, "YNX1" ---------
Kn(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) = Yn(x) if CF 01 ; f(x) = Kn(x) if SF 01 ; n = 0 , 1 , 2 , ........ , 69
Example: Evaluate Y2(3) and K2(3)
CF 01 2 ENTER^
3 XEQ "YNX1" >>>> Y2(3)
= -0.160400393 ( in 23 seconds )
SF 01 2 ENTER^
3 R/S
>>>> K2(3) = 0.061510458
--------------
Note: -If n is a positive integer, Y-n(x)
= (-1)n Yn(x) and K-n(x) =
Kn(x)
5°) Continued Fractions for Jn(x) &
Yn(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.52
- n2 )/( 2x + 2i + ( 1.52 - n2
)/( 2x + 4i + ..... ) ) ]
and gn = -1/(((2n
+ 2)/x) - 1/(((2n + 4)/x) - ..... ))
Jn(x) = sign(D) [ ( q2
+ ( p - gn - n/x )2 ) ( 2q/(x.Pi) ) ] -1/2
where D = the denominator of the second continued fraction
Yn(x) = [ ( p - gn - n/x )/q
] Jn(x)
Data Registers: R00 thru R13 are used
by "CFRZ" R01 = x , R14 = n
Flags: /
Subroutines: "CFR" & "CFRZ" ( cf "Continued
Fractions for the HP-41" )
01 LBL "JYNX"
02 STO 01
03 X<>Y
04 STO 14
05 "T"
( The same character as line 60 )
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"
( The same character as line 75 )
25 ASTO 00
26 CLX
27 CLA
28 RCL 01
29 XEQ "CFR"
30 CHS
31 STO 08
32 RCL 09
33 +
34 RCL 14
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 14
70 X^2
71 -
72 0
73 X<>Y
74 RTN
75 LBL "U"
76 RCL 02
77 RCL 14
78 +
79 ST+ X
80 RCL 01
81 /
82 1
83 CHS
84 END
( 121 bytes / SIZE 015 )
| STACK | INPUTS | OUTPUTS |
| Y | n >= 0 | Yn(x) |
| X | x > 0 | Jn(x) |
Examples:
10 ENTER^ XEQ "JYNX" >>>> J10(10) = 0.207486107 X<>Y Y10(10) = -0.359814151 ( in 2mn42s )
3.14 ENTER^
100 XEQ "JYNX" >>>> J3.14(100)
= 0.079535723 X<>Y
Y3.14(100) = 0.006582327
( in 4mn35s )
-The method doesn't work if n is a negative integer.
-If n < 0 we can use the relations Jn = J-n
cos n.Pi + Y-n sin n.Pi and Yn = -J-n
sin n.Pi + Y-n cos n.Pi
-If x < 0 the results are generally complex.
6°) An Asymptotic Expansion for Jn(x)
and Yn(x)
Formulae: Jn(x)
= (2/(pi*x))1/2 ( P.cos(x-(2n+1)pi/4) - Q.sin(x-(2n+1)pi/4)
) and Yn(x) = (2/(pi*x))1/2 ( P.sin(x-(2n+1)pi/4)
+ Q.cos(x-(2n+1)pi/4) )
where P ~ 1 - (4n2-1)(4n2-9)/(2!(8x)2) + (4n2-1)(4n2-9)(4n2-25)(4n2-49)/(4!(8x)4) - ......
and Q ~ (4n2-1)/(8x)
- (4n2-1)(4n2-9)(4n2-25)/(3!(8x)3)
+ ......
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 | Yn(x) |
| X | x | Jn(x) |
Example: 4
ENTER^ 100 XEQ "AEJY" >>>> J4(100)
= 0.026105809
X<>Y
Y4(100) = -0.075430120 ( in 10
seconds )
-The program will not work if x is too "small" or if n is too
"large"
7°) An Asymptotic Expansion for Kn(x)
Formula: Kn(x)
~ (pi/(2x))1/2 e-x ( 1 + (4n2-1)/(8x)
+ (4n2-1)(4n2-9)/(2!(8x)2) + (4n2-1)(4n2-9)(4n2-25)/(3!(8x)3)
+ ...... )
Data Registers: R00 = x ; R01 = 4n2
; 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 | Kn(x) |
Example: 2
ENTER^ 10 XEQ "AEK" >>>> K2(10)
= 0.000,021,509,817,05 ( in 15 seconds )
1.4 ENTER^ 19 R/S
>>>> K1.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 K2(7) = 0.0005545622
-Divergence is more typical than convergence for asymptotic series...
Remark: In(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 ( example:
I1.4(19) = 15,597,340 ( in 8 seconds ) )
8°) Integrals of Jn(x) & In(x)
-Here, we integrate the series expansions used by "JNX"
( n # -1 ; -2 ; -3 ; .... )
Data Registers: /
Flag: F01
CF 01 for §0x
Jn(t).dt
SF 01 for §0x In(t).dt
Subroutine: "1/G+" ( Reciprocal of
the Gamma Function )
-Synthetic registers M , N & O may be replaced by any unused data
registers.
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+" or XEQ
"1/G" but If you prefer XEQ "GAM"
or XEQ "GAM+"
47 RCL O
replace line 48 by X<>Y /
48 *
49 CLA
50 END
( 85 bytes / SIZE 000 )
| STACK | INPUTS | OUTPUTS |
| Y | n | / |
| X | x | §0x fn(t).dt |
where fn(t) = Jn(t) if
flag F01 is clear
and fn(t) = In(t)
if flag F01 is set
Examples:
CF 01
1.4 ENTER^
3 XEQ"ITIJ" >>>>
§03 J1.4(x).dx =
1.049262785 ( in 14s )
SF 01
1.4 ENTER^
3 XEQ"ITIJ" >>>>
§03 I1.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 Jn(x)
( integer order )
Formula: §0x
Jn(t).dt = 2 ( Jn+1(x) + Jn+3(x) + Jn+5(x)
+ ........ )
Data Registers: R00 thru R05 ( or R06
) are used by "JNX1" ( R01 = x ) R07: temp R08: partial sums and
eventually, (1/2) §0x
Jn(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"
XEQ "JNX" would work too, provided x is not too "large"
17 RCL 08
but in this case, "ITIJ" with CF 01 is much faster!
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 | §0x Jn(t).dt |
Examples:
1 ENTER^
3 XEQ "ITJ" >>>>
§03 J1(x).dx = 1.260051955
( in 2mn01s )
- "ITIJ" yields the same result in 14 seconds -
0 ENTER^
10 R/S
>>>> §010 J0(x).dx =
1.067011304 ( in 5mn54s )
- "ITIJ" yields 1.067011425 ( error ~ 10 -7 ) in
23 seconds -
50 ENTER^
30 R/S
>>>> §030 J50(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 §0100
J50(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
jn(x) = (Pi/(2x))1/2 Jn+1/2(x)
where n is an integer
-So, we can compute these functions by one of the programs above.
-However, another approach is used hereafter:
Formulae:
jn-1(x) = jn(x) (2n+1)/x - jn+1(x)
this recurrence relation is used , starting with jm(x)
= 1 , jm+1(x) = 0 for some large m
-Then, the results are normalized by the sum
1 j02(x) + 3 j12(x) + 5 j22(x)
+ ....... = 1
-We could also normalize by calculating j0(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: /
01 LBL "SB1"
02 X#0?
Lines 02 to 08 are only useful to yield j0(0) =
1 and jn(0) = 0 if n # 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^
The sum in register R04 may easily exceed 9.999999999 1099
36 -
So, in order to avoid an OUT OF RANGE,
37 ISG 00
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
38 STO 02
and add E40 STO 05 after line 23
39 ENTER^
40 X^2
Thus you'll get, for instance, j100(50) = 1.019012264
10 -22 ( in 3mn22s )
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 | jn(x) |
Examples:
2 PI XEQ "SB1" >>>> j2(Pi) = 0.303963551 ( execution time = 14 seconds )
10 ENTER^
2 R/S
>>>> j10(2) = 6.825300865 10 -8
( in 17 seconds )
100 ENTER^ R/S >>>> j100(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
yn(x) = (Pi/(2x))1/2 Yn+1/2(x)
where n is an integer
-But the following program uses another method:
Formulae:
yn+1(x) = yn(x) (2n+1)/x - yn-1(x)
; jn+1(x) = jn(x)
(2n+1)/x - jn-1(x)
y0(x) = -(cos x)/x , y1(x) = -(cos x)/x2
- (sin x)/x ;
j0(x) = (sin x)/x , j1(x) = (sin
x)/x2 - (cos x)/x
Data Registers: R00 = k.nnn
; R01 = x
Flag: F01 CF 01 to calculate
yn(x)
SF 01 ----------- jn(x)
Set flag F01 only if x is not significantly smaller than
n - the reccurence relation is unstable otherwise.
Subroutines: /
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
this line may be replaced by RTN and line 33 may be deleted
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 | fn(x) |
where fn(x) = yn(x) if
CF 01
and fn(x) = jn(x)
if SF 01 ( only if x is not significantly
smaller than n )
Examples:
a) CF 01
2 ENTER^
3.14 XEQ "SB2" >>>> y2(3.14)
= -0.222053752 ( in 2 seconds )
CF 01
30 ENTER^
5 XEQ "SB2" >>>>
y5(30) = -7.760717577 1018 ( in 15 seconds
)
b) SF 01
4 ENTER^
100 XEQ "SB2" >>>> j4(100)
= -4.179461836 10 -3 ( in 3 seconds )
SF
01
100 ENTER^ XEQ "SB2" >>>>
j100(100) = 1.088047702 10 -2 ( in 47
seconds )
-If n < 0 , we can use the relation yn(x)
= (-1)n+1 j -n-1(x) which actually holds
for all integers n ( n < 0 , n = 0 , n > 0 )
References:
Abramowitz and Stegun , "Handbook of Mathematical Functions" - Dover
Publications - ISBN 0-486-61272-4
Press , Vetterling , Teukolsky , Flannery , "Numerical Recipes in C++"
- Cambridge University Press - ISBN 0-521-75033-4
Go back to the HP-41 software library
Go back to the general software library
Go
back to the main exhibit hall