Hypersurface Area of a Hyperellipsoid for the HP-41

Overview
 

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

    A = § ...... §  [ 1 + ( ¶f / ¶x₁ )² + ................. + ( ¶f / ¶x_N-1 )² ] ^1/2 dx₁  ............. dx_N-1                  ( § = Integral )

-Applying this formula to an hyperellipsoid    x₁² / a₁² + ............... +  x_N² / a_N² = 1     leads to  ( after several changes of variables )
 

   A =  2^N a₁ ........ a_N-1   §₀^PI/2 ..... §₀^PI/2 Cos^N-2 µ₁ Cos^N-3 µ₂ .......... Cos µ_N-1  [ g( µ₁ , ....... , µ_N-1 ) ] ^1/2  dµ₁  .......... dµ_N-1

    where   g( µ₁ , ....... , µ_N-1 ) = 1 + (a_N²/a₁² - 1) Sin² µ₁ + (a_N²/a₂² - 1) Cos² µ₁ Sin² µ₂ + ......  + (a_N²/a_N-1² - 1) Cos² µ₁ Cos² µ₂ .... Sin² µ_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 = (PI)^(N-1)/2 [ a₂ ........ a_N / Gamma((N+1)/2) ]   §₀¹  f(u) du

       with     f(u) = (1-u)^(N-1)/2 u ^-1/2  [ 1 + SUM_i=2,...,N (a₁²/a_i²/b_i) ] ( b₂ .....  b_N ) ^-1/2        where        b_i  =  1 - u + u (a₁²/a_i²)
 

-To remove the singularities for u = 0 & u = 1 , we can use the change of variable  u = Sin² 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:

     §₀¹  [ (1-x) / x ]^1/2  g(x) dx  =  SUM_i=1,2,...,n  w_i g(x_i)

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

    x_i =  Sin² ([ (2i-1) / (2n+1) ] PI/2)    &   w_i = [ 2.PI / (2n+1) ] Cos² ([ (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 = a₁         •  R12 = a₂    .......................      •  R10+NN = a_N

      When the program stops:           R01 = A_n                                                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    A_n = corresponding approximation of the area

Examples:

  •  N = 4      a₁ = 4 ,  a₂ = 5 ,  a₃ = 6 ,  a₄ = 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 = a₁         •  R02 = a₂    .......................      •  RNN = a_N

      When the program stops:           R10 = A_n  ,  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    A_n = 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