Airy Functions & Related Functions for the HP-41

Overview
 

 1°)  Airy Functions Ai & Bi

  1-a)  Ascending Series ( 2 programs )
  1-b)  Asymptotic Expansion for Ai(x) & Bi(x), large positive arguments
  1-c)  Asymptotic Expansion for Ai(x) & Bi(x) , large negative arguments

 2°)  Scorer Functions Gi & Hi

  1-a)  Ascending Series
  1-b)  Asymptotic Expansion for Gi(x)
 
 

1°)  Airy Functions Ai & Bi
 

     a)  Ascending Series
 

    Ai(x) =   f(x) - g(x)                             with            f(x) = [ 3 ^-2/3 / Gamma(2/3) ] [ 1 + x³/3! + ( 1 . 4 x⁶ )/6! + ( 1 . 4 . 7 x⁹ )/9! + .....  ]
    Bi(x) = [ f(x) + g(x) ] sqrt(3)               and            g(x) = [ 3 ^-1/3 / Gamma(1/3) ] [ x + 2 x⁴/4! + ( 2 . 5 x⁷ )/7! + ( 2 . 5 . 8 x¹⁰ )/10! + .....  ]
 

Data Registers:   R00 to R05: temp
Flags: /
Subroutine:  "GAM"  ( cf "Gamma Function for the HP-41" )
 
 
 

01  LBL "AIRY" 02  STO 02 03  STO 05 04  3 05  Y^X 06  STO 00 07  CLX 08  STO 03 09  SIGN 10  STO 01 11  STO 04 12  LBL 01 13  RCL 00 14  RCL 03 15  1

16  + 17  STO 03        18  3 19  * 20  / 21  ST* 04 22  ST* 05 23  LASTX 24  DSE X 25  ST/ 04 26  LASTX 27  1 28  + 29  ST/ 05 30  RCL 04

31  RCL 01        32  + 33  STO 01 34  LASTX 35  RCL 05 36  RCL 02 37  + 38  STO 02 39  RDN 40  X#Y? 41  GTO 01 42  R^ 43  LASTX 44  X#Y? 45  GTO 01

46  2 47  3 48  STO Z 49  / 50  Y^X 51  ST/ 01 52  LASTX 53  XEQ "GAM" 54  ST/ 01 55  3 56  ENTER^ 57  1/X 58  Y^X 59  ST/ 02 60  LASTX

61  XEQ "GAM" 62  ST/ 02 63  RCL 01 64  RCL 02 65  + 66  3 67  SQRT 68  * 69  RCL 01 70  RCL 02 71  - 72  END

 
    ( 102 bytes / SIZE 006 )
 
 

STACK	INPUTS	OUTPUTS
Y	/	Bi(x)
X	x	Ai(x)

 
Example:

  0.4  XEQ "AIRY"  >>>>   Ai(0.4) = 0.254742355           ---Execution time = 9s---
                               X<>Y   Bi(0.4) = 0.801773001
 

   3   R/S   >>>>  Ai(3) = 0.006591141   X<>Y   Bi(3) = 14.03732897   ( in 15 seconds )

Note:

-Several bytes may be saved if you have loaded "0F1" in your HP-41 ( cf "Hypergeometric Functions for the HP-41" §3 )
 
 

01  LBL "AIRY" 02  STO 03 03  3 04  Y^X 05  9 06  / 07  2 08  3 09  /

10  STO 01 11  X<>Y 12  XEQ "0F1"   13  3 14  RCL 01 15  Y^X 16  / 17  STO 02 18  RCL 01

19  XEQ "GAM" 20  ST/ 02 21  2 22  ST* 01 23  RCL 00 24  XEQ "0F1" 25  3 26  ENTER^ 27  1/X

28  Y^X 29  / 30  ST* 03 31  3 32  1/X 33  XEQ "GAM" 34  ST/ 03 35  RCL 02 36  RCL 03

37  + 38  3 39  SQRT 40  * 41  RCL 02        42  RCL 03 43  - 44  END

 
   ( 74 bytes / SIZE 004 )
 
 

STACK	INPUTS	OUTPUTS
Y	/	Bi(x)
X	x	Ai(x)

 
Example:

  0.4  XEQ "AIRY"  >>>>   Ai(0.4) = 0.254742355           ---Execution time = 11s---
                               X<>Y   Bi(0.4) = 0.801773001
 

     b)  Asymptotic Expansion for Ai(x) & Bi(x) , Large Positive Arguments
 

    Ai(x) ~ ( x ^-1/4 / [ 2.sqrt(PI) exp(z) ] )  Sum_k=0,1,2,.... (-1)^k c_k z ^-k

  where  z = (2/3) x^3/2  ;  c₀ = 1 ,  c_k = c_k-1 ( 6k -5 )( 6k - 1 ) / (72 k )
 
 

Data Registers:      R00 = x , R01 = Sum , R02 to R04: temp

Flag:  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 term begins to increase

Subroutines: /
 
 

01  LBL "AI+" 02  SF 03 03  STO 00     04  1.5 05  Y^X 06  LASTX 07  / 08  CHS 09  STO 04 10  CLX 11  STO 02 12  SIGN

13  STO 01  14  STO 03     15  LBL 01 16  RCL 03 17  RCL 02 18  6 19  * 20  1 21  ST+ 02 22  + 23  ST* Y 24  4

25  + 26  * 27  72 28  RCL 02  29  * 30  RCL 04     31  * 32  / 33  ABS 34  LASTX 35  X<> 03 36  ABS

37  X<=Y? 38  GTO 02 39  RCL 03 40  RCL 01  41  + 42  STO 01     43  LASTX 44  X#Y? 45  GTO 01 46  CF 03 47  LBL 02 48  RCL 03

49  RCL 01 50  RCL 04 51  E^X 52  * 53  RCL 00     54  SQRT 55  PI 56  * 57  SQRT 58  ST+ X 59  / 60  END

 
   ( 79 bytes / SIZE 005 )
 
 

STACK	INPUTS	OUTPUTS
Y	/	eps
X	x > 0	Ai(x)

  where  eps is the first neglected term in the series

Examples:

   10  XEQ "AI+"   >>>>   Ai(10) = 1.104753252 10 ^-10     F03 is clear     ---Execution time = 11s---

    5       R/S           >>>>    Ai(5)  = 0.0001083444261       F03 is set        ---Execution time = 19s---
                             X<>Y     eps   = 3.5 10 ^-8

   R03 = eps is an estimation of the relative error, whence the absolute error is about 3.8 10 ^-12  for Ai(5)

Notes:

 "AIRY" returns accurate results for Bi(x)  ( x > 0 , x large ), but if you want to use an asymptotic series for Bi,

  Delete lines 58-36 , replace lines 33-34 by  ENTER^  and delete line 08

-It will yield for example Bi(10) = 455641154.0  ( "AIRY"  gives  455641153.9 )

-These asymptotic series may also be expressed in terms of hypergeometric functions, namely:

     Ai(x) ~ { 1/2sqrt[ pi sqrt(x) ]  }  exp ( -2/3 x^3/2 )  ₂F₀( 1/6 , 5/6 ; -3/(4x^3/2) )
     Bi(x) ~ { 1/sqrt[ pi sqrt(x) ]  }  exp ( 2/3 x^3/2 )  ₂F₀( 1/6 , 5/6 ; 3/(4x^3/2) )

-The following variant uses the M-Code routine HGF+  ( cf "Hypergeometric functions for the HP-41" )
 
 
 

01  LBL "AIBI+" 02  STO 00        03  6 04  1/X 05  STO 01 06  5 07  LASTX 08  / 09  STO 02

10  0 11  RCL 00        12  1.5 13  Y^X 14  LASTX 15  / 16  STO 03 17  1/X 18  2

19  STO T 20  / 21  CHS 22  HGF+ 23  STO 04        24  2 25  0 26  LASTX 27  CHS

28  HGF+ 29  RCL 03 30  E^X 31  * 32  RCL 04        33  RCL 03 34  CHS 35  E^X 36  *

37  PI 38  RCL 00        39  SQRT 40  * 41  SQRT 42  ST/ Z 43  ST+ X 44  / 45  END

 
    ( 63 bytes )
 
 

STACK	INPUTS	OUTPUTS
Y	/	Bi(x)
X	x > 0	Ai(x)

 
Example:

  6.4   XEQ "AIBI+"  >>>>   Ai(6.4) = 0.000003617762307                    ---Execution time = 11s---
                                  RDN   Bi(6.4) = 17400.13559

-With x = 6.3 ,  the series diverge too soon:  press any key to stop the infinite loop.
 

     c)  Asymtotic Expansion for Ai(x) & Bi(x) , Large Negative Arguments
 

 With  x > 0    Ai(-x) ~ ( x ^-1/4 / sqrt(PI) ) [ S sin ( z + 45° ) - T cos ( z + 45° ) ]
                      Bi(-x) ~ ( x ^-1/4 / sqrt(PI) ) [ S cos ( z + 45° ) + T sin ( z + 45° ) ]

    where   z = (2/3) x^3/2    ;     S = Sum_k=0,1,2,.... (-1)^k c_2k z ^-2k    ,    T = Sum_k=0,1,2,.... (-1)^k c_2k+1 z ^-2k-1

       and       c₀ = 1 ,  c_k = c_k-1 ( 6k -5 )( 6k - 1 ) / (72 k )
 

Data Registers:      R00 = -x , R01 = S , R02 = T , R03 to R05: temp

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 term begins to increase

Subroutines: /
 
 

01  LBL "AIBI-"   02  CHS   03  STO 00   04  1.5   05  Y^X   06  LASTX   07  /   08  STO 05   09  CLX   10  STO 02   11  STO 03   12  SIGN   13  STO 01   14  STO 04   15  SF 03   16  SF 10   17  GTO 01   18  LBL 00   19  RCL 04   20  RCL 03   21  6

22  *   23  1   24  ST+ 03   25  +   26  ST* Y   27  4   28  +   29  *   30  72   31  RCL 03          32  *   33  RCL 05   34  *   35  /   36  RTN   37  LBL 01   38  XEQ 00   39  ABS   40  LASTX   41  X<> 04   42  ABS

43  X<=Y?   44  GTO 02   45  RCL 04   46  RCL 02   47  +   48  STO 02   49  LASTX   50  X=Y?   51  CF 10   52  XEQ 00    53  CHS   54  ABS   55  LASTX   56  X<> 04   57  ABS   58  X<=Y?   59  GTO 02   60  RCL 04         61  RCL 01   62  +   63  STO 01

64  LASTX   65  X=Y?   66  FS? 10   67  GTO 01   68  CF 03   69  LBL 02   70  DEG   71  RCL 05    72  R-D   73  45   74  +   75  1   76  P-R   77  RCL 01   78  X<>Y   79  *   80  RCL 02         81  LASTX   82  *   83  CHS   84  RCL 01

85  R^   86  *   87  ST+ Y   88  X<> Z   89  RCL 02    90  LASTX   91  *   92  +   93  RCL 00         94  SQRT   95  PI   96  *   97  SQRT   98  ST/ Z   99  / 100  RCL04 101  X<> Z 102  END

 
     ( 137 bytes / SIZE 006 )
 
 

STACK	INPUTS	OUTPUTS
Z	/	eps
Y	/	Bi(x)
X	x < 0	Ai(x)

    where  eps is the first neglected term in the series

Examples:

   -10  XEQ "AIBI-"  >>>>   Ai(-10) = 0.04024123578                                 ---Execution time = 19s---
                                 X<>Y  Bi(-10) =  -0.314679830      ( F03 is clear )

    -5       R/S       >>>>    Ai(-5) =  0.350761000                                           ---Execution time = 21s---
                            RDN    Bi(-5) = -0.138369139        ( F03 is set )
                           X<>Y     eps   =  3.5 E-8

Note:

-As usual, the M-Code routine HGF+ may simplify the program:
 
 

01  LBL "AIBI-" 02  DEG 03  CHS 04  STO 00       05  1 06  STO 01 07  5 08  STO 02 09  7 10  STO 03 11  11 12  STO 04 13  .5 14  STO 05 15  1/X

16  X^2 17  12 18  ST/ 01 19  ST/ 02 20  ST/ 03 21  ST/ 04 22  SIGN 23  RCL 00       24  3 25  Y^X 26  4 27  * 28  9 29  / 30  STO 06

31  1/X 32  CHS 33  HGF+ 34  STO 07 35  1 36  ST+ 01 37  ST+ 02 38  ST+ 05 39  4 40  X<>Y 41  LASTX  42  HGF+ 43  5 44  * 45  RCL 00

46  RCL 05 47  Y^X 48  48 49  * 50  / 51  RCL 06       52  SQRT 53  R-D 54  45 55  + 56  STO 06  57  X<>Y 58  P-R 59  RCL 06 60  RCL 07

61  P-R 62  ST+ T 63  RDN 64  X<>Y 65  - 66  PI 67  RCL 00       68  SQRT 69  * 70  SQRT 71  ST/ Z 72  / 73  END

 
    ( 99 bytes )
 
 

STACK	INPUTS	OUTPUTS
Y	/	Bi(x)
X	x < 0	Ai(x)

 
Example:

     -7.4   XEQ "AIBI-"  >>>>   Ai(-7.4) =   0.341323752                   ---Execution time = 12s---
                                     RDN    Bi(-7.4) =  -0.021596519

-With x = -7.3 ,  the series diverge too soon:  press any key to stop the infinite loop ( or  ENTER^  ON ).
 

2°)  Scorer Functions Gi & Hi
 

     a)  Ascending Series
 

  Gi(x) = ( 1/Pi ) §₀^infinity  Sin ( t³ / 3  + t x ) dt   and   Hi(x) = ( 1/Pi ) §₀^infinity  Exp ( - t³ / 3  + t x ) dt

  Gi  is also the solution of the differential equation  w''(x) - x w = -1 / Pi     with Gi (0) = 3 ^-7/6 Gam(2/3)  &  Gi'(0) = 3 ^-5/6 Gam(1/3)
  Hi  --------------------------------------------   w''(x) - x w = + 1 / Pi   with Hi(0) = 2 Gi(0) & Hi'(0) = 2 Gi'(0)

"SCO"  calculates Gi(x) & Hi(x) by a series expansion
 

Data Registers:  R00 thru R05:  temp
Flag:   F01    CF 01 <-> Gi(x)
                      SF 01 <-> Hi(x)

Subroutines: /

-Line 32 = X+1  may be replaced by  ISG X  CLX
 
 

01  LBL "SCO"  02  STO 01  03  .1494294525  04  FS? 01  05  ST+ X  06  STO 03  07  *  08  .2049755425  09  FS? 01

10  ST+ X  11  STO 02  12  +  13  PI  14  ST+ X  15  1/X  16  FC? 01  17  CHS  18  STO 04

19  RCL 01  20  2  21  STO 00  22  Y^X  23  STO 05  24  *  25  +  26  LBL 01  27  RCL 04

28  X<> 03  29  X<> 02  30  RCL 00  31  ENTER^  32  X+1  33  STO 00  34  *  35  /  36  STO 04

37  RCL 01  38  RCL 05  39  *  40  STO 05          41  *  42  +  43  X#Y?  44  GTO 01  45  END

 
   ( 83 bytes / SIZE 006 )
 
 

STACK	INPUT	OUTPUT
X	x	Gi(x) or Hi(x)

    Where  X' = Gi(x)  if  flag F01 is clear
                X' = Hi(x)  if  flag F01 is set

Example:    Calculate  Gi(PI) & Hi(PI)

  •  CF 01

     PI  XEQ "SCO"  >>>>   Gi(PI) = 0.108572692                          ---Execution time = 21s---
 
  •  SF 01

     PI  XEQ "SCO"  >>>>   Hi(PI) = 17.63876165                          ---Execution time = 19s---
 

Note:

-We have  for all x ,  Gi(x) +Hi(x) = Bi(x)
 

     b)  Asymptotic Expansion for Gi(x) , Large Positive Arguments
 

-If x is a large positive integer, "SCO" produces inaccurate or even meaningless results
-The following asymptotic series may be used instead:

          Gi(x) ~ ( 1 / (x.Pi) )  Sum_{k=0,1,2,.....}  (3k)! / k! / ( 3x³ )^k
 

Data Registers:  R00 thru R03: temp
Flags: /
Subroutines: /
 
 

01  LBL "AEGI"  02  STO 01  03  3  04  Y^X  05  STO 03  06  CLX  07  STO 00  08  SIGN  09  STO 02  10  LBL 01  11  RCL 00  12  1  13  +  14  STO 00  15  3  16  *  17  DSE X  18  ENTER^  19  DSE X  20  RCL 02  21  *  22  *  23  RCL 03  24  /  25  STO 02  26  +  27  X#Y?  28  GTO 01  29  PI  30  RCL 01  31  *  32  /  33  END

 
   ( 45 bytes / SIZE 004 )
 
 

STACK	INPUT	OUTPUT
X	x > 0	Gi(x)

 
Example:    Calculate  Gi(10.3) & Gi(100)

    10.3  XEQ "AEGI"  >>>>  Gi(10.3)  =  0.03096153019                     ---Execution time = 8s---
     100       R/S           >>>>   Gi(100)  =  0.003183105228                    ---Execution time = 2s---

Note:

-The series already diverges too soon if  x = 10.2
 

References:

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