Bessel

# Bessel Functions for the HP-41

Overview

1°) Bessel & modified Bessel Functions of the first kind  Jn(x) &  In(x)
2°) Bessel Functions Jn(x) of integer order ( n >= 0 )
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)
6°) An Asymptotic Expansion for  Jn(x) and Yn(x)
7°) An Asymptotic Expansion for  Kn(x)
8°) Integrals of  Jn(x) &  In(x)
9°) Integral of  Jn(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  Jn(x) & Yn(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  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" 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) = 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 )
-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) = Jn(x)  if flag F01 is clear  and  f(x) = In(x)  if flag F01 is set.

Example:

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

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

-However, with or without "HGF+",  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? 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  XY 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 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) [ ( 2q/(x.Pi) ) / ( q2 + ( p - gn - n/x )2 ) ] 1/2    where  D = the denominator of the second continued fraction
Yn(x) =  [ ( p - gn - n/x )/q ] Jn(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 Yn(x) X x > 0 Jn(x)

Examples:

10   ENTER^  XEQ "JYNX"   >>>>   J10(10) =  0.207486107   X<>Y   Y10(10) =  -0.359814151    ( in 2mn27s )

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

-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 > 0 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"
-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 Yn(x) X x > 0 Jn(x)

Examples:

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

-Likewise, we find

JPI(11.6)  =  0.238578119
YPI(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  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)

Examples:

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   ( for instance,   I1.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 Kn(x)

Example:

PI   10.1    XEQ "AEK"   >>>>   KPi(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   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.
-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 §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" 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 §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: /

-The sum in register R04 may easily exceed  9.999999999 1099
-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,   j100(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  XY 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 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 to calculate   jn(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 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 )

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

a) Zeros of  Jn(x) & Yn(x)

-The s-th positive zero of  Jn(x) is denoted  jn,s and the s-th positive zero of  Yn(x) is denoted  yn,s

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

jn,s ~  b - (µ-1) / (8.b) - (4/3) (µ-1) (7µ-31) / (8.b)3 - (32/15) (µ-1) (83µ2-982µ+3779) / (8.b)5
- (64/105) (µ-1) (6949µ3-153855µ2+1585743µ-6277237) / (8.b)7 - ................

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

yn,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)2      R07 = last term of the series for J
R02 = b'        R04 = 8.b'      R06 = (8.b')2      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 yn,s X s jn,s

Examples:

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

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

0    ENTER^
3       R/S            >>>>    j0,3  = 8.653728            R07 =  -508 E-9
RDN   y0,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µ2+82µ-9) / (8.b)3 - (32/15) (83µ3+2075µ2-3039µ+3537) / (8.b)5
- (64/105) (6949µ4+296492µ3-1248002µ2+7414380µ-5853627) / (8.b)7 - ................

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

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)2      R07 = last term of the series for J'
R02 = b'        R04 = 8.b'      R06 = (8.b')2      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:

-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, delete lines 108 to 117 and lines  41-40.

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  Jn(z)
SF 01  to compute   In(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) = Jn(z)  if flag F01 is clear
z = x + i y                f(z) = In(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"   >>>>   J1+4i ( 2+3.i ) = 0.342267397 - 0.424349408 i        ---Execution time = 33s---

•   SF 01

4   ENTER^
1   ENTER^
3   ENTER^
2   XEQ "JNZ"   >>>>   I1+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  Yn(z)
SF 01  to compute   Kn(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) = Jn(z)  if flag F01 is clear
z = x + i y                f(z) = In(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"   >>>>   Y1+4i ( 2+3.i ) = -0.552969731 - 0.212458380 i        ---Execution time = 73s---

•   SF 01

4   ENTER^
1   ENTER^
3   ENTER^
2   XEQ "YNZ"   >>>>   K1+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  Yn(z)
SF 01  to compute   Kn(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) = Yn(z)  if flag F01 is clear     ;      f(z) = Kn(z)  if flag F01 is set.

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

•   CF 01

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

•   SF 01

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

-Likewise,   Y0 ( 1+2.i ) = 1.367418719 + 1.521506580 i    &   K0 ( 1+2.i ) = -0.242345103 - 0.176267190 i

Note:

-If n is a negative integer, one can use:   Yn(x) = (-1)n Y -n(x)  and  Kn(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