Kelvin

# Kelvin Functions for the HP-41

Overview

1°)  Kelvin Functions of the first kind:     bern(x) , bein(x)

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

2°)  Kelvin Functions of the second kind:     kern(x) , kein(x)

a)  Non-integer order
b)  Integer order

3°)  Asymptotic Series

1°)  Kelvin Functions of the first kind

a)  Program#1

Formulae:

bern(x)  =  (x/2)n  Sumk=0,1,2,...... cos ((3n/4+k/2)180°) (x2/4)k / ( k! Gam(n+k+1) )            where  Gam = Euler's Gamma function

bein(x)  =  (x/2)n  Sumk=0,1,2,...... sin ((3n/4+k/2)180°) (x2/4)k / ( k! Gam(n+k+1) )

Data Registers:   R00  R03  R06  R07:  temp

R01 = bern(x)   R04 = x/2
R02 = bein(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 bein(x) X x bern(x)

where n is not a negative integer.

Example:

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

b)  Program#2   ( faster )

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

Formulae:

bern(x) = [ cos(135°n) / Gam(n+1) ] (x/2)n S  -  [ sin(135°n) / Gam(n+2) ] (x/2)n+2 T

bein(x) = [ cos(135°n) / Gam(n+2) ] (x/2)n+2 S  +  [ sin(135°n) / Gam(n+1) ] (x/2)n T

where      S = 0F3(;(n+1)/2,1+n/2,1/2;-x4/256)       T = 0F3(;1+n/2,(n+3)/2,3/2;-x4/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 bein(x) X x bern(x)

where n is not a negative integer.

Example:

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

c)  Modulus

-The following routine returns   Mn(x) = [ bern2(x) + bein2(x) ]1/2

Formula:       Mn(x) = (x/2)n [ Sumk=0,1,2,...... (x2/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 Mn(x)

where n is not a negative integer.

Example:

2   SQRT
PI   XEQ "KLM"  >>>>    Msqrt(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:

kern(x) = - (PI/2) [ bein(x) - (ber-n(x)/sin(n.180°)) + (bern(x)/tan(n.180°)) ]

kein(x)  =   (PI/2) [ bern(x) + (bei-n(x)/sin(n.180°)) - (bein(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 kein(x) X x kern(x)

where n is not an integer.

Example:

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

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

kersqrt(2)(PI) =  0.025901894                                  ---Execution time = 33s---
keisqrt(2)(PI) =  0.089242867

b)  Integer order

Formulae:

kern(x) = (1/2) (x/2)-n SUMk=0,1,....,n-1 cos ((3n/4+k/2)180°) (n-k-1)!/(k!) (x2/4)k - ln(x/2) bern(x) + (PI/4) bein(x)
+ (1/2) (x/2)n SUMk=0,1,.....   cos ((3n/4+k/2)180°)  ( psi(k+1) + psi(n+k+1) ) (x2/4)k / (k!(n+k)!)

kein(x) = - (1/2) (x/2)-n SUMk=0,1,....,n-1 sin ((3n/4+k/2)180°) (n-k-1)!/(k!) (x2/4)k - ln(x/2) bein(x) - (PI/4) bern(x)
+ (1/2) (x/2)n SUMk=0,1,.....   sin ((3n/4+k/2)180°)  ( psi(k+1) + psi(n+k+1) ) (x2/4)k / (k!(n+k)!)

where  psi = digamma function

Data Registers:   R00  R03  R06 thru R10:  temp

R01 = bern(x)   R04 = x/2
R02 = bein(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 kein(x) X x kern(x)

where n is a non-negative integer.

Examples:

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

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

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

3°)  Asymptotic Series

Formulae:     if  x > 0

bern(x) = exp(x/sqrt(2)) (2.PI.x) -1/2 [ fn(x) cos A + gn(x) sin A ] - (1/PI) [ sin (360°.n) kern(x) + cos (360°.n) kein(x) ]
bein(x) = exp(x/sqrt(2)) (2.PI.x) -1/2 [ fn(x) sin A - gn(x) cos A ] + (1/PI) [ cos (360°.n) kern(x) - sin (360°.n) kein(x) ]

kern(x) = exp(-x/sqrt(2)) (PI/(2x))1/2 [ fn(-x) cos B - gn(-x) sin B ]       with    A = (180°/PI) x/sqrt(2) + ( n/2 - 1/8 ) 180°
kein(x) = exp(-x/sqrt(2)) (PI/(2x))1/2 [ -fn(-x) sin B - gn(-x) cos B ]                 B = A + 45°

fn(x)  ~  1 + Sumk=1,2,....  [ ( 4n2 - 1 ) ( 4n2 - 9 ) ...... ( 4n2 - ( 2k-1)2 ) ] / [ k! (-8x)k ]   cos (45°k)
fn(-x)  ~  1 + Sumk=1,2,....  [ ( 4n2 - 1 ) ( 4n2 - 9 ) ...... ( 4n2 - ( 2k-1)2 ) ] / [ k! (8x)k ]   cos (45°k)

gn(x)  ~   Sumk=1,2,....  [ ( 4n2 - 1 ) ( 4n2 - 9 ) ...... ( 4n2 - ( 2k-1)2 ) ] / [ k! (-8x)k ]   sin (45°k)
gn(-x)  ~  Sumk=1,2,....  [ ( 4n2 - 1 ) ( 4n2 - 9 ) ...... ( 4n2 - ( 2k-1)2 ) ] / [ k! (8x)k ]   sin (45°k)

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

R03 = fn(x)      R05 = fn(-x)    R07 = kern(x)         R01-R02 and R09-R10-R11: temp
R04 = gn(x)     R06 = gn(-x)   R08 = kein(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 / kein(x) Z / kern(x) Y n bein(x) X x > 0 bern(x)

Examples:

3.14   ENTER^
10     XEQ "AEKV"  >>>>    ber3.14(10) =  87.53643938                                ---Execution time = 72s---
RDN     bei3.14(10) = -58.97202184
RDN     ker3.14(10) =  0.0004680362504   = R07
RDN     kei3.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       >>>>    ber7.28(25) = -634767.8976                                ---Execution time = 36s---
RDN    bei7.28(25) =  -1673991.514
RDN    ker7.28(25) = 0.000000004145355465
RDN    kei7.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/