The Museum of HP Calculators

Spheroidal Harmonics 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°) Eigenvalues
 2°) Spheroidal Wave Functions ( | x | <= 1 )

-We want to solve the differential equation   ( 1 - x² ) d²S/dx² - 2.x  dS/dx + ( L_mn - c² x² - m²/( 1 - x² ) S = 0

  where  m , n are non-negative integers:   m = 0 , 1 , 2 , 3 , .......    n = m , m + 1 , m + 2 , ...... label the successive eigenvalues L_mn
  and  c² may be positive or negative ( i-e  c may be purely imaginary )
  c² > 0  corresponds to the prolate spheroidal wave functions
  c² < 0  ------------------  oblate   -------------------------           ( c = 0  leads to the associated Legendre functions )

-But first, we have to find the numbers  L_mn(c²)  so that the solutions are regular at  x = +/- 1
  ( if c = 0 , the eigenvalues are simply  n ( n + 1 ) )
 

1°) Eigenvalues
 

-The program below finds L_mn by solving the transcendental equation   U₁(L_mn) + U₂(L_mn) = 0   where    U₁ & U₂  are 2 continued fractions:

                                                 -b_r^m                      -b_r-2^m
   U₁(L_mn) = g_r^m - L_mn +    ---------------       ---------------   ..............
                                           g_r-2^m - L_mn +         g_r-4^m - L_mn +

                                                                                                                                r = n - m
                                     -b_r+2^m                    -b_r+4^m
   U₂(L_mn) =           ---------------       ----------------   ..............
                               g_r+2^m - L_mn +         g_r+4^m - L_mn +
 

   with      b_r^m  =  [ r.(r-1).(2m+r).(2m+r-1).c⁴ ] / [ (2m+2r-1)².(2m+2r+1).(2m+2r-3) ]
   and      g_r^m  =  (m+r).(m+r+1) + (c²/2).[ 1 - (4m²-1)/((2m+2r-1).(2m+2r+3)) ]
 

Data Registers:             R00 to R08 are used by "CFR"
                                        R09 = n-m   R11 = m   R15 thru R19: unused
                                        R10 = c²     R12 to R14: temp   R20 thru R23 are used by "SOLVE"
Flag:   F01
Subroutines:  "CFR"       ( cf "Continued Fractions for the HP-41" )
                        "SOLVE"  ( cf "Linear and non-linear Systems for the HP-41" ) after replacing registers R00 thru R03 by R20 thru R23 ( for instance )

-In other words:     LBL "SOLVE"     STO 23               ENTER^        /                       -                  GTO 01
                             STO 21                LBL 01               ENTER^        RCL 22            *                  END
                             X<>Y                  VIEW 21            X<> 23          RCL 21            +
                             STO 22                RCL 21              -                     STO 22            STO 21        ( 27 lines / 51 bytes )
                             XEQ IND 20       XEQ IND 20      X#0?              STO T              X#Y?
 

  01  LBL "LMN"
  02  STO 10
  03  RDN
  04  STO 09
  05  X<>Y
  06  STO 11
  07  -
  08  X<> 09
  09  ENTER^
  10  "T"                               ( the same character as the LBL line 27 )
  11  1
  12  +
  13  *
  14  RCL 10
  15  2
  16  /
  17  +
  18  ENTER^
  19  ASTO 20
  20  1
  21  +
  22  "U"                               ( the same character as the LBL line 43 )
  23  ASTO 00
  24  CLA
  25  XEQ "SOLVE"
  26  RTN
  27  LBL "T"
  28  STO 01
  29  RCL 09
  30  STO 12
  31  XEQ 01
  32  X<> L
  33  SF 01
  34  XEQ "CFR"
  35  STO 14
  36  CLX
  37  CF 01
  38  RCL 01
  39  XEQ "CFR"
  40  RCL 14
  41  +
  42  RTN
  43  LBL "U"
  44  RCL 09
  45  RCL 02
  46  1
  47  FC? 01
  48  CLX
  49  -
  50  ST+ X
  51  FS? 01
  52  CHS
  53  +
  54  STO 12
  55  1
  56  RCL 12
  57  -
  58  *
  59  RCL 11
  60  ST+ X
  61  RCL 12
  62  +
  63  ST* Y
  64  1
  65  -
  66  *
  67  RCL 10
  68  X^2
  69  *
  70  RCL 11
  71  RCL 12
  72  +
  73  ST+ X
  74  3
  75  -
  76  ST/ Y
  77  4
  78  +
  79  ST/ Y
  80  2
  81  FS? 01
  82  ST- 12
  83  -
  84  X^2
  85  /
  86  STO 13
  87  LBL 01
  88  1
  89  RCL 11
  90  ST+ 12
  91  ST+ X
  92  X^2
  93  -
  94  RCL 12
  95  ST+ X
  96  1
  97  -
  98  ST/ Y
  99  4
100  +
101  /
102  1
103  +
104  RCL 10
105  *
106  2
107  /
108  RCL 12
109  RCL 12
110  1
111  +
112  *
113  +
114  RCL 01
115  -
116  RCL 13
117  RTN
118  END

( 172 bytes / SIZE 024 )
 
 

STACK	INPUTS	OUTPUTS
Z	m	Lmn
Y	n	Lmn
X	c2	Lmn

Example1:      m = 4 ,  n = 11 ,  c² = -1

    4   ENTER^
   11  ENTER^
   -1  XEQ "LMN"   >>>>   L_mn = 131.5600809    ( in 1mn27s )     ( the HP-41 displays the successive approximations )

Example2:      m = 0 ,  n = 0 ,  c² = 3

   0   ENTER^  ENTER^
   3   XEQ "LMN"  >>>>   L_mn = 0.879933945    ( in 2mn30s )

Notes:

-This program uses 2 initial guesses:   L_i = n(n+1) + c²/2   and   L'_i = 1 + L_i   ( lines 02 to 21 )
-A better approximation is    L_i = n(n+1) + ( c²/2 ) [ 1 - (2m-1).(2m+1)/((2n-1).(2n+3)) ]
-More accurate L-approximations may be found in the "Handbook of Mathematical Functions" but the formulas are quite complex.

-Too bad guesses can produce correct eigenvalues, but corresponding to another value of n.
-For example, if  m = n = 0 and c² = -16 ,  L_i = 0 and  L'_i = 1  would yield  0.221407906  which is actually  L_mn for  m = 0 , n = 2 , c² = -16
-So, it could be preferable to improve the 2 guesses, especially for large c-values.

-Here are a few results obtained by "LMN" and Warren Furlow's excellent V41 emulator ( in "turbo" mode )
 ( of course - due to roundoff errors - the last decimal may be doubtful )
 
 

m     n      c2	Lmn	m     n      c2	Lmn	m     n      c2	Lmn	m     n      c2	Lmn
0     0      1	0.319000055	0     0     -1	-0.348602400	0     1      1	2.593084580	0     1      -1	1.393206310
4	1.127734066	-4	-1.594493214	4	4.287128544	-4	-0.505243981
9	2.136732226	-9	-4.343292705	9	6.820888328	-9	-3.899400290
16	3.172067425	-16	-9.150793382	16	9.805943842	-16	-9.025710674
25	4.195128875	-25	-16.07904274	25	12.91170325	-25	-16.05041268
0     2     1	6.533471801	0     2      -1	5.486800054	0     3     1	12.51446214	0     3     -1	11.49212090
4	8.225713002	-4	4.091509101	4	14.10020387	-4	10.00386381
9	11.19293865	-9	2.251269209	9	16.88903022	-9	7.611465213
16	15.30630000	-16	0.221407908	16	21.04896065	-16	4.339082933
25	20.17691472	-25	-2.448598906	25	26.58735961	-25	0.060929887
1     1      1	2.195548355	1     1     -1	1.795304587	1     2      1	6.424699144	1     2     -1	5.567527453
4	2.734111026	-4	1.118553391	4	7.653149562	-4	4.222747334
9	3.504795867	-9	-0.269420805	9	9.555998398	-9	1.821541436
16	4.399593067	-16	-2.909200759	16	11.94871939	-16	-1.869861296
25	5.350422298	-25	-7.493388287	25	14.64295625	-25	-7.127837527
1     3     1	12.46791533	1     3    -1	11.53481845	1     4     1	20.48169601	1     4     -1	19.52068334
4	13.88149342	-4	10.16324505	4	21.94014372	-4	18.09765523
9	16.23846623	-9	8.007074153	9	24.40831218	-9	15.77725227
16	19.46526054	-16	5.405903145	16	27.91188168	-16	12.62871852
25	23.39761313	-25	2.750367212	25	32.42194360	-25	8.694959247
2     2      1	6.140948992	2     2     -1	5.855162574	2     3      1	12.33110151	2     3     -1	11.66440927
4	6.542495275	-4	5.395010784	4	13.29825047	-4	10.62995129
9	7.151100524	-9	4.526460462	9	14.82777821	-9	8.809393921
16	7.903860949	-16	3.027624006	16	16.81295852	-16	6.046076609
25	8.747674251	-25	0.428957106	25	19.13598192	-25	2.109805814
2     4      1	20.40235305	2     4     -1	19.59722019	2     5      1	30.43614539	2     5     -1	29.56437148
4	21.60513304	-4	18.38832874	4	31.74704320	-4	28.26120113
9	23.58709333	-9	16.38610602	9	33.93615107	-9	26.10501394
16	26.29348662	-16	13.67312110	16	36.99626751	-16	23.12875359
25	29.62878614	-25	10.51236733	25	40.89293269	-25	19.38439052

         ........ and so on ........
 

2°) Spheroidal Wave Functions    ( -1 <= x <= +1 )
 

-The following program computes the angular spheroidal wave function of the first kind by a series expansion in powers of x
-But first, given  m , n and c² , the corresponding eigenvalue L_mn must be calculated by "LMN"
 

Formulae:       S_mn(x)  =  a₀ + a₂ x² + a₄ x⁴ + ........      with  a₀ = [ (-1)^(n-m)/2 (n+m)! ] / [ 2ⁿ ((n-m)/2)! ((n+m)/2)! ]                       if  n - m is even
                        S_mn(x)  =  a₁x + a₃ x³ + a₅ x⁵ + ........    with  a₁ = [ (-1)^(n-m-1)/2 (n+m+1)! ] / [ 2ⁿ ((n-m-1)/2)! ((n+m+1)/2)! ]          if  n - m is odd

-In both cases,   (k+1)(k+2) a_k+2 + ( L_mn - 2.k² - m² ) a_k + ( k² - 3.k + 2 - L_mn - c² ) a_k-2 + c² a_k-4 = 0
 

Data Registers:            R00 = x          R06 thru R10: temp                         ( Registers R01 thru R04 are to be initialized before executing "SMN" )

                                   •  R01 = m      •  R03 = c²        R05 = partial sums and finally  S_mn(x)
                                   •  R02 = n       •  R04 = L_mn
Flag:   F01
Subroutines: /
 

  01  LBL "SMN"
  02  STO 00
  03  SF 01
  04  RCL 02
  05  RCL 01
  06  -
  07  2
  08  MOD
  09  X=0?
  10  CF 01
  11  ST+ 01                                R01 is modified if n-m is odd, but its original content is restored line 51
  12  STO 06
  13  RCL 00
  14  1
  15  FS? 01
  16  X<>Y
  17  STO 05
  18  STO 07
  19  1
  20  CHS
  21  RCL 02
  22  RCL 01
  23  -
  24  2
  25  /
  26  STO 08
  27  Y^X
  28  RCL 01
  29  RCL 02
  30  +
  31  FACT
  32  *
  33  RCL 08
  34  FACT
  35  RCL 02
  36  LASTX
  37  -
  38  FACT
  39  *
  40  2
  41  RCL 02
  42  Y^X
  43  *
  44  /
  45  ST* 05
  46  STO 10
  47  CLX
  48  STO 08
  49  STO 09
  50  FS?C 01
  51  DSE 01
  52  LBL 01
  53  RCL 03
  54  RCL 08
  55  *
  56  RCL 06
  57  3
  58  -
  59  RCL 06
  60  *
  61  2
  62  +
  63  RCL 03
  64  -
  65  RCL 04
  66  -
  67  RCL 09
  68  STO 08
  69  *
  70  +
  71  RCL 04
  72  RCL 06
  73  X^2
  74  ST+ X
  75  -
  76  RCL 01
  77  X^2
  78  -
  79  RCL 10
  80  STO 09
  81  *
  82  +
  83  CHS
  84  RCL 06
  85  1
  86  +
  87  ST/ Y
  88  LASTX
  89  +
  90  STO 06
  91  /
  92  STO 10
  93  RCL 07
  94  RCL 00
  95  X^2
  96  *
  97  STO 07
  98  *                                              The program stops when 2 consecutive sums are equal.
  99  RCL 05                                    To avoid a premature interruption, add  DSE 11   GTO 02  after line 104
100  +                                              and  3  STO 11   LBL 02  after line 52
101  STO 05                                    in case 1 or 2 successive a_k = 0  and the next one is not.
102  LASTX                                    Practically, this seems highly improbable if c#0
103  X#Y?
104  GTO 01
105  END

( 123 bytes / SIZE 011 )
 
 
 

STACK	INPUTS	OUTPUTS
Y	/	Smn(x)
X	x	Smn(x)

Example1:     m = n = 2  and  c² = -25    "LMN"  gives   0.4289571064       We store these 4 numbers into  R01 thru R04  and

    0.6   XEQ "SMN"  >>>>     S_mn(0.6) = 4.564797329   ( in 19 seconds )
    0.9         R/S          >>>>     S_mn(0.9) = 3.188333453   ( in 24 seconds )

Example2:     m = n = 0  and  c² = -16    "LMN"  gives   -9.150793382        We store these 4 numbers into  R01 thru R04  and

    0.7   XEQ "SMN"  >>>>     S_mn(0.7) = 4.557370593   ( in 19 seconds )
     1           R/S          >>>>       S_mn(1)  = 12.41705452   ( in 20 seconds )

Example3:     m = 2 , n = 5  and  c² = 16    "LMN"  gives   36.99626751       We store these 4 numbers into  R01 thru R04  and

    0.3   XEQ "SMN"  >>>>     S_mn(0.3) = -9.214845515   ( in 15 seconds )
    0.7         R/S          >>>>     S_mn(0.7) = 10.51929252     ( in 19 seconds )

Note:

-If m > 0   S_mn(1) = S_mn(-1) = 0  therefore, you could add  ABS  1  X#Y?  GTO 00  RCL 01  0  X<Y?  RTN  LBL 00   after line 02

Reference:      Abramowitz & 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