Kelvin

Kelvin Functions for the HP-41

Overview

1°) Kelvin Functions of the first kind: ber_n(x) , bei_n(x)

    a) Program#1
    b) Program#2 ( faster )
    c) Modulus

2°) Kelvin Functions of the second kind: ker_n(x) , kei_n(x)

a) Non-integer order
b) Integer order

3°) Asymptotic Series

1°) Kelvin Functions of the first kind

a) Program#1

Formulae:

ber_n(x) = (x/2)ⁿ Sum_{k=0,1,2,......} cos ((3n/4+k/2)180°) (x²/4)^k / ( k! Gam(n+k+1) ) where Gam = Euler's Gamma function

bei_n(x) = (x/2)ⁿ Sum_{k=0,1,2,......} sin ((3n/4+k/2)180°) (x²/4)^k / ( k! Gam(n+k+1) )

Data Registers: R00 R03 R06 R07: temp

R01 = ber_n(x) R04 = x/2
R02 = bei_n(x) R05 = n
Flags: /
Subroutine: "GAM" or an M-code routine GAM ( cf "Gamma Function for the HP-41" )

01 LBL "KLV"
02 DEG
03 2
04 /
05 STO 04
06 X<>Y
07 STO 05
08 135
09 *
10 1
11 STO 00
12 STO 03
13 P-R
14 STO 01
15 X<>Y

16 STO 02
17 RCL 05
18 .75
19 *
20 STO 06
21 180
22 STO 07
23 LBL 01
24 RCL 03
25 RCL 04
26 X^2
27 *
28 RCL 05
29 RCL 00
30 ST+ Y

31 *
32 /
33 STO 03
34 RCL 06
35 RCL 00
36 2
37 /
38 +
39 RCL 07
40 *
41 X<>Y
42 ISG 00
43 CLX
44 P-R
45 RCL 01

46 +
47 STO 01
48 LASTX
49 -
50 X<>Y
51 RCL 02
52 +
53 STO 02
54 LASTX
55 -
56 R-P
57 X#0?
58 GTO 01
59 RCL 05
60 1

61 +
62 XEQ "GAM"
63 RCL 04
64 RCL 05
65 Y^X
66 X<>Y
67 /
68 ST* 01
69 ST* 02
70 RCL 02
71 RCL 01
72 END

( 95 bytes / SIZE 008 )

STACK	INPUTS	OUTPUTS
Y	n	bei_n(x)
X	x	ber_n(x)

where n is not a negative integer.

Example:

   2   SQRT
PI   XEQ "KLV" >>>>    ber_sqrt(2)(PI) = -0.674095951                  ---Execution time = 24s---
                             X<>Y    bei_sqrt(2)(PI) = -1.597357210

b) Program#2 ( faster )

-The results will be obtained faster if we use the hypergeometric function ₀F₃(;a,b,c;x)

Formulae:

ber_n(x) = [ cos(135°n) / Gam(n+1) ] (x/2)ⁿ S - [ sin(135°n) / Gam(n+2) ] (x/2)ⁿ⁺² T

bei_n(x) = [ cos(135°n) / Gam(n+2) ] (x/2)ⁿ⁺² S + [ sin(135°n) / Gam(n+1) ] (x/2)ⁿ T

where S = ₀F₃(;(n+1)/2,1+n/2,1/2;-x⁴/256) T = ₀F₃(;1+n/2,(n+3)/2,3/2;-x⁴/256)

Data Registers: R00 thru R03 R06 R07: temp

R04 = x/2 R05 = n
Flags: /
Subroutines: "GAM" ( cf "Gamma Function for the HP-41" ) & "0F3" ( cf "Hypergeometric Functions" )

01 LBL "KLV"
02 DEG
03 2
04 /
05 STO 04
06 X<>Y
07 STO 05
08 LASTX
09 /
10 .5
11 STO 03
12 +
13 STO 01
14 LASTX

15 +
16 STO 02
17 X<> L
18 *
19 4
20 Y^X
21 CHS
22 XEQ "0F3"
23 RCL 04
24 RCL 05
25 Y^X
26 STO 07
27 *
28 STO 06

29 RCL 03
30 ST+ 01
31 ST+ 02
32 ISG 03
33 RCL 00
34 XEQ "0F3"
35 RCL 04
36 X^2
37 *
38 RCL 05
39 1
40 +
41 /
42 ST* 07

43 LASTX
44 XEQ "GAM"
45 ST/ 06
46 ST/ 07
47 RCL 05
48 135
49 *
50 1
51 P-R
52 RCL 06
53 X<>Y
54 *
55 RCL 07
56 LASTX

57 *
58 RCL 06
59 R^
60 *
61 ST+ Y
62 X<> Z
63 RCL 07
64 LASTX
65 *
66 -
67 END

( 99 bytes / SIZE 008 )

STACK	INPUTS	OUTPUTS
Y	n	bei_n(x)
X	x	ber_n(x)

where n is not a negative integer.

Example:

   2   SQRT
PI   XEQ "KLV" >>>>    ber_sqrt(2)(PI) = -0.674095953                ---Execution time = 14s---
                             X<>Y    bei_sqrt(2)(PI) = -1.597357212

c) Modulus

-The following routine returns M_n(x) = [ ber_n²(x) + bei_n²(x) ]^1/2

Formula: M_n(x) = (x/2)ⁿ [ Sum_{k=0,1,2,......} (x²/4)^4k / ( Gam(n+k+1) Gam(n+2k+1) k! ) ] ^1/2

Data Registers: R00 thru R04: temp
Flags: /
Subroutines: "GAM" ( cf "Gamma Function for the HP-41" )

01 LBL "KVM"
02 2
03 /
04 STO 00
05 X<>Y
06 STO 02
07 1
08 STO 01
09 STO 03
10 STO 04

11 LBL 01
12 RCL 00
13 X^2
14 X^2
15 RCL 02
16 RCL 04
17 ST+ Y
18 *
19 RCL 02
20 RCL 04

21 ST+ X
22 +
23 ST* Y
24 1
25 ST+ 04
26 -
27 *
28 /
29 ST* 03
30 RCL 03

31 RCL 01
32 +
33 STO 01
34 LASTX
35 X#Y?
36 GTO 01
37 SQRT
38 RCL 00
39 RCL 02
40 Y^X

41 *
42 STO 01
43 RCL 02
44 1
45 +
46 XEQ "GAM"
47 RCL 01
48 X<>Y
49 /
50 END

( 68 bytes / SIZE 005 )

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

where n is not a negative integer.

Example:

2 SQRT
PI XEQ "KLM" >>>> M_sqrt(2)(PI) = 1.733769138 ---Execution time = 8s---

2°) Kelvin Functions of the second kind

a) Non-integer order

-If n is not an integer, we have:

ker_n(x) = - (PI/2) [ bei_n(x) - (ber_-n(x)/sin(n.180°)) + (ber_n(x)/tan(n.180°)) ]

kei_n(x) = (PI/2) [ ber_n(x) + (bei_-n(x)/sin(n.180°)) - (bei_n(x)/tan(n.180°)) ]

Data Registers: R00 thru R10: temp
Flags: /
Subroutine: "KLV"

01 LBL "KLV2"
02 STO 10
03 XEQ "KLV"
04 STO 08
05 X<>Y
06 STO 09
07 RCL 05
08 CHS
09 RCL 10

10 XEQ "KLV"
11 180
12 RCL 05
13 *
14 STO 07
15 COS
16 STO 06
17 RCL 08
18 *

19 -
20 RCL 07
21 SIN
22 STO 07
23 /
24 RCL 09
25 +
26 X<>Y
27 LASTX

28 RCL 06
29 *
30 -
31 RCL 07
32 /
33 RCL 08
34 -
35 X<>Y
36 PI

37 2
38 /
39 CHS
40 ST* Z
41 *
42 END

( 62 bytes / SIZE 011 )

STACK	INPUTS	OUTPUTS
Y	n	kei_n(x)
X	x	ker_n(x)

where n is not an integer.

Example:

   2   SQRT
PI   XEQ "KLV2" >>>>    ker_sqrt(2)(PI) = 0.025901896          ---Execution time = 56s---        with the first version of "KLV"
                               X<>Y    kei_sqrt(2)(PI) = 0.089242870

-If you prefer the second version of "KLV", it yields:

ker_sqrt(2)(PI) = 0.025901894 ---Execution time = 33s---
kei_sqrt(2)(PI) = 0.089242867

b) Integer order

Formulae:

ker_n(x) = (1/2) (x/2)^-n SUM_{k=0,1,....,n-1} cos ((3n/4+k/2)180°) (n-k-1)!/(k!) (x²/4)^k - ln(x/2) ber_n(x) + (PI/4) bei_n(x)
+ (1/2) (x/2)ⁿ SUM_k=0,1,..... cos ((3n/4+k/2)180°) ( psi(k+1) + psi(n+k+1) ) (x²/4)^k / (k!(n+k)!)

kei_n(x) = - (1/2) (x/2)^-n SUM_{k=0,1,....,n-1} sin ((3n/4+k/2)180°) (n-k-1)!/(k!) (x²/4)^k - ln(x/2) bei_n(x) - (PI/4) ber_n(x)
+ (1/2) (x/2)ⁿ SUM_k=0,1,..... sin ((3n/4+k/2)180°) ( psi(k+1) + psi(n+k+1) ) (x²/4)^k / (k!(n+k)!)

where psi = digamma function

Data Registers: R00 R03 R06 thru R10: temp

                               R01 = ber_n(x)   R04 = x/2
                               R02 = bei_n(x)    R05 = n
Flags: /
Subroutine:   "KLV"

-Lines 03 to 05 are superfluous if "KLV" is the program listed in §1a)
-Line 72 = 2 x Euler's Constant

01 LBL "KLV2I"
02 XEQ "KLV"
03 STO 01
04 X<>Y
05 STO 02
06 CLX
07 STO 06
08 STO 07
09 RCL 05
10 STO 03
11 X=0?
12 GTO 00
13 LBL 01
14 RCL 05
15 3
16 *
17 RCL 03
18 1
19 -
20 STO 00
21 ST+ X
22 +
23 45
24 *
25 RCL 05
26 RCL 03
27 -
28 FACT
29 RCL 00
30 FACT
31 /

32 RCL 04
33 X^2
34 RCL 00
35 Y^X
36 *
37 P-R
38 ST+ 06
39 X<>Y
40 ST+ 07
41 DSE 03
42 GTO 01
43 LBL 00
44 RCL 04
45 RCL 05
46 Y^X
47 ST+ X
48 ST/ 06
49 CHS
50 ST/ 07
51 RCL 04
52 LN
53 RCL 01
54 *
55 RCL 02
56 PI
57 4
58 /
59 STO 00
60 *
61 -
62 ST- 06

63 RCL 04
64 LN
65 RCL 02
66 *
67 RCL 01
68 RCL 00
69 *
70 +
71 ST- 07
72 1.15443133
73 STO 03
74 RCL 05
75 STO 08
76 LBL 02
77 RCL 08
78 X#0?
79 1/X
80 ST- 03
81 DSE 08
82 GTO 02
83 RCL 05
84 135
85 *
86 RCL 03
87 P-R
88 STO 08
89 X<>Y
90 STO 09
91 1
92 STO 00
93 STO 10

94 LBL 03
95 RCL 10
96 RCL 04
97 X^2
98 *
99 RCL 00
100 ST/ Y
101 RCL 05
102 +
103 /
104 STO 10
105 RCL 03
106 LASTX
107 1/X
108 -
109 RCL 00
110 1/X
111 -
112 STO 03
113 *
114 RCL 05
115 3
116 *
117 RCL 00
118 ST+ X
119 +
120 45
121 *
122 ISG 00
123 CLX
124 X<>Y

125 P-R
126 RCL 08
127 +
128 STO 08
129 LASTX
130 -
131 X<>Y
132 RCL 09
133 +
134 STO 09
135 LASTX
136 -
137 R-P
138 X#0?
139 GTO 03
140 RCL 04
141 RCL 05
142 Y^X
143 LASTX
144 FACT
145 ST+ X
146 /
147 ST* 08
148 ST* 09
149 RCL 07
150 RCL 09
151 -
152 RCL 06
153 RCL 08
154 -
155 END

( 203 bytes / SIZE 011 )

STACK	INPUTS	OUTPUTS
Y	n	kei_n(x)
X	x	ker_n(x)

where n is a non-negative integer.

Examples:

   0    ENTER^
PI   XEQ "KLV2I" >>>>    ker₀(PI) = -0.063594733          ---Execution time = 43s---        with the first version of "KLV"
                                X<>Y    kei₀(PI) = -0.039038451

   3    ENTER^
PI       R/S        >>>>    ker₃(PI) = -0.045754846          ---Execution time = 60s---        with the first version of "KLV"
                         X<>Y    kei₃(PI) = -0.196928535

-With the second "KLV", execution times are 38 seconds and 49 seconds respectively.

3°) Asymptotic Series

Formulae: if x > 0

ber_n(x) = exp(x/sqrt(2)) (2.PI.x)^-1/2 [ f_n(x) cos A + g_n(x) sin A ] - (1/PI) [ sin (360°.n) ker_n(x) + cos (360°.n) kei_n(x) ]
bei_n(x) = exp(x/sqrt(2)) (2.PI.x)^-1/2 [ f_n(x) sin A - g_n(x) cos A ] + (1/PI) [ cos (360°.n) ker_n(x) - sin (360°.n) kei_n(x) ]

ker_n(x) = exp(-x/sqrt(2)) (PI/(2x))^1/2 [ f_n(-x) cos B - g_n(-x) sin B ] with A = (180°/PI) x/sqrt(2) + ( n/2 - 1/8 ) 180°
kei_n(x) = exp(-x/sqrt(2)) (PI/(2x))^1/2 [ -f_n(-x) sin B - g_n(-x) cos B ] B = A + 45°

f_n(x) ~ 1 + Sum_k=1,2,.... [ ( 4n² - 1 ) ( 4n² - 9 ) ...... ( 4n² - ( 2k-1)² ) ] / [ k! (-8x)^k ] cos (45°k)
f_n(-x) ~ 1 + Sum_k=1,2,.... [ ( 4n² - 1 ) ( 4n² - 9 ) ...... ( 4n² - ( 2k-1)² ) ] / [ k! (8x)^k ] cos (45°k)

g_n(x) ~ Sum_k=1,2,.... [ ( 4n² - 1 ) ( 4n² - 9 ) ...... ( 4n² - ( 2k-1)² ) ] / [ k! (-8x)^k ] sin (45°k)
g_n(-x) ~ Sum_k=1,2,.... [ ( 4n² - 1 ) ( 4n² - 9 ) ...... ( 4n² - ( 2k-1)² ) ] / [ k! (8x)^k ] sin (45°k)

Data Registers: R00 = eps = the first neglected term ( before multiplying by sin or cos )

R03 = f_n(x) R05 = f_n(-x) R07 = ker_n(x) R01-R02 and R09-R10-R11: temp
R04 = g_n(x) R06 = g_n(-x) R08 = kei_n(x)

Flags: F10 & F03 If F03 is clear at the end, the accuracy is maximum.
If F03 is set at the end, the series have been truncated just before the terms begin to increase.
Subroutines: /

01 LBL "AEKV"
02 DEG
03 STO 01
04 X<>Y
05 STO 02
06 CLX
07 STO 04
08 STO 06
09 STO 07
10 SIGN
11 STO 00
12 STO 03
13 STO 05
14 STO 08
15 SF 03
16 SF 10
17 LBL 01
18 ISG 07
19 CLX
20 RCL 00
21 RCL 02
22 ST+ X
23 X^2
24 RCL 07
25 ST+ X
26 DSE X
27 X^2

28 X>Y?
29 CF 10
30 -
31 RCL 07
32 8
33 *
34 RCL 01
35 *
36 /
37 ST* 08
38 *
39 CHS
40 ABS
41 LASTX
42 X<> 00
43 ABS
44 FC? 10
45 X>Y?
46 X<0?
47 GTO 02
48 RCL 07
49 45
50 *
51 STO 09
52 RCL 00
53 P-R
54 RCL 03

55 +
56 STO 03
57 CLX
58 RCL 04
59 +
60 STO 04
61 RCL 09
62 RCL 08
63 P-R
64 RCL 05
65 +
66 STO 05
67 LASTX
68 -
69 X<>Y
70 RCL 06
71 +
72 STO 06
73 LASTX
74 -
75 R-P
76 X#0?
77 GTO 01
78 CF 03
79 LBL 02
80 RCL 01
81 2

82 SQRT
83 /
84 STO 10
85 R-D
86 RCL 02
87 4
88 ST* 02
89 1/X
90 -
91 90
92 ST* 02
93 *
94 +
95 STO 09
96 45
97 +
98 STO 08
99 RCL 06
100 P-R
101 RCL 08
102 RCL 05
103 P-R
104 R^
105 -
106 STO 07
107 RDN
108 +

109 CHS
110 RCL 10
111 E^X
112 STO 11
113 RCL 01
114 ST+ X
115 PI
116 /
117 SQRT
118 STO 01
119 *
120 ST/ 07
121 /
122 STO 08
123 RCL 02
124 RCL 07
125 P-R
126 RCL 02
127 RCL 08
128 P-R
129 R^
130 +
131 STO 02
132 RDN
133 -
134 STO 10
135 RCL 09

136 RCL 04
137 P-R
138 RCL 09
139 RCL 03
140 P-R
141 R^
142 +
143 X<> Z
144 -
145 RCL 11
146 RCL 01
147 /
148 ST* Z
149 *
150 RCL 10
151 +
152 X<>Y
153 RCL 02
154 -
155 PI
156 ST/ Z
157 /
158 RCL 08
159 RCL 07
160 R^
161 R^
162 END

( 194 bytes / SIZE 012 )

STACK	INPUTS	OUTPUTS
T	/	kei_n(x)
Z	/	ker_n(x)
Y	n	bei_n(x)
X	x > 0	ber_n(x)

Examples:

3.14   ENTER^
10     XEQ "AEKV" >>>>    ber_3.14(10) = 87.53643938                                ---Execution time = 72s---
                                  RDN     bei_3.14(10) = -58.97202184
                                  RDN     ker_3.14(10) = 0.0004680362504   = R07
                                  RDN     kei_3.14(10) = -0.00006810172135 = R08

-Flag F03 is set and since R00 ~ -5 E-10 is very small, almost all the decimals are correct.

7.28 ENTER^
   25       R/S       >>>>    ber_7.28(25) = -634767.8976                                ---Execution time = 36s---
                           RDN    bei_7.28(25) = -1673991.514
                           RDN    ker_7.28(25) = 0.000000004145355465
                           RDN    kei_7.28(25) = 0.00000001035356581

-Flag F03 is clear. The accuracy is maximum ( R00 ~ -3 E-10 ): the iterations stop when consecutive sums are equal.

References:

[1] Abramowitz and Stegun - "Handbook of Mathematical Functions" - Dover Publications - ISBN 0-486-61272-4
[2] http://functions.wolfram.com/

hp41programs

Kelvin Functions for the HP-41