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 + x3/3! + ( 1
. 4 x6 )/6! + ( 1 . 4 . 7 x9 )/9! + ..... ]
Bi(x) = [ f(x) + g(x) ] sqrt(3)
and g(x)
= [ 3 -1/3 / Gamma(1/3) ] [ x + 2 x4/4! + ( 2 . 5
x7 )/7! + ( 2 . 5 . 8 x10 )/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) ] ) Sumk=0,1,2,.... (-1)k ck z -k
where z = (2/3) x3/2 ; c0
= 1 , ck = ck-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 x3/2 ) 2F0( 1/6 , 5/6
; -3/(4x3/2) )
Bi(x) ~ { 1/sqrt[ pi sqrt(x) ] }
exp ( 2/3 x3/2 ) 2F0( 1/6 , 5/6
; 3/(4x3/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) x3/2 ; S = Sumk=0,1,2,.... (-1)k c2k z -2k , T = Sumk=0,1,2,.... (-1)k c2k+1 z -2k-1
and
c0 = 1 , ck = ck-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 ) §0infinity Sin ( t3 / 3 + t x ) dt and Hi(x) = ( 1/Pi ) §0infinity Exp ( - t3 / 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) ) Sumk=0,1,2,..... (3k)! / k! / ( 3x3
)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/