Parabolic Cylinder Functions for the HP-41
Overview
1°) U(a,x) & V(a,x)
a) Using Kummer's Functions
b) Asymptotic Expansions
c) Recurrence Relation
2°) W(a,x)
a) Series Expansion
b) Asymptotic Expansion
3°) Dn(x)
a) Series Expansion
b) Asymptotic Expansion
1°) U(a;x) & V(a;x)
a) Using
Kummer's Functions
-The standard solutions of the differential equation d2y/dx2 - ( x2/4 + a ).y = 0 may be expressed with the Kummer's Functions
U(a;x) = pi1/22-1/4-a/2e-x^2/4 (Gam(3/4+a/2))-1 M(1/4+a/2;1/2;x2/2) - pi1/221/4-a/2 x.e-x^2/4 (Gam(1/4+a/2))-1 M(3/4+a/2;3/2;x2/2)
V(a;x) = Gam(a+1/2) [ sin(a.pi) U(a;x) + U(a;-x) ].(pi)-1
-However, if a = -1/2 ; -3/2 ; .... the Gamma function will be infinite. Therefore, other formulas are used if a < 0 , namely:
U(a;x) = Y1 cos(pi(1/4+a/2)) - Y2 sin(pi(1/4+a/2)) and V(a;x) = (Gam(1/2-a))-1 [ Y1 sin(pi(1/4+a/2)) + Y2 cos(pi(1/4+a/2)) ]
where Y1 = pi-1/22-1/4-a/2e-x^2/4
(Gam(1/4-a/2)) M(1/4+a/2;1/2;x2/2) and Y2
= pi-1/221/4-a/2 x.e-x^2/4 (Gam(3/4-a/2))
M(3/4+a/2;3/2;x2/2)
Data Registers: R00 thru R06 ( R03
= x ; R04 = a )
Flags: /
Subroutines: "GAM" or "GAM+" ...
( Gamma Function ) , "KUM" ( Kummer's Functions )
| 01 LBL "PCF1"
02 STO 03 03 X^2 04 X<>Y 05 STO 04 06 .5 07 STO 02 08 ST* Y 09 ST* Z 10 X^2 11 LASTX 12 + 13 + 14 STO 01 15 X<>Y 16 ISG 02 17 XEQ "KUM" 18 STO 05 19 .5 20 ST- 01 21 STO 02 22 RCL 00 23 XEQ "KUM" 24 RCL 00 |
25 RCL 02
26 * 27 CHS 28 E^X 29 ST* 05 30 * 31 2 32 RCL 01 33 Y^X 34 / 35 STO 06 36 RCL 03 37 LASTX 38 RCL 02 39 SQRT 40 * 41 / 42 ST* 05 43 RCL 02 44 RCL 01 45 RCL 04 46 SIGN 47 * 48 + |
49 STO 00
50 XEQ "GAM" 51 PI 52 SQRT 53 / 54 RCL 04 55 SIGN 56 Y^X 57 ST/ 06 58 RCL 00 59 RCL 02 60 LASTX 61 * 62 - 63 XEQ "GAM" 64 PI 65 SQRT 66 / 67 RCL 04 68 SIGN 69 Y^X 70 ST/ 05 71 RCL 02 72 RCL 04 |
73 ABS
74 + 75 XEQ "GAM" 76 STO 02 77 1 78 CHS 79 ACOS 80 STO 00 81 ST* 01 82 RCL 04 83 X<0? 84 GTO 01 85 RCL 06 86 RCL 05 87 - 88 STO Y 89 RCL 00 90 RCL 04 91 * 92 SIN 93 * 94 RCL 05 95 RCL 06 96 + |
97 +
98 RCL 02 99 * 100 PI 101 / 102 GTO 02 103 LBL 01 104 RCL 01 105 RCL 05 106 P-R 107 RCL 01 108 RCL 06 109 P-R 110 R^ 111 - 112 RDN 113 + 114 RCL 02 115 / 116 LBL 02 117 R^ 118 END |
( 161 bytes / SIZE 007 )
| STACK | INPUTS | OUTPUTS |
| Y | a | V(a;x) |
| X | x | U(a;x) |
Example1: Calculate U(0.4;1.9) & V(0.4;1.9)
0.4 ENTER^
1.9 XEQ "PCF1" >>>>
U(0.4;1.9) = 0.194020564 X<>Y
V(0.4;1.9) = 1.882850363 (
in 42 seconds )
Example2: Calculate U(-0.4;1.9) & V(-0.4;1.9)
-0.4 ENTER^
1.9 XEQ "PCF1" >>>>
U(-0.4;1.9) = 0.376027811
X<>Y V(-0.4;1.9) = 1.376169516
( in 41 seconds )
b) Asymptotic
Expansions for U(a;x) & V(a;x) ( x large , a moderate
)
Formulas: U(a;x) ~ e-x^2/4 x-a-1/2 [ 1 - (a+1/2)(a+3/2)/(1.(2x2)) + (a+1/2)(a+3/2)(a+5/2)(a+7/2)/(2.(2x2)2) - ....... ]
V(a;x) ~ (2/pi)1/2 ex^2/4
xa-1/2 [ 1 + (a-1/2)(a-3/2)/(1.(2x2)) + (a-1/2)(a-3/2)(a-5/2)(a-7/2)/(2.(2x2)2)
+ ....... ]
Data Registers: R00 = x ; R01 = a ; R02 ,
R03: temp
Flags: F01 if F01 is clear,
AEUV computes U(a;x)
if F01 is set, AEUV computes V(a;x)
Subroutines: /
| 01 LBL "AEUV"
02 STO 00 03 X<>Y 04 STO 01 05 1 06 STO 02 07 STO 03 08 LBL 01 09 RCL 03 10 ST+ X 11 ENTER^ 12 DSE X 13 .5 |
14 ST- Z
15 - 16 RCL 01 17 FS? 01 18 CHS 19 ST+ Z 20 + 21 * 22 RCL 02 23 FC? 01 24 CHS 25 * 26 RCL 03 |
27 /
28 RCL 00 29 X^2 30 ST+ X 31 / 32 ISG 03 33 CLX 34 STO 02 35 + 36 X#Y? 37 GTO 01 38 RCL 00 39 X^2 |
40 FC? 01
41 CHS 42 4 43 / 44 E^X 45 * 46 RCL 00 47 RCL 01 48 FC? 01 49 CHS 50 .5 51 - 52 Y^X |
53 *
54 2 55 PI 56 / 57 SQRT 58 FC? 01 59 SIGN 60 * 61 END |
( 84 bytes / SIZE 004 )
| STACK | INPUTS | OUTPUTS |
| Y | a | / |
| X | x | f(x) |
where f(x) = U(a;x) if flag F01 is clear
f(x) =
V(a;x) if flag F01 is set.
Example:
CF 01 2 ENTER^ 10
XEQ "AEUV" >>>> U(2;10) = 4.210624069 10-14
( in 14 seconds )
SF 01 2 ENTER^ 10
R/S >>>>
V(2;10) = 1.823604920 1012
( in 7 seconds )
Notes:
-V(2;10) is accurately computed by "PCF1" ( it yields 1.823604925
1012 in 2mn26s )
but "PCF1" gives U(2;10) = -2000 ( a completely wrong
result! )
This meaningless value results from the subtraction of 2 large
numbers, and the leading digits are cancelled.
-If x is too small, the series will diverge too soon.
-However, U(-5;5) = 1.879976816 is correctly calculated
( the negative a-value helps convergence).
c) A Recurrence
Relation for U(a;x)
-If a and x are large, the recurrence relation (a+1/2)
U(a+1;x) + x.U(a;x) = U(a-1;x) may be used straightforward
from 2 values.
-However, this can lead to low accuracy by cancellation of leading
digits.
-Better accuracy is attainable if we use this relation in the reverse
direction, starting with two arbitrary values ( 1 and 0 ) for a+21 and
a+20 ( line 05 )
-Then, "AEUV" is called as a subroutine and finally, a normalization
factor gives U(a;x)
Data Registers: R00 = x ; R04 = a ;
R02 , R03 , R07: temp
Flag: CF01
Subroutine: "AEUV"
-Line 44: the recurrence relation is applied until 2x+a'
<= 0 but other choices may be better ( for instance | a | < 1 )
especially if line 54 is replaced by XEQ "PCF1"
| 01 LBL "UAX"
02 STO 00 03 X<>Y 04 STO 04 05 20 06 + 07 STO 01 08 CLX 09 STO 02 10 SIGN 11 STO 03 12 LBL 01 |
13 RCL 02
14 RCL 01 15 .5 16 + 17 * 18 RCL 03 19 STO 02 20 RCL 00 21 * 22 + 23 STO 03 24 RCL 01 |
25 1
26 - 27 STO 01 28 RCL 04 29 X#Y? 30 GTO 02 31 RCL 03 32 STO 07 33 LBL 02 34 RCL 03 35 ABS 36 E10 |
37 X>Y?
38 GTO 01 39 RCL 00 40 ST+ X 41 RCL 01 42 + 43 X>0? 44 GTO 01 45 RCL 04 46 RCL 01 47 X>Y? 48 GTO 01 |
49 RCL 03
50 ST/ 07 51 RCL 01 52 RCL 00 53 CF 01 54 XEQ "AEUV" 55 RCL 07 56 * 57 END |
( 81 bytes / SIZE 008 )
| STACK | INPUTS | OUTPUTS |
| Y | a | / |
| X | x | U(a;x) |
Examples:
5 ENTER^ XEQ "UAX"
>>>> U(5;5) = 1.552271290 10-7
( execution time = 45 seconds )
12 ENTER^ 7 R/S
>>>> U(12;7) = 3.282492495 10-17
( execution time = 55 seconds )
-In these examples, "PCF1" would yield meaningless results for U(a;x)
because of cancellation of the leading digits.
-If x is too "small" , the asymptotic expansion used in "AEUV"
may diverge too soon.
-In this case, replace line 54 with XEQ "PCF1" ( that's
why I've used register R07 instead of R05 )
2°) W(a;x)
a) Series
Expansion
-The differential equation d2y/dx2
+ ( x2/4 - a ).y = 0 is solved by a series expansion
y = a0 + a1.x + a2.x2 + ......
+ ak.xk + ......
-The standard solution W(a;x) is defined by
a0 = 2-3/4 | Gam(1/4 + i.a/2) / Gam(3/4 + i.a/2) |1/2 and a1 = -2-1/4 | Gam(3/4 + i.a/2) / Gam(1/4 + i.a/2) |1/2
-The relation | Gam(1/4 + i.a/2) |.| Gam(3/4 + i.a/2) | = pi.(2/cosh(a.pi))1/2 is also used to call "GAMZ" only once.
Data Registers: R00 thru R09 ( temp )
When the program stops, R01 = W(a;x)
Flags: /
Subroutines: "GAMZ" ( Gamma function
in the complex plane )
| 01 LBL "PCF2"
02 STO 06 03 X<>Y 04 STO 07 05 .5 06 ST* Y 07 X^2 08 XEQ "GAMZ" 09 LASTX 10 PI 11 RCL 07 12 * 13 E^X 14 ENTER^ 15 1/X 16 + |
17 .5
18 STO 04 19 * 20 SQRT 21 PI 22 / 23 SQRT 24 * 25 ST* 04 26 1/X 27 CHS 28 STO 05 29 RCL 06 30 STO 08 31 * 32 RCL 04 |
33 +
34 STO 01 35 CLX 36 STO 02 37 STO 03 38 SIGN 39 STO 00 40 LBL 01 41 3 42 STO 09 43 LBL 02 44 RCL 04 45 RCL 07 46 * 47 RCL 02 48 4 |
49 /
50 - 51 RCL 00 52 RCL 00 53 1 54 + 55 STO 00 56 * 57 / 58 X<> 05 59 X<> 04 60 X<> 03 61 STO 02 62 RCL 06 63 RCL 08 64 * |
65 STO 08
66 RCL 05 67 * 68 RCL 01 69 + 70 STO 01 71 LASTX 72 X#Y? 73 GTO 01 74 DSE 09 75 GTO 02 76 END |
( 100 bytes / SIZE 010 )
| STACK | INPUTS | OUTPUTS |
| Y | a | W(a;x) |
| X | x | W(a;x) |
Example: Calculate W(0.4;1.9)
0.4 ENTER^
1.9 XEQ "PCF2" >>>>
0.219336459 ( in 52 seconds )
b) An Asymptotic
Expansion for W(a;x)
( x large , a moderate )
Formulae: W(a;x) ~ (2k/x)1/2 ( s1(a;x) cos A - s2(a;x) sin A ) and W(a;-x) ~ (2/(kx))1/2 ( s1(a;x) sin A + s2(a;x) cos A ) ( x > 0 )
where
A = x2/2 - a.Ln(x) + pi/4 + B/2 with B = arg
( Gam(1/2+i.a) )
k = 1/(epi.a+(1+e2pi.a)1/2)
s1(a;x)+i.s2(a;x)
= Sumn=0,1,2,...... (-i)n Gam(2n+1/2+i.a)/Gam(1/2+i.a)
1/(n!(2x2)n)
Data Registers: R00 thru R08: temp
when the program stops, R02 = a ; R08 = x
Flags: /
Subroutine: "GAMZ" ( gamma function
in the complex plane )
| 01 LBL "AEW"
02 STO 08 03 X<>Y 04 STO 02 05 CLX 06 STO 01 07 STO 04 08 STO 07 09 SIGN 10 STO 03 11 STO 05 12 STO 06 13 LBL 01 14 ISG 01 15 CLX 16 RCL 02 17 RCL 01 18 ST+ X 19 .5 20 - 21 R-P 22 1 23 LASTX 24 - |
25 RCL 02
26 R-P 27 ST* Z 28 X<> T 29 + 30 RCL 04 31 + 32 STO 04 33 X<>Y 34 RCL 03 35 * 36 STO 03 37 RCL 05 38 RCL 08 39 X^2 40 ST+ X 41 * 42 RCL 01 43 * 44 STO 05 45 / 46 P-R 47 RCL 06 48 + |
49 STO 06
50 LASTX 51 - 52 ABS 53 X<>Y 54 RCL 07 55 + 56 STO 07 57 LASTX 58 - 59 ABS 60 + 61 X#0? 62 GTO 01 63 RCL 02 64 .5 65 XEQ "GAMZ" 66 R-P 67 CLX 68 2 69 / 70 RCL 08 71 X^2 72 PI |
73 +
74 4 75 / 76 RCL 02 77 RCL 08 78 LN 79 * 80 - 81 R-D 82 + 83 RCL 08 84 SIGN 85 X<0? 86 CLX 87 ASIN 88 + 89 1 90 P-R 91 RCL 07 92 * 93 X<>Y 94 RCL 06 95 * 96 + |
97 RCL 02
98 PI 99 * 100 E^X 101 ENTER^ 102 X^2 103 1 104 + 105 SQRT 106 + 107 RCL 08 108 SIGN 109 CHS 110 Y^X 111 ST+ X 112 RCL 08 113 / 114 SQRT 115 * 116 END |
( 138 bytes / SIZE 009 )
| STACK | INPUTS | OUTPUTS |
| Y | a | / |
| X | x | W(a;x) |
Example:
0.5 ENTER^ 10
XEQ "AEW" >>>> 0.092208658 (
in 59 seconds )
-0.5 ENTER^ 10
R/S >>>> -0.228640282
--------------
-The series diverge too soon if x is too small, but if you add RND after line 60
FIX 5 1 ENTER^ 5
XEQ "AEW" yields W(1;5) = 0.022808
and similarly W(-1;5) = -0.570255 ( in 57
seconds )
3°) Dn(x)
a) Series
Expansion - via Kummer's Function
-In fact we have: D -a-1/2(x) = U(a,x) therefore:
Dn(x) = 2n/2 Pi1/2
exp(-x2/4) [ 1/Gam((1-n)/2) M( -n/2 , 1/2 , x2/2
) - 21/2 ( x / Gam(-n/2) ) M [ (1-n)/2 , 3/2 , x2/2
]
Data Registers: R00 thru R05: temp
Flags: /
Subroutines: "1/G" or "1/G+" or a similar
M-Code function ( cf "Gamma Function for the HP-41" )
"KUM" ( cf "Kummer's Functions for the HP-41" )
| 01 LBL "DNX"
02 STO 03 03 X^2 04 X<>Y 05 .5 06 STO 02 07 ST* Z 08 * 09 STO 04 10 CHS |
11 STO 01
12 X<>Y 13 XEQ "KUM" 14 STO 05 15 RCL 02 16 ST+ 01 17 ISG 02 18 RCL 01 19 XEQ "1/G" 20 ST* 05 |
21 RCL 00
22 XEQ "KUM" 23 2 24 SQRT 25 * 26 ST* 03 27 RCL 04 28 CHS 29 XEQ "1/G" 30 ST* 03 |
31 RCL 05
32 RCL 03 33 - 34 RCL 00 35 2 36 / 37 E^X 38 / 39 2 40 RCL 04 |
41 Y^X
42 * 43 PI 44 SQRT 45 * 46 END |
( 77 bytes / SIZE 006 )
| STACK | INPUTS | OUTPUTS |
| Y | n | / |
| X | x | Dn(x) |
Example:
0.4 ENTER^
1.8 XEQ "DNX" >>>>
D0.4(1.8) = 0.579579485
---Execution time = 27s---
b) Asymptotic
Expansion
Formula: Dn(x) ~
xn exp(-x2/4) [ 1 - n(n-1) / ( 2 x2 )
+ n(n-1)(n-2)(n-3) / ( 2 ( 4 x4 ) ) - ....... ]
Data Registers: R00 thru R04: temp
Flags: /
Subroutines: /
| 01 LBL "AEDNX"
02 STO 00 03 X<>Y 04 STO 01 05 CLX 06 STO 03 07 SIGN 08 STO 02 09 STO 04 10 LBL 01 |
11 RCL 01
12 RCL 03 13 ST+ X 14 - 15 1 16 ST+ 03 17 LASTX 18 + 19 RCL 01 20 - |
21 *
22 RCL 02 23 * 24 RCL 03 25 / 26 RCL 00 27 X^2 28 ST+ X 29 / 30 STO 02 |
31 RCL 04
32 + 33 STO 04 34 LASTX 35 X#Y? 36 GTO 01 37 RCL 00 38 RCL 01 39 Y^X 40 * |
41 RCL 00
42 2 43 / 44 X^2 45 CHS 46 E^X 47 * 48 END |
( 62 bytes / SIZE 005 )
| STACK | INPUTS | OUTPUTS |
| Y | n | / |
| X | x | Dn(x) |
Examples:
4.5 ENTER^
5 XEQ "AEDNX" >>>>
D4.5(5) = 1.879976816
---Execution time = 9s---
PI 4.7 R/S >>>> Dpi(4.7) = 0.437982402 ---Execution time = 12s---
PI CHS 10 R/S >>>> D -pi(10) = 9.418973196 E-15 ---Execution time = 12s---
Notes:
-Of course, we'll fall into an infinite loop if x is too small.
-If you have the M-Code routine HGF+ ( cf "Hypergeometric Functions
for the HP-41" )
this program may be simplified as follows:
| 01 LBL "AEDN"
02 STO 00 03 X<>Y 04 STO 03 05 .5 06 STO 01 |
07 CHS
08 * 09 ST+ 01 10 STO 02 11 CLST 12 2 |
13 STO Z
14 RCL 00 15 X^2 16 CHS 17 STO 04 18 / |
19 HGF+
20 RCL 00 21 RCL 03 22 Y^X 23 * 24 RCL 04 |
25 4
26 / 27 E^X 28 * 29 END |
( 42 bytes / SIZE 005 )
| STACK | INPUTS | OUTPUTS |
| Y | n | / |
| X | x | Dn(x) |
Examples:
4.5 ENTER^
5 XEQ "AEDN" >>>>
D4.5(5) = 1.879976816
---Execution time = 4s---
PI 4.7 R/S >>>> Dpi(4.7) = 0.437982402 ---Execution time = 5s---
PI CHS 10 R/S >>>> D -pi(10) = 9.418973196 E-15 ---Execution time = 5s---
-The results are the same but they are obtained much faster.
-If an infinite loop occurs because x is too small, press any key to
stop it.
-Here, we have used the relation:
Dn(x) ~ xn exp(-x2/4)
2F0 [ (1-n)/2 , -n/2 ;; -2/x2 ]
References:
[1] Abramowitz and Stegun , "Handbook of Mathematical Functions"
- Dover Publications - ISBN 0-486-61272-4
[2] http://functions.wolfram.com/