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°)  D_n(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   d²y/dx² - ( x²/4 + a ).y = 0  may be expressed with the Kummer's Functions

     U(a;x) = pi^1/22^-1/4-a/2e^{-x^2/4} (Gam(3/4+a/2))^-1 M(1/4+a/2;1/2;x²/2) - pi^1/22^1/4-a/2 x.e^{-x^2/4} (Gam(1/4+a/2))^-1 M(3/4+a/2;3/2;x²/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) = Y₁ cos(pi(1/4+a/2)) - Y₂ sin(pi(1/4+a/2))     and     V(a;x) = (Gam(1/2-a))^-1 [ Y₁ sin(pi(1/4+a/2)) + Y₂ cos(pi(1/4+a/2)) ]

  where  Y₁ = pi^-1/22^-1/4-a/2e^{-x^2/4} (Gam(1/4-a/2)) M(1/4+a/2;1/2;x²/2)  and   Y₂ = pi^-1/22^1/4-a/2 x.e^{-x^2/4} (Gam(3/4-a/2)) M(3/4+a/2;3/2;x²/2)
 
 

Data Registers:  R00 thru R06  ( R03 = x ; R04 = a )
Flags: /
Subroutines:  "GAM" or  "GAM+" ...  ( Gamma Function )  ,  "KUM"  ( Kummer's Functions )
 
 

 
   ( 161 bytes / SIZE 007 )
 
 

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.(2x²)) + (a+1/2)(a+3/2)(a+5/2)(a+7/2)/(2.(2x²)²) -  ....... ]

                   V(a;x)  ~   (2/pi)^1/2  e^{x^2/4}  x^a-1/2 [ 1 + (a-1/2)(a-3/2)/(1.(2x²)) + (a-1/2)(a-3/2)(a-5/2)(a-7/2)/(2.(2x²)²) + ....... ]
 

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

 
   ( 84 bytes / SIZE 004 )
 
 

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 10¹²      ( in 7 seconds )

Notes:

-V(2;10) is accurately computed by "PCF1" ( it yields  1.823604925 10¹² 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"
 
 

 
   ( 81 bytes / SIZE 008 )
 
 

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   d²y/dx² + ( x²/4 - a ).y = 0  is solved by a series expansion   y = a₀ + a₁.x + a₂.x² + ...... + a_k.x^k + ......
-The standard solution  W(a;x) is defined by

       a₀ = 2^-3/4 | Gam(1/4 + i.a/2) / Gam(3/4 + i.a/2) |^1/2   and   a₁ = -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 )
 
 

 
   ( 100 bytes / SIZE 010 )
 
 

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  ( s₁(a;x) cos A - s₂(a;x) sin A )      and     W(a;-x)  ~  (2/(kx))^1/2 ( s₁(a;x) sin A + s₂(a;x) cos A )      ( x > 0 )

          where        A = x²/2 - a.Ln(x) + pi/4 + B/2   with  B = arg ( Gam(1/2+i.a) )
                            k = 1/(e^pi.a+(1+e^2pi.a)^1/2)

     s₁(a;x)+i.s₂(a;x)  =  Sum_{n=0,1,2,......}  (-i)ⁿ Gam(2n+1/2+i.a)/Gam(1/2+i.a)  1/(n!(2x²)ⁿ)
 

Data Registers:  R00 thru R08: temp  when the program stops, R02 = a ; R08 = x
Flags: /
Subroutine:  "GAMZ"  ( gamma function in the complex plane )
 
 

 
  ( 138 bytes / SIZE 009 )
 
 

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°)  D_n(x)
 

        a) Series Expansion - via Kummer's Function
 

-In fact we have:  D _-a-1/2(x) = U(a,x)   therefore:

    D_n(x) = 2^n/2 Pi^1/2 exp(-x²/4) [ 1/Gam((1-n)/2) M( -n/2 , 1/2 , x²/2 ) - 2^1/2 ( x / Gam(-n/2) ) M [ (1-n)/2 , 3/2 , x²/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" )
 
 

 
  ( 77 bytes / SIZE 006 )
 
 

Example:

   0.4   ENTER^
   1.8   XEQ "DNX"  >>>>   D_0.4(1.8) = 0.579579485        ---Execution time = 27s---
 

       b) Asymptotic Expansion
 

Formula:    D_n(x)  ~  xⁿ exp(-x²/4) [ 1 - n(n-1) / ( 2 x² ) + n(n-1)(n-2)(n-3) / ( 2 ( 4 x⁴ ) ) - ....... ]
 

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

 
   ( 62 bytes / SIZE 005 )
 
 

Examples:

   4.5   ENTER^
    5     XEQ "AEDNX"  >>>>   D_4.5(5) = 1.879976816        ---Execution time = 9s---

   PI   4.7    R/S    >>>>   D_pi(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:
 
 

 
   ( 42 bytes / SIZE 005 )
 
 

Examples:

   4.5   ENTER^
    5     XEQ "AEDN"  >>>>   D_4.5(5) = 1.879976816        ---Execution time = 4s---

   PI   4.7    R/S    >>>>   D_pi(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:

   D_n(x)  ~  xⁿ exp(-x²/4) ₂F₀ [ (1-n)/2 , -n/2 ;; -2/x² ]
 

References:

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