hp41programs

Hyperellipsoid

Hypersurface Area of a Hyperellipsoid for the HP-41


Overview
 

-The area A of an hypersurface defined by an equation of the form  xN = f( x1 , ............. , xN-1 )  may be calculated by the (N-1) multiple integral:

    A = § ...... §  [ 1 + ( f / x1 )2 + ................. + ( f / xN-1 )2 ] 1/2 dx1  ............. dxN-1                  ( § = Integral )

-Applying this formula to an hyperellipsoid    x12 / a12 + ............... +  xN2 / aN2 = 1     leads to  ( after several changes of variables )
 

   A =  2N a1 ........ aN-1   §0PI/2 ..... §0PI/2 CosN-2 µ1 CosN-3 µ2 .......... Cos µN-1  [ g( µ1 , ....... , µN-1 ) ] 1/2  dµ1  .......... dµN-1

    where   g( µ1 , ....... , µN-1 ) = 1 + (aN2/a12 - 1) Sin2 µ1 + (aN2/a22 - 1) Cos2 µ1 Sin2 µ2 + ......  + (aN2/aN-12 - 1) Cos2 µ1 Cos2 µ2 .... Sin2 µN-1
 

-For N = 4 , we have to compute a triple integral which remains possible in a "reasonnable" time
-For larger N-values, the execution time rapidly becomes prohibitive.

-Fortunately, in reference [1], a method is described to transform the (N-1) multiple integral into a univariate integral.
-With minor modifications and a few simplifications, it yields:

  A =  [ a2 ........ aN / Gamma((N+1)/2) ]   §01  f(u) du

       with     f(u) = (1-u)(N-1)/2 u -1/2  [ 1 + SUMi=2,...,N (a12/ai2/bi) ] ( b2 .....  bN ) -1/2        where        bi  =  1 - u + u (a12/ai2)
 

-To remove the singularities for u = 0 & u = 1 , we can use the change of variable  u = Sin2 v  and apply a Gauss-Legendre formula.
-In reference [1] , the authors employ another change of variable and the Romberg method.

-The following programs use a Gauss-Chebyshev formula, namely:

     §01  [ (1-x) / x ]1/2  g(x) dx  =  SUMi=1,2,...,n  wi g(xi)

-Unlike many Gaussian formulae, the abscissas & weights may be easily calculated:

    xi =  Sin2 [ (2i-1) / (2n+1) ] PI/2    &   wi = [ 2.PI / (2n+1) ] Cos2 [ (2i-1) / (2n+1) ] PI/2

-So, we can choose n large enough to get an almost perfect accuracy.
 

Program Listing
 
 

Data Registers:           •  R10 = N           ( Registers R10 thru R10+NN are to be initialized before executing "SAHE" )

                                      •  R11 = a1         •  R12 = a2    .......................      •  R10+NN = aN

      When the program stops:           R01 = An                                                R00 & R02 thru R09:  temp

Flags: /
Subroutine:  "1/G+" or "GAM+"  ( cf "Gamma Function for the HP-41 )
 
 
 

 01  LBL "SAHE"
 02  DEG
 03  STO 02
 04  RCL 10
 05  10.011
 06  +
 07  STO 09
 08  PI
 09  RCL 10 
 10  1
 11  +
 12  STO 08        
 13  Y^X
 14  SQRT
 15  ST+ X
 16  LBL 00
 17  RCL IND Y
 18  *
 19  DSE Y
 20  GTO 00
 21  STO 00
 22  RCL 08
 23  2
 24  /
 25  XEQ "1/G+"
 26  ST* 00
 27  LBL 01
 28  RCL 02        
 29  ST+ X
 30  1
 31  +
 32  STO 07
 33  90
 34  /
 35  STO 08
 36  CLX
 37  STO 01
 38  LBL 02
 39  RCL 02
 40  ST+ X
 41  DSE X
 42  RCL 08 
 43  /
 44  ENTER^
 45  SIGN
 46  P-R
 47  X^2
 48  STO 04        
 49  X<>Y
 50  X^2
 51  STO 05
 52  RCL 09
 53  STO 03
 54  SIGN
 55  STO 06
 56  LBL 03
 57  RCL 11
 58  RCL IND 03
 59  / 
 60  X^2
 61  RCL X
 62  RCL 05 
 63  *
 64  RCL 04        
 65  +
 66  ST* 06
 67  /
 68  +
 69  DSE 03
 70  GTO 03
 71  RCL 04
 72  RCL 10
 73  2
 74  /
 75  Y^X
 76  *
 77  RCL 06
 78  SQRT
 79  /
 80  ST+ 01
 81  DSE 02
 82  GTO 02
 83  RCL 00        
 84  RCL 07 
 85  /
 86  ST* 01
 87  RCL 01
 88  RTN
 89  STO 02
 90  GTO 01
 91  END

 
    ( 129 bytes / SIZE NNN+11 )
 
 

      STACK        INPUT      OUTPUT
           X             n            An

  where   n = number of points in the Gauss-Chebyshev formula
   and    An = corresponding approximation of the area

Examples:

  •  N = 4      a1 = 4 ,  a2 = 5 ,  a3 = 6 ,  a4 = 7

  N = 4  STO 10

     4  STO 11
     5  STO 12
     6  STO 13
     7  STO 14

-With  n = 4    4  XEQ "SAHE"  >>>>   3194.724649           ---Execution time = 19s---
                       8        R/S            >>>>   3194.584863                       in  32s
                      16       R/S            >>>>   3194.584860                       in  64s

-So    A = 3194.584860   is probably very close to the exact result.

  •  N = 5  ,  semi-axes =  1  2  4  8  16          Store these 6 numbers into R10 to R15

-With  n = 50    50  XEQ "SAHE"  >>>>   12926.73486        ---Execution time = 3m40s---
                        100        R/S          >>>>    12926.73511

-The exact value is given in reference [1]     A = 12926.73510  ( rounded to 5 decimals )

  •  N = 10  ,  semi-axes =  1  2  4  8  16  32  64  128  256  512          Store these 11 numbers into R10 to R20

-With  n = 50    50  XEQ "SAHE"  >>>>   2.971355318 E14        ---Execution time = 6m15s---
                        100        R/S          >>>>   2.971355398 E14

-The exact result is also given in reference [1]      A = 2.971355398 E14  ( rounded to 10 significant digits )

Notes:

-If you replace line 25 by  XEQ "GAM+" , replace line 26 by  ST/ 00

-Lines 27-88-89-90 just avoid to recalculate R00 at each new n-value
-Otherwise, they could be deleted, thus saving a few bytes.

-It seems preferable to store the semi-axes in increasing order into R11  R12  .........
-Gauss-Chebyshev formula performs well, even compared to Romberg method !

Variant:

-If  N < 10 , we can store the semi-axes into R01 thru Rnn , which is probably more "natural"
-The variant hereunder is also slightly shorter:
 

Data Registers:           •  R00 = N < 10          ( Registers R00 thru RNN are to be initialized before executing "SAHE" )

                                      •  R01 = a1         •  R02 = a2    .......................      •  RNN = aN

      When the program stops:           R10 = An  ,  R11 to R19:  temp         ( R12 = 2n+1 )

Flags: /
Subroutine:  "1/G+" or "GAM+"  ( cf "Gamma Function for the HP-41 )
 
 

 01  LBL "SAHE"
 02  DEG
 03  STO 11
 04  ST+ X
 05  1
 06  +
 07  STO 12        
 08  90
 09  /
 10  STO 13
 11  RCL 00
 12   E-3
 13  +
 14  STO 16
 15  CLX
 16  STO 10
 17  LBL 01
 18  RCL 11
 19  ST+ X
 20  DSE X
 21  RCL 13        
 22  /
 23  ENTER^
 24  SIGN
 25  P-R
 26  X^2
 27  STO 14
 28  X<>Y
 29  X^2
 30  STO 15
 31  RCL 16
 32  STO 17
 33  SIGN
 34  STO 18
 35  LBL 02
 36  RCL 01
 37  RCL IND 17
 38  /
 39  X^2
 40  RCL X
 41  RCL 15        
 42  *
 43  RCL 14
 44  +
 45  ST* 18
 46  /
 47  +
 48  DSE 17
 49  GTO 02
 50  RCL 14
 51  RCL 00
 52  2
 53  /
 54  Y^X
 55  *
 56  RCL 18
 57  SQRT
 58  /
 59  ST+ 10
 60  DSE 11
 61  GTO 01
 62  RCL 10        
 63  PI
 64  RCL 00
 65  1
 66  +
 67  2
 68  /
 69  STO 13
 70  Y^X
 71  *
 72  ST+ X
 73  LBL 03
 74  RCL IND 16
 75  *
 76  DSE 16
 77  GTO 03
 78  STO 10        
 79  RCL 13
 80  XEQ "1/G+"
 81  RCL 12
 82  /
 83  ST* 10
 84  RCL 10
 85  END

 
    ( 123 bytes / SIZE 019 )
 
 

      STACK        INPUT      OUTPUT
           X             n            An

  where   n = number of points in the Gauss-Chebyshev formula
   and    An = corresponding approximation of the area

Examples:

-The same examples give the same results, with perhaps a small difference in the last decimal.
-For the 3rd example ( with N = 10 ) , simply replace R10 by R19 in the listing above.

-Another one:  N = 5  ,   semi-axes =  4  5  6  7  8

   5  STO 00

   4  STO 01
   5  STO 02
   6  STO 03
   7  STO 04
   8  STO 05

-With  n = 30  ,   30   XEQ "SAHE"  >>>>  A =  31971.28857           ---Execution time = 137s---
-With  n = 40  ,   40           R/S          >>>>  A =  31971.28858           ---Execution time = 181s---
 

Notes:

-If you prefer  XEQ "GAM+"  ( line 80 ) , replace lines 81-82-83 by  RCL 12   *   ST/ 10
-Here, the loop LBL 03 is re-executed for each new n-value, which makes the program a little slower.
 

References:

[1]   Dunkl & Ramirez - "Computing Hyperelliptic Integrals for Surface Measure of Ellipsoids"
[2]   Abramowitz and Stegun - "Handbook of Mathematical Functions" -  Dover Publications -  ISBN  0-486-61272-4