The Museum of HP Calculators

HP Forum Archive 21

[ Return to Index | Top of Index ]

Perimeter of an Ellipse (HP-42S)
Message #1 Posted by Gerson W. Barbosa on 24 May 2013, 1:22 a.m.

00 { 48-Byte Prgm }
01>LBL "P"
02 RCL- ST Y
03 X<>Y
04 RCL+ ST L
05 /
06 LASTX
07 PI
08 *
09 X<>Y
10 X^2
11 STO ST Z
12 4
13 +
14 RCL* ST Z
15 RCL* ST Z
16 -3
17 *
18 256
19 /
20 R^
21 4
22 /
23 1
24 -
25 1/X
26 -
27 *
28 .END.

Exact formula:

  p = 4*a*E(1-b^2/a^2)
where
  E(x) is the complete elliptic integral of the second kind 
  a is the semi-major axis and
  b is the semi-minor axis
Example:
  a = 2, b = 1
p = 4*2*E(3/4) = 9.68844822054767619842850319639...
http://www.wolframalpha.com/input/?i=4*2*E%281-1%2F4%29

Approximate formula:

  p ~ pi*(a + b)*(4/(4-h) - 3/256*h^2*(4 + h))
where
  h = ((a - b)/(a + b))^2

The approximate formula is a rework of the Infinite Series 2 in this reference, whose terms are a subset of the infinite series SUM(k=0,inf,(h/4)^k), which converges to 4/(4-h). The resulting expression is

   4/(4-h) - 3/64*h^2 - 3/256*h^3 - 39/16384*h^4 - 15/65536*h^5 + 185/1048576*h^6 ...
This was done by hand and haven't been doublechecked. Anyway, the approximation formula above uses only terms up to h^3. The percent error ranges from 0% (circle) to about 0.1% in the worst case (two coincidental lines).
+-----+-----+---------------+---------------+
|  a  |  b  |  approximate  |     exact     |
+-----+-----+---------------+---------------+
|  1  |  0  | 4.00471251023 | 4.00000000000 |
|  2  |  1  | 9.68845167284 | 9.68844822055 |
|  3  |  2  | 15.8654396854 | 15.8654395893 |
|  4  |  3  | 22.1034921699 | 22.1034921607 |
|  5  |  4  | 28.3616678905 | 28.3616678890 |
|  5  |  5  | 31.4159265359 | 31.4159265359 |
+-----+-----+---------------+---------------+

Thanks Eduardo Duenez for his recent post below.

Edited to fix a typo per Ernie's observation below

Edited again to correct an error in the parameter of the elliptic function per Eduardo's observatin below

Edited: 24 May 2013, 4:44 p.m. after one or more responses were posted

      
Re: Perimeter of an Ellipse (HP-42S)
Message #2 Posted by Ernie Brock on 24 May 2013, 2:24 a.m.,
in response to message #1 by Gerson W. Barbosa

I don't have an HP-42S, but you may be able to improve on this by computing the line integral from 0 to pi/2 using the numeric integrator on the 42S

p = 4 * a * integral(sqrt(1-(a*a-b*b)/(a*a)*sin^2(t), t, 0, pi/2)

http://www.wolframalpha.com/input/?i=4+*+3+*+integrate+sqrt%281-%283*3-2*2%29%2F%283*3%29*sin%5E2%28t%29%29%2Ct%2C0%2Cpi%2F2

BTW quoted "exact" answer for a=3,b=2 should be 15.8654395893 not 15.8654396893

            
Re: Perimeter of an Ellipse (HP-42S)
Message #3 Posted by Gerson W. Barbosa on 24 May 2013, 2:43 a.m.,
in response to message #2 by Ernie Brock

Problem is numerical integration is rather slow to full accuracy on the HP-42S, but I'll check this later.

Quote:
BTW quoted "exact" answer for a=3,b=2 should be 15.8654395893 not 15.8654396893

Fixed, thanks!

            
Re: Perimeter of an Ellipse (HP-42S)
Message #4 Posted by Gerson W. Barbosa on 24 May 2013, 7:48 a.m.,
in response to message #2 by Ernie Brock

Quote:
I don't have an HP-42S, but you may be able to improve on this by computing the line integral from 0 to pi/2 using the numeric integrator on the 42S

p = 4 * a * integral(sqrt(1-(a*a-b*b)/(a*a)*sin^2(t), t, 0, pi/2)


The drawback is a somewhat long running time, but that's an option.

Regards,

Gerson.

00 { 41-Byte Prgm } 01>LBL "E2" 02 MVAR "A" 03 MVAR "B" 04 MVAR "T" 05 4 06 RCL* "A" 07 1 08 1 09 RCL "B" 10 RCL/ "A" 11 X^2 12 - 13 RCL "T" 14 SIN 15 X^2 16 * 17 - 18 SQRT 19 * 20 .END.

LLIM = 0 ULIM = 1.5707963268 (pi/2) ACC = 0.0000000010
+-----+-----+---------------+------+ | a | b | integral | t(s) | +-----+-----+---------------+------+ | 1 | 0 | 3.99999999999 | 44.9 | | 2 | 1 | 9.68844822056 | 45.6 | | 3 | 2 | 15.8654395893 | 46.2 | | 4 | 3 | 22.1034921608 | 45.7 | | 5 | 4 | 28.3616678889 | 45.8 | | 5 | 5 | 31.4159265360 | 2.9 | +-----+-----+---------------+------+

                  
Re: Perimeter of an Ellipse (HP-42S)
Message #5 Posted by hugh steers on 24 May 2013, 9:21 a.m.,
in response to message #4 by Gerson W. Barbosa

Hi. Good stuff!

Here's another way i found a while back, using a variation of arithmetic and geometric means. The iteration is fairly rapid and small to code.

This program is in C, but i hope you can see the method:

#include <math.h>
#include <stdio.h>

double agm2(double b) { double s = (1 + b*b)/2; double a = 1; double t = 1;

for (;;) { double a1 = (a + b)/2; if (a == a1) break; double c = (a - b)/2; b = sqrt(a*b); a = a1; s = s - t*c*c; t = t * 2; } return s/a; }

#define _PI 3.141592653589793238462643

double ellipse(double a, double b) { // perimeter of an ellipse // a semi-major axis // b semi-minor axis // e = eccentricity = sqrt(1-(b/a)^2) // p = 4*a*E(e)

// use AGM variant return 2*_PI*a*agm2(b/a); }

                        
Re: Perimeter of an Ellipse (HP-42S)
Message #6 Posted by Eduardo Duenez on 24 May 2013, 10:49 a.m.,
in response to message #5 by hugh steers

The algorithm used in Hugh Steers' C program is Semjon Adlaj's, upon which my earlier post for the WP34s is based.

(As an academic I feel the urge to always credit original ideas to their source.)

Eduardo

                              
Re: Perimeter of an Ellipse (HP-42S)
Message #7 Posted by Thomas Klemm on 26 May 2013, 5:42 a.m.,
in response to message #6 by Eduardo Duenez

In your program ELC you use both, the built-in routine AGM and your own program MAG which implements N(x), the modified arithmetic-geometric mean of 1 and x from the paper you mentioned. However Hugh Steers' program uses just one method agm2 to calculate the perimeter of the ellipse. I see that numerically both methods produce the same result but I don't understand how Hugh's method is related to the paper you cited. Could you explain that?

Kind regards
Thomas

                        
Re: Perimeter of an Ellipse (HP-42S)
Message #8 Posted by Gerson W. Barbosa on 24 May 2013, 2:28 p.m.,
in response to message #5 by hugh steers

Thanks for the C code, Hugh!

I am not proficient in C, but this appears to be easier to follow than RPL code. Anyway, my goal is a simple and short approximation formula that is accurate enough for practical cases, something that would give an error of a few meters in the length of the Earth's orbit for instance and takes about 20 or so steps on the HP-42S and no more than 1 second to run. Perhaps the approximate formula above meets this goal, but I haven't checked yet. Other approximate formulas are welcome, in case anyone knows about them (I've found only a few -- it appears measuring ellipses is not a popular sport :-)

Gerson.

      
Re: Perimeter of an Ellipse (HP-42S)
Message #9 Posted by Eduardo Duenez on 24 May 2013, 9:45 a.m.,
in response to message #1 by Gerson W. Barbosa

Quote:
Exact formula:
  p = 4*a*E((a + b)/(a + b*(b + 1)))
where
  E(x) is the complete elliptic integral of the second kind 
  a is the semi-major axis and
  b is the semi-minor axis

The above is not the correct relation between the complete elliptic integral E(x) and the perimeter p(a,b) of the ellipse. The argument "x" (usually denoted k, "the modulus") of E(k) is dimensionless, whereas (a+b)/(a+b(b+1)) is not. The correct relation is

p(a,b) = 4*a*E(c/a)
where a is the major semiaxis and
c = sqrt(a^2-b^2)
is the focal semidistance (so c/a is the eccentricity of the ellipse). Note that c/a is dimensionless.

Eduardo

EDIT: Coincidentally enough, Hugh Steers' C program above (posted as I was writing this, a few minutes earlier) documents the correct relation as well.

Edited: 24 May 2013, 11:25 a.m.

            
Re: Perimeter of an Ellipse (HP-42S)
Message #10 Posted by Gerson W. Barbosa on 24 May 2013, 1:45 p.m.,
in response to message #9 by Eduardo Duenez

Quote:
Quote:
Exact formula:
  p = 4*a*E((a + b)/(a + b*(b + 1)))
where
  E(x) is the complete elliptic integral of the second kind 
  a is the semi-major axis and
  b is the semi-minor axis

The above is not the correct relation between the complete elliptic integral E(x) and the perimeter p(a,b) of the ellipse. The argument "x" (usually denoted k, "the modulus") of E(k) is dimensionless, whereas (a+b)/(a+b(b+1)) is not. The correct relation is

p(a,b) = 4*a*E(c/a)

I'd found the formula

p(a,b) = 4*a*E(pi/2,e), where e is the eccentricity, sqrt(1 - (b/a)^2)
However, trying this on WolframAlpha for a = 3 and b =2 I get 14.567698609052450048.. instead of 15.865439589290589791..

http://www.wolframalpha.com/input/?i=4*3*E%28pi%2F2%2Csqrt%281-%282%2F3%29%5E2%29+%29

Oddly enough my wrong parameter works when the difference between a and b is 1 (at least for integer values) and when a = b or one of the semi-axis is 0, as in all of my examples. I should have tried a = 5 and b = 3, for instance, for which it fails. I wasn't able to find the right WolframAlpha syntax for this function. As soon as I find it, I will correct my original post. Thanks again.

Gerson.

                  
Re: Perimeter of an Ellipse (HP-42S)
Message #11 Posted by 聲gel Martin on 24 May 2013, 3:35 p.m.,
in response to message #10 by Gerson W. Barbosa

To add a log to the fire, here愀 the program on the 41 + SandMath:

Assumes a in Y and b in X, where ELI2 is the Incomplete Elliptic Integral of the second kind.-

01	LBL "ELIPER"
02	X^2
03	CHS
04	X<>Y
05	STO 05
06	X^2
07	+
08	SQRT
09	RCL 05
10	 /
11	90
12	DEG
13	ELI2
14	RCL 05
15	*
16	4
17	*
18	RCL 05
19	*
20	END

For a=3, b=2 I get: p = 14.56769861

in 1.8 seconds on the CL (gotta love this machine :-)

Cheers, 'AM

                        
Re: Perimeter of an Ellipse (HP-42S)
Message #12 Posted by Gerson W. Barbosa on 24 May 2013, 3:52 p.m.,
in response to message #11 by 聲gel Martin

Hi 聲gel,

What is the physical meaning of this result? Not the perimeter, I presume. That should be about 15.9, according to the best approximations I've seen.

Regards,

Gerson.

                              
Re: Perimeter of an Ellipse (HP-42S)
Message #13 Posted by 聲gel Martin on 24 May 2013, 4:20 p.m.,
in response to message #12 by Gerson W. Barbosa

This is what happens when rushing thru not paying attention to the details:-

I checked the formula and corrected the program accordingly, the result is now as yours:

http://www.wolframalpha.com/input/?i=perimeter+of+ellipse

01	LBL "ELIPER"
02	X<>Y
03	STO 05
04	 /
05	X^2
06	CHS
07	1
08	+
09	90
10	DEG
11	LEI2
12	RCL 05
13	ST+ X
14	*
15	ST+ X
16	END

3, ENTER^, 2, XEQ"ELIPER" --> 15,86543959

Cheers, 'AM

Edited: 24 May 2013, 4:24 p.m.

                                    
Re: Perimeter of an Ellipse (HP-42S)
Message #14 Posted by Gerson W. Barbosa on 24 May 2013, 4:48 p.m.,
in response to message #13 by 聲gel Martin

Thanks, 聲gel!

The parameter required by WolframAlpha is the square of the eccentricity, not the eccentricity. That's what I was doing wrong.

Cheers,

Gerson.

                  
Re: Perimeter of an Ellipse (HP-42S)
Message #15 Posted by Eduardo Duenez on 24 May 2013, 4:49 p.m.,
in response to message #10 by Gerson W. Barbosa

Quote:
I'd found the formula
p(a,b) = 4*a*E(pi/2,e), where e is the eccentricity, sqrt(1 - (b/a)^2)
However, trying this on WolframAlpha for a = 3 and b =2 I get 14.567698609052450048.. instead of 15.865439589290589791..

http://www.wolframalpha.com/input/?i=4*3*E%28pi%2F2%2Csqrt%281-%282%2F3%29%5E2%29+%29


I know why. The documentation linked to from Wolfram Alpha is inconsistent with Mathematica's definition.

What Wolfram Alpha's EllipticE(.) takes as an argument is the modulus k *squared*, and *not* k (in other words the eccentricity *squared*).

The convention used by W.Alpha is the one used by the function elliptic_ec(m) in venerable open-source program Maxima as correctly documented at:

http://maxima.sourceforge.net/docs/manual/en/maxima_16.html

Try:

http://maxima-online.org/#?in=4*3*elliptic_ec%20(5.0%2F9)%0A%20%09

Eduardo

                        
Re: Perimeter of an Ellipse (HP-42S)
Message #16 Posted by Gerson W. Barbosa on 24 May 2013, 4:54 p.m.,
in response to message #15 by Eduardo Duenez

Quote:
What Wolfram Alpha's EllipticE(.) takes as an argument is the modulus k *squared*, and *not* k (in other words the eccentricity *squared*).

Yes, we'd just found that out too (see my reply to 聲gel above). Thanks again!

Gerson.

                              
Re: Perimeter of an Ellipse (HP-42S)
Message #17 Posted by 聲gel Martin on 25 May 2013, 1:39 a.m.,
in response to message #16 by Gerson W. Barbosa

Exactly, so beating this dead horse just once more, this is the correct syntax for WolframAlpha:

4*3*E(pi/2,(1-(2/3)^2) )

[url:http://www.wolframalpha.com/input/?i=4*3*E%28pi%2F2%2C%281-%282%2F3%29^2%29+%29]

Best,

Edited: 25 May 2013, 1:41 a.m.


[ Return to Index | Top of Index ]

Go back to the main exhibit hall