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

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