# hp41programs

ParabolicCylinder

# 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/