The Museum of HP Calculators

Parabolic Cylinder Functions for the HP-41

This program is Copyright © 2006 by Jean-Marc Baillard and is used here by permission.

This program is supplied without representation or warranty of any kind. Jean-Marc Baillard and The Museum of HP Calculators therefore assume no responsibility and shall have no liability, consequential or otherwise, of any kind arising from the use of this program material or any part thereof.


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
 
 

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"  ( Gamma Function )
                        "KUM"  ( Kummer's Functions )
 

001  LBL "PCF1"
002  STO 03
003  X^2
004  X<>Y
005  STO 04
006  .5
007  STO 02
008  ST* Y
009  ST* Z
010  X^2
011  LASTX
012  +
013  +
014  STO 01
015  X<>Y
016  ISG 02
017  XEQ "KUM"
018  STO 05
019  .5
020  ST- 01
021  STO 02
022  RCL 00
023  XEQ "KUM"
024  RCL 00
025  RCL 02
026  *
027  CHS
028  E^X
029  ST* 05
030  *
031  2
032  RCL 01
033  Y^X
034  /
035  STO 06
036  RCL 03
037  LASTX
038  RCL 02
039  SQRT
040  *
041  /
042  ST* 05
043  RCL 02
044  RCL 01
045  RCL 04
046  SIGN
047  *
048  +
049  STO 00
050  XEQ "GAM"
051  PI
052  SQRT
053  /
054  RCL 04
055  SIGN
056  Y^X
057  ST/ 06
058  RCL 00
059  RCL 02
060  LASTX
061  *
062  -
063  XEQ "GAM"
064  PI
065  SQRT
066  /
067  RCL 04
068  SIGN
069  Y^X
070  ST/ 05
071  RCL 02
072  RCL 04
073  ABS
074  +
075  XEQ "GAM"
076  STO 02
077  1
078  CHS
079  ACOS
080  STO 00
081  ST* 01
082  RCL 04
083  X<0?
084  GTO 01
085  RCL 06
086  RCL 05
087  -
088  STO Y
089  RCL 00
090  RCL 04
091  *
092  SIN
093  *
094  RCL 05
095  RCL 06
096  +
097  +
098  RCL 02
099  *
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,    -----------------  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"
 

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                 the recurrence relation is applied until  2x+a' <= 0 but other choices may be better ( for instance | a | < 1 )
45  RCL 04                 especially if line 54 is replaced by  XEQ "PCF1"
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 )
 

Reference:        Abramowitz and Stegun , "Handbook of Mathematical Functions" -  Dover Publications -  ISBN  0-486-61272-4
 

Go back to the HP-41 software library
Go back to the general software library
Go back to the main exhibit hall